Solvent Fluctuations Induce Non-Markovian ... - ACS Publications

Subscriber access provided by UNIV OF CALIFORNIA SAN DIEGO LIBRARIES

Article

Solvent Fluctuations Induce Non-Markovian Kinetics in Hydrophobic Pocket-ligand Binding R. Gregor Weiss, Piotr Setny, and Joachim Dzubiella J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b01219 • Publication Date (Web): 23 Mar 2016 Downloaded from http://pubs.acs.org on March 24, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Solvent Fluctuations Induce Non-Markovian Kinetics in Hydrophobic Pocket-ligand Binding R. Gregor Weiß,†,‡ Piotr Setny,¶ and Joachim Dzubiella∗,†,‡ Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, D-12489 Berlin, Germany, Institut für Weiche Materie and Funktionale Materialen, Helmholtz-Zentrum Berlin, Hahn-Meitner Platz 1, D-14109 Berlin, Germany, and Centre of New Technologies, University of Warsaw, 00-927 Warsaw, Poland E-mail: [email protected]

Phone: +49 (30) 8062-42902.

∗ To

whom correspondence should be addressed für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, D-12489 Berlin, Germany ‡ Institut für Weiche Materie and Funktionale Materialen, Helmholtz-Zentrum Berlin, Hahn-Meitner Platz 1, D14109 Berlin, Germany ¶ Centre of New Technologies, University of Warsaw, 00-927 Warsaw, Poland † Institut

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract We investigate the impact of water fluctuations on the key-lock association kinetics of a hydrophobic ligand (key) binding to a hydrophobic pocket (lock) by means of a minimalistic stochastic model system. The latter describes the collective hydration behavior of the pocket by bimodal fluctuations of the water-pocket interface that dynamically couples to the diffusive motion of the approaching ligand via the hydrophobic interaction. This leads to a Markovian set of overdamped Langevin equations in 2D-coordinate-space spanned by the interface position and the ligand position. Numerical simulations demonstrate locally increased friction of the ligand, decelerated binding kinetics, and local non-Markovian (memory) effects in a reduced 1D-description along the ligand’s reaction (distance) coordinate as found previously in explicit-water simulations. Our minimalistic model elucidates the origin of locally enhanced friction that can be traced back to long-time decays in the force-autocorrelation function induced by the spatially fluctuating pocket-ligand interaction. Furthermore, we construct a generalized 1D-Langevin description of ligand binding including a spatially local memory function that reflects the dominant frequencies of the pocket wetting/dewetting process, enabling further interpretation and a semi-analytical quantification of our results.

1

Introduction

Nature expresses a strong versatility in its creation of substrates as ligands and binding sites as receptors, thereby utilizing the complex properties of water as natural solvent environment. This evolutionary framework facilitates the multifaceted kinetics of biomolecular recognition and association and has led to substantial interdisciplinary research in the last decades towards fundamentally comprehending the natural mechanisms of ligand-receptor (or key-lock) binding as a part of life’s cycle. Consequently, one substantive objective that recurs eminently in science is the detailed molecular understanding of ligand binding processes for the design and development of pharmaceutical substances. Many experimental as well as theoretical studies on the thermodynamics of an increasing num2


Page 2 of 30

Page 3 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ber of ligand-enzyme complexes have provided insight about binding free energies, namely binding affinities, of the individual systems. 1–5 Taken alone, however, thermodynamics cannot predict exact kinetic properties. Yet, rates of binding and unbinding events are crucial factors determining drugs’ efficiency. 6–8 Pioneering research recognized dynamic couplings as important component for estimating the time scale of molecular biological processes. 9–13 Beece et al. 9 discussed the impact of structural fluctuations of enzymes on the migration kinetics of substrates. They observed large effects of protein fluctuations on binding rates in the thoroughly explored process of carbon monoxide migration to myoglobin. Here, especially the fluctuations of opening-closing conformations of protein channels largely contribute. Different solvent viscosities, thus different environments to the ligand-enzyme complex, were argued to change the protein’s internal fluctuations which couple to the kinetics of ligand migration. In general, internal barriers governing conformational fluctuations in a protein can be comparable to the thermal energy, 14 facilitating their time scales to be similar to those of ligand diffusion and binding kinetics. 15,16 Specific work on inactive-active, e.g., open-closed, conformational transitions of biomolecular receptors proposed kinetic models based on the induced fit and conformational selection paradigms. 17,18 Within these models the conformational transitions are treated as distinct states taken with given probability and transition rates fulfilling detailed balance. Hence the ratios of transition rates and state-probabilities determine whether ligand migration induces the active conformation for binding, or whether binding occurs predominantly when the pocket conformation is active long before ligand association. Extending this picture, Zhou and co-workers 17,18 allowed for coupling of the conformational kinetics to ligand migration where they utilized Markovian kinetics in a two-state model. A more general discussion 10–12 describes dynamic coupling of substrates and enzymes by an underdamped kinetic description. It models ligand migration by a generalized Langevin equation (GLE) including memory on random velocity changes. The time scales of the memory kernel are incorporated in an additional multiplication factor to conventional transition state theory for

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

rate calculation over a barrier. Hence calculations for individual ligand-pocket systems estimate relative retardation or acceleration to a Markovian crossing rate. This extension to reaction rate theory is also known as Grote-Hynes theory. 10 Direct coupling of water dynamics to hydrophobic key-lock binding kinetics was recently observed by Setny et al. 19 By means of explicit-water molecular dynamics (MD) simulations of a model hydrophobic pocket-ligand system they found long-time correlation effects, i.e., the hint to memory, in position and force correlations when the ligand was situated in the immediate vicinity of the pocket. It was argued that their origin were pocket water occupancy fluctuations that occur due to capillary evaporation in the small confinement between hydrophobic ligand and hydrophobic pocket, 20 yielding dry states, without water inside the pocket, and wet states, with a maximally solvated pocket. Also the ligand position was shown to sensitively affect the bimodal (dry-wet) pocket hydration distribution and time scale. A Markovian description of mean binding times utilizing potential of mean force and a spatially dependent friction could not reproduce binding times directly calculated in the MD simulation. Hence, a non-Markovian treatment of ligand migration within hydrophobic key-lock binding processes was proposed. In a subsequent study by Mondal et al. 21 on MD simulations of a very related pocket-ligand model system the two hydration states were utilized in a reaction-diffusion model similar to the descriptions for inactive-active confirmation kinetics. 17,18 This two-state model then improved the rate predictions compared to direct MD simulations. 21 In this work we introduce a minimalistic stochastic model for the kinetic binding of a ligand to a hydrophobic pocket that exhibits bimodal wet-dry transitions. Here, one stochastic coordinate is the bimodally fluctuating pocket-water interface and the second the position of a ligand that travels in one spatial dimension. Both coordinates are coupled via the hydrophobic interaction between ligand and the fluctuating interface. Mathematically, we describe that by two coupled stochastic equations. In Section 2 the details of the model are described. Section 3 presents a numerical evaluation of the model system analyzing ligand binding kinetics. Comparison of mean binding times from numerical simulation to a corresponding memoryless stochastic process demonstrates

