Principles for Tuning Hydrophobic LigandâReceptor Binding Kinetics

Article pubs.acs.org/JCTC

Principles for Tuning Hydrophobic Ligand−Receptor Binding Kinetics R. Gregor Weiß,*,†,‡ Piotr Setny,*,¶ and Joachim Dzubiella*,†,‡ †

Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, D-12489 Berlin, Germany Institut für Weiche Materie and Funktionale Materialen, Helmholtz-Zentrum Berlin, Hahn-Meitner Platz 1, D-14109 Berlin, Germany ¶ Centre of New Technologies, University of Warsaw, Banacha 2c, 02-097 Warsaw, Poland ‡

S Supporting Information *

ABSTRACT: We investigate how to tune the rate of hydrophobic ligand−receptor association due to the role of solvent in adjustable receptor pockets by explicit-water molecular dynamics (MD) simulations. Our model considers the binding of a spherical ligand (key/guest) to a concave surface recess in a nonpolar wall as receptor (lock/host). We systematically modify the receptor’s physicochemical properties in terms of geometry and dispersion attraction which, in turn, alter the water occupancy and fluctuations within the pocket. We demonstrate that even minor pocket modifications can lead to a significant acceleration of the watermediated association. For example, the binding switches from comparably slow to fast if the binding pocket becomes only slightly deeper. We find that the degree of hydrophobicity, characterized by hydration occupancy and its fluctuations, clearly correlates with the binding times and, for instance, links the sudden acceleration to an abrupt increase in hydrophobicity. For a deeper analysis based on passage time theory, we quantify the intimate coupling between solvent fluctuations and the ligand’s local dynamics and friction. The coupling exhibits substantial nonequilibrium effects and maximizes shortly before binding, which slows down the binding kinetics in all cases. In summary, we rationalize how the physicochemical properties of a nonpolar, concave binding site tune key-lock binding kinetics due to water-mediated forces and fluctuations. Our study thus complements the profound understanding of the solvent’s influence in host−guest binding, which is essential for tailored solutions in catalysis and pharmaceutical applications.

1. INTRODUCTION The principal solvent of nature is water, which is the protagonist of hydrophobic effects.1−6 Their wide momentousness is exemplified in a variety of phenomena in natural, engineering, and pharmaceutical sciences, such as self-assembly of amphiphilic molecules to biological membranes,1,7 receptor− ligand binding,6,8−14 catalysis using cavitands,15−17 transport through nanopores,18−20 as well as protein folding, stability, and function.21−24 Having a common foundation in the incompatibility of water and nonpolar solutes, the hydrophobic effects are quite diverse in their physical nature. In particular, they exhibit a characteristic length scale dependence,1 which is a consequence of the different water structure, dynamics, and fluctuations around small and large hydrophobic objects25−35 and leads to an exchange between enthalpic and entropic thermodynamic driving forces.36 Hydrophobic solvation and association thus uniquely depends on the complex geometrical and chemical topology of a solute’s surface.14,37 This is expressed by nature by creating characteristic binding motifs (e.g., specifically shaped binding sites and pockets that steer protein association or ligand binding by means of the so-called key-lock principle).38,39 In particular, deep hydrophobic binding sites perturb the water © XXXX American Chemical Society

hydrogen bond network leaving many of these intrasolvent bonds unsatisfied,8,9,40 which can trigger wet−dry oscillations or even complete water expulsion from sterically accessible regions.3,39,41−43 Such an effect of dewetting transitions on binding thermodynamics has been examined recently.36 Turning the focus to affinities in hydrophobic ligand binding, early experimental studies intended to tune the binding affinity of enzyme pockets by chemical modifications and reported enhanced binding to increasingly hydrophobic binding sites.13 Similarly, synthetic host−guest systems, which can be utilized for catalysis, proofed to be highly adaptable and amenable to hands-on, guided mutations.16,17 Recent simulation studies indicated that tubularly confining protein pockets comprise stronger contributions to a ligand’s binding affinity than a shallower, thus hydrophobically less-confining binding site.12,39 Hence the influence of hydrophobicity on ligand binding to pockets is, in sum, tunable by chemical and topographical modifications of the binding partners.14,44,45 Our growing understanding of physical phenomena in living organisms leads to the recognition that efficacy and safety of a Received: March 1, 2017 Published: May 11, 2017 A

DOI: 10.1021/acs.jctc.7b00216 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 1. Tuning the ligand binding kinetics by pocket geometry and dispersion. (A) The simulation setup comprises a hydrophobic wall with a concave surface recess as a binding site, TIP4P water, and a spherical model ligand. The ligand is constrained to move only along the z axis. The pocket entrance is at z = 0. The hydrophobicity is tuned by modifying the dispersion interaction ϵ of the water and the wall, or changing the pocket geometry, as illustrated in (B) the tabulated illustration shows the different geometry settings that are used. The gray pocket illustrations differ in curvature and depth as indicated by the tuples (R, Δ). We choose the geometry (9.5, 7)Å as a reference system. (C) Definition of pocket geometry: when increasing the radius R of the blue spherical probe volume, the curvature of the axisymmetric pocket changes. If the position of the probe volume is shifted, the depth Δ of the pocket changes. (D) The mean binding time ; from a distance of 8 Å increases if the water−wall dispersion interaction ϵ is enhanced (i.e., hydrophobicity lowered) in the (9.5, 7)Å geometry. (E) If the depth of the pocket with R = 9.5 Å increases, the binding time abruptly accelerates, if deepening the pocket from Δ = 7 Å to Δ = 8 Å. (F) Decreasing pocket radius is slowing the binding kinetics while Δ = 7 Å. (G) Changing the depth in pockets with radii R = 11 and 9 Å also exhibits strongly accelerating binding times for deep pockets. The symbol and color coding of panels (D−G) is adopted throughout the paper.

