Hydrodynamic Steering in Protein Association Revisited: Surprisingly

Aug 18, 2017 - We investigate the previously postulated hydrodynamic steering phenomenon, resulting from complication of molecular shapes, its magnitu...
3 downloads 10 Views 4MB Size
Article pubs.acs.org/JPCB

Hydrodynamic Steering in Protein Association Revisited: Surprisingly Minuscule Effects of Considerable Torques Jan M. Antosiewicz,† Kamil Kamiński,‡ and Maciej Długosz*,§ Department of Biophysics, Faculty of Physics, University of Warsaw, Ż wirki i Wigury 93, Warsaw 02-089, Poland Faculty of Physics, University of Warsaw, Pasteura 5, Warsaw 02-093, Poland § Centre of New Technologies, University of Warsaw, Stefana Banacha 2c, Warsaw 02-097, Poland † ‡

ABSTRACT: We investigate the previously postulated hydrodynamic steering phenomenon, resulting from complication of molecular shapes, its magnitude and possible relevance for protein−ligand and protein−protein diffusional encounters, and the kinetics of diffusion-controlled association. We consider effects of hydrodynamic interactions in a prototypical model system consisting of a cleft enzyme and an elongated substrate, and real protein−protein complexes, that of barnase and barstar, and human growth hormone and its binding protein. The kinetics of diffusional encounters is evaluated on the basis of rigid-body Brownian dynamics simulations in which hydrodynamic interactions between molecules are modeled using the bead-shell method for detailed description of molecular surfaces, and the first-passage-time approach. We show that magnitudes of steering torques resulting from the hydrodynamic coupling of associating molecules, evaluated for the studied systems on the basis of the Stokes−Einstein type relations for arbitrarily shaped rigid bodies, are comparable with magnitudes of torques resulting from electrostatic interactions of binding partners. Surprisingly, however, unlike in the case of electrostatic torques that strongly affect the diffusional encounter, overall effects of hydrodynamic steering torques on the association kinetics, while clearly discernible in Brownian dynamics simulations, are rather minute. We explain this result as a consequence of the thermal agitation of the binding partners. Our finding is relevant for the general understanding of a wide spectrum of molecular processes in solution but there is also a more practical aspect to it if one considers the low level of shape detail of models that are usually employed to evaluate hydrodynamic interactions in particle-based Stokesian and Brownian dynamics simulations of multicomponent biomolecular systems. Results described in the current work justify, in part at least, such a low-resolution description.



INTRODUCTION Virtually all molecular association processes that occur in solution and involve proteins, nucleic acids, or small molecules (generally we will refer to such processes using the term receptor−ligand interactions), can be divided in at least two steps. First, a diffusive encounter takes place resulting in formation of intermolecular van der Waals contacts, that is eventually followed by a series of configurational and conformational changes of the binding partners that lead to the final complex.1,2 For many receptor−ligand interactions, the molecules not only must be brought into close proximity, but also must be oriented correctly in space for a fruitful (i.e., leading to a functional complex) encounter to occur. Apart from the possibility that after diffusional encounter in a random orientation binding partners are held together by short-range van der Waals forces that allow them to assume a proper orientation, sometimes entertained in the case of large receptors binding small ligands,3 the role of long-range interactions, most prominently electrostatic, that act to orient the binding partners prior to their encounter is usually accepted. From the electrostatic standpoint molecules can be treated as certain distributions of electric charges, and their © XXXX American Chemical Society

long-range interactions can be analyzed in terms of the electric multipole expansion. Opposite net charges (monopoles) result in attraction of associating molecules, whereas the existence of molecular electric dipoles gives rise to electrostatic torques that will orient approaching molecules so that their dipole moments are antiparallel. More than three decades ago in order to name these fundamental effects, a concept of electrostatic steering was introduced.4−6 There is a vast literature documenting the existence of electrostatic steering and its role in the kinetics of molecular association.4−19 Another long-range force that should be taken into account as we consider mechanisms of molecular association results from solvent-mediated hydrodynamic interactions (HI), each moving solute molecule sets the solvent medium in motion and the resulting flux in the solvent medium then tends to move all of the other solute molecules. Since the pioneering works of Friedman,20 Deutch and Felderhof,21 and Wolynes and Deutch,22 it has been known that the rates of diffusional Received: June 20, 2017 Revised: July 25, 2017 Published: August 18, 2017 A

DOI: 10.1021/acs.jpcb.7b06053 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

evaluated for the bead-shell models of molecular surfaces31−37 using an accurate approach described and validated previously.36,38 Electrostatic interactions between binding partners, screened by dissolved ions, are evaluated on the basis of the DLVO approach.39 In the case of the prototypical enzyme− substrate system arbitrary distributions of charges are assigned to binding partners, whereas for protein systems effective charges40,41 derived from molecular potentials resulting from numerical solutions of the Poisson−Boltzmann equation are employed.42 Kinetics of diffusional encounters is evaluated on the basis of rigid-body Brownian dynamics simulations with hydrodynamic interactions37 and the first-passage-time approach.43 Brownian dynamics29 is a method which (unlike molecular dynamics-type techniques for fluid simulations such as dissipative particle dynamics,44 stochastic rotation dynamics,45 multiparticle collision dynamics,46 lattice Boltzmann techniques,47 that utilize some explicit representations of fluid particles, although on a larger scale than a single solvent molecule) follows only motions of solute molecules immersed in the fluid. Hydrodynamic coupling of solutes is accounted for by the mobility matrix, which is a function of solutes relative configuration. Mobility matrix constitutes the relation between translational forces and torques that act on solutes and their translational and rotational velocities and determines the thermal agitation of solutes. Brownian dynamics is a tool suitable for the task at hand, as switching hydrodynamic torques either on or off in Brownian dynamics simulations is straightforward and allows us to quantify directly their effects on the association kinetics. With the mobility matrix being the central element of BD there is also a direct connection between the methodology applied in the current work and the approach adopted earlier by Brune and Kim.23 This article is structured as follows. We begin by presenting the prototypical cleft enzyme−substrate and the two protein systems that we consider in the current work. We describe their hydrodynamic and electrostatic models and justify the choice of study subjects. Techniques for the calculation of mobility tensors and the evaluation of hydrodynamic torques on the basis of the Stokes−Einstein relations for arbitrarily shaped rigid bodies are next presented, followed by a short presentation of the first-passage-time approach and the Brownian dynamics simulation technique. Finally, numerical results regarding the effects of torques resulting from the hydrodynamic coupling between the binding partners on their association kinetics are presented and discussed.

encounters of molecules are influenced by hydrodynamic interactions and that hydrodynamic interactions are expected to slow the formation of encounter complexes. However, one should also allow for a possibility that hydrodynamic interactions can act to orient receptors and ligands in solution. This idea was originally introduced by Brune and Kim23 who coined the term hydrodynamic steering. They postulated that the complementarity of shapes of binding partners, usually observed in molecular complexes, results in their better orientation upon closer approach. If binding partners are being pushed together by random thermal fluctuations they will be more favorable oriented for reaction as translations and rotations of molecules are coupled by hydrodynamic interactions which is a consequence of their complementary shapes. Conversely, if the binding partners are driven apart by fluctuations, their relative orientation is distorted by hydrodynamic interactions. According to Brune and Kim,23 such a behavior, i.e., a ligand that moves toward the receptor experiences a torque which corrects its misalignment, while a ligand that moves away from the receptor becomes even further misaligned, should increase the rate of the encounter. Using an effectively two-dimensional system consisting of models of a cleft enzyme and its substrate approaching each other with translational velocities calculated based on the equipartition theorem, for which they numerically solve the resistance problem,24 these authors indeed demonstrated the existence of hydrodynamic torques that act to orient the approaching molecules. They, however, pointed out, that detailed dynamic simulations are needed to quantify hydrodynamic coupling effects on reaction rate. To the best of our knowledge, since the seminal work of Brune and Kim,23 the problem of the hydrodynamic steering in molecular association has been raised only a few times.25−28 Brownian dynamics (BD) simulations29 that were performed for rather crude model systems and an approximate, far-field description of hydrodynamic interactions, and which do not allow for a direct evaluation of the role of torques resulting from the hydrodynamic coupling26−28 suggest that their effect on rates of diffusional encounters is only moderate, and diminishes with the increasing strength of electrostatic interactions between the binding partners. Hydrodynamic steering effects in the barnase−barstar system were studied through the analysis of the relative rotational velocity of the proteins approaching each other under the influence of electrostatic interactions using a rigorous hydrodynamic approach.25 It was shown that hydrodynamic interactions favor the proper (i.e., the one observed in their crystallographic complex) relative orientation of barnase and barstar purely as a consequence of the complementary shapes of the proteins’ interfaces, however, no quantification of hydrodynamic effects on the association kinetics was possible with the methodology employed. On a related note, the authors of the recently published discrete element simulation study postulated, using the work of Brune and Kim23 as a reference, that hydrodynamic steering torques may enhance gelation in colloidal suspensions by aligning clusters of aggregated particles.30 Picking up where Brune and Kim left off, in the current work we quantify effects of hydrodynamic steering on the association kinetics in the same prototypical cleft enzyme−substrate system that was considered by these authors and in two real protein− protein complexes, that of barnase and barstar, and human growth hormone and its binding protein. Hydrodynamic properties (mobility and resistance tensors, and hydrodynamic forces and torques acting on molecules) of these systems are



THEORETICAL AND COMPUTATIONAL METHODS Prototypical Cleft Enzyme−Substrate System. As described in the previous section, we choose to study hydrodynamic effects in the same model cleft enzyme− substrate system that was previously described by Brune and Kim,23 for the reason that its geometry results in significant hydrodynamic coupling between translational and rotational motions (see below). As it is shown in Figure 1, the enzyme is an assembly of two spheres, of radius 20 Å, with a center-tocenter distance of 50 Å, with a space between spheres corresponding to the enzyme cleft. The substrate is an elongated capsule (a cylinder capped on both ends by a hemisphere with a radius equal to the radius of the cylinder) of 40 Å length and 5 Å radius, so that when the substrate is oriented parallel to the cleft of the enzyme (i.e., the long axis of the substrate is perpendicular to the line connecting the centers of enzyme spheres, Figure 1A) it fits precisely into the cleft. B

DOI: 10.1021/acs.jpcb.7b06053 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

substrate (red line in Figure 2) and the z-axis of the laboratory coordinate frame equals either 0 (perpendicular encounter) or 90 (parallel encounter) degrees (see Figure 1). Binding partners in the prototypical cleft enzyme−substrate system are either neutral or charged. In the case of charged partners, we consider the following three scenarios: (1) central charges of similar magnitude but opposite signs are assigned to the molecules, (2) binding partners are assigned dipoles of similar magnitude; the enzyme dipole consists of two opposite charges, positioned on the line connecting centers of enzyme spheres, each charge at a distance of 15 Å from the center of the geometry of the enzyme; the substrate dipole consists of two opposite charges positioned on the long axis of the substrate, each charge at a distance of 15 Å from the substrate geometry center; as a consequence, the minimum of the enzyme− substrate electrostatic interaction is observed for their perpendicular encounter (ϕ = 0°, Figure 1), (3) binding partners are assigned dipoles of similar magnitude; the enzyme dipole consists of two opposite charges, positioned on the line connecting centers of enzyme spheres, each charge at a distance of 15 Å from the center of the geometry of the enzyme; the substrate dipole consists of two opposite charges positioned on the line perpendicular to the long axis of the substrate and to the x of the laboratory coordinate frame, each charge at a distance of 15 Å from the substrate geometry center, so that the minimum of the enzyme−substrate electrostatic interaction is observed for their parallel encounter (ϕ = 90°, Figure 1). Enzyme−substrate electrostatic interactions, screened by dissolved ions, and resulting electrostatic forces and torques acting on the substrate at a particular position and in a particular orientation, are evaluated using the following form of the potential energy:39