4


Page 4 of 30

Page 5 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the break-down of Markovian behavior in a reduced description with the single reaction (distance) coordinate of the ligand. Friction calculations indicate that additional damping in hydrophobic key-lock association originates from the fluctuating potential on slow time scales. In this way, our analysis that closely follows the procedure in Ref. 19 answers the previously open questions on the origin and nature of local friction and memory effects in hydrophobic pocket-ligand binding. To further corroborate these findings, Section 4 deals with a complementary theory describing an effective 1D-reaction coordinate ligand system in terms of a generalized Langevin equation including a local memory function, which enables further interpretation and a semi-analytical quantification of the results. The local memory kernel is constructed such that it relates to locally confined conditions of the coupled stochastic model. It then enables us to illustrate that the dominant memory timescale reflects the timescale of dry/wet hydration fluctuations in the stochastic 2D-model from Section 2. We conclude our study in Section 5.

2

Stochastic model

Our minimalistic stochastic model assumes that the ligand is a particle diffusing in one spatial direction z close to a binding pocket as depicted in Fig. 1. Water can reside in front of the pocket or penetrate it, giving rise to ’dry’ or ’wet’ pocket states, respectively. 20 This behavior we model by a Brownian pseudo-particle that effectively describes the water interface position and its motion in the pocket region. The two schematic setups in Fig. 1 illustrate the interface-ligand system in its ’dry’ and ’wet’-state, respectively. They show the interface position as a thick blue line, representing a sharp water-vapor surface at zs , and an orange, spherical ligand of radius R at zl . The ligand diffuses with the properties of a spherical particle in water utilizing the Einstein relation D = kB T /6πηR, where kB is the Boltzmann constant, η the viscosity of the bulk solvent and T the temperature. 22 We assume that all effects of pocket hydration fluctuations on the random forces on the ligand and the related frictional and viscous changes are explicitly considered in the fluctuating orthogonal reaction coordinate, i.e., the position of the interface. Dynamic coupling

5



double well

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

dry

|

desolvation 0 potential zs

Page 6 of 30

wet

z

z

R zl

| 0

zl

zs

Fig. 1: Illustration of the stochastic model: The interface position (thick blue ’surface’ line) at position zs fluctuates in the pocket in a double-well potential (Eq. (1), blue curve) in z-direction leaving the pocket either dehydrated/dry (left) or hydrated/wet (right). The ligand (orange sphere) at position zl diffuses freely in one-dimension on the z-axis perpendicular to the pocketed wall. It interacts with the interface by an attractive interaction potential (Eq. (3), red line). As illustrated here the interaction potential moves with interface. Note the ligand diffusion is spatially purely 1D and the 2D-schematics of the pocketed wall (gray) and bulk water area (blue) are just shown for illustration. may occur when interface fluctuations and ligand diffusion are on a similar time scale. Hence, for simplicity we choose equal diffusivity for both ligand and the Brownian pseudo-particle that mimics the interface, facilitating the condition having both time scales with comparable magnitude. In our model the energy scale kB T , the length scale R, and the diffusion constant D set the natural scales in our system. This introduces a Brownian time scale τB = R2 /D. A detailed comment on the relation of units and physical constants to the previous explicit-water MD simulations is provided in the Appendix A. As motivated from previous MD studies, 20 we assume that the interface fluctuates bimodally between positions inside and outside the pocket. This is a model for capillary fluctuations of the water confined in the pocket. We assume the interface moves as a Brownian pseudo-particle in an external double-well potential Vdw (zs ) =

h 2 2 2 z − R + b · zs s R4

(1)

which is drawn as blue curve in Fig. 1. For b = 0, we assume that the positions of the two wells are situated at ±R, and h is the height of the barrier which lies at z = 0. To further enable changes in relative depths of ’dry’ and ’wet’ wells we introduce a bias given by the linearity constant b in 6


Page 7 of 30

ζ (zs , zl )

z

zl zs

(a)

1.5 potential energy

1 0.5 zs /R

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0 −0.5 −1 −1.5

(b)

−1

0

1

2

3

4

zl /R

Fig. 2: (a) If the distance ζ = zl − zs < R then the water-vapor interface (vertical blue line) and the ligand (orange sphere) interact by the desolvation potential (Eq. (3), red curve). The interaction scales by the solvated volume of the ligand, which is the portion on the right hand side of the interface. (b) The key-lock model essentially describes Brownian motion in two dimensions (zs , zl ). Here the plot draws an example trajectory (orange) on the potential energy surface expressed by Eq. (5). units of kB T /R. A pair potential acting between the interface and the ligand accounts for the free energy gain of desolvation as the hydrophobic ligand passes through the water interface (see Fig. 2 (a)). The resulting desolvation potential is designed such that it drags the ligand out of the solvent into the pocket (zl < 0). At the same time, following the principle of action-reaction, the interaction pulls the interfacial water out of the pocket (zs > 0), which conceptually corresponds to the ligand-induced drying transition. For small solutes, the desolvation free energy (essentially entropy) approximately scales linearly with the solvent-excluded volume ∆G ∝ Vol, whereas after the transition at a crossover length-scale lc it is proportional to the solvent-accessible surface area A, ∆G = γ · A with γ as surface tension. 23 Modeling microscopic key-lock binding with a small-sized ligand, we choose the desolvation potential to scale linearly with the solvent-excluded volume.

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 30

We demand a reasonable proportionality constant Γ to fulfill ∆G(lc ) = Γ ·Vol(lc ) ≡ γ · A(lc ) at the crossover length-scale, which thus yields Γ = 3γ/lc

(2)

For pure water surface tension we calculate Γ = 2.95 kB T /R3 (see Appendix A) such that the effective desolvation free energy is roughly 12.36 kB T . This is comparable to the results of explicit water simulations. 24 Once the ligand penetrates the water-vapor interface, the solvent-excluded volume starts to decrease as a function of ligand-interface distance, ζ (zs , zl ) = zl − zs , as it is illustrated in Fig. 2 (a), leading to favorable attraction. The desolvation potential can then be written as 4π 3 π 2 R − (R − ζ ) (2R + ζ ) ∆G(zs , zl ) = Γ 3 3

(3)

which gives a parabolic pair force, acting on particle i = s, l (solvent or ligand) via i Fsol (zs , zl ) = −

d∆G(ζ ) dζ dζ = πΓ(ζ 2 − R2 ) d ζ dzi dzi

(4)

valid for |ζ (zs , zl )| = |zl − zs | < R We emphasize that the gray wall in Fig. 1 embedding the cavity is only drawn representatively. A steric repulsion and van der Waals (vdW) attraction to the wall are omitted. The former is not needed since we have fully adsorbing boundary conditions in the pocket (at zl /R = −1.25). Inclusion of additional vdW forces with the pocket only adds an energetic contribution that is constant in time. It would only serve as a constant energy shift inside the pocket that would slightly tilt the PMFs. Thus, the conclusions about ligand kinetics, friction and memory of our study would be unchanged. Also, we fix a reflective boundary to a given distance zmax to the pocket in order to avoid the ligand diffusing far away from the pocket. Throughout the main body of the paper zmax /R = 5, whereas in the Supporting Information (SI) the impact of varying zmax is discussed. In summary, our model essentially describes Brownian motion in two dimensions spanned by