that the additional friction can be quantitatively explained by the (solvent) fluctuations. Given the length scale dependence of the hydrophobic effects, and their not yet fully resolved influence on host−guest binding kinetics, we address in this paper the question of how such kinetics vary in response to changes in different aspects of the binding pocket hydrophobicity. Our in silico represented hydrophobic model comprises a spherical ligand, of the size and nature of a methane molecule, and a concave binding site in a nonpolar wall, as was introduced by Setny and co-workers.8−10,41,50 It enables systematic modifications of the host− guest system’s topography and its water attraction strength. In particular, we tune the curvature and depth of the binding pocket on an Ångström length scale which is likely to fall within the range of natural conformational fluctuations of single, flexible host molecules. We also scale receptor hydrophobicity by modifying water−wall dispersion interaction, thus reflecting changes likely to occur due to packing fluctuations in proteins. In order to obtain qualitative, generic trends, we then analyze the variety of obtained binding times by correlating them with phenomenological hydrophobicity scales (e.g., as signified by water fluctuations).25−32 Finally, in a subsequent in depth analysis for the purpose of a more quantitative theoretical assessment, we uncover further the intricate details of nonequilibrium hydrophobic association processes. In particular, we further corroborate our previously derived relation between solvent fluctuations and ligand friction52 and apply a nonequilibrium correction to quantitatively reproduce binding times derived from atomistic MD simulations. Our study helps understanding the molecular determinants of hydrophobic binding kinetics and shall thus be helpful in guiding the molecular engineering of proteins, drugs, or host−guest systems for biotechnological applications.

chemical agent in open, in vivo systems is highly influenced not only by binding equilibrium to its target receptors but also by its kinetic rates.46−49 Accordingly, a drug’s association and dissociation rates determine how quickly a drug can bind and how long it will remain bound to the receptor (i.e., its efficacious time frame). Thus, the binding affinity, described by the binding constant Ka, alone does not suffice as sole parameter for tailored drug development. One rather has to keep in mind that the binding constant relates to the ratio of the association rate kon and the dissociation rate koff Ka =

kon koff

(1)

This seemingly simple relation, however, constitutes a nontrivial principle for tailored key-lock solutions: even though the ratio of kon/koff might stay the same, absolute changes in the respective rates might be affected by modification of transition barriers.46 Very recently, solvent fluctuations were found to create a kinetic barrier in hydrophobic key-lock model simulations.50−52 Especially in the study by Setny et al.,50 it was demonstrated that enhanced wetting transitions in a hydrophobic model pocket can lead to locally increased friction encountered by a vicinal hydrophobic ligand. These findings were also consistent with observations from Berne and co-workers on a similar pocket-ligand system51 as well as for the association of two hydrophobic plates53 or two fullerene spheres.54 Within these studies, Berne et al. additionally demonstrated that an enhanced friction can be tuned while changing the dispersion interaction of water with the hydrophobic objects. In that context, we recently proposed a description of the friction peak by a minimalistic model of a set of two coupled stochastic differential equations which reproduce the kinetic barrier of the key-lock binding process.52 With that, we demonstrated B


Article


2.2. Constrained Ligand Simulations. In order to probe selected observables as functions of the ligand−pocket separation, we utilized umbrella sampling simulations. In each umbrella setup, the ligand was restrained to a given position zi using an external harmonic potential with spring constant K = 835 kJ mol−1 nm−2. The first umbrella potential was placed 3.14 Å from the pocket’s bottom, namely the inner crystal layer of the pocket, such that in terms of the ligand−pocket separation it was situated at, for example, z = −4 Å in a pocket with Δ = 7 Å and shifted respectively for different pocket depths. Up to 48 additional umbrella windows with 0.5 Å spacing were introduced to cover increasing ligand−pocket distances. Individual umbrella simulations were carried out for 10 ns with a step size of 2 fs, while the ligand coordinate was stored for every time step and the water coordinates were stored every 20 fs. The potential of mean force (PMF) was obtained by the weighted histogram analysis method (WHAM).56,57 Additionally, we evaluated the average pocket water occupancy ⟨N⟩ and its fluctuation autocorrelation as a function of zi:

2. METHODS The principal MD setup is akin to the hydrophobic pocket-ligand model of ref 50. It comprises a concave surface recess in a hydrophobic wall and a spherical hydrophobic ligand, which are solvated by a slab of TIP4P water molecules, as illustrated in Figure 1A. The wall is mirrored along the z axis, on both sides of the ∼4 nm wide water slab, technically implementing two identical, independent binding pockets which effectively correspond to a single binding site with a reflective boundary at a distance of z ∼ 2 nm.50 The spherical ligand is placed on the symmetry axis of the axisymmetric pocket in the z direction and is only free to move along that. The pocket entrance defined as the plane containing the outermost layer of wall particles sets the zero point (z = 0) of the reaction coordinate. The hydrophobic model ligand interacts with a conventional Lennard-Jones potential with the parameters (σL = 3.73 Å, ϵL = 1230.09 J/mol). The molded wall is also composed of Lennard-Jones spheres with (σP = 4.15 Å, ϵP = 2.40 J mol) which are arranged in a hexagonal closed-packed crystal structure. Details on the underlying model parameters can be found in ref 50 and citations therein. We systematically adjust the geometry of the pocket as shown in Figure 1B, illustrated by the bluish spherical probe volume in Figure 1C: the probe volume’s inverse radius R determines the inner pocket curvature, and its position along z determines the pocket depth Δ. We define a given geometric setup by the tuple (R, Δ) in units of Ångstrom. The pocket depth varies in 1 Å steps between 4 and 10 Å, while efficacious variations of the pocket radius scan the values R = 9, 9.5, 10, and 11 Å. As reference geometry, we choose the (9.5, 7)Å setup which represents the well-studied system in, for example, refs 41 and 50, for which the pocket hydration bimodally fluctuates between wet and dry states. Finally, the hydrophobicity of the wall material is tuned by the wall−water oxygen interaction with a reference effective value of ϵ = 39 J/mol.50 This parameter was initially chosen in a previous work8 to construct a finely grained pocket structure with high particle density which still resembles the interaction of a paraffin wall. For the reference geometry (9.5, 7)Å, we consider an additional parameter range from ϵ = 29 J/mol to ϵ = 49 J/mol in steps of 5 J/mol. Thus, we change the wall−water Hamaker constant by approximately ±13% and ±25%, respectively.55 Note that the interaction of the ligand neither to the water nor to the wall is modified by this procedure. This is also why we emphasize only the wall−water interaction parameter. All simulations were prepared as NVT ensembles at T = 298K with a Nosé Hoover thermostat, with a coupling frequency of 1 ps for equilibration runs and 10 ps for production runs. Note that only the solvent, for example, the water, was coupled to the thermostat throughout the production runs. Although for technical reasons we used constant volume simulations, and the density of the water slab was individually tuned for each system to (997.0 ± 0.3) g/dm3 by slight adjustments of wall separation. All simulations were performed with the GROMACS-4.6.3 simulation package. 2.1. Unconstrained Ligand Simulations. For each setup, we stored a production run of 20 ns in steps of 0.2 ps in which the ligand was fixed at the reflective boundary. This primary trajectory served as a source for randomly seeded initial configurations for subsequent binding event simulations. We then sampled >1000 binding events starting from independent initial configurations by randomly picking one from the previously generated production run. To ensure a selection’s randomness, upon possible reselection, we applied an additional annealing step: within a short simulation of 50 ps, the configuration was heated up to 350 K in a stochastic integrator scheme. After this, the heated configuration was equilibrated for 100 ps to 298 K using the Velocity Verlet algorithm. In the final production run, the ligand was released, free to move along the z direction. The runs were terminated once the ligand was bound 2 Å within the pocket, namely at zf = −2 Å. The time for binding, that is, the first passage time (FPT) was then sorted to calculate FPT distributions and the mean FPT, ; .

CδNδN (t |zi) =

⟨δN (t )δN (0)⟩i ≈ χN (zi)e−t / τN(zi) ⟨N ⟩i

(2)

where δN(t) = N(t) − ⟨N⟩i is the deviation from the mean and ⟨...⟩i denotes the average taken at a given umbrella position zi. The approximate relation on the r.h.s. of eq 2 defines the fluctuation strength χN(zi) and its time scale τN(zi) at zi. Throughout the paper, we will absorb an explicit spatial dependence into the indices [e.g., χN(zi) = χN,zi and especially χN(z → ∞) = χN,∞]. Note that the fluctuation time scale τN,zi quantifies slow (capillary) dewetting modes which we obtain by fitting the exponential tail of CδNδN(t|zi). The local friction ξ(zi) of the ligand was evaluated by58 ∞

βξ(zi) =

∫0 ⟨δz(t )δz(0)⟩i dt ⟨δz 2⟩i2

(3)

where ⟨δz(t) δz(0)⟩i is the position autocorrelation function (PACF) with δz(t) = z(t) − ⟨z⟩i. The thermal energy kBT = β−1 is a constant throughout. If the ligand is far away, the friction in the bulk water slab measures βξ(z → ∞) = βξ∞ = (0.40 ± 0.02) ns nm−2. All the obtained profiles are summarized in the Supporting Information.

3. RESULTS AND DISCUSSION 3.1. Mean Binding Time and Pocket Hydration. Since the process of the ligand approaching from far away to the pocket is standard diffusion, we focus in our following discussion on the mean binding time ; from z = 8 Å, expressed by ; (8 Å). This refers to a little more than one ligand diameter and is the location at which the kinetic coupling of pocket water occupancy fluctuations and ligand friction starts to occur (see ref 50 and the Supporting Information). The plots in Figure 1 (panels D−G) show how the mean binding time ; (8 Å) is affected by modulated pocket properties. It is evident that an increasing binding site affinity to water by increasing pocket−water dispersion interaction, as well as decreasing the level of confinement by increasing the pocket radius, both lead to gradual slowing down of the binding kinetics, as shown in Figure 1 (panels D and F, respectively). The opposite effect is achieved while increasing the pocket depth (Figure 1, panels E and G). Interestingly, here we observe an abrupt reduction of the binding time once a certain critical pocket depth is reached, which is related to a rapid increase in the amplitude of the solvent density fluctuations, as is discussed in the following. In a unifying view, it is the pocket’s general hydrophobicity resulting from the combination of physicochemical properties (affinity to water) and topography (the level of confinement) C