Figure 1. Bead-shell models of the enzyme (orange) and the substrate (green). (A) An encounter complex in which the substrate is oriented parallel to the cleft of the enzyme. (B) An encounter complex in which the substrate is oriented perpendicular to the cleft of the enzyme. ϕ denotes the angle between the line connecting centers of enzyme spheres and the long axis of the substrate. Drawings were done using the VMD package.48 Electric dipoles are assigned to the enzyme and the substrate in such a way that their orientation in an encounter complex is the same (as shown in the inset), regardless the orientation (either parallel or perpendicular) assumed by the binding partners.

Similarly as Brune and Kim,23 we assume that the motions of the substrate are limited to translations along a single axis and rotations around this axis. We define the laboratory coordinate frame so that its origin coincides with the center of geometry of the enzyme (whose position and orientation is fixed in the laboratory frame), and its z-axis is parallel to the line that connects centers of enzyme spheres (Figure 2). Substrate

Nenzyme Nsubstrate

,(x , ϕ) =

∑ ∑ i=1

j=1

Q iQ j ϵrij(x , ϕ)

e−κrij(x , ϕ) (1)

where Qi (Qj) symbol denotes enzyme (substrate) charges, rij denotes the distance between the ith charge of the enzyme and jth charge of the substrate that depends on the position (x) and orientation (ϕ) of the substrate in the laboratory coordinate frame, ϵ is the dielectric constant of the solvent (78.54), and κ is the inverse of Debye length.39 We consider different magnitudes of enzyme and substrate central charges that result in electrostatic interaction energies upon encounter varying between 1kBT and 10 kBT when κ is set to infinity (which corresponds to zero ionic strength). Electric dipoles assigned to molecules vary between 250 D and 750 D which is a reasonable range for protein-like molecules.49 Calculations and simulations for the prototypical system were performed for ionic strength values of 0, 50, and 150 mM. Apart from Brownian dynamics simulations, we evaluate (as described below) and analyze the torque that acts on the substrate due to its hydrodynamic interactions with the enzyme, as a function of its position in the laboratory coordinate frame x and the rotation angle ϕ. Evaluation of hydrodynamic interactions in the studied system requires calculations of mobility matrices for a ligand molecule moving in the vicinity of the fixed receptor (see below). These matrices are evaluated for bead-shell models31,32,34,35,37 of molecules. Enzyme model consists of 5992 spherical subunits while the substrate bead-shell model consists of 824 spherical subunits, of equal radii of 0.55 Å, equidistributed on their surfaces (Figure 1). Bead-shell models were

Figure 2. Diffusional association in the prototypical cleft enzyme− substrate system. Substrate (green) undergoes translational (translations along the x-axis of the laboratory coordinate frame) and rotational (rotations about an angle ϕ around the x-axis of the laboratory coordinate frame) diffusion, under the influence of the potential resulting from its electrostatic interactions with the enzyme (orange). Substrate diffuses in a finite domain, between planar boundaries located at x = xo (absorbing boundary) and x = L (reflective boundary).

geometry center translates along the x-axis of the laboratory frame and the substrate rotates around this axis (about the angle ϕ) while its long axis stays in the yz-plane (Figure 2). For the purpose of Brownian dynamics simulations on the basis of which we derive the enzyme−substrate association kinetics and analyze the role of hydrodynamic coupling, we assume that the substrate diffuses in a finite domain, enclosed between two planar boundaries located at x = xo and x = L. xo is chosen as a minimal distance between the substrate and the enzyme that do not result in the overlap of their surfaces, regardless the orientation of the substrate. Boundary located at x = L is reflective. For the association to take place, the substrate must cross the boundary located at x = xo, oriented in such a way that the angle ϕ between the long axis of the C

DOI: 10.1021/acs.jpcb.7b06053 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

second system) due to hydrodynamic interactions with its binding partner. Hydrodynamic torques are calculated for protein−protein configurations generated using the procedure depicted schematically in Figure 4.

created using algorithms and methods described and applied previously.37,38 Protein−Protein Systems. Two protein−protein complexes were considered for the purpose of the current work. The first complex is that of barnase and barstar50 (Figure 3A).

Figure 3. Protein systems considered in the current work. (A) The complex of barnase (blue) and barstar (green) (PDB ID 1BRS50). (B) The complex of the human growth hormone (orange) and a fragment of the hormone-binding protein (cyan) (PDB ID 1A2255). Atomic structures of proteins and their molecular surfaces are shown at the top of each panel. Hydrodynamic bead-shell models of proteins are depicted at the bottom of each panel. Drawings were done using UCSF Chimera56 and VMD packages.48 Figure 4. (A) Generation of protein−protein configurations for which hydrodynamics torques are evaluated. Hemispheres represent the binding partners. One of the proteins forming the complex (either barstar or the growth hormone-binding protein) is translated along the line connecting the geometric centers of proteins in their crystallographic complex and rotated around this line by the angle ϕ in a direction indicated with an arrow. (B) Example barnase−barstar configurations. In all of the configurations, barstar (green) is translated by 12 Å from its position in the crystallographic complex so that regardless the orientation of proteins there is no overlap between their bead-shell models. (C) Example configurations of human growth hormone and its binding protein. In all of the configurations, growth hormone-binding protein (cyan) is translated by 27 Å from its position in the crystallographic complex so that regardless the orientation of proteins there is no overlap between their bead-shell models.