8


Page 9 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the coordinates zs and zl . One may write the sum of the double-well and desolvation potential as Vtot (zs , zl ) = ∆G(zs , zl ) · Θ [R − |ζ (zs , zl )|] +Vdw (zs )

(5)

which makes evident that the coupling via the desolvation potential simply yields the associated two-dimensional potential surface. An example for arbitrary potential parameters is illustrated in Fig. 2 (b) together with an example trajectory. Then, our nonlinearly coupled stochastic system of equations describes the key-lock model by h ξs z˙s (t) = +πΓ R2 − ζ (zs , zl )2 · Θ [R − |ζ (zs , zl )|] − 4 4 zs (z2s − R2 ) − b + Fs (t) R ξl z˙l (t) = −πΓ R2 − ζ (zs , zl )2 · Θ [R − |ζ (zs , zl )|] + Fl (t)

(6a) (6b)

with ξi as friction coefficients, i = s, l, and Fi (t) denoting δ -correlated random forces fulfilling the fluctuation-dissipation theorem hFi (t)Fi (t0 )i = 2 kB T ξi δ (t − t0 ). Note that both ξs = ξl = 1, since we set the diffusivities of both ligand and pseudo-particle (that models the interface position) equal for simplicity.

3

Numerical simulations

In the following we discuss the binding kinetics obtained from integrating Eq. (6) numerically, where we use the numeric scheme proposed by Ermak and McCammon. 25 We focus on how the stochastic coupling affects the ligand’s motion and kinetics along its (distance) reaction coordinate. This effectively projects the Markovian two dimensional treatment onto the reaction coordinate zl . We thus omit the associated subscript such that zl ≡ z throughout the following Sections 3.1 and 3.2. In general, pocket hydration can be affected by changing its hydrophobicity, geometry or size. Such changes, however, simultaneously affect pocket occupancy and the solvent fluctuation time scale. In our model, we are able to disentangle both effects and their influence on ligand binding. Water occupancy can be tuned by changes in the biasing parameter b of the double well poten9


-1

0

1

2

-2

zs /R

-1

0

1

Vdw /kB T

.4 −4

-2

Page 10 of 30

6 5 4 3 2 1 0

(b) h

b=

(a)

4.4

6 4 2 0 -2 -4 -6

b=

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Vdw /kB T


2

zs /R

Fig. 3: The double-well potential which governs hydration fluctuations in the pocket: The fluctuations are tuned by (a) biasing from b = −4.4 kB T /R to 4.4 kB T /R in steps of 0.8 kB T /R, hence breaking rate symmetry, or (b) barrier height from h = 1 kB T to 5 kB T in steps of 0.5 kB T , thus changing the magnitude of the interface’s characteristic time scales. Black dashed lines draw the reference double well setting (h = 1 kB T , b = 0). Example biasing b or changes in magnitude h are indicated by example guide lines (a) or arrow (b), respectively. (Color coding is maintained throughout the paper.) tial, Vdw from Eq. (1). Potentials with biasing ranging from b = −4.4 kB T /R to b = 4.4 kB T /R in steps of 0.8 kB T /R are drawn in Fig. 3 (a). The black dashed line refers to the reference with barrier height h = 1 kB T and no biasing, which is motivated from explicit and implicit solvent simulation studies. 19,20 The barrier height can separately be tuned from h = 1 kB T to h = 5 kB T in steps of 0.5 kB T as plotted in Fig. 3 (b). Changes in barrier height directly influence the wet-dry transition time and thus the effective interface fluctuation time scale. From Kramer’s rate theory one knows that the rate r of crossing the double well barrier scales exponentially with barrier height r ∝ e−h/kB T . 26 The color coding from both plots in Fig. 3 is consistently adopted to other plots throughout this paper. Further we note that the equilibrium distribution of the water interface depends on the ligand position due to the nonlinear coupling evident from contributions of Eq. (4) in Eq. (6), if |ζ | < R. A schematic plot in Fig. 4 illustrates how the bimodal distribution changes upon the approach of the ligand. When the ligand and the interface interact, their pair potential adds to the effective equilibrium free energy landscape of the interface itself. This can tilt the bimodal distribution of the interface and thus influences the equilibrium wetting behavior of the pocket. We observe that the pocket drying is enhanced for close ligand positions which coarsely mimics also how pocket hydration couples to ligand position in all-atom and implicit solvent simulations. 19,20

10


Page 11 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


zl /R = −1.5 zl /R = 1.5 zl /R = 2.0 zl /R = 2.5

z/R 1.5 2.0 2.5

−1.5

Fig. 4: The equilibrium distribution (blue lines) of the water interface is affected by the ligand position due to the hydrophobic ligand-interface interaction in Eq. (3) and Eq. (4). As the ligand (orange circle) approaches the pocket, the bimodality of the water interface distribution is lost for intermediate states, e.g. zl /R = 1.5, but is recovered if the ligand is fully bound to the pocket, e.g., zl /R = −1.5.

3.1

Mean first passage time and memory

As first measure of ligand binding kinetics the mean first passage time (MFPT) is sampled from each point z to a bound configuration at z f /R = −1.25. Therefore, for each setup with given

biasing b and barrier height h around 2 × 105 trajectories are simulated and analyzed for the ligand

starting at zmax /R = 5 until it is zl /R < −1.25 inside the pocket. Hence, resulting MFPT curves describe the mean first passage time T1 (z, z f ) of the ligand crossing z f /R = −1.25, given it started at z with a reflective boundary at zmax /R = 5. Fig. 5 (a.1) shows the MFPT curves corresponding to simulation setups with varying bias in the double-well potential the interface coordinate is subject to. Regarding setups with a negative bias (greenish lines), hence a preferentially dry pocket, the ligand’s mean binding time is faster than without biasing (black). With growing values for the bias the MFPT slows down by a factor of two. Also, as the barrier height increases, a deceleration of the MFPT is observed in Fig. 5 (b.1). There the MFPT curves exhibit a hunch around z ≈ 1 which is enhanced with growing barrier whereas at z . 0 all curves T1 (z, z f ) in panel (b.1) coincide. Further, we analyze for all considered parameter values (h, b) the potential of mean force (PMF) of the ligand using the weighted histogram analysis method (WHAM) 27–29 along z. Fig. 5 (a.2) shows a strong dependence of the PMFs on the biasing parameter b. Besides small changes in

11


-1 0

1

(a.1)

(b.1)

(a.2)

(b.2)

2 3

4

5 -1 0

z/R

1

2

3

4

18 15 12 9 6 3 0 -1 -3 -5 -7 -9 -11

T1 (z, z f )/τB

