Effects of Spatially Dependent Mobilities on the Kinetics of the

Article pubs.acs.org/JPCB

Effects of Spatially Dependent Mobilities on the Kinetics of the Diffusion-Controlled Association Derived from the First-PassageTime Approach Maciej Długosz*,† and Jan M. Antosiewicz‡ †

Centre of New Technologies, University of Warsaw, Stefana Banacha 2c, Warsaw 02-097, Poland Department of Biophysics, Faculty of Physics, University of Warsaw, Ż wirki i Wigury 93, Warsaw 02-089, Poland

‡

ABSTRACT: Brownian dynamics (BD) simulations and the first-passage-time approach are applied to investigate diffusion-controlled association in a biologically relevant model system consisting of a fixed receptor with an elongated cavity and a capsule-like ligand that fits this cavity precisely. Before the binding at the receptor cavity, the ligand undergoes translational and rotational diffusion, either free or under the influence of electrostatic interactions with the receptor. The spatial dependence of the translational and rotational mobilities of the ligand resulting from its hydrodynamic interactions (HIs) with the receptor is accounted for in BD simulations, and an accurate numerical approach is applied for the evaluation of the spatially dependent mobility tensor of the ligand. Different magnitudes of electrostatic interactions (either attraction or repulsion) between the ligand and receptor are considered. The effective range of receptor−ligand electrostatic interactions is varied to account for their screening under different conditions of ionic strength. The effects of HIs on the kinetics of the diffusion-controlled association, evaluated for different electrostatic properties of binding partners, are thoroughly analyzed and discussed.

■

INTRODUCTION Biochemical reactions are usually preceded by a diffusive encounter step, in which the ligand and receptor (both terms may correspond to a diverse range of biomolecules) approach one another through translational and rotational diffusion. This first step usually takes place in a complex force field resulting from electromagnetic properties of biomolecules and their environments.1,2 Diffusion of molecules toward the encounter is often the rate-limiting step of the reaction with the relative mobility of binding partners controlling the kinetics of the association.3−5 Under low Reynolds number conditions that are typical for systems of biological relevance,6 the translational and rotational diffusive motions of the molecules suspended in a fluid are impacted by electromagnetic intermolecular interactions and by fluid-mediated hydrodynamic interactions (HIs) with their surroundings. Both types of interactions are of longrange. Electromagnetic interactions of molecules arise as a consequence of their electronic structure. HIs result from the fact that each moving molecule sets the solvent medium in motion and the resulting flux in the solvent medium then tends to move all of the other molecules. As a consequence, the mobility of a particular molecule is spatially dependent, as a complex function of the position and the orientation of the molecule evaluated relative to its surroundings. The general effects of HIs, their interplay with electromagnetic interactions, and their significance in biological soft matter have been of long-standing interest.7−26 In the particular case of a bimolecular association, the primary effect of HI seems to be the slowing of the diffusional © XXXX American Chemical Society

encounter step. HI reduces the relative mobility of binding partners, and as a result, the association rate constant is decreased.7−10 Friedman7 considered association of spherical molecules, either neutral or charged, and estimated, using a method originally developed by Smoluchowski,27−29 that the decrease in the association rate constant resulting from HIs between molecules is roughly 15% when compared to that from classical Smoluchowski’s theory, which neglects HI. Larger reduction was predicted by Deutch and Felderhof8 on the basis of an approximate solution of the Fokker−Planck equation. These authors predicted the reduction in the association rate constant due to HIs to be 46% for the association of hard spheres. They also predicted the reduction between 25 and 60% for association of ionic species, depending on the extent of attraction or repulsion between particles. Wolynes and Deutch10 showed that when slip boundary conditions are applied instead of stick boundary conditions (that were applied in the work of Deutch and Felderhof8), the reduction in the rate of association of neutral particles is smaller, roughly 30% relative to that from classical Smoluchowski’s theory without HIs. Whereas the aforementioned works dealt with the association of spherical, isotropically reactive objects, Shushin14 considered the association between model molecules with highly anisotropic reactivities. He predicted, on the basis of the Received: May 25, 2016 Revised: July 1, 2016

A

DOI: 10.1021/acs.jpcb.6b05281 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B method of local analysis,30 that HIs can reduce association rate constants in such systems, even about 3−5 times. Apart from the approaches described above that can be generally termed theoretical or analytical, there are also studies12,13,20,23,26 that investigate the role of HI in molecular association by means of direct computer simulations. Such studies usually employ the Brownian dynamics (BD) technique.31,32 Antosiewicz and coauthors12,13 studied diffusional encounter in a simple system consisting of a dumbbell model of an enzyme and a dumbbell model of a ligand, taking into account screened electrostatic interactions and HIs between molecules. Whereas the focus of their work was primarily on the hydrodynamic steering effects that were postulated earlier by Brune and Kim,11 these authors also concluded that association rate constants predicted using the approach of Northrup et al. 33 are overestimated by approximately 15−20% when HIs between associating molecules are neglected. Długosz et al.20 performed BD simulations of the diffusional encounter in the periodic systems of model isotropic (spheres) and anisotropic (dumbbell) molecules. The maximal magnitude of the HI effect on the association rates (i.e., the decrease in association rate constants) that the authors observed was about 30%. These authors also concluded that the relative magnitude of the HI effect does not depend on the ionic strength for the association of isotropic charged objects and for nonspecific association (i.e., without restrictions regarding the relative orientation of binding partners in an encounter complex) of anisotropic objects. However, they observed such a dependence for orientationspecific association of anisotropic objects. Frembgen-Kesner and Elcock23 investigated more elaborate systems, namely, barnase and barstar and their mutants, using coarse-grained Gõlike models of these proteins (their models consisted of 33 pseudoatoms for barnase and 27 pseudoatoms for barstar). These authors performed BD simulations either with or without barnase−barstar HI under various ionic strength conditions. They observed that simulations with HI turned on yield association rates (calculated using the approach of Northrup et al.33) that are lower up to 80% relative to those obtained from BD simulations without intermolecular HI. They also observed that effects of HI are more pronounced when electrostatic interactions between molecules are weakened. Saglam and Chong26 used the same models of proteins as those used by Frembgen-Kesner and Elcock to estimate the influence of HI on the basal rate (i.e., the rate constant for association in the absence of electrostatic interactions) of barnase−barstar association. The value of the basal rate that they obtained from BD simulations (based on the approach of Northrup et al.33) with HI turned on is roughly 2 times smaller than the value obtained from simulations with HI turned off.26 The most serious shortcoming of almost all studies, both theoretical/analytical7,8,10 and those utilizing BD simulations12,13,20,23,26 that we have briefly described above, is the employment of far-field approximations for modeling of HIs between associating molecules. The only exception is the work of Shushin14 in which long-range HIs are actually assumed to be irrelevant for the association and only their short-range behavior is considered. However, apart from Shushin, all other authors mentioned above model HI using approximate analytical formulas given by Oseen,34 Rotne−Prager,35 or Rotne−Prager−Yamakawa.36 Strictly speaking, applications of these formulas should be limited to systems consisting of wellseparated objects. The Oseen formula is applicable for widely

separated point-like particles.34 The Rotne−Prager35 and Rotne−Prager−Yamakawa36 formulas account for finite sizes of spherical particles, but even for the simplest system, consisting of two spheres, they fail to correctly describe its dynamics when the distance between the surfaces of spheres is comparable to their sizes. Moreover, predictions regarding the role and magnitude of the effects of HIs on the kinetics of diffusional encounter of molecules described above are made for systems consisting of rather crude molecular models (perhaps with the exception of barnase and barstar models constructed by Frembgen-Kesner and Elcock23) that are unable to accommodate complex molecular shapes. Thus, it appears that despite a number of attempts at evaluating the effects of HIs on the kinetics of diffusioncontrolled association of biomolecules the knowledge regarding their role and significance in such processes is still limited. Thus, we believe that the current work provides a valuable contribution to this subject. We present an in-depth analysis of hydrodynamic effects in a biologically relevant model system consisting of a fixed spherical receptor with an elongated cavity located at its surface and a capsule-like ligand that fits this cavity precisely. We performed BD simulations in which the ligand molecule undergoes translational and rotational diffusion, either free or under the influence of electrostatic interactions with the receptor, and eventually associates with the receptor. We consider different magnitudes of receptor−ligand electrostatic interactions and vary their effective range, accounting for their screening by ions dissolved in solution. The mobility of the ligand in the studied system is spatially dependent because of its HIs with the receptor. HIs in our BD simulations are evaluated using the bead−shell approach22,37−40 that is capable of accommodating arbitrarily complex molecular shapes and numerical calculations of mobility tensors for systems consisting of rigid clusters of spherical elements that result in the correct description of short- and long-range behaviors of HIs between clusters.39 Hence, we overcome the shortcomings of far-field approximations mentioned above. The kinetics of receptor−ligand encounters, either nonspecific or specific, is evaluated on the basis of BD simulations using the first-passagetime approach.41,42