Article


pockets increase the binding time on average and vice versa. The blue line illustrates the aforementioned abrupt speedup for deepening pocket geometries. Apparently, increasing the pocket depth slightly past a certain critical value can significantly enhance the occupancy fluctuations, manifesting the transition to complete dewetting in an increasingly confining geometry. Intuitively, in such cases, water is more easily replaced by the ligand, thus allowing faster binding. In the extreme case of an all-time dry binding site and χN → ∞, the ligand is expected to bind most quickly. On the contrary, as the pocket occupancy fluctuations decrease, the pocket dehydration is less likely. Hence for χN → 0, the binding time substantially increases, maybe even diverges, ; → ∞, given that the pocket is permanently wet. The above conclusions, clearly indicate the link between the level of hydrophobicity and the association speed, similar to the previously observed relation between hydrophobicity and binding affinity.12,13,39 In order to investigate the directness of such correspondence, here we employ the equilibrium profiles of the potential energy landscape V(z) and the friction profile ξ(z) (see also the Supporting Information), with the most accurate approach based on the Smoluchowski equation to estimate the mean binding times59 for the sake of a subsequent comparison with direct MD-based observations:

that directly governs binding kinetics. Adopting the solventrelated perspective, water fluctuations χN,∞ have been recently established as an index for a surface’s degree of hydrophobicity.25−32 In the thermodynamic limit (i.e., for large N), the fluctuation strength χ N,∞ relates to a isothermal compressibility.25 It is proportional to the second cumulant of the pocket water occupancy distribution p(N), which in turn quantifies the free energy of desolvating the volume of interest. Namely, large values of χN,∞ indicate tendencies for dewetting and thus express stronger hydrophobicity. Both quantities, ⟨N⟩ and χN,∞ are expected to respond significantly to changes in pocket properties due to the sensitivity of dewetting (or capillary evaporation) to the geometry of the nanoconfinement and its affinity to water.2,25

;theory(z) = β

∫z

z f

dz′ξ(z)e βV (z′)

z max

∫z′

dz″e−βV (z″) (4)

As a primary test to the above conclusions from Figures 1 and 2A, we neglect the spatial dependence of ξ because it is hardly accessible in experiments. Thus, if we use eq 4 with a constant ligand resistance (i.e., βξ∞ = 0.4 ns nm−2), we can estimate a lower boundary to the observed binding time.50,52 In Figure 2D, we draw the observable ; (8 Å) against this theoretical estimate ;theory (8 Å). The comparison supports the general, qualitative conclusion that the binding time remains to be dominated by the binding affinity, whereas the latter is known to vary due to solvent repulsion effects (i.e., hydrophobicity).12,39 Importantly, however, as indicated in Figure 2D, we observe a systematic shift of the current theoretical predictions, ;theory (8 Å), with respect to the simulated observables, ; (8 Å), toward shorter binding times. This error, highly unsatisfying for quantitative rate predictions, apparently arises from neglecting the friction’s variations with the ligand−pocket separation. Unfortunately, the inclusion of spatially dependent friction into the analysis of hydrophobic association is highly nontrivial,50,52−54 being also inaccessible in experiments.51 In the following section, we elaborate on the nature of the ligand friction and propose a model allowing its inclusion in a simple, yet quantitatively satisfactory manner. 3.2. Kinetic Coupling. As reported previously,50−54 solvent fluctuations intimately couple to ligand kinetics, such that the effective ligand friction along the reaction coordinate, ξ(z), becomes spatially dependent and is locally enhanced prior to association. For example, in our reference setup (Figure 3), the density fluctuations χN (red ○), their time scale τN (gray ○), and the ligand friction ξ (black ○) exhibit a clear maximum around z = zp ∼ (3−4) Å, while the solvent within the pocket region becomes additionally confined by the approaching ligand. The solvent effects can be explained by the occurrence of a free energy barrier for pocket capillary evaporation,41,50 which clearly separates mostly solvated (wet) and mostly

Figure 2. Binding times ; from z = 8 Å correlate with (A) the pocket water occupancy fluctuations χN,∞ (characterizing the pocket’s solvent repulsion) and (B) the mean occupancy normalized by the maximal water occupancy, ⟨N⟩/Nmax. Low hydration and large water fluctuations give rise to faster binding of the ligand. The blue line connects the data points from Figure 1E and illustrates that the abrupt binding speedup is associated with an abrupt increase of the solvent fluctuations. Solid black lines are empirical fitting laws of the form A + Bxα with A = 0.31, B = 0.27, and α = −1.26 in (A) and A = 0.26, B = 0.90, and α = 1 in (B). In plot C, the sampled mean binding time ; (8 Å) is compared to the theoretical values ;theory (8 Å) using eq 4 with constant friction ξ∞ = 0.4 ns nm−2 and PMF V(z). The horizontal distance to the black linear slope of one quantifies the discrepancy Δ; = ; (8 Å)− ;theory (8 Å). Symbol types and colors correspond to those of Figure 1 (panels D−G).

In Figure 2A, we show the correlation of the mean binding times ; (8 Å) from Figure 1 (panels D−G) with fluctuation strength, χN,∞, whereas the correlation with average pocket occupancy, ⟨N⟩/Nmax (here Nmax is the individual maximal occupation number), is shown in Figure 2B. The comparisons suggest that as a general trend, favorably wet, less hydrophobic D