18 15 12 9 6 3 0 -1 -3 -5 -7 -9 -11

Page 12 of 30

V (z)/kB T

V (z)/kB T

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

T1 (z, z f )/τB


5

z/R

Fig. 5: MFPT curves T1 (z) and PMF V (z) of ligand binding to the pocket depend on changes of interface’s fluctuation behavior by double-well variations from Fig. 3. Panel (a.1)/(a.2) shows ligand’s MFPT/PMF dependent on double-well biasing b. Panel (b.1)/(b.2) draws ligand’s MFPT/PMF dependent on double-well barrier height h. Blue dash-dotted line in panel (b.1) draws an example MFPT curve T1M (z) using Eq. (7) for the reference setting (h = 1 kB T , no bias) with respective PMF and constant diffusivity. (Color coding is adopted from Fig. 3.) shape, the attractive part of the PMFs essentially shifts towards smaller values of z, if the interface bias shifts towards an increasingly wetted pocket. On the other hand, the ligand PMF negligibly depends on the barrier height, h, as revealed in Fig. 5 (b.2). These PMFs, as equilibrium quantities, are essentially unaffected since the interface kinetics mainly change with h. This is especially noteworthy since the corresponding MFPT curves in panel (b.1) alter relatively strongly with h, suggesting that the effect on ligand binding times originates from modified interface kinetics. In the case of a Markovian process the PMF, V (z), together with the possibly spatially dependent diffusivity, D(z), determine the n-th moment of the first passage time distribution 30

TnM (z, z f ) = n

Zz zf

z βV (z0 ) Zmax 00 M 0e dz00 e−βV (z ) · Tn−1 dz 0 D(z ) z0

(7)

where the zeroth moment T0M = 1 determines normalization and β = kB T −1 = 1. The blue, dashdotted line in Fig. 5 (b.1) is a numerically integrated solution of Eq. (7) using constant diffusivity and a spatially dependent V (z) of the reference case b = 0 and h = 1 kB T . It should be compared to the black, dashed MFPT curve in the same plot. Only for negative z-values both coincide. Effects that can no longer be treated by Markovian kinetics occur around z ≈ 1, where the hunch in 12


(a.1)

(b.1)

(a.2)

(b.2)

6 5 4 3 2 1 0 12 10 8 6 4 2 0

(T1 − T1M )/τB

6 5 4 3 2 1 0 12 10 8 6 4 2 0

σ2 /τB

σ2 /10τB

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(T1 − T1M )/τB

Page 13 of 30

-1 0 1 2 3 4 5 -1 0 1 2 3 4 5

z/R

z/R

Fig. 6: Markovian assumption breaks down as manifested in differences of MFPT T1 (z) from simulation to T1M (z) from Eq. (7) and in memory index σ2 (z) from Eq. (8). Panels (a.1) and (a.2) show the respective measures dependent on the double well bias b. Panels (b.1) and (b.2) draw both measures dependent on double well barrier height h. Note that the effect is much bigger as the biasing towards a wet pocket increases. Therefore the scale of (a.2) is σ2 /10τB . (Color coding is adopted from Fig. 3.) simulated MFPT curves qualitatively deviates to a dent in the solution of Eq. (7). For even bigger values of the reaction coordinate, shapes of the MFPT curves of both methods only conform, but the Markovian solution represents overall faster association. Together, the ligand dynamics can only be modeled by a pure Markovian description when the ligand is inside the pocket. As a general measure, calculating MFPT curves T1M in a Markovian picture for all considered cases of bias strengths and barrier heights, and using the respective PMFs, enables direct observations where the deviations in simulation occur. Accordingly, the difference T1 (z, z f ) − T1M (z, z f ) is plotted in Fig. 6 (a.1) and (b.1). In all cases the difference vanishes inside the pocket and increases towards a maximum situated just in front of the pocket mouth. It then plateaus to a constant positive value for large z, which indicates slowed ligand kinetics in all considered cases of pocket water fluctuations. For the cases of biased wetting, the difference T1 − T1M in Fig. 6 (a.1) is very small, if the pocket is preferably dry, namely with a strong negative bias. As the biasing parameter b increases, and thus the interface’s distribution tends towards mainly hydrated pocket states, the difference in MFPT accumulates to a peak when b = 2.8 kB T /R, and alleviates for even higher bias. In Fig. 6 (b.1) deviations to the Markovian picture are enhanced by growing barrier height, hence slowed-down wetting fluctuations. 13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 30

In order to investigate further the break-down of Markovian dynamics and possibly accompanied memory, we additionally determine the so-called memory index 31 σ2 =

T2 (z, z f ) − T2M (z, z f ) T1 (z, z f )

(8)

introduced by Hänggi et al. 31 who noted its possible value to ligand migration studies. It provides an additional spatially resolved measure of the character of a random process. It indicates a process to be non-Markovian if the difference of the second moment of a random process to the second moment of the corresponding process with Markovian assumption does not vanish. Here, the memory index vanishes for ligand positions far from the pocket in all considered cases of pocket wetting behavior in Fig. 6 (a.2) and (b.2). Note, however, that finite σ2 values are reached already at intermediate positions, that is z/R ∼ 3, where the ligand is still out of reach of the desolvation potential (ζ (zs , z) > R). Finally, the memory index peaks at the position z/R . 1, at which initial deviations in the first moments occur in Fig. 6 (a.1) and (b.1). Subsequently, for z/R < 1, it steeply recedes to zero. Inside the pocket σ2 diverges once more, majorly due to numerical inaccuracies close to the target distance as these are amplified by division in Eq. (8). Also, actual non-Markovian effects reoccur inside the pocket as it becomes evident from results and discussions in following sections and in the SI. All together it is inept to assume constant ligand diffusivity within a Markovian description of our ligand-interface system in order to estimate correct mean binding time. The process rather indicates non-Markovian contributions, slowing the ligand binding which predominantly arise in the region where the ligand and the interface start interacting. Overall, the first and second moments are influenced far outside where the ligand is driven by delta correlated random noise by definition in Eq. (6). In addition, strongly dissimilar MFPTs but almost similar PMFs in Fig. 5 (b.1) and (b.2), respectively, suggest that local coupling of interface time scales alone is enough to impact ligand binding kinetics.

14


Page 15 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.2

Spatially dependent friction

In order to illuminate the kinetic effect we restrict further analysis to systems considering an unbiased interface distribution, and only vary the barrier height h in the double-well potential. For these particular cases, a purely kinetic nature of the observed effects was evident from the perturbation of average ligand binding times occurring with unchanged PMFs, as shown in Fig. 5 b. For clear comparison of ligand kinetics and free energy we replot the MFPT and PMF data for varying double well barrier height in Fig. 7 together with the following results on spatially dependent ligand friction. Spatially resolved friction ξl (z) can be obtained from the position auto-correlation function (PACF), hδ z(t)δ z(0)ir , with δ z(t) = z(t) − zr , derived from simulations with a harmonically restrained ligand at zr . A detailed description of the ligand umbrella simulation is provided in SI. Integration of the auto-correlation function provides the time scale which is divided by the square of the position fluctuations, yielding the local friction 32