■

THEORETICAL AND COMPUTATIONAL METHODS First-Passage-Time Approach to the Diffusion-Controlled Association. The formulation of the first-passage-time problem that is given below follows that presented in ref 41. We consider a ligand molecule that diffuses in a finite domain Ω. The boundary of the Ω domain that we denote δΩ is reflective. Ω encloses a subdomain ω whose boundary δω is absorbing. We identify ω with a fixed receptor. The dissociated state of the ligand corresponds to the situation in which the ligand stays in Ω\ω (i.e., between boundaries δω and δΩ). The ligand molecule, generated initially at some point in Ω\ω, diffuses (possibly under the influence of the potential resulting from its interactions with the receptor) and eventually reacts upon reaching the boundary δω (the ligand associates with the receptor). The probability of finding the ligand molecule at time t at the position Q⃗ (in general, Q⃗ can be a generalized vector that describes not only the position of the diffusing molecule but also its orientation in a given coordinate frame), P(Q⃗ ,t), obeys the Smoluchowski equation43 δtP(Q⃗ , t ) = −∇Q⃗ j ⃗ (Q⃗ , t ) B

(1) DOI: 10.1021/acs.jpcb.6b05281 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B ⃗ ⃗ ,t) being the flux with j(Q

obtained from BD simulations in which the configuration space trajectories of the ligand molecule in Ω\ω, Q⃗ (t), are composed of successive displacements ΔQ⃗ , each taken over a short time step Δt. Such a displacement is evaluated as31,49

⎛ 1 j ⃗ (Q⃗ , t ) = −kBTμ(Q⃗ )⎜∇Q⃗ P(Q⃗ , t ) + P(Q⃗ , t )∇Q⃗ kBT ⎝ ⎞ E(Q⃗ )⎟ ⎠

ΔQ⃗ = kBT ∇Q⃗ μ(Q⃗ )Δt + μ(Q⃗ )∇Q⃗ E(Q⃗ )Δt + R⃗ Q⃗ (Δt ) (2)

(10)

kB is the Boltzmann constant, T is the temperature, μ(Q⃗ ) is the spatially dependent mobility of the diffusing molecule, and E(Q⃗ ) denotes either the potential affecting the ligand or potential energy of the system in a configuration given with Q⃗ . ∇Q⃗ denotes the gradient with respect to the components of Q⃗ . We attempt to solve the above equation with the following boundary conditions

where displacement R⃗ Q⃗ (Δt) is a random displacement with a Gaussian distribution function whose average value is 0 and whose variance−covariance is given by (the † symbol denotes the transpose operation) ⟨R⃗ Q⃗ (Δt )(R⃗ Q⃗ (Δt ))† ⟩ = 2kBTμ(Q⃗ )Δt

(11)

First, a nonreactive (before the activation or triggering) system is considered, subjected to the following boundary conditions

j ⃗ (Q⃗ , t ) = 0

on

δΩ

(3)

P(Q⃗ , t ) = 0

on

δω

(4)

j ⃗ (Q⃗ , t ) = 0

on

δΩ

(12)

(5)

j ⃗ (Q⃗ , t ) = 0

on

δω

(13)

and for the δ-function initial condition P(Q⃗ , 0) = δ(Q⃗ − Q⃗ o)

Q⃗ o ∈ Ω\ω

The desired solution is the conditional probability, P(Q⃗ ,t|Q⃗ o,0), that if the ligand molecule was at Q⃗ o at time 0 then it will be at Q⃗ at time t. The probability that the ligand molecule located initially at Q⃗ o has not yet reacted at time t is denoted S(Q⃗ o,t) (S(Q⃗ o,t) is the survival probability function) and is given by S(Q⃗ o , t ) =

∫Ω\ω P(Q⃗ , t|Q⃗ o , 0) dQ⃗

that is, both boundaries of the space in which the ligand diffuses, outer δΩ and inner δω, are reflective. BD simulations are performed for such a system, and positions, Q⃗ , from BD trajectories of a molecule diffusing in Ω\ω are collected forming an equilibrium or a steady-state ensemble. This ensemble corresponds to the probability distribution, Peq(Q⃗ o), in eqs 8 and 9. We note that the divergence (or drift) term, ∇Q⃗ μ(Q⃗ )Δt, that is present in eq 10 is required for the probability distribution, Peq(Q⃗ o), resulting from BD simulations with a spatially dependent mobility to approach the Boltzmann distribution.49 As in the studied system the mobility of the ligand decreases near the acceptor (see further in the text), the neglect of the drift term would result in an artificially increased equilibrium probability of finding the ligand near the δω boundary. Next, the studied system is activated (or triggered), that is, its inner boundary δω is made to be absorbing. Multiple BD trajectories are generated. These trajectories start from the configurations of the system taken from the equilibrium ensemble (in accordance with the δ-function initial condition, eq 5). Each trajectory is propagated until the ligand molecule associates with the receptor (i.e., the ligand reaches δω). Times needed for the receptor−ligand association to take place (mean first-passage times, eq 7) are evaluated for each trajectory, and the global first-passage time is calculated as their average according to eq 9. Mobility of the Ligand Moving in the Presence of a Fixed Receptor. Below, we present the construction of the spatially dependent mobility tensor, μ(Q⃗ ), for a ligand molecule that moves in the vicinity of an unmoving receptor. A similar description can be found in refs 22, 50. Initially, we consider a general form of the mobility tensor, M, for two arbitrarily shaped objects in a particular configuration: the ligand particle (1) and receptor (2), suspended in a viscous fluid. Such a tensor is represented by a 12 × 12 matrix

(6)

As −δtS(Q⃗ o,t) dt is the probability that the ligand reaches δω in the time interval (t,t + dt), the average time τ(Q⃗ o) required for the receptor−ligand association to take place (the mean firstpassage time41,42,44) can be evaluated as τ(Q⃗ o) = −

∫0

∞

tδtS(Q⃗ o , t ) dt =

∫0

∞

S(Q⃗ o , t ) dt

(7)

We may now consider a situation in which δω is reflective for t < 0 and at t = 0 the studied system needs to be activated (or triggered) for a reaction to occur (as, e.g., in the case of experiments involving various pump−probe techniques45−48). Before the activation, the probability of finding the ligand at the position Q⃗ o ∈ Ω\ω follows the equilibrium distribution Peq(Q⃗ o). The probability that the ligand starting from the equilibrium distribution Peq(Q⃗ o) stays in Ω\ω at time t after activation is given by S(t ) = =

∫Ω\ω Peq(Q⃗ o)S(Q⃗ o , t ) dQ⃗ o ∫Ω\ω ∫Ω\ω Peq(Q⃗ o)P(Q⃗ , t|Q⃗ o , 0) dQ⃗ o dQ⃗

and the average time for association passage time44) is τ=

∫Ω\ω Peq(Q⃗ o)τ(Q⃗ o) dQ⃗ o = ∫0

41

(8)