Barnase is an extracellular ribonuclease and barstar is its intracellular inhibitor. Barnase binds strongly to negatively charged RNA molecules, and its interface is very polar and contains a number of positively charged lysine and arginine residues. Barstar, on the other hand, bears a strong negative overall charge, and its interface is negatively charged. Barnase− barstar complex is one of the fastest binding protein−protein complexes known, with an association rate constant on the order of 108M−1s−1.51,52 In the case of the barnase−barstar system, their strong electrostatic complementarity enhances the binding speed over the Smoluchowski rate that applies to association by simple random search and diffusion.52−54 While both barnase and barstar can be classified as globular proteins (Figure 3A), a local complementarity of shapes of their interfaces can be observed.25 The second complex that we consider here is that of the human growth hormone and a fragment of the hormonebinding protein (Figure 3B).55 Association in this system is roughly 3 orders of magnitude slower than in the case of the barnase−barstar system,57 due to the lack of significant electrostatic interactions. Hormone-binding protein possesses a distinct, boomerang-like shape (Figure 3B), which makes this particular complex interesting from the hydrodynamic standpoint. Both for the barnase−barstar system, and for human growth hormone and its binding protein, the rigid-body description of binding partners works reasonably well, when applied to evaluate the kinetic aspect of their association.52,58−60 In the case of the protein systems, we evaluate and analyze torques acting on one of the molecules forming a complex (barstar in the first system and growth hormone-binding protein in the

First, coordinates of proteins taken from the crystallographic structure of a given complex are transformed, so that the line connecting centers of geometry of proteins in the complex coincides with the x-axis of the laboratory coordinate frame. The origin of the laboratory coordinate frame is chosen as the center of geometry of the fixed molecule. Next, the selected protein (barstar in the first system and growth hormonebinding protein in the second system) is translated away from its fixed partner along the x-axis, in predefined increments. For each new position x, the chosen protein is rotated around the xaxis, counterclockwise (Figure 4), by the angle ϕ that covers the range from 0° to 360°. Some example configurations of barnase−barstar and human growth hormone−hormone-binding protein are given in Figure 4. Hydrodynamic torques are evaluated only for configurations in which hydrodynamic models (see below) of binding partners do not overlap. Only translations along the x-axis and rotations around this axis are D

DOI: 10.1021/acs.jpcb.7b06053 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

the molecule by small beads.34,71,72 Resulting hydrodynamic representations of proteins are presented in Figures 3 and 4. Shells of generated models consist of 5726 beads in the case of barnase, 4623 in the case of barstar, 9162 in the case of growth hormone, and 10193 in the case of growth hormone-binding protein. Radius of a single bead in the shell is the same for all proteins and equals 0.41 Å. Mobility Tensor of a Molecule Moving in the Presence of the Fixed Receptor. Detailed description of the procedure that allows one to calculate the spatially dependent mobility tensor, μ(q⃗), for a molecule that moves in the vicinity of an unmoving receptor can be found in refs 34,37,73, hence we present it here only briefly. Initially, we consider a general form of the mobility tensor, 4 , for a system consisting of two arbitrarily shaped rigid objects (molecules) represented with bead-shell models, in a particular configuration. To derive 4 hydrodynamic interactions between particular beads are first evaluated at the two-body Rotne− Prager level38,74,75 and then the projection operation is applied that effectively constrains velocities of beads within a given object so that their relative motions are not possible.38,76 We note that formally only the analytical expressions for translational self-mobilities of spherical particles and couplings between translations of different particles are referred to as the Rotne−Prager tensor.74 However, we employ an extended set of expressions that accounts also for rotational selfmobilities, and couplings between rotations and rotations, and translations and rotations of different spheres; these expressions are in the usual form of power series of the inverse of distances between spherical elements (r) up to the order of r−3.38,75 Resulting 4 is represented with a 12 × 12 matrix:

considered. Thus, we choose configurations of proteins that are possible to occur in solution as they are sterically allowed, and that are geometrically close to the final crystallographic complex. Electrostatic interactions within a pair of proteins in a given configuration were evaluated with the effective charges approach.40,41 We created atomistic proteins’ models based on X-ray structures of their complexes.50,55 As crystallographic complexes contain only heavy atoms of proteins, we used the PROPKA utility,61−64 a well established tool for empirical prediction of protein pKa values, to add hydrogens to the heavy atoms of proteins with the solution pH set to 8.0 in the case of the barnase−barstar complex,51 and 7.4 in the case of the human growth hormone-growth hormone-binding protein complex.57 Partial charges and radii, that are needed for calculation o proteins’ electrostatic potentials were assigned to proteins’ atoms with the PDB 2PQR tool65 on the basis of the AMBER99 force-field parameters66,67 provided with the tool. The resulting total protein net charges were +2e for barnase and −5e for barstar, and −3e for human growth hormone and −4e for its binding protein. We derived effective charges positioned at the heavy atoms of proteins by fitting the electrostatic potential resulting from the Debye−Hückel approximation39 to the molecular potential obtained as a numerical solution of the linearized Poisson−Boltzmann equation.42 The electrostatic potentials of proteins were computed using APBS68 at ionic strength of 150 mM with the dielectric constants of the proteins and the solvent set respectively to 1 and 78.54, and with dielectric boundary determined as the van der Waals surface of a given protein. For each protein, fitting of the electrostatic potential was performed in the volume outside of a protein, defined as a 8 Å-thick skin around the protein with the lower boundary of the skin determined as the van der Waals surface of the protein inflated by 2 Å, using methods and tools described previously.26,69 Protein−protein electrostatic interactions, screened by dissolved ions, and resulting forces and torques acting on proteins in a given configuration, are conveniently evaluated on the basis of effective charges, using the form of potential energy defined above (eq 1). As effective charges are derived by fitting Poisson−Boltzmann potentials, they account for electrostatic effects resulting from immersing a particular protein in the high-dielectric solvent (i.e., solvation or hydration effects). Eq 1 lacks a term accounting for the so-called desolvation penalty40,70 resulting from interactions of charges of a particular protein with the low-dielectric cavity occupied by its binding partner, which are neglected here due to their short-range character40,70 and the fact that we further consider electrostatic forces and torques in encounter complexes in which van der Waals surfaces of proteins are rather well separated. Short-range nonpolar (i.e., hydrophobic) effects are also neglected in evaluation of protein−protein interactions. We also note that desolvation penalties and nonpolar effects are rather difficult to calculate with a reasonable accuracy using analytical expressions. Similarly as in the prototypical enzyme−substrate system described above, torques resulting from hydrodynamic interactions in protein systems are evaluated using bead-shell models of molecules. These are created using an approach described previously, in which a coarse-grained representation of the protein is first generated by replacing each amino acid residue with a single sphere and subsequently covering areas of surfaces of these spheres that are accessible from the outside of

⎛ μ TT ⎜ 11 ⎜ μ RT ⎜ 11 4=⎜ TT ⎜ μ21 ⎜⎜ RT ⎝ μ21

μ11TR μ12TT μ12TR ⎞ ⎟ μ11RR μ12RT μ12RR ⎟⎟ TR TT TR ⎟ μ21 μ22 μ22 ⎟ RR RT RR ⎟ μ21 μ22 μ22 ⎟⎠

(2)

consisting of 3 × 3 blocks μklij . Indices i ∈ (1,2) and j ∈ (1,2) correspond to molecules. Blocks μklii (k ∈ (T, R) and l ∈ (T, R)) correspond to translations (TT) and rotations (RR) of the ith molecule and their couplings (TR/RT). Blocks μkli ≠ j (k ∈ (T, R) and l ∈ (T, R)) describe how the motions of the ith molecule affect motions of the jth molecule. Again, couplings between translations (TT), rotations (RR), translations and rotations (TR/RT) are accounted for. Motions of each molecule are evaluated in a coordinate frame whose origin coincides with the geometric center of that molecule; the corresponding axes of these molecule-centered frames are parallel to each other. Having calculated 4 we further assume that the first molecule (the ligand) diffuses freely while the position and orientation of the second molecule (the receptor) are fixed by some external forces and torques, so that its translational and rotational velocities vanish (that may correspond for instance to a situation where receptor is much larger than the ligand and thus its movement occurs on a time scale that is significantly longer than a time scale characteristic for ligand′s motion, or the receptor is embedded in some biological structure, such as a membrane; we note that such an assumption does not limit the generality of conclusions drawn from the results of the current work). In this scenario, the mobility tensor of the ligand, E

DOI: 10.1021/acs.jpcb.7b06053 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B moving relative to the fixed receptor, is represented with a 6 × 6 matrix that can be derived from the components of 4 as34,37,73

Generally, μ is represented by a matrix composed of four 3 × 3 blocks that vary with the configuration of the system (described with the vector q⃗) due to HIs between the ligand and the receptor. In the case of the prototypical cleft enzyme−substrate system

described above, in which motions of the substrate are limited to the x-axis, q⃗ = (x, ϕ)†, μ takes the form:

where components μii, i ∈ (x, y, z) correspond to translations in directions x, y, and z, components μii, i ∈ (ϕ, β, γ) correspond to rotations around the axes that are parallel to the axes x, y, z of the laboratory coordinate frame and pass through the geometric center of the substrate, and components μij, i ≠ j, i, j ∈ (x, y, z, ϕ, β, γ) correspond to translation-translation, translation-rotation, and rotation-rotation couplings. As translations in the x-direction and rotations around the x-axis are not hydrodynamically coupled to translations in directions y and z, and rotations around axes y and z, the μ matrix for the prototypical cleft enzyme−substrate system that we consider here reduces simply to ⎛ μ xx (x , ϕ) μ xϕ (x , ϕ) ⎞ ⎟ μ(x , ϕ) = ⎜⎜ ⎟ ϕx ϕϕ ( , ) ( , ) μ x ϕ μ x ϕ ⎝ ⎠

Due to the lack of symmetry, coupling between all translations and rotations should now be considered. Similarly as in the case of the cleft enzyme−substrate system, we define the μo tensor evaluated for the moving protein in the absence of the binding partner. For an arbitrarily shaped rigid body, such as a protein, μo is not necessarily diagonal, and while components of μo do not depend on the position of the moving protein in the laboratory coordinate frame, they vary (with the exception of ϕϕ μxx o and μo in the case considered here) with the varying orientation of the moving molecule in the laboratory coordinate system, which is a consequence of general properties of mobility tensors.77−79 Using the μxx component of μo we define the translational hydrodynamic radius of the moving molecule (be that the model substrate or a protein) as

(5)

For this system, we define also two additional mobility tensors: ⎛ μ xx 0 ⎞ o ⎟ μo = ⎜ ⎜ 0 μϕϕ ⎟ ⎝ o ⎠

(6)

and ⎛ μ xx (x , ϕ) ⎞ 0 ⎟ μ μxϕ = 0, μϕx = 0 (x , ϕ) = ⎜⎜ ⎟ μϕϕ (x , ϕ)⎠ 0 ⎝

σH = (6πημoxx )−1

(9)

with η being the viscosity of he solvent. σH is 11.7 Å in the case of the prototypical substrate, 16.7 Å in the case of barsta,r and 25.3 Å in the case of growth hormone-binding protein. All hydrodynamic calculations were conducted at the temperature of 298.15 K with the solvent viscosity set to 0.089 poise using an in-house software. Fixing the position and orientation of one of the molecules that may for example correspond to the situation in which its diffusion occurs on a much longer time scale than the diffusion of its partner, do not affect the analysis and conclusions that we present further as mobility matrix given with eq 3 correctly accounts for the hydrodynamic coupling between the binding partners. However, such an approach simplifies evaluation of hydrodynamic interactions in BD simulations that, as described

(7)

Mobility tensor μo describes the substrate in the absence of the receptor. Similarly, as for the tensor μ, μo is evaluated in the coordinate system with an origin located at the geometric center of the substrate and axes that are parallel to the corresponding axes of the laboratory coordinate system. Mobility tensor μμxϕ=0,μϕx=0 is derived from μ by zeroing the off-diagonal components μxϕ/ϕx that correspond to the coupling of translations and rotations induced by the presence of the fixed receptor. In the case of protein systems, the μ tensor of the protein moving along the x-axis of the laboratory coordinate frame and rotating around this axis, is of the form: F

DOI: 10.1021/acs.jpcb.7b06053 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B below, are performed using precomputed mobility tensors stored on a grid whose dimensions correspond to the number of degrees of freedom of the simulated system. Hydrodynamic Torques. In what follows, we consider forces and torques evaluated relative to the center of geometry of a moving body. Hydrodynamic forces and torques that act on a molecule moving in a potential resulting from its interactions with a fixed receptor are evaluated on the basis of the generalized Stokes’ law that relates forces, -⃗el , and torques, ;el⃗ , acting on an arbitrarily shaped rigid body translating and rotating in an incompressible viscous solvent (we use here the el index as the forces and torques considered in the current work are of electrostatic origin), and its translational, =⃗ , and rotational, ω⃗ , velocities. The moving molecule is assigned a resistance tensor ξ that can be derived from its mobility tensor using the inverse operation, i.e., ξ = (μ)−1. The resistance tensor of a molecule rotating and translating in three dimensions, ξ is represented with a 6 × 6 matrix of form: −1

ξ = (μ)

⎛ ξ TT ξ TR ⎞ ⎟ = ⎜⎜ ⎟ RT ξ RR ⎠ ⎝ξ

where -⃗H and ;⃗H correspond respectively to (additional) hydrodynamic forces and torques exerted by the fluid on the moving molecule in the presence of the receptor. Considering the prototypical cleft enzyme−substrate system described above, in which translations in the x-direction and rotations around the x-axis are not hydrodynamically coupled to translations in directions y and z, and rotations around axes y and z, the expression for the hydrodynamic torque reduces simply to

where we have distinguished the term resulting from the hydrodynamic coupling of translations and rotations of the moving molecule resulting from its hydrodynamic interactions with the receptor. Action of this hydrodynamic coupling torque may result in hydrodynamic steering and its effects on the association kinetics are the focus of our work. Components ;H, y and ;H, z are zero due to distributions of charges that we have chosen to assign to binding partners that result in vanishing translational and rotational velocities in directions orthogonal to x. In the case of protein−protein systems, however, the x component of the hydrodynamic torque must be evaluated based on the general form of eq 14, accounting for possible couplings of translational and rotational motions in different directions, as well as y and z components of translational and rotational velocities that are not zero. The recipe for calculation of hydrodynamic torques acting on a molecule is thus the following. For a given enzyme−substrate (or protein−protein) configuration we calculate electrostatic forces, and mobility tensors μ and μo. Having calculated electrostatic forces and torques, and mobilities, we evaluate translational and rotational velocities of the moving molecule in the absence of hydrodynamic interactions between the binding partners, on the basis on the relation:

(10)

where 3 × 3 blocks TT, RR, and TR/RT correspond respectively to translations, rotations, and their couplings, and depend on the position and orientation of the molecule relative to the fixed receptor. In the case of a creeping motion the following linear relation holds true: ⎛ -⃗ ⎞ ⎛ TT TR ⎞⎛ ⎞ ξ ⎜ el ⎟ = ⎜ ξ ⎟⎜ =⃗ ⎟ ⎜ ⎟ ⎜ ⃗ ⎟ RT ξ RR ⎠⎝ ω⃗ ⎠ ⎝ ;el ⎠ ⎝ ξ

(11)

Similar relation can be written if we consider a hypothetical situation in which the molecule moves in exactly the same electrostatic force field, resulting in forces -⃗el and torques ;el⃗ , but with the receptor absent, i.e., its hydrodynamic interactions with the receptor are neglected (the friction tensor ξo is now evaluated as the inverse of the μo tensor): ⎛ -⃗ ⎞ ⎛ ξ TT ξ TR ⎞⎛ ⃗ ⎞ o ⎟ =o ⎜ el ⎟ = ⎜ o ⎜ ⎟ ⎜ ⃗ ⎟ ⎜ ξ RT ξ RR ⎟⎜ ω⃗ ⎟ ⎝ ;el ⎠ ⎝ o o ⎠⎝ o ⎠

⎛ =⃗ ⎞ ⎛ μTT μTR ⎞⎛ -⃗ ⎞ o ⎟⎜ el ⎟ ⎜ o⎟ = ⎜ o ⎜ ⎟ ⎜ RT RR ⎟⎜ ⃗ ⎟ μo ⎠⎝ ;el ⎠ ⎝ ωo⃗ ⎠ ⎝ μo

We then invert μ, μo to obtain resistance tensors ξ, ξo and use eq 14 to obtain ;⃗H . Hydrodynamic torques in the studied systems were evaluated using an in-house software. First-Passage-Time Approach to the Diffusion Controlled Association. Only a short description of the firstpassage-time approach that we apply to evaluate the kinetics of diffusional association is given below; for a more detailed formulation we point the reader to the ref 43. We consider a fixed receptor and a ligand that diffuses in a finite domain, of volume V, enclosed between two boundaries, external δΩ and internal δω. We assume that the δΩ boundary is reflective, while the δω boundary is absorbing. Absorption of the ligand on the δω boundary corresponds to its association with the fixed receptor. The ligand molecule, generated initially between δω and δΩ diffuses in the potential resulting from its interactions with the receptor and eventually reaches δω where it is absorbed. At any time, position of the ligand molecule in the laboratory frame (and, in general, also its orientation) is

(12)

Denoting changes in translational and rotational velocities resulting from “turning on” hydrodynamic interactions in the system respectively by Δ=⃗ and Δω⃗ , we may write eq 11 as ⎛ -⃗ ⎞ ⎛ TT TR ⎞⎛ ⃗ ξ = − Δ=⃗ ⎞ ⎜ el ⎟ = ⎜ ξ ⎟ ⎟⎜ o ⎜ ⃗ ⎟ ⎜ ξ RT ξ RR ⎟⎜ ω⃗ − Δω⃗ ⎟ ⎠⎝ o ⎠ ⎝ ;el ⎠ ⎝

(13)

Combining eqs 12 and 13 leads to ⎛ ξ TT ξ TR ⎞⎛ ⃗ ⎞ ⎟⎜ Δ= ⎟ −⎜⎜ RT RR ⎟⎝ ξ ⎠ Δω⃗ ⎠ ⎝ξ ⎛⎛ ξ TT ξ TR ⎞ ⎛ TT TR ⎞⎞⎛ ⃗ ⎞ ⎛ ⃗ ⎞ ξ ξ ⎟ =o o ⎟ ⎜ o ⎟ ⎜ ⎟ ≡ ⎜ H⎟ = ⎜⎜ −⎜ ⎜ RT RR ⎟ ⎜ RT RR ⎟⎟⎜ ω⃗ ⎟ ⎜ ⃗ ⎟ ξo ⎠ ⎝ ξ ξ ⎠⎠⎝ o ⎠ ⎝ ;H ⎠ ⎝⎝ ξo

(16)

(14) G

DOI: 10.1021/acs.jpcb.7b06053 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B given with the vector q⃗. Probability of finding the ligand at time t at q⃗, P(q⃗, t), can be obtained as the solution of the Smoluchowski equation:80 ⎛ ⎛ ⎞⎞ 1 δtP(q ⃗ , t ) = kBT ∇q ⃗ ⎜⎜μ(q ⃗)⎜∇q ⃗ P(q ⃗ , t ) + P(q ⃗ , t )∇q ⃗ ,(q ⃗)⎟⎟⎟ k T ⎝ ⎠⎠ ⎝ B

⟨R⃗ q (⃗ Δt )(R⃗ q (⃗ Δt ))† ⟩ = 2kBTμ(q ⃗)

First, a nonreactive system is considered, in which reflective boundaries enclose the diffusion domain. BD simulations are performed for such a system and generated positions, q⃗, of the ligand molecule are collected forming an equilibrium ensemble. This ensemble corresponds to the probability distribution Peq(q⃗o) Next, the studied system is activated and its internal boundary is made absorbing. Multiple BD trajectories are generated. These trajectories start from the configurations of the system taken from the equilibrium ensemble. Each trajectory is propagated until the ligand associates with the receptor (i.e., the ligand is absorbed on the internal boundary). Times elapsed between the beginning of each trajectory and the association are collected and the global first passage time is calculated as their average, according to eq 20. All Brownian dynamics simulations conducted for the purpose of the current work utilize the algorithm given above (eq 21) and an in-house software.37 In the prototypical enzyme−substrate system considered here, the substrate is allowed to diffuse between two planes (see Figure 2), located at x = 25 Å and x = 125 Å (thus L = 100 Å and considering the diameter of the substrate equal 10 Å the studied system is rather dilute). The plane located at x = 125 Å is reflective, whereas the plane located at x = 25 Å is either reflective or absorbing for a particular orientation of the substrate relative to the enzyme (i.e., either parallel or perpendicular). L is so chosen that μ converges to μo when x approaches the reflective boundary (μo is the mobility tensor for the substrate evaluated neglecting its HI with the enzyme). Reflective boundaries are implemented in BD simulations as described previously.37,88 Three kinds of BD simulations were performed: simulations with full enzyme−substrate hydrodynamic interactions (in which the mobility tensor of the substrate is described as in eq 5), simulations without enzyme−substrate hydrodynamic interactions (the substrate mobility tensor described as in eq 6), and simulations in which the hydrodynamic coupling between translations and rotations of the substrate that results from the presence of the enzyme is turned off (eq 7). To avoid performing hydrodynamic calculations on-the-fly during BD simulations, all mobility tensors are precalculated on twodimensional curvilinear grids (x, ϕ) with x ∈ ⟨25 Å, 125 Å⟩ and ϕ ∈ ⟨0°,360°), corresponding to different positions and orientations of the substrate in the laboratory coordinate frame. Resolution of grids in x dimension is set to 1 Å, and in the ϕ dimension it is set to 5°. Precalculated tensors are stored in memory. During a BD simulation, for a current position and orientation of the substrate in the laboratory coordinate frame, its mobility tensor is evaluated on the basis of stored tensors by means of a bilinear interpolation. BD simulations were conducted at a temperature of 298.15 K. The integration step (Δt in eq 21) was set to 5 ps. BD simulations were conducted for different charges or dipole moments assigned to the enzyme and the substrate and for different values of the ionic strength (see above). Electrostatic interactions were derived from the potential energy given with eq 1. For a given set of charges or dipoles and for a particular value of ionic strength we first generate an ensemble of equilibrium configurations of the enzyme−substrate system. This ensemble consists of 106 configurations collected from BD trajectories of a system in which both boundaries are reflective. Then, for a given set of charges or dipoles assigned to the binding partners, and for a given value of ionic strength, 105 trajectories of a reactive/