ξ (zr ) =

R∞ 0

hδ z(t)δ z(0)ir dt hδ z2 i2r

(9)

Solid curves in Fig. 7 (c) show Gaussian fits to ξ (zr ), gathered from PACF evaluation (see also SI). Simulations restraining the ligand far away from the pocket yield the preset friction of value one. However, ξ (z) strongly peaks in front of the pocket mouth, where the ligand is subject to interaction with the interface. Growing barrier height, and thus exponentially slowed double-well transition rates of the interfacial motion, increases the peak up to a factor of approximately 85. It indicates that the effect arises due to the ligand interacting with the bimodally fluctuating interface. While the interface penetrates the pocket, the ligand, still remaining around z ∼ 2.0, is only subject to the δ -correlated random force of the Brownian model (Eq. (6)). Whereas when the interface is in the outer well, in front of the pocket, the desolvation potential acts on the ligand. Hence the desolvation force acts as additional fluctuating force which introduces additional friction. Peaking friction occurs in regions in which fluctuations are most pronounced, where interface and ligand

15


V (z)/kB T ξ M (z), ξ (z)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

T1 (z)/τB


15 13 11 9 7 5 3 1 0 -2 -4 -6 -8 -10 -12 90 80 70 60 50 40 30 20 10 0

Page 16 of 30

(a)

h

(b)

ξ M (z) ξ (z) (c)

-1.0 0.0

1.0

2.0

3.0

4.0

5.0

z/R

Fig. 7: Panel (a) draws again the MFPT curves from simulation already shown in Fig. 5 (b.1). Here it also compares to a MFPT curve T1M (z) for the reference setting calculated from Eq. (7) with PMF and spatially dependent friction from PACF plotted as green dash-dotted line. Panel (b) draws again the PMFs already shown in Fig. 5 (b.2). Panel (c) plots spatially resolved friction ξ (z) from fits to PACF data (Eq. (9), solid) and from MFPT data using (Eq. (10), dashed) with Einstein relation ξ M (z) = kB T /DM (z). (Color coding is adopted from Fig. 3.) might interact or not, and thus the desolvation potential can essentially be on or off. Together, with the observations throughout the previous MFPT analysis it seems that the additional force fluctuations serve as a source of memory. Their time scale is proportional to the Brownian time, τB , and scales exponentially with double well barrier height: τdw ∝ τB eh/kB T (see also Eq. (18)). Thus the time scale of the desolvation force fluctuations does not separate from the time scale of ligand migration, which leaves memory to the association process. As it will be shown in the following section a generalized Langevin model derives an exponential growth of the friction peak values with barrier height, which directly relates to exponentially growing time scales of the interface fluctuations. Again we obtain the MFPT curve T1M (z, z f ) using Eq. (7), with the PMF V (z) of the reference case (h = 1 kB T, b = 0), and now additionally with a spatially resolved diffusion D(z) = kB T · ξ −1 (z) 16


Page 17 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


from Einstein’s relation, 22 with previously evaluated ξ (z). Only in the interval z ∈ (−1.25, 0.5) the result, plotted as green dash dotted line in Fig. 7 (a), coincides with MFPT curves from simulation (black dashed), and moreover, with evaluation of Eq. (7) without spatially resolved friction (blue dash dotted). Subsequently, a steep edge in the curve yields values which overestimate the actually simulated results far outside, z 2. So, on one hand, the solution of Eq. (7) overestimates the results using both spatially resolved profiles V (z) and ξ (z). On the other hand, it is underestimated using only spatially resolved PMF, but constant friction of value one. For comparison we also calculate spatially dependent profiles ξ M (z) = kB T / DM (z) by solving Eq. (7) for D(z) ≡ DM (z) as primarily introduced by Hinczewski et al. 33 eβV (z) D (z) = ∂ T1 (z)/∂ z M

Z zmax z

0

dz0 e−βV (z )

(10)

Note that ξ M (z) uses the Markovian assumption, and thus, is certainly not the proper friction profile fulfilling fluctuation-dissipation theorem for our non-Markovian ligand migration process. In detail, it does not measure the quantity friction/dissipation which can be proportionally related to the system’s fluctuations. However, it will trivially reproduce the correct MFPT T1 from simulation when using it in Eq. (7). Dashed lines in Fig. 7 (c) show curves for ξ M (z) which similarly peak in front of the pocket mouth. In exact comparison to ξ (z) from PACF, the results from Eq. (10) show a positional shift closer to the pocket and differ in peaking value as well as peak width. (Additional information on how ξ M (z) depends on the choice of reflective boundary is discussed in the SI.) Again, observing essential discrepancy between local friction ξ (z) from simulation and a spatial profile ξ M (z) from a memoryless picture in Eq. (10) illustrates the location and strength of non-Markovian effects within the stochastic system (6). Also, whether calculating MFPTs from Eq. (7) either with constant friction or spatially resolved friction ξ (z), yields under- and overestimated results, respectively. In neither case the Markovian property is fulfilled, which is the requirement for obtaining a proper system description with Eq. (7). We emphasize that our system exhibits similar non-Markovian effects as those resolved by explicit water simulations from

17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 30

Setny et al. 19 We note that in principle rate calculations on multi-dimensional energy landscapes are accessible by Langer theory 34,35 which could be utilized in cases of high barriers in our simplified model with two degrees of freedom. Our goal, however, is to characterize the relations between pocket water fluctuations and ligand diffusion in order to mark off methods that can be used in more complex techniques like MD simulations when the full energy landscape is more complicated to extract. Hence our work focuses on the reduction on the one-dimensional description with ligand position as reaction coordinate for which we want to rationalize the yet found impact of memory and friction.

4

Generalized Langevin model

Having identified fluctuations of the desolvation potential as the origin of local memory and friction in the ligand’s reaction coordinate we show in the following how their impact can be quantified. Further, this section shall formulate the proper stochastic characteristics when dealing with the ligand coordinate alone. For simplicity it focuses on local conditions, namely, constrained ligand position. Yet it elucidates a proper non-Markovian formulation to classify possible treatment with conventional theory. Generally in the case of a known memory kernel η(t) one can directly investigate the corresponding one-dimensional generalized Langevin equation (GLE) ∂Veq (q) − mq(t) ¨ =− ∂q

Z t

dt 0 η(t − t 0 )q(t ˙ 0 ) + F (t)

(11)

with mass m, equilibrium potential Veq (q), and a random force fulfilling fluctuation dissipation hF (t)F (t 0 )i = 2 kB T η(t − t 0 ). Simple systems of two coupled Langevin equations can be analyt-

ically contracted onto a one-dimensional GLE 36,37 and vice versa. A prominent example is that of an underdamped Brownian particle in a harmonic potential. For the coupled system described by Eq. (6) analytic contraction from 2D to 1D is not feasible due to higher than harmonic coupling and 18


Page 19 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