(or the global first-

∞

S (t ) d t

(9)

where we have effectively integrated out the dependence of the mean first-passage time on the initial position of the ligand Q⃗ o. BD Simulations. The average reaction time (or the global first-passage time) for the above-described diffusion-controlled process in a force field, governed by the Smoluchowski equation with a space-dependent mobility (eqs 1 and 2), can be

⎛ μ11 μ12 ⎞ M=⎜ ⎟ ⎝ μ21 μ22 ⎠

(14)

that consists of four 6 × 6 blocks, μij, of form C


Article

The Journal of Physical Chemistry B ⎛ μ TT μ TR ⎞ ij ⎟ ⎜ ij μij = ⎜ ⎟ RT ⎜μ μijRR ⎟⎠ ⎝ ij

(15)

where indices i ∈ (1,2) and j ∈ (1,2) correspond to molecules. Subtensors μαβ (α ∈ (T,R), β ∈ (T,R)) correspond to ii translations (TT) and rotations (RR) of the ith molecule and their couplings (TR/RT). Subtensors μαβ i≠j (α ∈ (T,R), β ∈ (T,R)) describe how the motions of one molecule (i) affect the motions of the second molecule (j). Thus, we have couplings between translations of molecules (TT), rotations of molecules (RR), and their translations and rotations (TR/RT). The motions of each molecule are evaluated in a coordinate system whose origin coincides with the geometric center of that molecule. The corresponding axes of these molecule-centered coordinate systems are parallel to each other. We may now introduce the translational and rotational velocities of molecules, represented by generalized vector v ⃗ = (v†1⃗ ,v†2⃗ ), and forces and torques acting on molecules, represented by generalized vector F⃗ = (F⃗†1,F⃗†2). A linear relation, applicable at a low Reynolds number, holds v ⃗ = MF ⃗

(16)

v1⃗ = μ11 F1⃗ + μ12 F2⃗

(17)

v2⃗ = μ21 F1⃗ + μ22 F2⃗

(18)

Figure 1. Bead−shell models of the receptor (orange) and the capsulelike ligand (blue). Drawings were done with the Visual Molecular Dynamics package.51

number of spherical elements distributed on its surface is maximal, whereas at the same time, there are no overlaps between elements. Overlapping spherical elements belonging to surfaces of neighboring primary beads are removed. Moreover, spherical elements that are not accessible from the outside of the capsule-like ligand are also removed. The total number of spherical elements in the final hydrodynamic model of the ligand is 1442 (Figure 1). To create the bead−shell model of the receptor, first, we cover the surface of a single primary bead with the radius (Σ) equal to 15 Å with a maximal number of equidistributed, nonoverlapping spherical elements of radii 0.24 Å. Next, we prepare a bead−shell model of the ligand particle in which the radii of primary beads are inflated by 1 Å (the arrangement of primary beads of the ligand remains unchanged); the radii of spherical subunits covering the surface of the inflated ligand model again are set to 0.24 Å. We position this bead−shell model of the inflated ligand so that its center of geometry lies on the surface of the receptor primary bead and its long axis is perpendicular to the line connecting the center of the receptor bead and the geometric center of the ligand. From the surface of the receptor primary bead, we remove all spherical subunits that are enclosed by the surface of the ligand, and at the same time, from surfaces of inflated ligand beads, we remove all spherical subunits located above the surface of the receptor. As a result, we obtain a bead−shell model of the receptor with an elongated cavity able to accommodate the bead−shell model of the original capsule-like ligand (Figure 1). The total number of spherical subunits in this receptor model is 8827. The details of the procedure for evaluating the mobility tensor, M, for a system consisting of an arbitrary number of bead−shell models are given elsewhere;21,39 thus, we describe it here only briefly. First, for the receptor−ligand system in a particular configuration, the general mobility tensor for N spherical subunits is evaluated, with N being the sum of the number of spherical subunits in the bead−shell model of the capsule-like ligand and the number of spherical subunits in the bead−shell model of the receptor. At this stage, these spherical subunits are treated as independent entities. The elements of this tensor are evaluated at the two-body level using the analytical Rotne−Prager approach35,39 (translational and rotational self-mobilities and couplings between translations and

or

In the considered system, the ligand molecule (1) moves freely, whereas the position and orientation of the receptor (2) are fixed by some external forces and torques. As a consequence, v2⃗ vanishes v2⃗ = μ21 F1⃗ + μ22 F2⃗ = 0

(19)

and we may eliminate F⃗2 from eqs 17 and 18 and write v1⃗ = (μ11 − μ12 (μ22 )−1μ21 )F1⃗ = μF1⃗

(20)

where we have introduced the 6 × 6 mobility tensor, μ = μ11 − μ12(μ22)−1μ21, of the ligand that moves in the presence of the fixed receptor. Generally, μ is represented by a matrix composed of four 3 × 3 blocks that vary with the configuration of the system (i.e., the position and orientation of the ligand relative to the fixed receptor), described with the vector Q⃗ , due to HIs between the ligand and receptor ⎛ μTT μTR ⎞ ⎛ μTT (Q⃗ ) μTR (Q⃗ )⎞ ⎟ = μ(Q⃗ ) ⎟≡⎜ μ = ⎜⎜ ⎟ ⎜ RT ⃗ RT RR ⎟ RR ⃗ μ ⎠ ⎝ μ (Q ) μ (Q ) ⎠ ⎝μ

(21)

Calculation of Mobility Tensors for the Receptor− Ligand System. Both molecules, ligand and receptor, are represented in hydrodynamic calculations of mobility tensors using bead−shell models22,37−40 (Figure 1). To create the bead−shell model of the capsule-like ligand, we arrange three spheres (we will refer to these spheres using the term primary beads or simply beads) of radii 5 Å so that their centers are collinear. The distances between the centers of neighboring primary beads are equal to 3 Å. The surfaces of primary beads are covered with equidistributed spherical elements with equal radii of 0.24 Å. The positions of spherical elements on the surfaces of primary beads are calculated using an algorithm described elsewhere.39 For each primary bead, the D


Article

The Journal of Physical Chemistry B

(Figure 2). This line coincides with the x axis of the laboratory coordinate system. Rotational diffusion is also one-dimensional, the ligand molecule rotates around the x axis of the laboratory coordinate system, and this motion corresponds to changes in the angle, ϕ, measured in the yz plane between the long axis of the ligand (denoted z′ in Figure 2) and the z axis of the laboratory coordinate system (angle ϕ samples values in the range between 0 and 360°). The position and orientation of the ligand in the laboratory coordinate frame are thus described with the vector Q⃗ = (x,ϕ). As described above, the mobility tensor of the ligand depends on Q⃗ (both on x and ϕ) as a result of HIs of the ligand with the receptor (we should also note here that regardless of HIs, the components of mobility tensors depend on the orientation of the ligand in the coordinate frame in which they are evaluated54−56). This tensor, evaluated in a coordinate system with the origin located at the geometric center of the ligand (which is also the diffusion center of the ligand54) and the axes that are always parallel to the axes of the laboratory coordinate system, x, y, z (Figure 2), can be represented by a 6 × 6 matrix that for the system studied here has a general form

translations, rotations and rotations, and translations and rotations of different spherical elements are evaluated using expressions in the form of power series of the inverse of the interelement distances, r, up to the order of 1/r3). With the sufficient number of spherical subunits in the hydrodynamic models of interacting bodies, this level of theory is able to provide a description of HIs between rigid bodies composed of these subunits that is in a reasonable quantitative agreement39 with the description resulting from exact numerical schemes such as the induced force method.52 The resulting general mobility tensor matrix is of size 6N × 6N. Having calculated the 6N × 6N general mobility tensor, we apply the projection operation39,53 to obtain from this tensor the mobility tensor M that is given with a 12 × 12 matrix (eq 14), which describes a system consisting of two rigid conglomerates of spherical subunits (i.e., molecules).21,39 Once M is known, the 6 × 6 mobility tensor μ (eq 21) for a ligand moving in the presence of the fixed receptor can be evaluated as described in the previous section. A similar procedure is applied to evaluate the mobility tensor for the bead−shell model of the ligand in the absence of the receptor. Hydrodynamic bead−shell models of molecules were prepared using an in-house software. All hydrodynamic calculations were conducted at the temperature of 293.15 K with the solvent viscosity set to 0.102 poise using an in-house software. Diffusion-Controlled Association in the Model Receptor−Ligand System. Our model system consists of a spherical receptor with an elongated cavity located at its surface and an axially symmetric, capsule-like ligand (Figures 1 and 2). The ligand molecule undergoes translational and rotational Brownian motions, whereas the position and orientation of the receptor are fixed. The translational diffusion of the geometric center of the ligand is limited to a line passing through the center of the receptor and the center of the receptor cavity

⎛ μ xx ⎜ ⎜0 ⎜ ⎜0 μ=⎜ ⎜ μϕx ⎜ ⎜0 ⎜ ⎜ ⎝0

0

0

μ xϕ 0

μ yy μ yz 0

μ yβ

μ zy μ zz 0

μ zβ

0

0

μϕϕ 0

μ βy μ βz 0

μ ββ

μγy μγz 0

μγβ

0 ⎞ ⎟ μ yγ ⎟ ⎟ zγ ⎟ μ ⎟ 0 ⎟ ⎟ βγ ⎟ μ ⎟ ⎟ μγγ ⎠

(22)

Components μ , i ∈ (x,y,z), correspond to translations in directions x, y, z. Components μii, i ∈ (α,β,γ), correspond to rotations around the axes that are parallel to directions x, y, z and pass through the geometric center of the ligand. Components μij, i ≠ j i,j ∈ (x,y,z,α,β,γ), correspond to translation−translation, translation−rotation, and rotation− rotation couplings. As described above, in our model system, we consider only the translational motions of the ligand in the x direction and its rotations around the x axis. In the considered system, these motions are uncoupled from other motions, that is, translations along y and z and rotations about angles β and γ (eq 22). Thus, the μ matrix that we use here simply reduces to ii

⎛ μ xx μ xϕ ⎞ ⎟ μ = ⎜⎜ ⎟ ϕx μϕϕ ⎠ ⎝μ

(23)

Additionally, we define mobility matrix μo ⎛ μ xx μ xϕ ⎞ o o ⎟ μo = ⎜⎜ ϕx ϕϕ ⎟ ⎝ μo μo ⎠

Figure 2. Scheme of the receptor−ligand system considered in this work. The ligand undergoes translational (translations along the x axis of the laboratory coordinate system) and rotational (rotations around the x axis of the laboratory coordinate system described with angle ϕ) diffusion in the presence of the fixed receptor. Ligand diffusion occurs in a finite domain with planar boundaries located at x = 0 (absorbing boundary) and x = L (reflective boundary). Nonspecific association takes place when the center of the ligand crosses the boundary located at x = 0. For a specific association to take place, the ligand must cross the boundary located at x = 0, oriented parallel to the receptor’s cavity (i.e., the long axis of the ligand, z′, must be either parallel or antiparallel to the z axis of the laboratory coordinate frame). Σ denotes the radius of the receptor, and σT is the translational hydrodynamic radius of the ligand.

(24)

for the capsule-like ligand in the absence of the receptor (i.e., evaluated neglecting receptor−ligand HIs). Similarly, as for the μ tensor, μo is evaluated in the coordinate system with an origin located at the geometric center of the ligand and axes that are always parallel to the axes of the laboratory coordinate system. Unlike the components of the μ tensor, the components of μo do not depend on the position (x) of the ligand in the laboratory coordinate frame; however, they vary with the varying orientation of the ligand in the laboratory coordinate E


Article

The Journal of Physical Chemistry B frame. This results from general properties of mobility tensors already mentioned above.54−56 We define the hydrodynamic translational radius of the ligand, σT, as σ T = (6πημoxx )−1 ,

Table 1. Energies (E) of the Receptor−Ligand Complex, Calculated on the Basis of Equation 27 (with x Set to 0; Figure 2) at Different Ionic Strengths (SDH ) and for Different Magnitudes of Electrostatic Interactions (α)a

(25)

E(x = 0)/kBT

where η is the viscosity of the solvent and the symbol ∥ denotes that μxx o is evaluated for the parallel orientation of the ligand in the laboratory coordinate frame (i.e., the angle, ϕ, between the z axis of the laboratory frame and the long axis, z′, of the ligand equals 0; Figure 2). In the case of the capsule-like ligand considered in this work, σT equals roughly 6.6 Å. Similarly, the rotational hydrodynamic radius of the ligand, σR, can be defined as

σ R = (8πημoϕϕ )−1/3 ,

a

(26)

In the studied case, σR equals roughly 6.8 Å. The ligand is allowed to diffuse in the finite domain between planar boundaries located at the origin of the laboratory coordinate system x = 0 (the receptor is located at x = −(Σ + σT), with Σ = 15 Å being the radius of the receptor sphere) and x = L (Figure 2). The plane located at x = 0 is absorbing, whereas the plane located at x = L is reflective. The value of L, 150 Å, is so chosen that (for ϕ ∈ ⟨0°,360°)) μ ≃μo when x = L, that is, at x = L, the ligand almost does not feel the flow disturbance resulting from the presence of the receptor. Consequently, with the L/σT ratio equal to roughly 22.6, the studied system is rather dilute. In the case of the nonspecific association, the ligand is absorbed when its center crosses the x = 0 plane, regardless of the orientation of the ligand in the laboratory coordinate frame. In the case of the specific association, the angle ϕ between the long axis of the ligand and the z axis of the laboratory coordinate frame must be in the ⟨350°,10°⟩ range, for the ligand to be absorbed on the x = 0 plane. However, we note that the position of the absorbing plane is so chosen that all orientations of the ligand are sterically allowed on this plane; in our simulations, the ligand actually never enters the cavity of the receptor. Electrostatic interactions between the ligand and receptor, screened by dissolved ions, are modeled using the following form of the potential energy function57 E(Q⃗ ) ≡ E(x) =