(17)

evaluated assuming that P(q⃗, t) equals zero on δω, the flux of ligand molecules vanishes on δΩ, and that at t = 0 the ligand is located between δω an δΩ at some qo⃗ . In the above equation kB is the Boltzmann constant, T is the temperature, μ(q⃗) is the spatially dependent mobility of the ligand molecule, and ,(q ⃗) denotes either the potential affecting the ligand or the potential energy of the system in configuration given with q⃗;∇q⃗ denotes the gradient with respect to the components of q⃗. Introducing the survival probability function, S(qo⃗ , t), i.e., the probability that the ligand molecule located initially at q⃗o has not yet reacted at time t: S(qo⃗ , t ) =

∫V P(q ⃗ , t|qo⃗ , 0)dq ⃗

(18)

where P(q⃗, t|qo⃗ ,0) is the conditional probability that if the ligand molecule was at q⃗o at time 0 then it will be at q⃗ at time t, resulting from the solution of the Smoluchowski eq 17, and the integration is over the whole volume V of the domain enclosed between δω and δΩ, one may evaluate the time τ(qo⃗ ) of the receptor−ligand association (the mean first-passage time43,81,82) as τ(qo⃗ ) =

∫0



S(qo⃗ , t )dt

(19)

We assume that for t < 0 the studied system is nonreactive, i.e., both δω and δΩ are reflective, and only an activation (or triggering) at t = 0 makes the association between the ligand and the receptor possible, such a situation is usually encountered in experiments involving pump−probe techniques.83−86 Before the activation, the probability of finding the ligand at qo⃗ follows the equilibrium distribution Peq(q⃗o) and thus the average time for association, τ, (or the global firstpassage time81) can be evaluated by integrating out the dependence of the mean first-passage time (eq 19) on the initial position of the ligand: τ=

∫V Peq(qo⃗ )τ(qo⃗ )dqo⃗

(20)

Again, the integration is over the whole volume V of the domain enclosed between δω and δΩ. Brownian Dynamics Simulations. Global first-passage times are evaluated from BD simulations, in which configuration space trajectories q⃗(t) of the ligand molecule that undergoes translational and rotational diffusion, either free or under the influence of electrostatic and hydrodynamic interactions with the fixed receptor, and eventually associates with the receptor are generated. BD trajectories are composed of successive displacements Δq, taken over a time step Δt, according to an algorithm given with29,87 Δq ⃗ = kBT ∇q ⃗ μ(q ⃗)Δt + μ(q ⃗)∇q ⃗ ,(q ⃗)Δt + R⃗ q (⃗ Δt )

(22)

(21)

where R⃗ q⃗(Δt) is a random displacement with a Gaussian distribution function whose average value is 0 and whose variance-covariance is given by († denotes the transpose operation): H

DOI: 10.1021/acs.jpcb.7b06053 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B activated system are generated starting from different configurations, chosen randomly from the equilibrium ensemble. Each trajectory is propagated until the substrate associates with the enzyme. Durations of these trajectories are collected for analysis. Parallel and perpendicular encounters are simulated separately.



RESULTS Substrate Mobility in the Cleft Enzyme−Substrate System. To illustrate the extent to which hydrodynamic interactions with the cleft enzyme modify the mobility of the substrate, in Figures 5 and 6 we show different components of

Figure 6. Components of the mobility tensor of the model ligand, evaluated in the presence of the receptor, as functions of the position and orientation of the ligand in the laboratory coordinate frame (xx, translations; ϕϕ, rotations; xϕ, translation-rotation coupling). kB is the Boltzmann constant, T is the temperature, and σH is the translational hydrodynamic radius of the ligand.