nonlinearity in the double-well potential. Therefore we reinterpret a method which is usually used to expand a one-dimensional GLE to a set of two coupled equations without memory. With that we are able to approximate friction from local conditions of the pocket-ligand system. We restrict the analysis to the location of the friction peaks discussed above and can predict the peaking value max(ξ (z)) as function of barrier height in a bimodally fluctuating force. To this end we reverse the approach from Pollak et al. 36,38,39 It originally extends a onedimensional GLE of reaction coordinate q such as Eq. (11) by an auxiliary variable x to receive two coupled equations. Each of the resulting equations then omits memory and only the auxiliary variable x is driven by a temporally delta correlated random force, N(t). Taking unit mass m = 1, the GLE (11) is mapped on the two dimensional, underdamped system ∂V (q, x) =0 ∂q ∂V (q, x) ξ x˙ + = N(t) ∂x q¨ +

(12a) (12b)

The driving noise N(t) is a delta correlated Gaussian noise hN(t)i = 0, hN(t)N(t 0 )i = 2kB T ξ δ (t − t 0 )

(13)

There are two further requirements that memory η(t) and coupling potential V (q, x) must fulfill for proper mapping: 38,39 (a) The kernel η(t) may be represented by a sum of exponentials, and for this very example even η(t) =

ξ −t/τ e ≡ Ωe−t/τ τ

(14)

(b) The coupling between auxiliary and reaction coordinate should be harmonic such that ∂V (q, x) dVeq (q) d f (q) = − Ω[x − f (q)] ∂q dq dq

19


(15)


h/kB T

1

2

4

5 200

(a)

Eq. (20) Fig. 7 (c) Eq. (21)

Eq. (18) simul.

0 2

3

180 160 140

(b)

(c)

120 100

1

ξ

1

log10 (ξ ) log10 (τdw /τB )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 30

80 60

0

40 20 0 1

2

3

4

5

h/kB T

Fig. 8: Panel (a) plots time scale τdw (h) of the interfacial motion in the double-well without ligand. Black line refers to Eq. (18) and orange squares are values obtained by interface’s position autocorrelation from simulation. Panel (b) shows ξ (h) (Eq. (20)) from constructed GLE using ε = 0.36 (black line) and blueish shaded area delimited by gray dashed lines indicates the range covered by Eq. (20) with standard deviation ∆ε = ±0.24. The green lines shows ξ (h) from Eq. (21). the The peaking values from ξ (z) in Fig. 7 (c) are here shown as circular symbols with corresponding color coding from Fig. 7. Panel (c) plots the same data as in (b) on a log-scale.

In our case, let us focus on the situation at the position of the friction peak, max(ξ (z)) in Fig. 7 (c). In that case, the expansion T1ε of the solvation force in Eq. (4) at fixed ligand-interface distance ε up to first order, with respect to a perturbation δ ζ = ζ − ε, gives the harmonic contribution of our solvation coupling T1ε Fsol (ζ ) = πΓ[ε 2 − R2 ] + 2πΓε(ζ − ε) + O(ζ 2 )

(16)

This identifies the memory kernel constant, Ω = ξ /τ ≡ 2πΓε, by comparison with Eq. (15). The value of ε ≈ 0.36 ± 0.24 is estimated from simulations constraining the ligand at the position of the friction peaks from Fig. 7. It is the mean and standard deviation of the distribution of distance ζ between constrained ligand and bimodally fluctuating interface. Detailed evaluations of ε are discussed in the SI. The set of coupled equations of motion of a free ligand (here ql ) coupled to an

20


Page 21 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


auxiliary variable xs are adopted to the requirements (a) and (b) described above such that q¨l − 2πΓε(ql − xs − R) = 0 ξ x˙s + 2πΓε(ql − xs − R) = N(t)

(17a) (17b)

Note that the system of the two above equations is not fully equivalent to the original coupled stochastic system (6). With the aim to formulate the influence of interface fluctuations on local friction encountered by the ligand, it describes only a single system configuration, which determines the ε. In general, ε depends on ligand position in our original coupled key-lock system. A striking difference is that Eq. (17) does not implement the double-well itself. Rather the time scale determining the memory is chosen to be that of a Brownian particle in a double-well. A compact approximate solution of that time scale is given by 40,41 √ τB (eβ h − 1) p (π β h + 21− β h ) τdw (h) = p 2 2β h 2β h

(18)

To confirm this approximation for the setups which we previously considered, we probe the time scale of interface fluctuations in the double-well within simulations without ligand. Tuning the barrier height from h = 1 kB T to h = 5 kB T reveals that the approximate formula (18) is in very good agreement within the range of interest in h (Fig. 8 (a)). The GLE corresponding to Eq. (17) with memory from Eq. (18) is followingly given by q¨l (t) = −2πΓε

Z t

0

dt 0 e−(t /τdw ) q˙l (t 0 ) + F (t)

(19)

Comparison of its memory kernel with Eq. (14) determines the friction for the constructed system such that ξ (h) = 2πΓε · τdw (h) ≡ Ω · τ

(20)

Fig. 8 (b) and (c) demonstrate the strong resemblance of both systems, the fully coupled keylock binding model (6) and the non-Markovian model (19). The circular symbols with error bars 21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 30

from Gaussian fits draw the maxima of the friction peaks max(ξ (z)) (from Fig. 7) from PACF calculations of the original key-lock model (6). The black line draws expression (20) found for ξ (h). The blue shade indicates the error from variance calculations to ε. In a more rigorous approach one can calculate the average interface-ligand distance utilizing the expression of the energy landscape Vtot (zs , zl ) from Eq. (5). If we fix the ligand to zl /R = 1.7, which is roughly the peaking friction position in Fig. 7, the Boltzmann factor exp[−βVtot (zs , 1.7)] serves as weighting function to calculate the average distance ε = hζ (zs , 1.7)i = h1.7 − zs i = 1.7 − hzs i. In that case the expression for the barrier dependent friction constant in the GLE system reads dzs zs eVtot (zs ,1.7,h) ξ (h) = 2πΓ 1.7 − R dzs eVtot (zs ,1.7,h) R

! · τdw (h)

(21)

where the dependence on barrier height of Vtot (zs , 1.7, h) is emphasized. In Fig. 8 this expression for ξ (h) is drawn as green curve. This also suitably approximates the barrier height dependence of the peaking friction from simulation. The assumption of fixed zl /R = 1.7, however, suits less the conditions of minor ligand position fluctuations within simulations. Hence extracting ε from simulations for Eq. (20) yields a slightly better approximations here.

5

Concluding Remarks

Our investigations presented here reveal the origin of increased friction and additional memory in hydrophobic pocket-ligand binding as it was observed in previous work using all-atom simulations. 19,21 We employ a simple stochastic model of two nonlinearly coupled equations, each driven with memoryless Gaussian noise. One equation models pocket hydration in terms of a continuously diffusing pseudo-particle mimicking the water-pocket interface fluctuating between dry and wet states. Another one describes a ligand, freely diffusing in one-dimension and subjected to a (de)solvation force when in contact with the interface, motivated from the well-known result that the (de)solvation free energy of a microscopic hydrophobic ligand scales with the solvent-excluded volume. 23 The coupling is bi-directional as it affects not only the ligand diffusion but also the mean 22