⎛ σT + α exp ⎜− SDH ε( σ T + x ) ⎝

x⎞ ⎟ ⎠

SDH /σ

α/e2 = ±1

α/e2 = ±4

0.50 0.75 1.00 1.50 2.00 3.00 ∞

±0.14 ±0.28 ±0.39 ±0.55 ±0.65 ±0.77 ±1.07

±0.58 ±1.13 ±1.58 ±2.21 ±2.61 ±3.08 ±4.30

T

e denotes the unit charge.

∈ ⟨0°,360°), corresponding to different positions and orientations of the ligand in the laboratory coordinate frame. Resolution of grids in the 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 ligand in the laboratory coordinate frame, Q⃗ = (x,ϕ), its mobility tensor is evaluated on the basis of tensors precalculated on the grid by means of a bilinear interpolation. Reflective boundaries are implemented in BD simulations as follows.58 From a particular configuration of the system (Q⃗ ), a BD step is made, that is, translational and rotational displacements of the ligand are calculated according to eq 10, and its current position and orientation are changed accordingly. If, as a result of the BD step, the ligand crosses the reflective boundary, then the system is brought back to its initial configuration (Q⃗ ) and a BD step is attempted with different random vectors R⃗ Q⃗ (eq 10) until either the new position of the ligand falls within the finite domain or the ligand moves away from the boundary. A similar procedure is applied in the case of the specific association, when the boundary positioned at x = 0 becomes reflective for particular orientations of the ligand (i.e., for ϕ ∈ (10°,350°)). BD simulations were conducted at a temperature of 293.15 K. The integration step (Δt in eq 10) was set to a rather conservative value of 2.5 ps. BD simulations were conducted for different values of parameters α and SDH (eq 27). For a given pair of parameters α, SDH , we first generate an ensemble of equilibrium configurations of the ligand/receptor system. This ensemble consists of 2 × 106 configurations collected from BD trajectories of a nonreactive system (in these simulations, both boundaries of the system are reflective). For a given pair (α, SDH ), a total number of 5 × 105 BD trajectories of a reactive (i.e., activated) system are then generated starting from different configurations, chosen randomly from the equilibrium ensemble. Each trajectory is propagated until the ligand associates with the receptor. Durations of these trajectories are collected for analysis. BD simulations with and without HIs between the ligand and receptor are conducted using the same protocol. All BD simulations were performed using an in-house software. Further in the text, we use the translational hydrodynamic radius of the ligand molecule, σT (eq 25), as the length unit. As

(27)

where ε is the dielectric constant of the solvent (set here to 80), SDH denotes the Debye−Hückel length,57 σT is the translational hydrodynamic radius of the ligand, and α is the magnitude of receptor−ligand electrostatic interactions. The energies of the receptor−ligand complex calculated for the values of parameters α and SDH considered in this work are given in Table 1. Diffusion of the ligand in the finite domain limited by planes x = 0 and L is simulated according to the BD algorithm described above (eq 10). Simulations are conducted either with or without HIs between the receptor and ligand. In the former case, the μ tensor given in eq 23 is used in the evaluation of ligand translational and rotational motions. In the latter case, the μ tensor is replaced with the μo tensor given in eq 24. To avoid performing hydrodynamic calculations on-the-fly during BD simulations, we precalculate mobility tensors μ and μo on two-dimensional curvilinear grids (x,ϕ), with x ∈ ⟨0,L⟩ and ϕ

the time unit, we use the quantity

3πησ TL2 , kBT

that is, time, for

which an average square displacement of a particle with the translational mobility coefficient 1 T that undergoes Brownian 6πησ

F


Article


changes are much smaller, as at x = 0, μϕϕ drops only by approximately 10% when compared with the value calculated at x = L. Thus, the rotational dynamics of the ligand is less affected by HIs with the receptor than its translational dynamics in the considered finite domain. Moreover, the effects of HI on the translational mobility of the ligand are of much 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 x = 0 and for distances of the ligand from the absorbing boundary that are above σT, the coupling between translations and rotations of the ligand molecule practically vanishes (Figure 3). Contribution of HI to the Slowing Down of the Receptor−Ligand Association. As we described above, the primary effect of HIs is the slowing down of the first, diffusive step of the association.7−10 This effect is clearly visible in Figure 4, depicting survival probability functions (S(t), eq 8) for the

motion in one dimension in a fluid of viscosity η at the temperature T, equals L2. Unless stated otherwise, all results presented below are obtained for L equal to 150 Å (22.6 σT).

■

RESULTS Spatially Dependent Ligand Mobility. To accentuate the impact of HI on the translational and rotational dynamics of the ligand, in Figure 3 we show the components μxx, μϕϕ, and μxϕ