perpendicular to parallel. Rotational component of μ reaches maximum for ϕ = 90° while minima are observed for ϕ = 0° and ϕ = 180°. However, changes in μϕϕ observed with the varying ϕ are not as significant as in the case of μxx, an approximately 15% increase in μϕϕ is observed when orientation of the substrate relative to the cleft of the enzyme is changed from perpendicular to parallel. Hydrodynamic coupling of substrate translations and rotations (the μxϕ = μϕx component of μ), induced by the presence of the enzyme, vanishes both for parallel and perpendicular encounters (Figure 5), while its magnitude is maximal for ϕ values between 30° and 40° (140° and 150°). Notably, in the case of the prototypical cleft-enzyme substrate system, coupling components of μ, obtained for some ϕ values, are of comparable order of magnitude to its translational and rotational components. In Figure 6 we show the dependence of the substrate mobility tensor on the distance from the absorbing boundary xo (Figure 2). All components of the μ tensor are monotonic functions of the distance between the substrate and the fixed enzyme. The largest variation with the distance from the absorbing boundary is observed in the case of the translational mobility, μxx, which at the absorbing boundary drops to approximately 30% (60%) of the value calculated at the reflective boundary for parallel (perpendicular) encounters. μxx curves obtained for different values of ϕ increase with the increasing distance from the absorbing boundary and all converge to μxxo at the external boundary of the considered finite domain. In the case of the rotational mobility, changes are much smaller, as at the absorbing boundary, μϕϕ drops only by approximately 10% in the case of parallel encounters (and

Figure 5. Components of the mobility matrix of the model ligand molecule, evaluated in the presence of the receptor, (xx, translations; ϕϕ, rotations; xϕ, translation-rotation coupling), μ, as functions of the ligand’s orientation in the receptor−ligand encounter complex (given with the angle ϕ, see Figure 2). kB is the Boltzmann constant, T is the temperature, and σH is the translational hydrodynamic radius of the ligand.

the substrate mobility tensor μ (eq 5) as functions of the position and the orientation of the substrate relative to the enzyme. Components μxx, μxϕ (which in the considered case is equal to μϕx), and μϕϕ are scaled respectively by σ2H, σ1H, and σ0H to allow for their direct comparison. Figure 5 shows the rotation angle dependence of the substrate’s mobility upon encounter with the enzyme wherein the center of the substrate is located at the absorbing boundary xo (see Figure 2). All components of μ vary with the rotation angle ϕ. Translational component μxx is maximal for ϕ = 90°, for a parallel orientation of the substrate relative to the cleft of the enzyme (Figure 1A); μxx is minimal for ϕ = 0° and ϕ = 180° (when the substrate is oriented perpendicular to the cleft of the enzyme, Figure 1B). Variation in the value of the translational component with ϕ is quite significant as about 50% increase in μxx is observed when orientation of the substrate relative to the cleft of the enzyme is changed from I

DOI: 10.1021/acs.jpcb.7b06053 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B about 25% in the case of perpendicular encounters) when compared with the value calculated at the reflective boundary. Dependence of μϕϕ on the distance from the absorbing boundary reaches plateau (that corresponds to μoϕϕ) for distances above 2σH. Thus, the rotational dynamics of the substrate is less affected by HI with the enzyme than its translational dynamics. Moreover, the effects of HI on the translational mobility of the substrate are of longer range than the effects of HI on its rotational mobility. As for the coupling between translations and rotations, maximal absolute values of μxϕ are observed at the absorbing boundary and for distances of the substrate from the absorbing boundary that are above 3σH (where the hydrodynamic radius σH is defined in eq 9), the coupling between translations and rotations of the ligand molecule practically vanishes. Hydrodynamic Torques in the Cleft Enzyme−Substrate System. We begin the analysis of torques resulting from the hydrodynamic coupling of motions of the model cleft enzyme and its substrate by considering the situation in which the binding partners are assigned central charges of equal magnitudes but opposite signs. As a result, an electrostatic attraction causes the substrate to approach the enzyme along the line connecting their centers (which is also the x axis of the laboratory coordinate frame) but there is no electrostatic torque that would act to rotate the substrate relative to the enzyme. Hydrodynamic coupling torque (eq 15) that act on the substrate approaching the enzyme has only one component, x, that we show in Figure 7 as a function of the orientation of the substrate (the ϕ angle, Figure 2), the strength of the enzyme− substrate attraction and the ionic strength; center of the substrate is located at the absorbing boundary xo (see Figure 2). According to data shown in Figure 7, for a given value of the rotation angle and ionic strength, the hydrodynamic coupling torque acting on the substrate increases monotonically with the increasing magnitude of the attractive force; for a given value of the rotation angle and the magnitude of the attractive force, the hydrodynamic coupling torque acting on the substrate decreases monotonically with the increasing ionic strength. Such behavior is expected considering the linear relation between forces (torques) and velocities defined by eqs 12−14. Much more interesting is, however, the dependence of the hydrodynamic coupling torque on the substrate’s rotation angle, ϕ. Torque vanishes when ϕ equals 0°, 90° or 180°, i.e., values for which the translation-rotation coupling terms in the substrate’s mobility tensor vanish (Figure 5 and eq 15). This means, in particular, that when the substrate approaching the enzyme is oriented properly (ϕ = 90° and the substrate fits the cleft of the enzyme), there is no hydrodynamic coupling torque that would act to disrupt their encounter. If, however, the substrate approaching the enzyme is not oriented properly (i.e., ϕ ≠ 90°), the torque resulting from the hydrodynamic coupling always act to orient the substrate so that it fits the cleft of the enzyme. Thus, as already described by Brune and Kim,23 the hydrodynamic coupling of translations and rotations of the binding partners may result in their proper orientation upon encounter if they fit together, i.e., if their shapes are complementary. A further step in our analysis of the possible role of hydrodynamics in the model system, is to consider binding partners possessing permanent electric dipole moments. As described above, we consider two types of enzyme−substrate encounters, parallel (when the substrate fits the cleft of the enzyme) and perpendicular (Figure 1). In both cases, the

Figure 7. Variations of the hydrodynamic coupling torque acting on the substrate in the model enzyme−substrate system, with the orientation of the substrate relative to the enzyme. Electrostatic force acting on the substrate results from central electric charges (Q) of similar magnitudes but opposite signs, assigned to the binding partners. There is no net electrostatic torque acting on the substrate. I denotes ionic strength, kB is the Boltzmann constant, T is the temperature, and e is an elementary charge.

dipole moment assigned to the enzyme is aligned along the line connecting centers of enzyme’s spheres while the substrate dipole moment is either aligned along the long axis of the substrate (the perpendicular encounter) or along the line perpendicular to the long axis of the substrate in the yz-plane of the laboratory coordinate frame (the parallel encounter). In Figure 8 we show variations of electrostatic and hydrodynamic (the total hydrodynamic torque and the contribution to the hydrodynamic torque resulting from the translation-rotation coupling, see eq 15) torques acting on the substrate as a function of the orientation of the substrate relative to the enzyme (angle ϕ, Figure 2), the ionic strength and the orientation of its permanent dipole; center of the substrate is located at the absorbing boundary xo (see Figure 2). Overall, magnitudes of torques resulting from hydrodynamic interactions are comparable to magnitudes of electrostatic torques, and as electrostatic torques decrease with the increasing ionic strength so do hydrodynamic torques, which again is a conseqeunce of linear relations given with eqs 12−14. Electrostatic torques act to orient dipoles of the binding partners in an antiparallel fashion, and thus, depending on the orientation of the substrate’s dipole, to restore enzyme− substrate orientation corresponding to either the parallel or the perpendicular encounter (Figure 1). However, as we have shown above, the complementarity of shapes of the binding partners results in hydrodynamic coupling that favors only the J

DOI: 10.1021/acs.jpcb.7b06053 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 8. Variations of hydrodynamic and electrostatic torques acting on the substrate in the model enzyme−substrate system, with the orientation of the substrate relative to the enzyme. Electrostatic forces and torques acting on the substrate result from electric dipoles of similar magnitudes (750 D) assigned to the binding partners. Left panel: substrate’s electrostatic dipole oriented as in the case of a parallel encounter. Right panel: substrate’s electrostatic dipole oriented as in the case of a perpendicular encounter (see text for details). I denotes ionic strength, kB is the Boltzmann constant, T is the temperature.

Table 1. Effects of Hydrodynamic Interactions on the Global First-Passage Time (τ) for the Association in the Model EnzymeSubstrate Systema I (mM) 0

50

150

|Q|(e) 0.0 1.86 2.62 3.71 4.55 5.25 5.87 1.86 2.62 3.71 4.55 5.25 5.87 1.86 2.62 3.71 4.55 5.25 5.87

o τHI ∥ /τ∥

1.42 1.40 1.39 1.38 1.38 1.35 1.31 1.40 1.40 1.40 1.40 1.40 1.40 1.42 1.41 1.41 1.40 1.41 1.40

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02

o τHI ⊥ /τ⊥

1.83 1.80 1.75 1.70 1.70 1.65 1.58 1.80 1.80 1.75 1.72 1.70 1.66 1.84 1.81 1.80 1.79 1.77 1.77

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02

=0,μϕx=0)



τHI(μ ∥

1.05 1.03 1.03 1.02 1.00 1.00 1.00 1.04 1.04 1.03 1.02 1.02 1.01 1.05 1.04 1.04 1.04 1.04 1.03

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

/τHI ∥

0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01

τHI(μ ⊥



=0,μϕx=0)

0.94 0.94 0.94 0.93 0.92 0.91 0.91 0.94 0.94 0.95 0.95 0.94 0.94 0.94 0.95 0.95 0.95 0.95 0.95

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

/τHI ⊥

0.01 0.01 0.01 0.02 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01

a Electric central charges (Q) of similar magnitudes but opposite signs are assigned to binding partners. Upper index o-simulations without hydrodynamic interactions; upper index HI-simulations with full enzyme-substrate hydrodynamic interactions, upper index HI(μxϕ = 0,μϕx = 0) simulations in which components of the substrate’s mobility matrix that correspond to hydrodynamic couplings between translations and rotations are set to zero. Lower indices ∥ and ⊥ correspond respectively to parallel and perpendicular encounters (see Figure 1). I denotes ionic strength, e is an elementary charge. Errors were evaluated based on standard deviations of means, using the uncertainty propagation method. In all cases τo∥ = τo⊥.

K

DOI: 10.1021/acs.jpcb.7b06053 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Table 2. Effects of Hydrodynamic Interactions on the Global First Passage Time (τ) for the Association in the Model Enzyme− Substrate Systema I (mM) 0

50

150

p (D) 250 500 750 250 500 750 250 500 750

o τHI ∥ /τ∥

1.36 1.34 1.32 1.37 1.35 1.33 1.37 1.37 1.36

o τHI ⊥ /τ⊥

± ± ± ± ± ± ± ± ±

0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02

1.73 1.55 1.54 1.77 1.68 1.55 1.78 1.75 1.70

± ± ± ± ± ± ± ± ±

0.02 0.02 0.03 0.02 0.02 0.02 0.02 0.02 0.02

=0,μϕx=0)



τHI(μ ∥

1.03 1.01 1.00 1.05 1.03 1.02 1.05 1.03 1.02

± ± ± ± ± ± ± ± ±

/τHI ∥

0.01 0.01 0.02 0.01 0.01 0.01 0.01 0.01 0.01

=0,μϕx=0)