Article


Figure 4. (A) The product on the r.h.s. of eq 5, using simulation results of τN,zp and δFzp at the peak position zp, is plotted against values of βξ(zp) from PACF sampling in eq 3. The black line draws a linear regression with a slope of one. Thus, the enhanced ligand friction is determined and tuned by changes in the pocket occupancy fluctuation time scale τN,zp and the rate of change of the average forces δFzp. The symbol and color coding is adopted from Figure 1 (panels D−G). (B) The occupancy fluctuations χN(z) are plotted in the example of the reference geometry (9.5,7)Å. Its spatial dependence from unconstrained simulations (solid line) differs to that from constraining simulations (dashed line). The red dashed and red solid lines indicate two different peak positions zp.

Figure 3. Potential of mean force (PMF), local friction, and local hydration fluctuations for the (9.5, 7)Å reference setup. The lower part plots the average PMF V(z) (black □) with black arrows, which indicate its fluctuations. The upper part of this figure shows how the ligand friction ξ (black ○; right axis), the pocket occupancy fluctuations χN (brown ○) and their time scale τN (gray ○) depend on the pocket-ligand separation z. The local occupancy fluctuations χN(z) and τN(z) inside the pocket are 1 order of magnitude higher than in a similarly sized bulk probe volume, namely χbulk and τbulk.

kinetic barrier.50 As a result, the fluctuation peak is shifted by roughly 3 Å toward the pocket interior (Figure 4B), while retaining a similar shape and amplitude to the one derived from umbrella simulations. Thus, the pocket’s solvation response in constraining simulations barely materializes in the same way during the enhanced nonequilibrium binding processes in unconstrained simulations. In other words, the nonequilibrium solvent effects (i.e., the solvent fluctuations) manifest differently in the two simulation schemes. Consequently, in reference to the aforementioned kinetic coupling, we assume that the increased friction must similarly shift to smaller values of z. Thus, as a straight solution to correct the error of the mean binding time predictions evidenced in Figure 2D, we suggest to use the friction peaks from umbrella simulations but to translate them by the individual shifts in χN(z) from constrained and unconstrained simulations (see the Supporting Information). In Figure 5A, we illustrate how the theoretical estimate for the FPT in our most hydrophilic system with the (9.5, 7)Å geometry and ϵ = 49 J/mol improves upon the correction for the existence of the kinetic barrier and its subsequent shift. First, utilizing constant (bulk) friction ξ(z) = ξ∞ within eq 4, as before for Figure 2D, yields a curve for ;theory(z) (black dash-dotted), which certainly underestimates the actual simulation sampled binding time curve (blue symbols with error bars). If we simply use the friction profile from umbrella simulations (reddish dashed line), the theoretical binding time (black dashed line) lies above the simulation sampled binding times. Yet, shifting the friction (reddish solid line) and using it within eq 4 yields a result that very well coincides with the sampled values. As one can see in Figure 5A, the estimates for ;theory(z) which utilize constant bulk friction or the simplified friction profile from umbrella sampling assess a lower and an upper boundary, respectively, for the mean binding time (gray area). A more reliable computation incorporates the nonequilibrium effects that are not sampled in the ligand constraining simulations. Here we have incorporated parts of these by a straightforward shift of the equilibrium friction profile. As one

desolvated (dry) pocket states, and governs the time scale of enhanced fluctuations. The strong fluctuations point to the fluctuation−dissipation theorem, which states that a system’s thermal fluctuations depicts its resistance and vice versa. This means that, in our case, the strong solvent fluctuations depict the enhanced ligand resistance before association. Concentrating on the region of locally perturbed ligand diffusion around zp, we applied a nonMarkovian stochastic theory to derive the relation between the enhanced solvent fluctuations time scale, τN(zp), and the friction peak, ξ(zp).52 Accordingly, we found that the magnitude of the ligand friction peak is linearly proportional to the occupancy fluctuation time scale: ξ(z p) = ξz p = δFz pτN , z p

(5)

with the scaling factor, δFzp, corresponding to the rate of change of the solvent mean force, determined by the second spatial derivative of the PMF, δFzp = (∂2V(z)/∂z2)zp.52 The relevance of this dependence for various, modified pocket−ligand systems considered in the current study is confirmed by reasonable agreement between estimates obtained with the use of eq 5 and friction peaks directly evaluated from simulations based on the ligand PACF (eq 3), as shown in Figure 4A. So far we concentrated on the presence of the kinetic coupling, the coupling of fluctuations and friction, by means of umbrella sampling simulations. Yet, while the ligand is freely diffusing and eventually binds to the pocket, its passage time occurs on a similar time scale as the solvent response. As evidenced by the ratio of τn/τbulk in Figure 3, the time scale of water fluctuations in the binding pocket is natively ∼10 times longer than that in a similarly sized bulk water probe volume and becomes even longer while the ligand passes through the E


Article