Figure 4. Exemplary survival probability functions (S(t)) for nonspecific association reactions, obtained on the basis of BD simulations performed either with or without HIs between the receptor and ligand. Survival probability functions corresponding to different magnitudes of electrostatic interactions (α) between the receptor and ligand and zero ionic strength (infinite value of the Debye−Hückel length SDH ) are shown. A base-10 logarithmic scale is used for the x axis in the plot.

Figure 3. Components of the mobility tensor of the 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). For the considered values of the x coordinate, translational and rotational mobilities practically do not depend on the orientation of the ligand in the laboratory coordinate system and thus their values are shown only for ϕ = 0°.

nonspecific receptor−ligand association. The curves presented in Figure 4 are obtained on the basis of BD simulations conducted either with or without HIs. At this point, we consider the results of BD simulations in which electrostatic interactions between the ligand and receptor are either switched off (by setting α in eq 27 to 0), maximally repulsive (α is set to 4e2 and there is no screening by dissolved ions, that is, SDH → ∞), or maximally attractive (α is set to −4e2 and there is no screening by dissolved ions). As described above, the area under a given S(t) curve corresponds to the value of the global first-passage time for the corresponding association process (eq 9). Clearly, the inclusion of HIs results in larger areas under survival probability curves and longer global first-passage times and thus slower association kinetics. Interestingly, whereas survival probability curves obtained for α = 0 and 4e2 and with HI either switched on or off can be reasonably approximated by single-exponential decays, the same is not true for the case α = −4e2 for which the behavior of S(t) is more complicated. Because of a strong electrostatic

of the mobility tensor μ, evaluated as described above, for the ligand in the presence of the fixed receptor (eq 23), as functions of the distance of the ligand from the absorbing boundary, and for different orientations of the ligand in the laboratory coordinate frame (angle ϕ, Figure 2). Because of the symmetry of the studied system and the choice of the rotation axis, the two components of μ that describe coupling between translations and rotations of the ligand molecule are equal, that is, μxϕ = μϕx. As a result of HIs, all components of the μ tensor are monotonic functions of the distance between the ligand and fixed receptor. The largest variation with the distance from the absorbing boundary is observed in the case of translational mobility, μxx, which at x = 0 drops to approximately 25% of the value calculated at x = L (as we have described above, L ≃ 22.6 σT is so chosen that μ converges to μo when x approaches L; μo is the mobility tensor for the ligand evaluated neglecting its HIs with the receptor). In the case of the rotational mobility, G


Article


Figure 5. Global first-passage times (τ) for the nonspecific association reaction as functions of the magnitude of receptor−ligand electrostatic interactions (α) at different values of the ionic strength (SDH ). Results of BD simulations (A) with and (B) without HIs. Error bars, evaluated as standard deviations of means, are comparable to symbol sizes. A base-10 logarithmic scale is used for the y axis in both plots.

Figure 6. Global first-passage times (τ) for the specific association reaction as functions of the magnitude of receptor−ligand electrostatic interactions (the magnitude of α) at different values of the ionic strength (SDH ). Results of BD simulations (A) with HIs and (B) without HIs. Error bars, evaluated as standard deviations of means, are comparable to symbol sizes. A base-10 logarithmic scale is used for the y axis in both plots.

the random Brownian forces. Thus, in the α = −4e2 case, only a relatively small fraction of BD trajectories, initiated using the equilibrium ensemble, samples far regions of the diffusion space in the vicinity of the reflective, outer boundary of the system. As a result, the association occurs predominantly on a short-tomoderate time scale (Figure 4), without the ligand feeling the presence of the boundary; thus, the resulting survival probability function deviates significantly from the singleexponential behavior. The absolute values of global first-passage times calculated on the basis of BD simulations conducted for different combinations of parameters α and SDH (magnitude of

attraction between the receptor and ligand, the diffusional dynamics of the ligand observed when α is set to −4 is much less affected by the finite size of the diffusion space than in cases when α = 0 and 4e2. First, in the α = −4e2 case, the equilibrium probability of finding the ligand increases toward the absorbing boundary, whereas for α = 0, the equilibrium probability of finding the ligand at a particular position is constant in the whole diffusion space, and for α = 4e2, the equilibrium probability of finding the ligand at a particular position increases with an increasing distance from the absorbing boundary. Second, the smaller the receptor−ligand distance, the greater the dominance of attractive electrostatic forces over H


Article


Figure 7. Effects of the spatial dependence of the ligand mobility, resulting from its HIs with the receptor, on the global first-passage time (τ) for nonspecific (A) and specific (B) association reactions as functions of the magnitude of receptor−ligand electrostatic interactions (α) and ionic strength of the solution (SDH ). Error bars, evaluated using the uncertainty propagation method, are comparable to symbol sizes. The values of 1 − τHI/τnoHI obtained for α/e2 = 0: −0.402 ± 0.003, nonspecific association; −0.769 ± 0.005, specific association. The values of 1 − τHI/τnoHI obtained for SDH /σT → ∞: −0.187 ± 0.005 (α/e2 = −4), −0.301 ± 0.005 (α/e2 = −1), −0.556 ± 0.006 (α/e2 = 1), and −1.508 ± 0.018 (α/e2 = 4) for nonspecific association; −0.267 ± 0.004 (α/e2 = −4), −0.548 ± 0.006 (α/e2 = −1), and −1.006 ± 0.011 (α/e2 = 1) for specific association.

Let us first consider the effects of the spatial dependence of the ligand mobility, resulting from its HIs with the receptor, in the absence of electrostatic interactions between the receptor and ligand. In the case of nonspecific association, the HI contribution is roughly −40% (i.e., HI increases the global firstpassage time by 40% when compared with the global firstpassage time obtained in the absence of HI). In the case of specific association, the HI contribution is almost two times larger, roughly −80% (Figure 7, caption). The contribution of HIs to the association kinetics depends on the magnitude and effective range of electrostatic interactions between binding partners. Let us consider a particular value of the ionic strength and examine the dependence of the HI contribution on α (i.e., the magnitude of electrostatic interactions). From Figure 7, one may conclude that the largest HI effects are observed in the case of the strongest repulsion between molecules (for α = 4e2), whereas the smallest HI effects are observed in the case of the strongest attraction between molecules (for α = −4e2). Overall, the absolute value of the contribution of HI to the kinetics of association systematically increases upon changing α from −4e2 to 4e2. This holds true for all values of the ionic strength considered here and for both nonspecific and specific associations (Figure 7). Additionally, for a particular value of α, we observe the dependence of the HI contribution on the ionic strength (Figure 7). For negative values of α (attractive electrostatic interactions), the absolute value of the HI contribution decreases with a decreasing ionic strength. For positive values of α (repulsive electrostatic interactions), the absolute value of the HI contribution increases with a decreasing ionic strength. Upon varying the ionic strength, largest changes in the HI contribution are seen for the strongest repulsion, that is, for α = 4e2. The above holds true for both nonspecific and specific

receptor−ligand electrostatic interactions and the ionicstrength-dependent Debye−Hückel length, eq 27) are shown in Figure 5 (nonspecific association) and Figure 6 (specific association). We were not able to obtain trustworthy values of global first-passage times for the specific association in the system with HIs switched on for α = 4e2 and either SDH = 3σT or SDH → ∞ because of unreasonable long durations of trajectories that need to be simulated for the association to be observed. For the same reason, we were not able to obtain the global first-passage time for the specific association in the system with HIs switched off for α = 4e2 and SDH → ∞. Inclusion of HIs result in longer association times for all combinations of parameters α and SDH considered and for both nonspecific and specific reactions (Figures 5 and 6). As expected, in systems in which electrostatic interactions are repulsive, the values of the global first-passage time increase with an increasing absolute value of α. Additionally, for a particular value of α, association becomes faster with an increasing ionic strength. The reverse behavior is observed in systems in which electrostatic interactions are attractive. Namely, the values of the global first-passage time decrease with an increasing absolute value of α. Additionally, for a particular value of α, the association becomes faster with a decreasing ionic strength. This holds true for both nonspecific and specific reactions (Figures 5 and 6). Furthermore, we examine the relative contributions of HIs to the association kinetics. This can be done using the ratio τ HI

1 − noHI (with indices HI and noHI denoting global firstτ passage times resulting from BD simulations with and without receptor−ligand HI taken into account, respectively). Its values are given in Figure 7 as functions of the ionic strength for a particular magnitude of receptor−ligand electrostatic interactions, for both nonspecific and specific associations. I


Article


Table 2. Effects of Electrostatic Interactions on the Global First-Passage Time (τ) for the Nonspecific Association Reaction in the Receptor−Ligand System either with (Upper Index HI) or without (Upper Index noHI) HIsa SDH /σT

α/e2 = −4

α/e2 = −1

α/e2 = 1