τHI(μ ⊥

0.94 0.95 0.88 0.94 0.95 0.95 0.94 0.94 0.94

± ± ± ± ± ± ± ± ±

/τHI ⊥

0.01 0.01 0.02 0.01 0.01 0.01 0.01 0.01 0.01

a Electric dipoles (p) of similar magnitudes are assigned to binding partners. Upper index o-simulations without hydrodynamic interactions; upper index HI-simulations with full enzyme-substrate hydrodynamic interactions, upper index HI(μxϕ = 0,μϕx = 0) simulations in which components of the substrate’s mobility matrix that correspond to hydrodynamic couplings between translations and rotations are set to zero. Lower indices ∥ and ⊥ correspond respectively to parallel and perpendicular encounters (see Figure 1). I denotes ionic strength, D denotes the unit of Debye. Errors were evaluated based on standard deviations of means, using the uncertainty propagation method. In all cases τo∥ = τo⊥.

parallel encounter. Thus, in the case of the parallel encounter, the torque resulting from the hydrodynamic coupling of translations and rotations will act in concert with electrostatic torque to accelerate enzyme−substrate association, while in the case of the perpendicular encounter electrostatic torques will act to properly orient the binding partners against the torque resulting from the hydrodynamic coupling. This interplay between electrostatics and hydrodynamics is clearly evidenced in Figure 8. In the next section, we quantify effects of the hydrodynamic coupling torque on the kinetics of parallel and perpendicular encounters. Hydrodynamic Steering Effects on the Enzyme− Substrate Association Kinetics. As described above, three setups of BD simulations were considered with regard to the treatment of enzyme−substrate HI. We performed simulations without HI (using the mobility tensor μo, eq 6), with full HI (the mobility tensor μ(x, ϕ), eq 5), and with zeroed coupling terms in the mobility tensor (the mobility tensor μμxϕ=0,μϕx=0(x, ϕ), eq 7), i.e., effectively neglecting the hydrodynamic coupling torque (eq 15). We evaluate the kinetics of parallel and perpendicular encounters (Figure 1) for different distributions of charges assigned to binding partners and different values of the ionic strength. Results of BD simulations are gathered in Tables 1 and 2. Let us first consider the situation where the binding partners are either neutral or central charges of similar magnitude but opposite signs are assigned to them (Table 1). As a result a net attractive force acts between the enzyme and the substrate, but there is no electrostatic torque that would act to rotate the substrate in the yz-plane of the laboratory coordinate frame (Figure 2). Overall, the general effect of HI, that we quantify in terms of the ratio of global first-passage times calculated based on BD simulations with and without enzyme−substrate HI, τHI/τo, is to slow down the encounter, regardless the relative orientation of the binding partners in the complex. However, the kinetics of the parallel encounter is less affected by HI than the kinetics of the perpendicular encounter; for instance, in the case of the neutral system, HI slow down the kinetics of the parallel encounter by roughly 40% when compared with the kinetics obtained in the absence of enzyme−substrate HI, while the kinetics of the perpendicular encounter is slowed down by roughly 80% (Table 1).

When the attractive force acting between the binding partners is increased (with the increasing magnitude of charges and for a particular magnitude of charges with the decreasing ionic strength) the magnitude of the HI effects on the encounter kinetics decreases, for both parallel and perpendicular encounters. Our focus is, however, on the effects of the hydrodynamic coupling torque (i.e., hydrodynamic steering effects) on the association kinetics, and thus, in Table 1 we compare global first-passage times obtained on the basis of BD simulations in which hydrodynamic coupling torque acting on the substrate is switched off with global first-passage time obtained on the basis of simulations with full enzyme−substrate HI. In the case of the parallel encounter, switching-off hydrodynamic coupling torques results in longer global first-passage times, hydrodynamic steering speeds up the association by helping the binding partners in assuming the correct orientation, in which the substrate fits the cleft of the enzyme (Figure 1). In the case of the perpendicular encounter, switching off coupling torques results in shorter association times as a consequence of the fact, that, as we have already described above, coupling torques act to prevent the substrate to assume orientation that is perpendicular to the cleft of the enzyme. What is however surprising is that the relative magnitude of hydrodynamic steering effects on the kinetic of association is quite small, with maximal values of roughly 5% in the case of the parallel encounter of neutral molecules, and roughly 10% in the case of the perpendicular encounter of charged molecules at zero ionic strength. (Table 1). Magnitude of the hydrodynamic coupling torque increases with the increasing strength of electrostatic attraction between the binding partners, as evidenced in Figure 7. Relative magnitude of hydrodynamic steering effects on the association kinetics decreases with the increasing strength of the enzyme−substrate attraction (resulting either from the increase in magnitudes of charges or the decrease in the ionic strength) in the case of the parallel encounter and for the largest charges considered, at zero ionic strength, effects of the hydrodynamic steering vanish completely. In the case of the perpendicular encounter opposite is observed at zero ionic strength with the relative magnitude of hydrodynamic steering effects increasing with the increasing attraction between molecules. At higher ionic strengths, in the case of the perpendicular encounter, the relative magnitude of hydrodynamic steering effects is essentially independent of the L

DOI: 10.1021/acs.jpcb.7b06053 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 9. Components of mobility tensors of barstar (left panel) and the growth hormone-binding protein (right panel), μ, as functions of their orientation in encounter complexes with their binding partners (given with the angle ϕ, see Figure 4). kB is the Boltzmann constant, T is the temperature, and σH denotes translational hydrodynamic radii of proteins.

magnitude of charges assigned to binding partners. We note that for both types of the encounter, dependence of the magnitude of hydrodynamic steering effects on the strength of electrostatic interactions is rather weak. When central charges are replaced with permanent dipole moments, and as a result, in addition to net forces (which are now either repulsive or attractive, depending on the relative enzyme−substrate orientation), electrostatic torques are introduced in the system, the overall picture does not change (Table 2). Hydrodynamic interactions slow down the encounter, the parallel encounter to the smaller extent than the perpendicular one. The hydrodynamic coupling torque accelerates the parallel encounter while slowing down the perpendicular encounter. These effects are, however, no larger than in the case described above. Again, only a weak dependence of the relative magnitude of hydrodynamic steering on the strength of electrostatic interactions is observed. Hydrodynamic Steering in Protein Association. Having analyzed hydrodynamic steering effects on the kinetics of association in the prototypical enzyme−substrate system, which, as evidenced above are surprisingly (considering the magnitudes of hydrodynamic coupling torques) small, we now proceed with the analysis of protein systems. We consider encounter complexes of barnase and barstar (Figure 3A), and human growth hormone and its binding protein (Figure 3B) that we generated using a procedure described above and depicted in Figure 4A. We consider here relative configurations of proteins in which all values of the rotation angle are sterically allowed. Some of these configurations are shown in Figure 4B,C. While hydrodynamic interactions with their binding partners introduce orientation dependent modifications to mobilities of moving proteins, the extent of these modifications, that can be quantified on the basis of the mobility data presented in Figure

9 is much smaller in comparison to the enzyme−substrate system described in the previous section. The maximal relative change in the translational mobility observed upon the variation of the rotation angle is roughly 15% in the case of barstar, and roughly 20% in the case of the hormone-binding protein. Even smaller changes, below 5% are observed for both proteins in the case of rotational mobilities. Thus, in both real protein systems that we consider here, the angular dependence of translational and rotational mobilities resulting from hydrodynamic interactions with their binding partners, is much less pronounced than in the case of the model enzyme−substrate system. Additionally, the coupling components of μ obtained for different values of ϕ are now order of magnitude smaller than its translational and rotational components. Hydrodynamic torques (similarly as in the case of the enzyme−substrate system we differentiate the (total) hydrodynamic torque and the hydrodynamic coupling torque) acting on barnase and hormone-binding protein are depicted in Figure 10 as functions of distance from and orientation relative to their respective binding partners. According to data shown in this figure, for both proteins, the (total) hydrodynamic torque consists almost exclusively of the coupling term, as both for barstar and hormone-binding protein, curves corresponding to the two torques are nearly the same. This results from the fact that for considered positions and orientations, rotational mobilities of moving proteins are only slightly affected by hydrodynamic interactions RR with their binding partners so that the (ξRR o −ξ ) ω⃗ o term in the expression for ;⃗H (eq 14) almost vanishes. In the case of barstar, magnitudes of the hydrodynamic coupling torque are comparable with magnitudes of the electrostatic torque. In the case of hormone-binding protein, the hydrodynamic coupling torque vanishes in configurations of the binding partners for which the protein−protein distance is so chosen that hormoneM

DOI: 10.1021/acs.jpcb.7b06053 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

performed for the model system, we may conclude that these effects are not likely to be significant. This is discussed in more detail in the next section.



DISCUSSION We investigated possible effects of the hydrodynamic steering on the kinetics of the diffusion-controlled association in receptor−ligand systems of biological relevance. Hydrodynamic steering may occur due to the hydrodynamic coupling of translations and rotations of the binding partners, provided that they possess complementary shapes,23 which is often true in biological systems. In such a case every time molecules are pushed together, a hydrodynamic torque arises that rotate molecules allowing a better fit of their shapes. It was postulated that this coupling torque should influence the association kinetics23 similarly as it has been observed in the case of the electrostatic steering.4−19 We consider the hydrodynamic steering phenomenon in three different molecular systems: a prototypical cleft enzyme− substrate system, and two protein systems that differ with regard to shapes of the binding partners and their electrostatic properties: the barnase and barstar systems in which the binding partners are of globular shapes and their association is driven by the strong electrostatic attraction and complementarity of charge distribution at their interfaces, and the system consisting of human growth hormone and its binding protein of a distinct boomerang-like shape, in which electrostatic interactions are rather weak. In the case of the enzyme− substrate system, a global lock and key shape complementarity of the binding partners results in a substantial hydrodynamic coupling of their translational and rotational motions (see mobility tensor’s components depicted in Figures 5 and 6). In protein systems, however, hydrodynamic coupling of translations and rotations of binding partners (Figure 9) is less significant. Assuming that molecules drift toward each other with translational and rotational velocities resulting from their electrostatic interactions we calculate and analyze hydrodynamic coupling torques that act on the substrate in the prototypical system, for different charge distributions assigned to the binding partners and for different values of the ionic strength. We show, that magnitudes of hydrodynamic torques are of similar order as magnitudes of electrostatic torques, and that hydrodynamic torques aid molecules in assuming a relative orientation in which their shapes fit together, regardless the action of electrostatic torques (see Figures 7 and 8). Further, we perform Brownian dynamics simulations and evaluate how switching off the hydrodynamic coupling torque affects enzyme−substrate association kinetics (Tables 1 and 2). Generally, we observe that the association is accelerated by hydrodynamic steering, if the shapes of binding partners forming the complex are complementary so that they fit together. In the absence of such a complementarity, the hydrodynamic steering would actually result in a slower association kinetics. While these results and the widely accepted consensus that the formation of receptor−ligand complexes usually involves some degree of shape complementarity allow one to conclude that molecular complexes are optimal from the hydrodynamic standpoint, hydrodynamic steering effects on the association kinetics are minuscule at best, especially, that the hydrodynamic coupling in the prototypical system is strongly exaggerated in comparison with real protein systems (as the barnase−barstar, and growth hormone−hormone binding

Figure 10. Magnitudes of hydrodynamic and electrostatic torques (evaluated at 150 mM ionic strength) acting on barstar (top) and the human growth hormone-binding protein (bottom) in encounter complexes with their binding partners, as functions of the distance from the binding partner (i.e., the displacement of the protein from its position in the X-ray complex) and the orientation of binding partners (given with the angle ϕ, see Figure 4). Note that for shorter distances (i.e., 10 Å in the case of the barnase−barstar system and 20 Å in the case of the growth hormone and its binding protein) some relative orientations of proteins in encounter complexes are sterically prohibited, hence no torques are given for some values of the rotation angle. kB is the Boltzmann constant, T is the temperature.

binding protein is able to rotate freely (Figure 4), which is a consequence of weaker electrostatic attraction, and thus lower velocity assumed upon approach, than in the barnase−barstar case. We do not perform a detailed analysis of directions of hydrodynamic torques and their role upon encounter in protein systems as configurations that we constructed using the procedure presented above do not correspond to actual paths of protein−protein association. Our aim is simply to illustrate that magnitudes of hydrodynamic coupling torques acting in real protein system are not likely to be larger than the magnitudes calculated in the prototypical enzyme−substrate system. Moreover, we argue that there is no need to perform Brownian dynamics simulations to quantify effects of the hydrodynamic steering on the kinetics of the association in the considered protein systems. On the basis of the analysis N

DOI: 10.1021/acs.jpcb.7b06053 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B ORCID

protein systems considered here). Even though magnitudes of hydrodynamic torques that arise if we consider molecules that move toward each other with net drift velocities resulting from their electrostatic interactions are comparable with magnitudes of electrostatic torques (which are known to determine rates of diffusional encounters in biological systems4−19), their effect on the association kinetics is almost completely destroyed by thermal agitation. In dynamic simulations that we have performed, at any given moment Brownian forces may either push associating molecules together, which, due to the hydrodynamic translational-rotational coupling, results in their better alignment, or move them apart which in turn distorts their relative orientation; as a consequence the net effect of the hydrodynamic steering on the association kinetics in biological receptor−ligand systems is, as we have shown, not likely to be significant. Brownian dynamics simulations were performed for a model enzyme−substrate system which is effectively two-dimensional, with one rotational and one translational degree of freedom, as shown in Figure 2. Such a simplification is introduced for a 2fold reason. First, our aim is to perform dynamic simulations of the same exact system that was investigated previously by Brune and Kim who performed single point hydrodynamic calculations assuming the same two degrees of freedom.23 Such a construction of the model system does not, however, limit the generality of our conclusions regarding the possible role of hydrodynamic steering in biological systems, with an additional benefit of making the analysis of results more comprehensible. Second, assuming two degrees of freedom simplifies dynamics simulations and lowers significantly their computational cost. As described above, in Brownian dynamics simulations we employ mobility matrices precomputed on a curvilinear (twodimensional) grid. The number of grid points is 7272. Calculation of the 6 × 6 mobility matrix corresponding to a particular grid point requires an inversion of a general mobility matrix of size 40896 × 40896 (as described above, this size results from the resolution of bead-shell models employed to model surfaces of molecules) - which requires significant amount of time and memory. Allowing the substrate to translate and rotate in three dimensions (i.e., considering its six degrees of freedom) would result in an enormous number of calculations that would need to be performed to generate the grid of mobility matrices. As an alternative one could consider calculations of mobility matrices on-the-fly, during BD simulations that would require fast algorithms for solving large linear systems of equations. Such algorithms have started to emerge only recently.35,36 Our work contributes to the general debate on the role and significance of HI in biological systems.20−23,25−28,34,37,89−97 Moreover, the main result of our work, namely the lack of significant hydrodynamic steering effects on the kinetics of the diffusional encounter of arbitrarily shaped molecules, may (at least to some extent) serve as a justification of different computational studies96,98−100 in which shapes of biomolecules are neglected and their hydrodynamic interactions are modeled on the basis of analytical expressions and numerical approaches, either exact96,100 or approximate,98,99 that are appropriate for spherical particles.



Maciej Długosz: 0000-0002-2534-0744 Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS Authors acknowledge support from National Science Centre, Poland (UMO-2014/13/B/NZ1/00004). REFERENCES

(1) Johnson, K. A. Advances in Transient-State Kinetics. Curr. Opin. Biotechnol. 1998, 9, 87−89. (2) Johnson, K. A. Transient-State Kinetic Analysis of Enzyme Reaction Pathways. Enzymes 1992, 20, 1−61. (3) Chou, K. C.; Zhou, G. P. Role of the Protein Outside Active Site on the Diffusion-Controlled Reaction of Enzymes. J. Am. Chem. Soc. 1982, 104, 1409−1413. (4) Allison, S. A.; Srinivasan, N.; McCammon, J. A.; Northrup, S. H. Diffusion-Controlled Reactions Between a Spherical Target and Dumbbell Dimer by Brownian Dynamics Simulation. J. Phys. Chem. 1984, 88, 6152−6157. (5) Northrup, S. H.; Boles, J. O.; Reynolds, J. Electrostatic Effects in the Brownian Dynamics of Association and Orientation of Heme Proteins. J. Phys. Chem. 1987, 91, 5991−5998. (6) Head-Gordon, T.; Brooks, C. L. The Role of Electrostatics in the Binding of Small Ligands to Enzymes. J. Phys. Chem. 1987, 91, 3342− 3349. (7) Sines, J. J.; Allison, S. A.; McCammon, J. A. Point Charge Distributions and Electrostatic Steering in Enzyme/Substrate Encounter: Brownian Dynamics of Modified Copper/Zinc Superoxide Dismutases. Biochemistry 1990, 29, 9403−9412. (8) Tan, R. C.; Truong, T. N.; McCammon, J. A.; Sussman, J. L. Acetylcholinesterase: Electrostatic Steering Increases the Rate of Ligand Binding. Biochemistry 1993, 32, 401−403. (9) Luty, B. A.; Wade, R. C.; Madura, J. D.; Davis, M. E.; Briggs, J. M.; McCammon, J. A. Brownian Dynamics Simulations of Diffusional Encounters Between Triose Phosphate Isomerase and Glyceraldehyde Phosphate: Electrostatic Steering of Glyceraldehyde Phosphate. J. Phys. Chem. 1993, 97, 233−237. (10) Kozack, R. E.; D’Mello, M. J.; Subramaniam, S. Computer Modeling of Electrostatic Steering and Orientational Effects in Antibody-Antigen Association. Biophys. J. 1995, 68, 807−814. (11) Wade, R. C.; Gabdoulline, R. R.; Lüdemann, S. K.; Lounnas, V. Electrostatic Steering and Ionic Tethering in Enzyme-Ligand Binding: Insights from Simulations. Proc. Natl. Acad. Sci. U. S. A. 1998, 95, 5942−5949. (12) Wlodek, S. T.; Shen, T.; McCammon, J. A. Electrostatic Steering of Substrate to Acetylcholinesterase: Analysis of Field Fluctuations. Biopolymers 2000, 53, 265−271. (13) Herzfeld, J.; Tounge, B. NMR Probes of Vectoriality in the Proton-Motive Photocycle of Bacteriorhodopsin: Evidence for an ’Electrostatic Steering’ Mechanism. Biochim. Biophys. Acta, Bioenerg. 2000, 1460, 95−105. (14) Myles, T.; Le Bonniec, B. F.; Betz, A.; Stone, S. R. Electrostatic Steering and Ionic Tethering in the Formation of Thrombin-Hirudin Complexes: The Role of the Thrombin Anion-Binding Exosite-I. Biochemistry 2001, 40, 4972−4979. (15) Huang, Y. M. M.; Huber, G.; McCammon, J. A. Electrostatic Steering Enhances the Rate of cAMP Binding to Phosphodiesterase: Brownian Dynamics Modeling. Protein Sci. 2015, 24, 1884−1889. (16) Blöchliger, N.; Xu, M.; Caflisch, A. Peptide Binding to a PDZ Domain by Electrostatic Steering via Nonnative Salt Bridges. Biophys. J. 2015, 108, 2362−2370. (17) Wang, X.; Putkey, J. A. PEP-19 Modulates Calcium Binding to Calmodulin by Electrostatic Steering. Nat. Commun. 2016, 7, 13583. (18) Mohan, R. R.; Huber, G. A.; Morikis, D. Electrostatic Steering Accelerates C3d:CR2 Association. J. Phys. Chem. B 2016, 120, 8416− 8423.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]; Phone: +48 22 5543 692; Fax: +48 22 5540 801. O