underlying physical effects, seem to offer appealing possibilities for controlled adjustments for instance within the scope of synthetic cavitands or drug discovery. As one of the first steps toward understanding the link between hydrophobicity and binding kinetics, we investigated how the rate of hydrophobic key-lock/host−guest association can be steered while systematically modifying physicochemical properties of a binding site and thus the resulting solventmediated effects. The spectrum of our geometric modifications is likely to reflect naturally occurring shape variations on an Ångström length scale of receptor binding sites. Similarly, the magnitude of scaling the water−wall interaction reflects fluctuations in a receptor’s interaction strength with water, resulting from reconfigurations of interacting chemical groups or modifications of the packing density. Also the properties of synthetic hosts (e.g., synthetic cavitands),15,17 generally are adjustable within the considered range such that the binding speed could be steered. These comprise current challenges in engineered (bio)chemical processes for synthetic cavitands15,17 and drug discovery46, for example. Our general, somehow intuitively expected, finding is that the association speed positively correlates with an increasing degree of hydrophobicity.12,13,39 However, our results less intuitively highlight that the binding kinetics are not directly governed by host−guest interactions but rather indirectly, substantially mediated by solvent effects. Specifically, we point out that the degree of pocket wetting and the amplitude of solvent density fluctuations can be considered as standalone descriptors of the binding times. Such a dependence on solvent behavior can sometimes lead to drastic effects, such as switching from slow to fast binding by increasing the pocket depth merely by 1 Å, which is enough to abruptly change the magnitude of solvent fluctuations. In a second part, we dissected the maximized ligand resistance prior to binding due to the high solvent fluctuations. We employed the fluctuation−dissipation relation (eq 5), which we have previously derived52 from a non-Markovian kinetic description. This relation expresses the intimate coupling of the enhanced pocket water occupancy fluctuation time scale and the accompanying increase in ligand friction and applies to all investigated physicochemical pocket variations. Thus, the slowing ligand friction can be directly tuned by means of the fluctuation time scale and is a potential obstacle that one must overcome if binding rates shall be optimized. Yet, as a remaining challenge, future investigations must fill the gap for accurate water solvation models that determine the solvent (interface) fluctuations due to a solute’s surface physicochemistry. Finally, we laid out that quantitative binding time predictions must account for the nonequilibrium nature of the kinetic coupling. We found that the process of dewetting occurs delayed if comparing unconstrained and constrained simulations. The fluctuations materialize spatially shifted if the ligand is or is not constrained in an umbrella potential. Thus, probing the friction under ligand constraining conditions must subsequently incorporate a delaying spatial shift, in order to advance the theoretical binding time estimates. In doing so, we essentially were able to empirically improve the binding time predictions. Future work along these lines will show for what kind of realistic scenarios these nonequilibrium effects are important when understanding and guiding the binding kinetics are desired.

Figure 5. (A) illustrates how the ligand friction influences the theoretical binding time estimation in the example of the most hydrophilic pocket (i.e., the (9.5, 7)Å-geometry with ϵ = 49 J/mol). The dotted black line draws eq 4 given ξ(z) = ξ∞. The dash-dotted black line plots the binding time while using ξ(z), shown as a dashed reddish curve. The black solid line represents the binding time curve using the friction profile (reddish solid line) that accounts for the nonequilibrium effect by a negative spatial shift. The latter well coincides with the simulation sampled ;(z), which is drawn as blue symbols with error bars. (B) The symbols plot ; (8 Å) from simulations against the respective estimate ;theory (8 Å) from eq 4 utilizing the accordingly shifted ξ(z). The black line has a slope of one and quantifies the remaining discrepancies. The symbol and color coding is adopted from Figure 1 (panels D−G).

consequence, the impact of the friction peak is actually diminished by the nonequilibrium effects, since the binding time would be slower without the shift (see again the black dashed line in Figure 5A). This and further detailed information on how the enhanced friction impacts the binding time is presented in the Supporting Information. Finally, we present how the mean binding time prediction improves for all/most systems, if eq 4 includes a positiondependent friction, that is yet shifted by the respective delay observed from occupancy fluctuations. In Figure 5B, we the show binding times measured from MD simulations, ; (8 Å), versus the new theoretical predictions, ;theory (8 Å). The black line with a slope of one now corroborates that the binding time estimates improve if the friction peaks are shifted accordingly. Note that possible peak height or width modulations due to nonequilibrium effects are not accessible.

4. SUMMARY AND CONCLUDING REMARKS It is becoming increasingly appreciated that the optimization of active compounds for their overall efficacy in open, for instance in vivo, systems must reach beyond a simple maximization of the affinity to a respective receptor. One important aspect, increasingly amenable to in silico studies, is the kinetics of assembly with the target (macro)molecule. In this respect, hydrophobicity-mediated interactions, with their wide variety of F


Article


■