1 − τ (α) /τ (α = 0) 0.023 ± 0.003 0.049 ± 0.003 0.073 ± 0.003 0.114 ± 0.002 0.143 ± 0.002 0.182 ± 0.002 0.284 ± 0.002 1 − τnoHI(α) /τnoHI(α = 0) 0.010 ± 0.002 0.026 ± 0.002 0.042 ± 0.002 0.071 ± 0.002 0.094 ± 0.002 0.128 ± 0.002 0.228 ± 0.002 HI

a

0.50 0.75 1.00 1.50 2.00 3.00 ∞

0.077 0.162 0.236 0.353 0.438 0.547 0.769

± ± ± ± ± ± ±

0.002 0.002 0.002 0.002 0.002 0.001 0.001

0.50 0.75 1.00 1.50 2.00 3.00 ∞

0.040 0.096 0.155 0.263 0.350 0.469 0.726

± ± ± ± ± ± ±

0.002 0.002 0.002 0.002 0.002 0.001 0.001

α/e2 = 4

HI

−0.022 −0.055 −0.090 −0.151 −0.198 −0.267 −0.454

± ± ± ± ± ± ±

0.002 0.003 0.003 0.003 0.003 0.003 0.004

−0.105 −0.295 −0.542 −1.112 −1.705 −2.744 −6.748

± ± ± ± ± ± ±

0.003 0.003 0.004 0.006 0.008 0.012 0.037

−0.011 −0.028 −0.047 −0.083 −0.114 −0.161 −0.310

± ± ± ± ± ± ±

0.002 0.002 0.002 0.003 0.003 0.003 0.003

−0.049 −0.138 −0.255 −0.532 −0.815 −1.316 −3.333

± ± ± ± ± ± ±

0.002 0.003 0.003 0.003 0.004 0.005 0.009

Errors were evaluated using the uncertainty propagation method.

Table 3. Effects of Electrostatic Interactions on the Global First-Passage Time (τ) for the Specific Association Reaction in the Receptor−Ligand System Either with (Upper Index HI) or without (Upper Index noHI) HIsa SDH /σT

α/e2 = −4

α/e2 = −1

α/e2 = 1

1 − τ (α) /τ (α = 0) 0.077 ± 0.003 0.142 ± 0.003 0.198 ± 0.003 0.262 ± 0.003 0.306 ± 0.002 0.352 ± 0.002 0.445 ± 0.002 1 − τnoHI(α) /τnoHI(α = 0) 0.047 ± 0.002 0.096 ± 0.002 0.137 ± 0.002 0.192 ± 0.002 0.227 ± 0.002 0.271 ± 0.002 0.367 ± 0.001 HI

a

0.50 0.75 1.00 1.50 2.00 3.00 ∞

0.247 0.413 0.516 0.631 0.691 0.757 0.869

± ± ± ± ± ± ±

0.003 0.002 0.001 0.001 0.001 0.001 0.001

0.50 0.75 1.00 1.50 2.00 3.00 ∞

0.164 0.298 0.390 0.506 0.578 0.661 0.817

± ± ± ± ± ± ±

0.002 0.002 0.001 0.001 0.001 0.001 0.001

α/e2 = 4

HI

−0.076 −0.170 −0.261 −0.405 −0.503 −0.645 −0.946

± ± ± ± ± ± ±

0.004 0.004 0.004 0.005 0.005 0.006 0.008

−0.381 −1.032 −1.875 −3.840 −5.645

± ± ± ± ±

0.005 0.008 0.012 0.025 0.038

−0.052 −0.114 −0.176 −0.278 −0.357 −0.461 −0.717

± ± ± ± ± ± ±

0.002 0.002 0.002 0.003 0.003 0.003 0.004

−0.245 −0.651 −1.187 −2.426 −3.688 −5.887

± ± ± ± ± ±

0.003 0.003 0.005 0.008 0.011 0.018

Errors were evaluated using the uncertainty propagation method.

Effects of Electrostatic Interactions on the Association Kinetics. Another issue that we investigate is whether the magnitude of the effects of electrostatic interactions on the association kinetics depends on HIs. For this, we evaluate the

associations (Figure 7). For the strongest attraction that we consider here with α set to −4e2 and the nonspecific association, the HI contribution changes monotonically from roughly −40% to roughly −20% upon going from the practically infinite ionic strength (no electrostatic interactions) to the zero ionic strength (for which SDH → ∞). In the case of specific association, the HI contribution changes monotonically from roughly −80% (neutral system or an infinite ionic strength) to roughly −30% (SDH → ∞). For the strongest repulsion that we consider here with α set to 4e2 and the nonspecific association, the HI contribution changes monotonically from roughly −40% to roughly −150% upon going from the practically infinite ionic strength (no electrostatic interactions) to the zero ionic strength. In the case of specific association, the HI contribution changes monotonically from roughly −80% (neutral system or an infinite ionic strength) to roughly −150% (at SDH = 2σT); as described above, we were not able to determine specific association times for SDH = 3σT and SDH → ∞ in the case of α = 4e2.

ratio 1 −

τ(α) , τ(α = 0)

which corresponds to the contribution of

electrostatic interactions to the association kinetics (τ(α) denotes the global first-passage time obtained for a particular (nonzero) value of α, and τ(α = 0) denotes the global firstpassage time obtained for a neutral system) for both nonspecific and specific reactions, simulated either with or without HIs between molecules, and for different ionic strengths. The resulting values of the electrostatic contribution are given in Tables 2 and 3. As expected, there is no qualitative difference between systems simulated with and without HIs regarding the influence of electrostatic interactions on the association kinetics. Electrostatic interactions either slow down (repulsion) or speed up (attraction) the association, and the absolute values of J


Article


Figure 8. Comparison of global first-passage times (τ) for specific (upper index S) and nonspecific (upper index NS) association reactions for different values of ionic strength (SDH ) and different magnitudes of receptor−ligand electrostatic interactions (α). (A) Results of simulations with HIs. (B) Results of simulations without HIs. Error bars, evaluated using the uncertainty propagation method, are comparable to symbol sizes. Values of 1 − τS/τNS obtained for α/e2 = 0: −1.447 ± 0.008, simulations with HIs; −0.940 ± 0.002, simulations without HIs. Values of 1 − τS/τNS obtained for SDH /σT → ∞: −0.383 ± 0.003 (α/e2 = −4), −0.893 ± 0.004 (α/e2 = −1), and −2.276 ± 0.011 (α/e2 = 1) for simulations with HIs and −0.295 ± 0.003 (α/e2 = −4), −0.592 ± 0.003 (α/e2 = −1), and −1.541 ± 0.004 (α/e2 = 1) for simulations without HIs.

kinetics decreases with a decreasing ionic strength; however, HI effects are rather considerable even at zero ionic strength. In systems with repulsive electrostatic interactions, the impact of HI on the electrostatic contribution to the association kinetics for a particular value of the ionic strength increases with an increasing absolute value of α as can be deduced from data in Tables 2 and 3 (the opposite is observed in systems in which molecules attract each other). For example, in the case of nonspecific association at zero ionic strength, the ratio of electrostatic contributions obtained from simulations with and without HIs is approximately 1.5 for α = 1e2 and approximately 2 for α = 4e2 (Table 2). Specific and Nonspecific Association. Regardless of the presence or absence of HI in studied systems, the kinetics of the specific association is much slower than the kinetics of the nonspecific association (Figures 5 and 6). One may ask whether the relative slowing down of the association that results from requesting that molecules need to be oriented in a particular way for the reaction to occur is affected by HI. We define the

magnitudes of these effects increase monotonically with a decreasing ionic strength; this is true for both nonspecific and specific associations (Tables 2 and 3). However, results obtained for systems with and without HIs are quantitatively different. Overall, the absolute value of the electrostatic contribution to the association kinetics is larger in systems in which HI between molecules is taken into account. Moreover, this apparent increase of electrostatic effects by HIs is overall larger for nonspecific reactions than for specific ones (Tables 2 and 3). One should keep in mind that there are no electrostatic torques in the studied systems. It is evident from Table 2 (nonspecific association) and Table 3 (specific association) that when the receptor and ligand attract each other (α = −4e2 or −1e2) the impact of HI on the electrostatic contribution to the association kinetics for a particular value of α decreases with a decreasing ionic strength. As an illustration, let us consider the following. The electrostatic contribution to the nonspecific association kinetics for α/e2 = −4e2 and SDH /σT = 0.5 is roughly 8% when HIs between molecules are accounted for and roughly 4% when we neglect HI (Table 2); the ratio of these values is 2. For zero ionic strength (when SDH → ∞), the electrostatic contribution is roughly 77% when HIs between molecules are accounted for and roughly 73% when we neglect HI; the ratio of these two values is slightly above 1. HI effects diminish at zero ionic strength when electrostatic interactions are the strongest because of the absence of screening by dissolved ions. Moreover, the impact of HI on the electrostatic contribution to the association kinetics for a particular value of the ionic strength decreases with a increasing absolute value of α. For systems in which interactions between molecules are repulsive, we also observe that for a particular value of α the impact of HI on the electrostatic contribution to the association