DOI: 10.1021/acs.jpcb.7b06053 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Algorithm and Its Application to the Nucleosome. Biopolymers 2001, 58, 106−115. (42) Gilson, M. K.; Honig, B. Energetics of Charge-Charge Interactions in Proteins. Proteins: Struct., Funct., Genet. 1988, 3, 32−52. (43) Szabo, A.; Schulten, K.; Schulten, Z. First Passage Time Approach to Diffusion Controlled Reactions. J. Chem. Phys. 1980, 72, 4350−4357. (44) Hoogerbrugge, P. J.; Koelman, J. M. V. A. Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics. Europhys. Lett. 1992, 19, 155−160. (45) Ihle, T.; Kroll, D. M. Stochastic Rotation Dynamics. I. Formalism, Galilean Invariance, and Green−Kubo Relations. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2003, 67, 066705. (46) Malevanets, A.; Kapral, R. Mesoscopic Model for Solvent Dynamics. J. Chem. Phys. 1999, 110, 8605−8613. (47) Ladd, A. J. C. Short-Time Motion of Colloidal Particles: Numerical Simulation via a Fluctuating Lattice-Boltzmann Equation. Phys. Rev. Lett. 1993, 70, 1339−1342. (48) Humphrey, W.; Dalke, A.; Schulten, K. VMD - Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (49) Felder, C. E.; Prilusky, J.; Silman, I.; Sussman, J. L. A Server and Database for Dipole Moments of Proteins. Nucleic Acids Res. 2007, 35, W512−W521. (50) Buckle, A. M.; Schreiber, G.; Fersht, A. R. Protein-Protein Recognition: Crystal Structural Analysis of a Barnase-Barstar Complex at 2.0Å Resolution. Biochemistry 1994, 33, 8878−8889. (51) Schreiber, G. A.; Fersht, A. R. Energetics of Protein-Protein Intercations: Analysis of the Barnase-Barstar Interface by Single Mutations and Double Mutant Cycles. J. Mol. Biol. 1995, 248, 478− 486. (52) Qin, S.; Pang, X.; Zhou, H. X. Automated Prediction of Protein Association Rate Constant. Structure 2011,.19174410.1016/ j.str.2011.10.015 (53) Dong, F.; Vijayakumar, M.; Zhou, H. X. Comparison of Calculation and Experiment Implicates Significant Electrostatic Contributions to the Binding Stability of Barnase and Barstar. Biophys. J. 2003, 85, 49−60. (54) Lee, L. P.; Tidor, B. Optimization of Binding Electrostatics: Charge Complementarity in the Barnase-Barstar Protein Complex. Protein Sci. 2001, 10, 362−377. (55) Clackson, T.; Ultsch, M. H.; Wells, J. A.; de Vos, A. M. Structural and Functional Analysis of the 1:1 Growth Hormone:Receptor Complex Reveals Molecular Basis for Receptor Affinity. J. Mol. Biol. 1998, 277, 1111−1128. (56) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. UCSF Chimera - a Visualization System for Exploratory Research and Analysis. J. Comput. Chem. 2004, 25, 1605−1612. (57) Cunningham, B. C.; Wells, J. A. Comparison of a Structural and a Functional Epitope. J. Mol. Biol. 1993, 234, 554−563. (58) Spaar, A.; Dammer, C.; Gabdoulline, R. R.; Wade, R. C.; Helms, V. Diffusional Encounter of Barnase and Barstar. Biophys. J. 2006, 90, 1913−1924. (59) Gabdoulline, R. R.; Wade, R. C. Simulation of the Diffusional Association of Barnase and Barstar. Biophys. J. 1997, 72, 1917−1929. (60) Schluttig, J.; Alamanova, D.; Helms, V.; Schwarz, U. S. Dynamics of Protein-Protein Encounter: A Langevin Equation Approach with Reaction Patches. J. Chem. Phys. 2008, 129, 155106. (61) Li, H.; Robertson, A. D.; Jensen, J. H. Very Fast Empirical Prediction and Interpretation of Protein pKa Values. Proteins: Struct., Funct., Genet. 2005, 61, 704−721. (62) Bas, D. C.; Rogers, D. M.; Jensen, J. H. Very Fast Prediction and Rationalization of pKa Values for Protein-Ligand Complexes. Proteins: Struct., Funct., Genet. 2008, 73, 765−783. (63) Olsson, M. H. M.; Søndergaard, C. R.; Rostkowski, M.; Jensen, J. H. PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa Predictions. J. Chem. Theory Comput. 2011, 7, 525−537.