(15) Shi, Q.; Mower, M. P.; Blackmond, D. G.; Rebek, J. J. Watersoluble cavitands promote hydrolyses of long-chain diesters. Proc. Natl. Acad. Sci. U. S. A. 2016, 113, 9199−9203. (16) Hillyer, M. B.; Gibb, B. C. Molecular shape and the hydrophobic effect. Annu. Rev. Phys. Chem. 2016, 67, 307−329. (17) Wanjari, P. P.; Gibb, B. C.; Ashbaugh, H. S. Simulation optimization of spherical non-polar guest recognition by deep-cavity cavitands. J. Chem. Phys. 2013, 139, 234502. (18) Hu, Y.; Devegowda, D.; Striolo, A.; Phan, A.; Ho, T. A.; Civan, F.; Sigal, R. The dynamics of hydraulic fracture water confined in nanopores in shale reservoirs. J. Unconv. Oil Gas Resour. 2015, 9, 31− 39. (19) Phan, A.; Cole, D. R.; Striolo, A. Aqueous methane in slitshaped silica nanopores: high solubility and traces of hydrates. J. Phys. Chem. C 2014, 118, 4860−4868. (20) Phan, A.; Cole, D. R.; Weiß, R. G.; Dzubiella, J.; Striolo, A. Confined water determines transport properties of guest molecules in narrow pores. ACS Nano 2016, 10, 7646−7656. (21) Rose, G. D.; Geselowitz, A. R.; Lesser, G. J.; Lee, R. H.; Zehfus, M. H. Hydrophobicity of amino acid residues in globular proteins. Science 1985, 229, 834−828. (22) Dill, K. A. Dominant forces in protein folding. Biochemistry 1990, 29, 7133−7155. (23) Kauzmann, W. Some factors in the interpretation of protein denaturation. Adv. Protein Chem. 1959, 14, 1−63. (24) Camilloni, C.; Bonetti, D.; Morrone, A.; Giri, R.; Dobson, C. M.; Brunori, M.; Gianni, S.; Vendruscolo, M. Towards a structural biology of the hydrophobic effect in protein folding. Sci. Rep. 2016, 6, 1−9. (25) Mittal, J.; Hummer, G. Static and dynamic correlations in water at hydrophobic interfaces. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 20130−20135. (26) Patel, A. J.; Varilly, P.; Jamadagni, S. N.; Hagan, M. F.; Chandler, D.; Garde, S. Sitting at the edge: How biomolecules use hydrophobicity to tune their interactions and function. J. Phys. Chem. B 2012, 116, 2498−2503. (27) Patel, A. J.; Varilly, P.; Chandler, D. Fluctuations of water near extended hydrophobic and hydrophilic surfaces. J. Phys. Chem. B 2010, 114, 1632−1637. (28) Jamadagni, S. N.; Godawat, R.; Garde, S. Hydrophobicity of proteins and interfaces: Insights from density fluctuations. Annu. Rev. Chem. Biomol. Eng. 2011, 2, 147−171. (29) Gianti, E.; Delemotte, L.; Klein, M. L.; Carnevale, V. On the role of water density fluctuations in the inhibition of a proton channel. Proc. Natl. Acad. Sci. U. S. A. 2016, 113, E8359−E8368. (30) Xi, E.; Remsing, R. C.; Patel, A. J. Sparse sampling of water density fluctuations in interfacial environments. J. Chem. Theory Comput. 2016, 12, 706−713. (31) Underwood, R.; Ben-Amotz, D. Communication: Length scale dependent oil-water energy fluctuations. J. Chem. Phys. 2011, 135, 201102. (32) Remsing, R.; Patel, A. J. Water density fluctuations relevant to hydrophobic hydration are unaltered by attractions. J. Chem. Phys. 2015, 142, 024502. (33) Xu, Y.; Havenith, M. Perspective: Watching low-frequency vibrations of water in biomolecular recognition by THz spectroscopy. J. Chem. Phys. 2015, 143, 170901. (34) Laage, D.; Stirnemann, G.; Hynes, J. T. Why water reorientation slows without iceberg formation around Hydrophobic solutes. J. Phys. Chem. B 2009, 113, 2428−2435. (35) Weiß, R. G.; Heyden, M.; Dzubiella, J. Curvature dependence of hydrophobic hydration dynamics. Phys. Rev. Lett. 2015, 114, 187802. (36) Setny, P.; Baron, R.; McCammon, J. A. How can hydrophobic association be enthalpy driven? J. Chem. Theory Comput. 2010, 6, 2866−2871. (37) Giovambattista, N.; Lopez, C. F.; Rossky, P. J.; Debenedetti, P. G. Hydrophobicity of protein surfaces: separating geometry from chemistry. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 2274−2279.

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00216. Details on potential of mean force, ligand friction, and water fluctuations; fits to the profiles; mean binding time and spatially shifted friction (PDF)

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. *E-mail: [email protected]. ORCID

R. Gregor Weiß: 0000-0002-4345-4091 Piotr Setny: 0000-0002-0769-0943 Joachim Dzubiella: 0000-0001-6751-1487 Funding

The authors thank the Deutsche Forschungsgemeinschaft (DFG) for financial support for this project. P.S. is supported by EMBO IG 3051/2015. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Chandler, D. Interfaces and the driving force of hydrophobic assembly. Nature 2005, 437, 640−647. (2) Pratt, L. R.; Chandler, D. Theory of the hydrophobic effect. J. Chem. Phys. 1977, 67, 3683−3704. (3) Lum, K.; Chandler, D.; Weeks, J. D. Hydrophobicity at small and large length scales. J. Phys. Chem. B 1999, 103, 4570−4577. (4) Djikaev, Y. S.; Ruckenstein, E. Recent developments in the theoretical, simulational, and experimental studies of the role of water hydrogen bonding in hydrophobic phenomena. Adv. Colloid Interface Sci. 2016, 235, 23−45. (5) Southall, N. T.; Dill, K. A.; Haymet, A. D. J. A view of the hydrophobic effect. J. Phys. Chem. B 2002, 106, 521−533. (6) Berne, B. J.; Weeks, J. D.; Zhou, R. Dewetting and hydrophobic interaction in physical and biological systems. Annu. Rev. Phys. Chem. 2009, 60, 85−103. (7) Nagarajan, R.; Ruckenstein, E. Theory of surfactant self-assembly: a predictive molecular thermodynamic approach. Langmuir 1991, 7, 2934−2969. (8) Setny, P.; Geller, M. Water properties inside nanoscopic hydrophobic pocket studied by computer simulations. J. Chem. Phys. 2006, 125, 144717. (9) Setny, P. Water properties and potential of mean force for hydrophobic interactions of methane and nanoscopic pockets studied by computer simulations. J. Chem. Phys. 2007, 127, 054505. (10) Setny, P. Hydrophobic interactions between methane and a nanoscopic pocket: Three dimensional distribution of potential of mean force revealed by computer simulations. J. Chem. Phys. 2008, 128, 125105. (11) Baron, R.; Setny, P.; Andrew McCammon, J. Water in cavityligand recognition. J. Am. Chem. Soc. 2010, 132, 12091−12097. (12) Young, T.; Abel, R.; Kim, B.; Berne, B. J.; Friesner, R. A. Motifs for molecular recognition exploiting hydrophobic enclosure in proteinligand binding. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 808−813. (13) Nair, S. K.; Calderone, T. L.; Christianson, D. M.; Fierke, C. A. Altering the mouth of a hydrophobic pocket - Structure and kinetics of human carbonic anhydrase-II mutants at residue VAL-121. J. Biol. Chem. 1991, 266, 17320−17325. (14) Carey, C.; Cheng, Y.-K.; Rossky, P. J. Hydration structure of the α-chymotrypsin substrate binding pocket: the impact of constrained geometry. Chem. Phys. 2000, 258, 415−425. G