τS

ratio 1 − NS , where τS denotes the global first-passage time for τ the specific association and τNS denotes the global first-passage time for the nonspecific association, and evaluate its value for systems either with or without HIs between the ligand and receptor and for different values of parameters α and SDH . The results are given in Figure 8. We observe that regardless of the presence or absence of HI the considered ratio depends on the magnitude of electrostatic interactions in the system (i.e., α) and on the ionic strength. Overall, smallest absolute values of τS

1 − NS are observed for α = −4e2 (attraction), and largest τ absolute values are observed for α = 4e2 (repulsion). The absolute value of the considered ratio decreases with a decreasing ionic strength when the receptor and ligand attract K


Article


molecules. In the current work, we apply an approach that is able to accurately describe both short- and long-range HI effects. Thus, we overcome the most serious shortcoming of the previous studies on the effects of HI on the diffusion-controlled association.7−10,12−14,20,23,26 Although in our work we evaluate the association kinetics in terms of global first-passage times that we derive from BD simulations, these are usually assumed to be inversely proportional to association rate constants. Thus, we may, similarly as our predecessors,7−10,12−14,20,23,26 estimate the overall effects of HI on association rates in the considered receptor−ligand system. Thus, in the case of nonspecific association between neutral objects, the association rate constant evaluated with HI between the binding partners accounted for is 30% smaller than that evaluated in the absence of HI. In the case of specific association between neutral objects, HI decreases the rate constant by 45%. For the strongest attraction considered in this work, we find the decrease in the rate constant to be 16% in the case of nonspecific association and 35% in the case of specific association. For the strongest repulsion considered in this work, the maximal decrease in the rate constant is 60% for the nonspecific association and even larger for the specific association. However, one should keep in mind that as the mobility of the ligand decreases with its decreasing distance from the receptor the above numbers depend on the size of the ligand diffusion space in the studied system. The smaller the diffusion space (or the less dilute the system), the smaller, on average, the receptor−ligand distance and the more pronounced the effects of HI on the association kinetics (see Figure 3). Our findings regarding the effects HI on the kinetics of diffusion-controlled association in the model receptor−ligand system can be summarized as follows. HI slows down both nonspecific and specific associations, the latter to a much larger extent than the former. The effects of HI on the association kinetics increase with an increasing ionic strength when there is an electrostatic attraction between binding partners. However, the opposite is observed when binding partners repel each other and effects of HI on the association kinetics decrease with an increasing ionic strength. Whereas it was previously described for dumbbell molecules that HI effects on the kinetics of their specific association are most pronounced at moderate ionic strengths,13,20 we do not confirm this observation, as in our systems, HI effects are monotonic functions of the ionic strength. This holds true for both nonspecific and specific associations. One possible reason for such a behavior could be that in the receptor−ligand system studied here electrostatic torques are absent as electrostatic forces (either attractive or repulsive) act along the line connecting the diffusion centers of binding partners. In the aforementioned studies,13,20 electrostatic interactions resulted from the distribution of charges, giving rise to nonzero electrostatic torques acting on binding partners. However, at this point, we prefer to refrain from speculations. Considering charge distributions of real biomolecules, one should expect a complicated interplay between hydrodynamic and electrostatic interactions upon their diffusion-controlled association. We are currently investigating this subject. Moreover, the less favorable the energy of electrostatic interactions between the receptor and ligand, the larger the overall HI effect on the association kinetics (a similar effect was described by Frembgen-Kesner and Elcock; these authors

each other. The opposite is observed in the case of repulsion between binding partners. This can be understood easily. When binding partners are kept close to each other by electrostatic attraction, the residence times of the ligand near the receptor are long compared to the characteristic time of ligand rotations. Moreover, the ligand moves toward rather than away from the receptor (in an average sense). As the rotational dynamics of the ligand is much faster than the dynamics of its escape (via translational diffusion) from the vicinity of the receptor, the values of the considered ratio are small. As the attraction becomes stronger (with an increasing α or a decreasing ionic strength), the residence times of the ligand near the receptor τS

become longer and thus the absolute value of the ratio 1 − NS τ becomes smaller. However, when receptor−ligand interactions are weakened or the receptor repels the ligand, the distance between the ligand and receptor increases. Now, in the average sense, the ligand moves away rather than toward the receptor. The translational dynamics of the ligand escaping from the vicinity of the receptor comes into play as the ligand must assume the proper orientation before diffusing away from the receptor. As a τS

result, the absolute values of 1 − NS become larger with τ increasing repulsion. Let us consider a situation in which electrostatic interactions between the ligand and receptor are switched off, that is, α in eq 27 is set to 0. The relative increase of the global first-passage time observed in this case is roughly 145% when HIs are accounted for and roughly 95% when we neglect HI (Figure 8, caption). Thus, one may conclude that the presence of HI between binding partners somehow strengthens the effect that the required stereospecificity has on the association kinetics. Switching the receptor−ligand electrostatic interactions on does not change this conclusion. It holds true for the whole range of ionic strengths considered and regardless of the character (attraction or repulsion) of electrostatic interactions. We interpret this predominantly as the effect of the fact that the translational mobility of the ligand decreases toward the receptor as a result of HI. After an unsuccessful encounter (i.e., an encounter upon which the ligand is not correctly aligned with the receptor cavity), the ligand diffuses away from the receptor. After some time, when the ligand again reaches the receptor via translational diffusion, another encounter (either successful or unsuccessful) takes place. Intervals between successive unsuccessful encounters, which contribute to global first-passage times for specific reactions, are overall longer in systems with enabled HI because of the varying translational mobility of the ligand with the distance to the τS

receptor. As a result, the absolute values of the 1 − NS ratio τ obtained for systems with enabled HI are always larger than those obtained for systems in which HIs are neglected.

■

DISCUSSION We investigated the diffusion-controlled association in the receptor−ligand system. The position and orientation of the receptor molecule are fixed, whereas the ligand molecule undergoes translational and rotational diffusion in the force field resulting from its electrostatic interactions with the receptor. The mobility of the ligand is a function of its position and orientation relative to the receptor due to HIs between the L



■

observed that the association kinetics of slowly associating mutants of barnase and barstar was more affected by HIs than the kinetics of fast-associating mutants of these proteins). Switching off HIs between associating molecules apparently decreases the electrostatic contribution to the association kinetics, that is, the difference between the global first-passage time for the association of charged species and the global firstpassage time for the association of neutral species is smaller when HIs are neglected than the corresponding difference obtained in the presence of HI. This holds true for both nonspecific and specific associations; however, in the considered receptor−ligand model system, larger effects are observed in the case of nonspecific association. The kinetics of the specific association, when evaluated relative to the nonspecific association kinetics, is slowed down by HI because of longer (than in the absence of HI) times needed for a re-encounter of binding partners that diffused away after an unsuccessful encounter. The bead−shell method is able to account for arbitrarily complicated molecular shapes,22,37−39 and the BD algorithm that we employed in our work can be generalized to two-body systems with any number of degrees of freedom. However, we should point out an important bottleneck of the described approach. Calculations of the mobility tensors of the ligand in the presence of the receptor involve operations on large matrices, whose sizes depend on the resolution of bead−shell models of molecules (in the current study, the general mobility tensor for the receptor−ligand system in a particular orientation is of size 61 614 × 61 614), and require significant amount of computer memory. Moreover, the complexity of the linear algebra operations involved makes the evaluation of mobility tensors time-consuming. BD simulations of association between molecules that undergo translational and rotational diffusion in three dimensions are not easily feasible when detailed bead−shell models are employed as they would require mobilities for an enormous number of possible configurations of binding partners to be evaluated. Simplifications are thus unavoidable. The studied system, with its rather idealized molecular shapes (i.e., a spherical receptor with a cavity and an elongated capsule-like ligand, models that were chosen on the basis of the widely accepted consensus that the formation of molecular complexes usually involves some degree of shape complementarity) and with the ligand diffusion limited effectively to 2 degrees of freedom can be considered as only a toy model, but its biological relevance should not be easily dismissed. We argue that our observations and conclusions can, for the most part, be generalized to actual, much more complicated biomolecular systems.

■

REFERENCES