Page 23 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


position and fluctuation time scale of the interface. The model then enables the investigation of the effects of tuneable pocket hydration by biasing the balance between ’wet’ and ’dry’ states, as well as by changing the time scale of fluctuations in a controlled fashion. Bimodal water occupancy distributions were observed in explicit-water simulations of related systems. 19 The first passage time analysis of ligand binding revealed non-Markovian contributions to binding kinetics, where ligand association kinetics from numerical simulation of the full 2D model was decelerated in comparison to an effective 1D picture in which the equilibrium PMF and a spatially constant friction was utilized. Detailed analysis of the first and the second moments of the first passage times indicated the strongest non-Markovian contributions to occur shortly before binding, at positions where the ligand is subject to intermittent interactions with the bimodally fluctuating interface. While increasing the fluctuation time scales leaves the ligand (equilibrium) PMFs unchanged, the deceleration of ligand binding is enhanced, evidencing that ligand kinetics couples to the time scale of water fluctuations due to the fluctuating desolvation potential, which in all-atom simulations refers to a fluctuating PMF. 21 We resolved the effective friction of the ligand in a reduced 1D description by calculating position auto-correlation functions and found spatially dependent friction in front of the pocket. The friction was found to peak at positions where coupling to the water interface occurs most prominently. In comparison, the friction peak determined from an inversion of the MFPT that assumed Markovian behavior was found to be shifted in position due to memory effects, as previously indicated already by the explicit-water simulations. 19 The enhanced local friction is essentially induced by additional (long-time) force fluctuations in the random force of the ligand due to the fluctuating desolvation force acting in a limited spatial range. We note that it was recently shown that in general for stochastic dynamics on multi-dimensional energy landscapes effective local friction and possibly memory effects will occur in a reduced, lower dimensional description if relevant time scales in the system’s dynamics are not separable. 42 We further corroborate the origin of memory by constructing a generalized Langevin model (19) restricted to local conditions of the original two coupled stochastic equations (6). The GLE uti-

23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

lized the coupling desolvation potential between the ligand and the interface, as well as the autocorrelation time of the interface fluctuation as time scale of the memory kernel. By construction the dominant frequency of colored noise in the GLE is given by the time scale of bimodal pocketdewetting. When comparing friction of the constructed GLE to local modulations in ligand friction from simulations of the coupled stochastic system (6) both models quantitatively coincide. Again this substantiates that non-Markovian kinetics within the coupled stochastic system is the outcome of projecting the binding kinetics onto the ligand reaction coordinate alone. Note that in systems like Eq. (6) the colored noise is spatially dependent due to the local interface-ligand coupling in a limited interaction range. For extreme cases with strictly bimodal fluctuations, a two-state approach has been successfully applied recently, 21 suggesting possible consideration in reaction-diffusion models, multistate models, 17,43 or even Markov-state models. 44 Other studies 45 make evident that ligand binding is not necessarily accompanied by bimodal pocket occupancy fluctuations. If in such cases memory remains, corrections to a Kramer’s rate of ligand binding and unbinding, can be determined within Grote-Hynes theory 10 and generalizations of it. 11 Beyond these known kinetic approaches to ligand-receptor binding it remains to future work to develop fundamental and worthwhile treatments that consider the general issue of position-dependent memory in ligand binding.

A

Units and constants

The stochastic system in Eq. (6) is strongly inspired by the all-atom simulation setup in Ref. 19 Therefore all units are directly related to the system of the corresponding MD setup. Those were conducted in ambient conditions which is why effective temperature during the numerical simulations of Eq. (6) relates to T = 300 K. The particle size of the ligand is that of a methane molecule which is roughly R = 0.4 nm which sets the Brownian length scale R. With viscosity η ≈ 10−3 Pa·s

of water 46 the corresponding diffusion constant relates to D = kB T /6πηR = 0.54 nm2 ns−1 .

For the parametrization of the desolvation potential we assume the crossover length-scale at

24


Page 24 of 30

Page 25 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


lc = 1 nm which relates to lc = 2.5 R. With surface tension γ ≈ 15.36 kB T /nm2 for water 47 we calculate the solvation volume scaling constant to be Γ = 3γ/lc = 46.1 kB T /nm3 = 2.95 kB T /R3 .

Supporting Information. Details on umbrella sampling setups; data and fits of PACF calculations; system size effects in the utilized Markovian description; numeric estimates of mean forcing for GLE

Acknowledgement R. Gregor Weiß, Piotr Setny, and Joe Dzubiella wish to express their sincere gratitude to their mentor J. Andrew ’Andy’ McCammon for all the exciting opportunities in science, his neverending support, and his kind hospitality in San Diego, California. The authors thank the Deutsche Forschungsgemeinschaft (DFG) for financial support of this project.

References (1) Gellman, S. H. Introduction: Molecular recognition. Chem. Rev. 1997, 97, 1231–1734. (2) Gilson, M. K.; Zhou, H.-X. Calculation of protein-ligand binding affinities. Annu. Rev. Biophys. Biomol. Struct 2007, 36, 21–42. (3) Woo, H.-J.; Roux, B. Calculation of absolute protein-ligand binding free energy from computer simulations. Proc. Natl. Acad. Sci. (USA) 2005, 102, 6825–6830. (4) Wang, J.; Zheng, X.; Yang, Y.; Drueckhammer, D.; Yang, W.; Verkhivker, G.; Wang, E. K. Quantifying intrinsic specificity: A potential complement to affinity in drug screening. Phys. Rev. Lett. 2007, 99, 198101. (5) Head, R. D.; Smythe, M. L.; Oprea, T. I.; Waller, C. L.; Green, S. M.; Marshall, G. R. VALIDATE: A new method for the receptor-based prediction of binding affinities of novel ligands. J. Am. Chem. Soc. 1996, 118, 3959–3969.