Article

Journal of Chemical Theory and Computation (38) Koshland, D. E., Jr. The key-lock theory and the induced fit theory. Angew. Chem., Int. Ed. Engl. 1995, 33, 2375−2378. (39) Young, T.; Hua, L.; Huang, X.; Abel, R.; Friesner, R. A.; Berne, B. J. Dewetting transitions in protein cavities. Proteins: Struct., Funct., Genet. 2010, 78, 1856−1869. (40) Vaitheeswaran, S.; Yin, H.; Rasaiah, J. C.; Hummer, G. Water clusters in nonpolar cavities. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 17002−17005. (41) Setny, P.; Wang, Z.; Cheng, L.-T.; Li, B.; McCammon, J. A.; Dzubiella, J. Dewetting-controlled binding of ligands to hydrophobic pockets. Phys. Rev. Lett. 2009, 103, 187801. (42) Qvist, J.; Davidovic, M.; Hamelberg, D.; Halle, B. A dry ligandbinding cavity in a solvated protein. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 6296−6301. (43) Remsing, R. C.; Xi, E.; Vembanur, S.; Sharma, S.; Debenedetti, P. G.; Garde, S.; Patel, A. Pathways to dewetting in hydrophobic confinement. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 8181−8186. (44) Cheng, Y.-K.; Rossky, P. J. Surface topography dependence of biomolecular hydrophobic hydration. Nature 1998, 392, 696−699. (45) Cheng, Y.-K.; Rossky, P. J. The effect of vicinal polar and charged groups on hydrophobic hydration. Biopolymers 1999, 50, 742−750. (46) Pan, A. C.; Borhani, D. W.; Dror, R. O.; Shaw, D. E. Molecular determinants of drug-receptor binding kinetics. Drug Discovery Today 2013, 18, 667−673. (47) Swinney, D. C. The role of binding kinetics in therapeutically useful drug action. Curr. Opin. Drug Discovery Dev. 2009, 12, 31−39. (48) Copeland, R. A.; Pompliano, D. L.; Meek, T. D. Drug-target residence time and tis implications for lead optimization. Nat. Rev. Drug Discovery 2006, 5, 730−739. (49) Tiwary, P.; Mondal, J.; Morrone, J. A.; Berne, B. J. Role of water and steric constraints in the kinetics of cavity−ligand unbinding. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 12015−12019. (50) Setny, P.; Baron, R.; Kekenes-Huskey, P. M.; McCammon, J. A.; Dzubiella, J. Solvent fluctuations in hydrophobic cavity-ligand binding kinetics. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 1197−1202. (51) Mondal, J.; Morrone, J. A.; Berne, B. J. How hydrophobic drying forces impact the kinetics of molecular recognition. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 13277−13282. (52) Weiß, R. G.; Setny, P.; Dzubiella, J. Solvent fluctuations induce non-Markovian kinetics in hydrophobic pocket-ligand binding. J. Phys. Chem. B 2016, 120, 8127−8136. (53) Li, J.; Berne, B. J.; Morrone, J. A. Are Hydrodynamic Interactions Important in the Kinetics of Hydrophobic Collapse? J. Phys. Chem. B 2012, 116, 11537−11544. (54) Morrone, J. A.; Li, J.; Berne, B. J. Interplay between hydrodynamics and the free energy surface in the assembly of nanoscale hydrophobes. J. Phys. Chem. B 2012, 116, 378−389. (55) Hamaker, H. C. The London − van der Waals attraction between spherical particles. Physica 1937, 4, 1058−1072. (56) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. The weighted histogram analysis method for freeenergy calculations on biomolecules. I. The method. J. Comput. Chem. 1992, 13, 1011−1021. (57) 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. (58) 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. (59) Weiss, G. H. First Passage Time Problems in Chemical Physics; Advances in Chemical Physics; John Wiley & Sons, Inc.: Hoboken, NJ, 1996; Vol. 13, pp 1−1810.1002/9780470140154.ch1.

H


Principles for Tuning Hydrophobic LigandâReceptor Binding Kinetics

Recommend Documents

Principles for Tuning Hydrophobic LigandâReceptor Binding Kinetics