(1) Schreiber, G.; Haran, G.; Zhou, H. X. Fundamental Aspects of Protein-Protein Association Kinetics. Chem. Rev. 2009, 109, 839−860. (2) Alsallaq, R.; Zhou, H. Electrostatic Rate Enhancement and Transient Complex of Protein-Protein Association. Proteins 2008, 71, 320−335. (3) Atkins, P. W. Physical Chemistry, 6th ed.; Freeman: New York, 1998. (4) Alberty, R. A.; Hammes, G. G. Application of the Theory of Diffusion-Controlled Reactions to Enzyme Kinetics. J. Phys. Chem. 1958, 62, 154−159. (5) Calef, D. F.; Deutch, J. M. Diffusion-Controlled Reactions. Annu. Rev. Phys. Chem. 1983, 34, 493−524. (6) Purcell, E. M. Life at Low Reynolds Number. Am. J. Phys. 1977, 45, 3−11. (7) Friedman, H. L. A Hydrodynamic Effect in the Rates of Diffusion-Controlled Reactions. J. Phys. Chem. 1966, 70, 3931−3933. (8) Deutch, J. M.; Felderhof, B. U. Hydrodynamic Effect in Diffusion-Controlled Reaction. J. Chem. Phys. 1973, 59, 1669−1671. (9) Honig, E. P.; Roebersen, G. J.; Wiersema, P. H. Effect of Hydrodynamic Interaction on the Coagulation Rate of Hydrophobic Colloids. J. Colloid Interface Sci. 1971, 36, 97−109. (10) Wolynes, P.; Deutch, J. M. Slip Boundary Conditions and the Hydrodynamic Effect on Diffusion Controlled Reactions. J. Chem. Phys. 1976, 65, 450−454. (11) Brune, D.; Kim, S. Hydrodynamic Steering Effects in Protein Association. Proc. Natl. Acad. Sci. U.S.A. 1994, 91, 2930−2934. (12) Antosiewicz, J.; McCammon, J. A. Electrostatic and Hydrodynamic Orientational Steering Effects in Enzyme-Substrate Association. Biophys. J. 1995, 69, 57−65. (13) 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. (14) 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. (15) Allison, S. Diffusion-Controlled Reactions: Hydrodynamic Interaction Between Charged, Uniformly Reactive Spherical Reactants. J. Phys. Chem. A 2006, 110, 13864−13867. (16) Szymczak, P.; Cieplak, M. Hydrodynamic Effects in Proteins. J. Phys.: Condens. Matter 2011, 23, No. 033102, DOI: 10.1088/09538984/23/3/033102. (17) Ando, T.; Skolnick, J. Crowding and Hydrodynamic Interactions Likely Dominate in Vivo Macromolecular Motion. Proc. Natl. Sci. Acad. U.S.A. 2010, 107, 18457−18462. (18) Ando, T.; Skolnick, J. On the Importance of Hydrodynamic Interactions in Lipid Membrane Formation. Biophys. J. 2013, 104, 96− 105. (19) 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, No. 046011, DOI: 10.1088/1478-3975/7/4/046011. (20) Długosz, M.; Antosiewicz, J. M.; Zieliń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. (21) Długosz, M.; Antosiewicz, J. M. Hydrodynamic Effects on the Relative Rotational Velocity of Associating Proteins. J. Phys. Chem. B 2013, 117, 6165−6174. (22) 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. (23) 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.

AUTHOR INFORMATION

Corresponding Author

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

The authors declare no competing financial interest.

■

Article

ACKNOWLEDGMENTS

Authors acknowledge support from National Science Centre (UMO-2014/13/B/NZ1/00004). M


Article

The Journal of Physical Chemistry B (24) 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. (25) 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. (26) Saglam, A. S.; Chong, L. T. Highly Efficient Computation of the Basal kon Using Direct Simulation of Protein-Protein Association with Flexible Molecular Models. J. Phys. Chem. B 2016, 120, 117−122. (27) Smoluchowski, M. Drei Vorträge über Diffusion, Brownsche Molekularbewegung und Koagulation von Kolloidteilchen. Physik. Z. 1916, 17, 557−571 , 585−599. (28) Smoluchowski, M. Versuch Einer Mathematischen Theorie der Koagulationskinetik Kolloider Lösungen. Z. Physik. Chem. 1917, 92, 129−168. (29) Chandrasekhar, S. Stochastic Problems in Physics and Astronomy. Rev. Mod. Phys. 1943, 15, 1−89. (30) Shushin, A. I. Recombination Kinetics of Molecules with Anisotropic Reactivity and Interaction Potential. Exact Solution for a Small Reactive Region. Chem. Phys. Lett. 1986, 130, 452−457. (31) Ermak, D. L.; McCammon, J. A. Brownian Dynamics with Hydrodynamic Interactions. J. Chem. Phys. 1978, 69, 1352−1360. (32) Dickinson, E.; Allison, S. A.; McCammon, J. A. Brownian Dynamics with Rotation-Translation Coupling. J. Chem. Soc., Faraday Trans. 2 1985, 81, 591−601. (33) Northrup, S. H.; Allison, S. A.; McCammon, J. A. Brownian Dynamics Simulation of Diffusion-Influenced Bimolecular Reactions. J. Chem. Phys. 1984, 80, 1517−1524. (34) Kirkwood, J. G.; Riseman, J. The Intrinsic Viscosities and Diffusion Constants of Flexible Macromolecules in Solution. J. Chem. Phys. 1948, 16, 565−574. (35) Rotne, J.; Prager, S. Variational Treatment of Hydrodynamic Interaction in Polymers. J. Chem. Phys. 1969, 50, 4831−4837. (36) Yamakawa, H. Transport Properties of Polymer Chains in Dilute Solution - Hydrodynamic Interaction. J. Chem. Phys. 1970, 53, 436−443. (37) Garcia de la Torre, J. Building Hydrodynamic Bead-shell Models for Rigid Bioparticles of Arbitrary Shape. Biophys. Chem. 2001, 94, 265−274. (38) 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. (39) 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. (40) Swan, J. W.; Wang, G. Rapid Calculations of Hydrodynamic and Transport Properties in Concentrated Solutions of Colloidal Particles and Macromolecules. Phys. Fluids 2016, 28, No. 011902, DOI: 10.1063/1.4939581. (41) Szabo, A.; Schulten, K.; Schulten, Z. First Passage Time Approach to Diffusion Controlled Reactions. J. Chem. Phys. 1980, 72, 4350−4357. (42) Hänggi, P.; Borkovec, M. Reaction-rate Theory: Fifty Years after Kramers. Rev. Mod. Phys. 1990, 62, 251−341. (43) Dhont, J. K. G. An Introductions to Dynamics of Colloids; Elsevier: Amsterdam, 1996. (44) Godec, A.; Metzler, R. First Passage Time Distribution in Heterogeneity Controlled Kinetics: Going Beyond the Mean First Passage Time. Sci. Rep. 2016, 6, No. 20349, DOI: 10.1038/srep20349. (45) Faas, G. C.; Mody, I. Measuring the Kinetics of Calcium Binding Proteins with Flash Photolysis. Biochim. Biophys. Acta 2012, 1820, 1195−1204. (46) 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.

(47) Boscá, F.; Sastre, G.; Andreu, J. M.; Jornet, D.; Tormos, R.; Miranda, M. A. Drug-Tubulin Interactions Interrogated by Transient Absorption Spectroscopy. RSC Adv. 2015, 5, No. 49451, DOI: 10.1039/C5RA05636E. (48) 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. (49) Lau, A. W. C.; Lubensky, T. C. State-Dependent Diffusion: Thermodynamic Consistency and its Path Integral Formulation. Phys. Rev. E 2007, 76, No. 011123, DOI: 10.1103/PhysRevE.76.011123. (50) 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. (51) Humphrey, W.; Dalke, A.; Schulten, K. VMD - Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (52) Cichocki, B.; Felderhof, B. U.; Hinsen, K.; Wajnryb, E.; Bławździewicz, J. Friction and Mobility of Many Spheres in Stokes Flow. J. Chem. Phys. 1994, 100, 3780−3790. (53) Cichocki, B.; Hinsen, K. Stokes Drag on Conglomerates of Spheres. Phys. Fluids 1995, 7, 285−291. (54) 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. (55) Brune, D.; Kim, S. Predicting Protein Diffusion Coefficients. Proc. Natl. Acad. Sci. U.S.A. 1993, 90, 3835−3839. (56) Wegener, W. A. Transient Electric Birefringence of Dilute RigidBody Suspensions at Low Field Strengths. J. Chem. Phys. 1986, 84, 5989−6004. (57) 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. (58) Zhou, H. X. Kinetics of Diffusion-Influenced Reactions Studied by Brownian Dynamics. J. Phys. Chem. 1990, 94, 8794−8800.

N


Effects of Spatially Dependent Mobilities on the Kinetics of the

Recommend Documents