(19) Xu, J.; Xie, Y.; Lu, B.; Zhang, L. Charged Substrate and Product Together Contribute Like a Nonreactive Species to the Overall Electrostatic Steering in Diffusion-Reaction Processes. J. Phys. Chem. B 2016, 120, 8147−8153. (20) Friedman, H. L. A Hydrodynamic Effect in the Rates of Diffusion Controlled Reactions. J. Phys. Chem. 1966, 70, 3931−3933. (21) Deutch, J. M.; Felderhof, B. U. Hydrodynamic Effect in Diffusion-Controlled Reaction. J. Chem. Phys. 1973, 59, 1669−1671. (22) Wolynes, P. G.; Deutch, J. M. Slip Boundary Conditions and the Hydrodynamic Effect on Diffusion Controlled Reactions. J. Chem. Phys. 1976, 65, 450−454. (23) Brune, D.; Kim, S. Hydrodynamic Steering Effects in Protein Association. Proc. Natl. Acad. Sci. U. S. A. 1994, 91, 2930−2934. (24) Kim, S.; Karrila, S. J. Microhydrodynamics. Principles and Selected Applications; Dover Publications, Inc.: Mineola, NY, 2005; pp 189− 196. (25) Długosz, M.; Antosiewicz, J. M. Hydrodynamic Effects on the Relative Rotational Velocity of Associating Proteins. J. Phys. Chem. B 2013, 117, 6165−6174. (26) Długosz, M.; Antosiewicz, J. M.; Zieliński, P.; Trylska, J. Contributions of Far-Field Hydrodynamic Interactions to the Kinetics of Electrostatically Driven Molecular Association. J. Phys. Chem. B 2012, 116, 5437−5447. (27) Antosiewicz, J.; McCammon, J. A. Electrostatic and Hydrodynamic Orientational Steering Effects in Enzyme-Substrate Association. Biophys. J. 1995, 69, 57−65. (28) Antosiewicz, J.; Briggs, J. M.; McCammon, J. A. Orientational Steering in Enzyme-Substrate Association: Ionic Strength Dependence of Hydrodynamic Torque Effects. Eur. Biophys. J. 1996, 24, 137−141. (29) Ermak, D. L.; McCammon, J. A. Brownian Dynamics with Hydrodynamic Interactions. J. Chem. Phys. 1978, 69, 1352−1360. (30) Varga, Z.; Swan, J. Hydrodynamic Interactions Enhance Gelation in Dispersions of Colloids with Short-Ranged Attraction and Long-Ranged Repulsion. Soft Matter 2016, 12, 7670−7681. (31) Carrasco, B.; Garcia de la Torre, J. Hydrodynamic Properties of Rigid Particles: Comparison of Different Modeling and Computational Procedures. Biophys. J. 1999, 76, 3044−3057. (32) Garcia de la Torre, J. Building Hydrodynamic Bead-shell Models for Rigid Bioparticles of Arbitrary Shape. Biophys. Chem. 2001, 94, 265−274. (33) Ortega, A.; Amorós, D.; Garcia de la Torre, J. Prediction of Hydrodynamic and Other Solution Properties of Rigid Proteins from Atomic- and Residue-Level Models. Biophys. J. 2011, 101, 892−898. (34) Długosz, M. Effects of Hydrodynamic Interactions on the Apparent 1D Mobility of a Nonspecifically Bound Protein Following a Helical Path around DNA. J. Phys. Chem. B 2015, 119, 14433−14440. (35) Swan, J. W.; Wang, G. Rapid Calculations of Hydrodynamic and Transport Properties in Concentrated Solutions of Colloidal Particles and Macromolecules. Phys. Fluids 2016, 28, 011902. (36) Usabiaga, F. B.; Kallemov, B.; Delmotte, B.; Bhalla, A. P. S.; Griffith, B. E.; Donev, A. Hydrodynamics of Suspensions of Passive and Active Rigid Particles: A Rigid Multiblob Approach. Comm. Appl. Math. Comp. Sci. 2016, 11, 217−296. (37) Długosz, M.; Antosiewicz, J. M. Effects of Spatially Dependent Mobilities on the Kinetics of the Diffusion-Controlled Association Derived from the First-Passage-Time Approach. J. Phys. Chem. B 2016, 120, 7114−7127. (38) Długosz, M.; Antosiewicz, J. M. Toward an Accurate Modeling of Hydrodynamic Effects on the Translational and Rotational Dynamics of Biomolecules in Many-Body Systems. J. Phys. Chem. B 2015, 119, 8425−8439. (39) Bell, G. M.; Levine, S.; McCartney, L. N. Approximate Methods of Determining the Double-layer Free Energy of Interaction Between Two Charged Colloidal Spheres. J. Colloid Interface Sci. 1970, 33, 335− 359. (40) Gabdoulline, R. R.; Wade, R. C. Effective Charges for Molecules in Solvent. J. Phys. Chem. 1996, 100, 3868−3878. (41) Beard, D. A.; Schlick, T. Modeling Salt-Mediated Electrostatics of Macromolecules: The Discrete Surface Charge Optimization P

DOI: 10.1021/acs.jpcb.7b06053 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (64) Olsson, M. H. M. Improving the Desolvation Penalty in Empirical Protein pKa Modeling. J. Mol. Model. 2012, 18, 1097−1106. (65) Dolinsky, T. J.; Czodrowski, P.; Li, H.; Nielsen, J. E.; Jensen, J. H.; Klebe, G.; Baker, N. A. PDB2PQR: Expanding and Upgrading Automated Preparation of Biomolecular Structures for Molecular Simulations. Nucleic Acids Res. 2007, 35, W522−W525. (66) Ponder, J.; Case, D. Force Fields for Protein Simulations. Adv. Protein Chem. 2003, 66, 27−85. (67) Wang, J.; Cieplak, P.; Kollman, P. A. How Well Does a Restrained Electrostatic Potential (RESP) Model Perform in Calculating Conformational Energies of Organic and Biological Molecules? J. Comput. Chem. 2000, 21, 1049−1074. (68) Baker, N. A.; Sept, J. D.; Joseph, S.; Holst, M. J.; McCammon, J. A. Electrostatics of Nanosystems: Application to Microtubules and the Ribosome. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 10037−10041. (69) Długosz, M.; Antosiewicz, J. M. Anisotropic Diffusion Effects on the Barnase-Barstar Encounter Kinetics. J. Chem. Theory Comput. 2013, 9, 1667−1677. (70) Arora, N.; Bashford, D. Solvation Energy Density Occlusion Approximation for Evaluation of Desolvation Penalties in Biomolecular Interactions. Proteins: Struct., Funct., Genet. 2001, 43, 12−27. (71) Długosz, M.; Antosiewicz, J. M. Transient Effects of Excluded Volume Interactions on the Translational Diffusion of Hydrodynamically Anisotropic Molecules. J. Chem. Theory Comput. 2014, 10, 2583− 2590. (72) Długosz, M.; Antosiewicz, J. M. Evaluation of Proteins’ Rotational Diffusion Coefficients from Simulations of Their Free Brownian Motion in Volume-Occupied Environments. J. Chem. Theory Comput. 2014, 10, 481−491. (73) Cichocki, B.; Jones, R. B.; Kutteh, R.; Wajnryb, E. Friction and Mobility for Colloidal Spheres in Stokes Flow Near a Boundary: The Multipole Method and Applications. J. Chem. Phys. 2000, 112, 2548− 2561. (74) Rotne, J.; Prager, S. Variational Treatment of Hydrodynamic Interaction in Polymers. J. Chem. Phys. 1969, 50, 4831−4837. (75) Wajnryb, E.; Mizerski, K. A.; Zuk, P. J.; Szymczak, P. Generalization of the Rotne-Prager-Yamakawa Mobility and Shear Disturbance Tensors. J. Fluid Mech. 2013, 731, R3. (76) Cichocki, B.; Hinsen, K. Stokes Drag on Conglomerates of Spheres. Phys. Fluids 1995, 7, 285−291. (77) Harvey, S. C.; Garcia de la Torre, J. Coordinate Systems for Modeling the Hydrodynamic Resistance and Diffusion Coefficients of Irregularly Shaped Rigid Macromolecules. Macromolecules 1980, 13, 960−964. (78) Brune, D.; Kim, S. Predicting Protein Diffusion Coefficients. Proc. Natl. Acad. Sci. U. S. A. 1993, 90, 3835−3839. (79) Wegener, W. A. Transient Electric Birefringence of Dilute RigidBody Suspensions at Low Field Strengths. J. Chem. Phys. 1986, 84, 5989−6004. (80) Dhont, J. K. G. An Introductions to Dynamics of Colloids; Elsevier: Amsterdam, 1996; pp 183−186. (81) Godec, A.; Metzler, R. First Passage Time Distribution in Heterogeneity Controlled Kinetics: Going Beyond the Mean First Passage Time. Sci. Rep. 2016, 6, 20349. (82) Hänggi, P.; Borkovec, M.; Talkner, P. Reaction-Rate Theory: Fifty Years after Kramers. Rev. Mod. Phys. 1990, 62, 251−341. (83) Faas, G. C.; Mody, I. Measuring the Kinetics of Calcium Binding Proteins with Flash Photolysis. Biochim. Biophys. Acta, Gen. Subj. 2012, 1820, 1195−1204. (84) Yanagisawa, Y.; Kageyama, T.; Wada, N.; Tanaka, M.; Ohno, S. Time Courses and Time-Resolved Spectra of Firefly Bioluminescence Initiated by Two Methods of ATP Injection and Photolysis of Caged ATP. Photochem. Photobiol. 2013, 89, 1490−1496. (85) Boscá, F.; Sastre, G.; Andreu, J. M.; Jornet, D.; Tormos, R.; Miranda, M. A. Drug-Tubulin Interactions Interrogated by Transient Absorption Apectroscopy. RSC Adv. 2015, 5, 49451. (86) Petrov, V.; Stanimirov, S.; Petrov, I. K.; Fernandes, A.; de Freitas, V.; Pina, F. Emptying the β-Cyclodextrin Cavity by Light:

Photochemical Removal of the trans-Chalcone of 4′,7-Dihydroxyflavylium. J. Phys. Chem. A 2013, 117, 10692−10701. (87) Lau, A. W. C.; Lubensky, T. C. State-Dependent Diffusion: Thermodynamic Consistency and its Path Integral Formulation. Phys. Rev. E 2007, 76, 011123. (88) Zhou, H. X. Kinetics of Diffusion-Influenced Reactions Studied by Brownian Dynamics. J. Phys. Chem. 1990, 94, 8794−8800. (89) Shushin, A. I. Influence of Hydrodynamic Interaction on the Diffusion-Controlled Reaction Kinetics of Molecules with Highly Anisotropic Reactivity. J. Chem. Phys. 2003, 118, 1301−1311. (90) Szymczak, P.; Cieplak, M. Hydrodynamic Effects in Proteins. J. Phys.: Condens. Matter 2011, 23, 033102. (91) Wojciechowski, M.; Szymczak, P.; Cieplak, M. The Influence of Hydrodynamic Interactions on Protein Dynamics in Confined and Crowded Spaces-Assessment in Simple Models. Phys. Biol. 2010, 7, 046011. (92) Allison, S. Diffusion-Controlled Reactions: Hydrodynamic Interaction Between Charged, Uniformly Reactive Spherical Reactants. J. Phys. Chem. A 2006, 110, 13864−13867. (93) Frembgen-Kesner, T.; Elcock, A. H. Striking Effects of Hydrodynamic Interactions on the Simulated Diffusion and Folding of Proteins. J. Chem. Theory Comput. 2009, 5, 242−256. (94) Frembgen-Kesner, T.; Elcock, A. H. Absolute Protein-Protein Association Rate Constants from Flexible, Coarse-Grained Brownian Dynamics Simulations: The Role of Intermolecular Hydrodynamic Interactions in Barnase-Barstar Association. Biophys. J. 2010, 99, L75− L77. (95) Ando, T.; Skolnick, J. On the Importance of Hydrodynamic Interactions in Lipid Membrane Formation. Biophys. J. 2013, 104, 96− 105. (96) Ando, T.; Skolnick, J. Crowding and Hydrodynamic Interactions Likely Dominate in Vivo Macromolecular Motion. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 18457−18462. (97) Skolnick, J. Perspective: On the Importance of Hydrodynamic Interactions in the Subcellular Dynamics of Macromolecules. J. Chem. Phys. 2016, 145, 100901. (98) Mereghetti, P.; Wade, R. C. Atomic Detail Brownian Dynamics Simulations of Concentrated Protein Solutions with a Mean Field Treatment of Hydrodynamic Interactions. J. Phys. Chem. B 2012, 116, 8523−8533. (99) Kondrat, S.; Zimmermann, O.; Wiechert, W.; von Lieres, E. The Effect of Composition on Diffusion of Macromolecules in a Crowded Environment. Phys. Biol. 2015, 12, 046003. (100) Chow, E.; Skolnick, J. Effects of Confinement on Models of Intracellular Macromolecular Dynamics. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 14846−14851.

Q

DOI: 10.1021/acs.jpcb.7b06053 J. Phys. Chem. B XXXX, XXX, XXX−XXX