25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(6) McCammon, J. A. Theory of biomolecular recognition. Curr. Opin. Struct. Biol. 1998, 8, 245–249. (7) Pan, A. C.; Borhani, D. W.; Dror, R. O.; Shaw, D. E. Molecular determinants of drug-receptor binding kinetics. Drug Discov. Today 2013, 18, 667–673. (8) Tiwary, P.; Limongello, V.; Salvalaglio, M.; Parrinello, M. Kinetics of protein-ligand unbinding: Predicting pathways, rates, and rate-limiting steps. Proc. Natl. Acad. Sci. 2014, 112, E386–E391. (9) Beece, D.; Eisenstein, L.; Frauenfelder, H.; D.Good,; Marden, M. C.; L.Reinisch,; Reynolds, A.; Sorensen, L.; Yue, K. Solvent viscosity and protein dynamics. Biochemistry 1980, 19, 5147–5157. (10) Grote, R. F.; Hynes, J. T. The stable states picture of chemical reactions. II. Rate constants for condensed and gas phase reaction models. J. Chem. Phys. 1980, 73, 2715–2732. (11) Hänggi, P. Physics of ligand migration in biomolecules. J. Stat. Phys. 1983, 30, 401–412. (12) Doster, W. Viscosity scaling and protein dynamics. Biophys. Chem. 1983, 17, 97–103. (13) Frauenfelder, H.; Wolynes, P. G. Rate theories and puzzles of hemeprotein kinetics. Science 1985, 229, 337–345. (14) Batista., P. R.; Pandey, G.; Pascutti, P. G.; Bisch, P. M.; Perahia, D.; Robert, C. H. Free energy profiles along consensus normal modes provide insight into HIV-1 protease flap opening. J. Chem. Theory Comput. 2011, 7, 2348–2352. (15) Ishikawa, H.; Kwak, K.; Chung, J. K.; Kim, S.; Fayer, M. D. Direct observation of fast protein conformational switching. Proc. Natl. Acad. Sci (USA) 2008, 105, 8619–8624. (16) Bernardi, R. C.; Cann, I.; Schulten, K. Molecular dynamics study of enhanced Man5B enzymatic activity. Biotech. for Biofuels 2014, 7, 83. 26


Page 26 of 30

Page 27 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(17) Zhou, H.-X. From induced fit to conformational selection: A continuum of binding mechanism controlled by the timescale of conformational transitions. Biophys. J. 2010, 98, L15– L17. (18) Cai, L.; Zhou, H.-X. Theory and simulation on the kinetics of protein-ligand binding coupled to conformational change. J. Chem. Phys. 2011, 134, 105101. (19) Setny, P.; Baron, R.; Kekenes-Huskey, P.; McCammon, J. A.; Dzubiella, J. Solvent fluctuations in hydrophobic cavity-ligand binding kinetics. Proc. Natl. Acad. Sci. (USA) 2013, 110, 1197–1202. (20) Setny, P.; Wang, Z.; Cheng, L.-T.; Li, B.; McCammon, J. A.; Dzubiella, J. Dewettingcontrolled binding of ligands to hydrophobic pockets. Phys. Rev. Lett. 2009, 103, 187801. (21) Mondal, J.; Morrone, J. A.; Berne, B. J. How hydrophobic drying forces impact the kinetics of molecular recognition. Proc. Natl. Acad. Sci. (USA) 2013, 110, 13277–13282. (22) Einstein, A. Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen. Ann. Phys. (Berlin, Ger.) 1905, 322, 549–560. (23) Chandler, D. Interfaces and the driving force of hydrophobic assembly. Nature 2005, 437, 640–647. (24) Setny, P.; Baron, R.; McCammon, J. A. How can hydrophobic association be enthalpy driven? J. Chem. Theory Comput. 2010, 6, 2866–2871. (25) Ermak, D. L.; McCammon, J. A. Brownian dynamics with hydrodynamic interactions. J. Chem. Phys. 1978, 69, 1352–1360. (26) Hänggi, P.; Talkner, P.; Borkovec, M. Reaction-rate theory: Fifty years after Kramers. Rev. Mod. Phys. 1990, 62, 251–341.

27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(27) Ferrenberg, A. M.; Swendsen, R. H. Optimized Monte Carlo data analysis. Phys. Rev. Lett. 1989, 63, 1195–1198. (28) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method. J. Comput. Chem. 1992, 13, 1011–1021. (29) Zhu, F.; Hummer, G. Convergence and error estimation in free energy calculations using the weighted histogram analysis method. J. Comput. Chem. 2012, 33, 453–465. (30) Siegert, A. J. F. On the first passage time probability problem. Phys. Rev. 1951, 81, 617–623. (31) Hänggi, P.; Talkner, P. Memory index of the first passage time: A simple measure of nonMarkovian character. Phys. Rev. Lett. 1983, 51, 2242–2245. (32) Hummer, G. Position-dependent diffusion coefficients and free energies from Bayesian analysis of equilibrium and replica molecular dynamics simulations. New J. Phys. 2005, 7, 34. (33) Hinczewski, M.; von Hansen, Y.; Dzubiella, J.; Netz, R. R. How the diffusivity profile reduces the arbitrariness of protein folding free energies. J. Chem. Phys. 2010, 132, 245103. (34) Langer, J. S. Statistical theory of decay of metastable states. Ann. Phys. 1969, 54, 258–275. (35) Zhou, H.-X. Rate theory for biologists. Q. Rev. Biophys. 2010, 93, 219–293. (36) Risken, H. The Fokker-Planck equation: Methods of solution and application; Springer, Berlin, 1996; pp 179–195. (37) van Kampen, N. G. Stochastic processes in physics and chemistry; North Holland, Amsterdam, 2007; pp 240–242. (38) Pollak, E.; Berezhkovsii, A. Fokker-Planck equation for nonlinear stochastic dynamics in the presence of space and time dependent friction. J. Chem. Phys. 1993, 99, 1344–1346.

28


Page 28 of 30

Page 29 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(39) Berezhkovsii, A. M.; Frishman, A. M.; Pollak, E. Variational transition state theory for multidimensional activated rate processes in the presence of anisotropic friction. J. Chem. Phys. 1994, 101, 4778–4789. (40) Cregg, P. J.; Crothers, D. S. F.; Wickstead, A. W. An approximate formula for the relaxation time of a single domain ferromagnetic particle with uniaxial anisotropy and collinear field. J. Appl. Phys. 1994, 76, 4900–4902. (41) Kalmykov, Y. P.; Coffey, W. T.; Waldron, J. T. Exact analytic solution for the correlation time of a Brownian particle in a double-well potential from the Langevin equation. J. Chem. Phys. 1996, 105, 2112–2118. (42) Berezhkovsii, A. M.; Szabo, A. Time scale separation leads to position-dependent diffusion along a slow coordinate. J. Chem. Phys. 2011, 135, 074108. (43) Szabo, A.; Shoup, D.; Northrup, S. H.; McCammon, J. A. Stochastically gated diffusioninfluenced reactions. J. Chem. Phys. 1982, 77, 4484–4493. (44) Hummer, G.; Szabo, A. Optimal dimensionality reduction of multistate kinetic and MarkovState Models. J. Phys. Chem. B 2015, 119, 9029–9037. (45) Mondal, J.; Friesner, R. A.; Berne, B. J. Role of desolvation in thermodynamics and kinetics of ligand binding to a kinase. J. Chem. Theory Comput. 2014, 10, 5696–5705. (46) Smith, P. E.; van Gunsteren, W. F. The viscosity of SPC and SPC/E water at 277 and 300 K. Chem. Phys. Lett. 1993, 215, 315–318. (47) Vega, C.; de Miquel, E. Surface tension of the most popular models of water by using the test-area simulation method. J. Chem. Phys. 2007, 126, 154707.

29



Graphical TOC Entry

local ligand friction

local ligand friction

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ligand distance to pocket

pocket-dewetting timescale

30


Page 30 of 30

Solvent Fluctuations Induce Non-Markovian ... - ACS Publications

Recommend Documents