Affinity, Kinetics, and Pathways of Anisotropic Ligands Binding to Hydrophobic Model Pockets
R. Gregor Weiß, Richard Chudoba, Piotr Setny, Joachim Dzubiella
AAffinity, kinetics, and pathways of anisotropic ligands binding tohydrophobic model pockets
R. Gregor Weiß,
1, 2, ∗ Richard Chudoba,
1, 3, 4
Piotr Setny, and Joachim Dzubiella
1, 3, 4, † Institut f¨ur Physik, Humboldt-Universit¨at zu Berlin,Newtonstrasse 15, D-12489 Berlin, Germany Laborartory of Physical Chemistry,ETH Z¨urich, CH-8093 Z¨urich, Switzerland Research Group Simulations of Energy Materials,Helmholtz-Zentrum Berlin, Hahn-Meitner-Platz 1, D-14109 Berlin, Germany Physikalisches Institut, Albert-Ludwigs-Universit¨at Freiburg,Hermann-Herder Strasse 3, D-79104 Freiburg, Germany Centre of New Technologies, University of Warsaw,Stefana Banacha 2c, 02-927 Warsaw, Poland
Abstract
Using explicit-water molecular dynamics (MD) simulations of a generic pocket-ligand model weinvestigate how chemical and shape anisotropy of small ligands influences the affinities, kinetic ratesand pathways for their association to hydrophobic binding sites. In particular, we investigate aro-matic compounds, all of similar molecular size, but distinct by various hydrophilic or hydrophobicresidues. We demonstrate that the most hydrophobic sections are in general desolvated primarilyupon binding to the cavity, suggesting that specific hydration of the different chemical units cansteer the orientation pathways via a ‘hydrophobic torque’. Moreover, we find that ligands with bi-modal orientation fluctuations have significantly increased kinetic barriers for binding compared tothe kinetic barriers previously observed for spherical ligands due to translational fluctuations. Weexemplify that these kinetic barriers, which are ligand specific, impact both binding and unbindingtimes for which we observe considerable differences between our studied ligands. ∗ [email protected] † [email protected] a r X i v : . [ c ond - m a t . s o f t ] S e p olecular recognition in aqueous solution is of fundamental importance in living andchemically engineered systems. As an example, enzymes bind a substrate to an often com-plementary, concavely shaped binding site [1]. Also, receptors are activated or inhibitedif their binding pockets take up a small molecule, such as a neurotransmitter [2], a hor-mone [3] or a pharmaceutical drug [4]. Moreover, this binding principle from nature iscopied in chemical engineering of supramolecular chemistry, where so-called cavitands [5, 6]or macrocycles [7] are designed as molecular containers. The superior principle is oftenpictured by a binding agent representing a ’key’ or ’guest’ selectively fitting into a com-plementary shaped ’lock’, or ’host’, respectively, giving the names of key-lock or host-guestprinciples [8, 9]. Modern drug screening and design is based on this principle.Key-lock binding in water exhibits strong solvent-mediated effects which in fact can bediverse [10] but are very often of hydrophobic nature. By today, it is well accepted fromcomputer simulations and experimental studies that hydrophobic protein pockets comprisestrong contributions to the binding affinity of ligands through non-trivial dehydration ef-fects [11–17]. Motivated by the recent recognition of the importance of drug-receptor bindingrates for the efficacy of the drug [18], most novel studies have now added on the role of wa-ter in hydrophobic association kinetics. In particular, Setny et al. [19] documented a directcoupling between water fluctuations in hydrophobic pockets and the ligand binding rate.Using explicit-water MD simulations, they studied the binding of a spherical ligand to ahydrophobic pocket represented by a hemispherical surface recess in a model wall [20–22].They demonstrated that hydrophobically driven wet-dry fluctuations inside the pocket couldlead to locally enhanced ligand friction and kinetic barriers in the vicinity of the bindingsite. These findings were consistent with observations of Berne and coworkers in a study ona similar hydrophobic key-lock model setup [23]. In a most recent work of ours [24], we alsoshowed that increased hydrophobicity by modulating shape or water affinity of the pocketinfluences the friction peak and can speed up the binding kinetics.All previous modeling work on hydrophobic pocket-ligand binding exclusively focusedon simple spherical ligands, mimicking a methane-like molecule [15–17, 19–22, 24] or otheridealized carbon-based assemblies [23, 25, 26]. Chemical and shape anisotropy of ligands,however, are of fundamental importance for biological function [27, 28]. For instance, anal-yses of a series of chemical derivatives demonstrated that altering the ensemble of ligandbinding orientations changes signaling output, providing a novel mechanism for titrating al-2osteric signaling activity [29]. Other work found that only the orientation of the substratescorrelated with the conjugation capacity in in-vitro experiments, where the conjugationreaction proceeded only when the hydroxyl group of the ligand is oriented towards the coen-zyme [30]. Hence, the widely used spherical models of ligands is most of the time inadequatefor the identification of potential binding pockets in computational methods [31]. A moresystematic investigation of the effects of chemical and shape anisotropy of a ligand to ahydrophobic pocket on affinity as well as kinetics is therefore of high fundamental inter-est. Since complex anisotropic ligands have more degrees of freedom than a simple sphere,interesting behavior in the coupling to hydration and in the association pathways can beexpected, possibly opening more opportunites in drug design.So far, fundamental studies on solvent-influenced binding pathways of anisotropic sub-strates can be only found in more coarse-grained descriptions discussing the role of solventdepletion in molecular pair association processes. Kinoshita [32], for instance, calculated thedepletion potential of two freely rotating plates which favors an association pathway withvery tilted orientations to each other. Moreover, Roth et al. [33] estimated the associationpotential of a rod associating to a planar wall. They found that solvent depletion generatesa torque that will favor an equivalently tilted association pathway. Also in the applicationto key-lock models [34] with a spheroidal ligand the depletion forces were calculated andfound to exhibit comparably high barriers if the ligand faces the pocket with its extendedside during the association process. In essence, they conclude that an aspherical ligand willnon-parallely associate to a concave binding site and only in the last step fold/lay downinto the pocket. Dlugosz et al. [35] studied the binding time in Brownian simulations of anspheroidal ligand that electrostatically associated to a concave pocket. If the hydrodynamicinteractions were turned off, the mean binding time was faster than in calculations whichconsidered hydrodynamic interactions. Hence the water-mediation, i.e., hydrodynamics,plus the reorientation process of the aspherical ligand increased the absolute binding time.In this paper, we investigate how chemical and shape anisotropy of small ligands influencethe affinities, kinetic rates and pathways for the association to hydrophobic binding sitesusing explicit-water MD simulations of the generic pocket model used previously [15–17, 19–22, 24]. We highlight in particular how water influences the binding and the unbinding ofvarious aromatic ligands to hydrophobic binding sites, and how the new rotational degreesof freedom couple to our previously reported water fluctuation effects [19, 24, 36]. First, we3 ll atom e t h y l b e n z e n e t o l u e n e b e n z e n e ph e n o l OH b e n z y l a l c o h o l OH Figure 1: All ligands studied in this work contain an aromatic ring. The reference benzene ringin the center is alkylated with a methyl and ethyl group, respectively, while stepping to the leftcreating toluene and ethylbenzene. Stepping to the right, a hydroxyl group and a hydroxymethylmoiety are respectively introduced to form phenol and benzyl alcohol. focus on the binding of benzene to a hydrophobic planar wall compared to a hydrophobicbinding pocket and present the ligand reorientation potential, pathway, and friction profile.We elaborate on an interpretation of enhanced friction as kinetic barrier in a rescaled energylandscape which was originally introduced by Hinczewski et al. [37] in the context of proteinfolding. Moreover, we discuss how the kinetic barrier varies upon our various ligands, whichare all of similar molecular size, but distinct by various hydrophilic or hydrophobic residues.We describe that the general concept of a kinetic barrier is not only influenced by thepreviously reported water fluctuations but also by the new rotational degree of freedom ofthe ligands. We analytically discuss the impact of the resulting slowdown and see that themost significant implications must follow for the unbinding process, whereas the unbindingtime can be influenced by hundreds of microseconds. Thus the kinetic barrier can be steeredby ligand shape and, in turn, steer the binding and unbinding rates such that the action ofa toxin and the efficacy of a drug can be optimized.
I. METHODSA. Varyious non-spherical ligands
In contrast to our previous study [24] in which we modified the physicochemical prop-erties of the binding site, we investigate here the binding of various aromatic compoundsto hydrophobic binding sites. Therefore we use ethylbenzene, toluene, benzene, phenol and4enzyl alcohol, which are illustrated in Fig. 1. The aromatic compounds are represented bythe OPLS-AA force field [38] whereas we employ the LINCS algorithm to constrain all bondlengths. Phenol can be considered to be the most hydrophilic compound based on its ratioof polar to non-polar solvent-accessible surface area while hydrophobicity increases for theligands left and right from it in Fig. 1. We use the TIP4P model for water [39, 40].
B. Constrained simulations for PMF and ligand orientation
The simulation setup is illustrated in Fig. 2 and 3 which is the same as in our previousstudy [24]. The ligand binding process is constrained to one-dimensional diffusion along z , the distance of the ligand center-of-mass perpendicular to the pocketed wall shown inFig. 2 and 3. Movement of the center of mass (COM) along x - and y -direction was stronglyrestrained with a harmonic potential with spring constant k x/y = 42000 kJ mol − nm − .To probe selected observables as functions of the ligand separation to the binding site weutilized umbrella sampling simulations along z . In each umbrella setup, the center of massof the ligand was constrained to a given position z i using an external harmonic potentialwith spring constant K = 835 kJ mol − nm − . We define the origin z = 0 by the firstlayer of the wall that is in contact with the water (Fig. 2) or the pocket’s bottom, namelythe inner crystal layer of the pocket (Fig. 3). The first umbrella potential was placed at z = 3 .
14 ˚A. Up to 48 additional umbrella windows with 0.5 ˚A spacing were introducedto cover increasing ligand-pocket distances. The individual umbrella simulations produced10 ns with a step size of 2 fs, while the ligand coordinate was stored for every time stepand the water coordinates were stored every 20 fs. The potential of mean force (PMF) wasobtained by the weighted histogram analysis method (WHAM) [41, 42].The molecular snapshots in Figs. 2 and 3 also illustrate how we define the ligand orien-tation by the angle θ . It is the angle between the normal vector of the aromatic ring andthe z -axis, and consequentially the angle between the ring’s plane spanned by its ring atomsand the x - y -plane (gray), which is the parallel plane to the wall. Note that θ runs from0 (parallel to the x - y plane) to π/ x - y plane) given the molecule’s ringsymmetry and thus its otherwise degenerate orientations, such as 0 and π . We sample thedistribution of θ in each umbrella window and hence its Boltzmann inversion, the angularpotential W ( θ, z ). Note, that the distribution must be normalized by the 2sin( θ ) scaling to5alculate W ( θ, z ).The additional residue breaks the ring symmetry which we assumed for the benzene ring.Therefore, we formally replace our angular coordinate θ by a new angle φ . It defines theangle between the respective residue R and the z -axis as illustrated in the upper sketch ofFig. 4. Note that the unique mapping onto φ ranges from 0 to π , which is the necessarydescriptor range to distinguish all orientations of the ligands other than benzene. If theangle is zero the residue points into the water, away from the pocket. If the angle is equalto π the residue points into the pocket, away from the bulk. These two orientations aredegenerate for the benzene ring and the definition of θ , as mentioned above. We sample thepotential U ( φ, z ) from the orientation distributions in all umbrella windows in the same waywe obtained the potential W ( θ, z ). Note, that the distributions must be normalized by thesin( φ ) scaling to calculate U ( φ, z ). C. Unconstrained simulations for mean first passage times
For each setup, we store a production run of 20 ns in steps of 0.2 ps in which the ligand isconstrained at the reflective boundary at z = 22 ˚A for the planar binding site and z = 29 ˚Afor the pocketed binding site. This initial trajectory served as a source for randomly seededinitial configurations for subsequent binding event simulations. We then sampled more thana thousand binding events starting from independent, initial frames by randomly pickingone from the previously generated production run. To ensure a selection’s randomness, uponpossible re-selection, we applied an additional annealing step: within a short simulation of50 ps, the configuration was heated up to 350 K in a stochastic integrator scheme. After this,the heated simulation was equilibrated for 100 ps at 298 K using the Berendsen thermostatand 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 bound tothe binding site at z = 4 ˚A. The time for binding, that is, the first passage time (FPT),was then averaged to calculate the mean first passage time (MFPT) T ( z, z f ), the time theligand takes to bind from z to the final bound position z f = 4 ˚A. The curves for T ( z, z f )are discussed in the SI.In simple cases, the MFPT curve can be theoretically calculated using a Markovianapproach for diffusion in one dimension. Given an energy landscape V ( z ) such as a PMF6 ( θ ) [ k B T ] θ [rad] ˚A ˚A z [˚A]0 π π π π θ [ r a d ] ≥ W ( θ ) [ k B T ] π π π π θθ x y z z = (a)(b) Figure 2: The simulation snapshot illustrates benzene, the planar wall and part of the water. Thebinding process is constrained to one-dimensional diffusion along z . The ligand’s orientation isquantified by the angle θ between the atomic plane of the ring atoms and the x - y plane (gray).It is the same angle like the one between the normal vector of benzene and the z -axis. (a) Thedependence of the angular potential W ( θ, z ) on z illustrates the pathway upon binding. Here thetransition from a favorably perpendicular ( θ = π/
2) to a lateral ( θ = 0) orientation for decreasingligand-wall separations occurs on a narrow range from roughly z = 5 ˚A to z = 7 ˚A. Panel (b)exemplifies sampled data for W ( θ, z ) as gray symbols including blue lined fits for z = 6 ˚A and z = 5 ˚A. Strikingly the perpendicular orientation is favored by over 2 k B T . and a possibly spatially dependent friction ξ ( z ), the MFPT can be calculated by [43, 44] T ( z, z f ) = β Z zz f d z ξ ( z )e βV ( z ) Z z max z d z e − βV ( z ) . (1)where z f and z max denote an absorbing and reflective boundary, respectively. We exploitand elaborate on the framework of Eq. (1) in section B, C, and D. II. RESULTSA. Ligand reorientation
We first study the reference ligand, i.e., benzene associating to a planar hydrophobic wall.Fig. 2 (a) and (b) draw the orientation potential W ( θ, z ) from umbrella sampling. Benzeneundergoes a clear orientation pathway which depends on the ligand separation to the wall.7 ( θ ) [ k B T ] θ [rad] . ˚A z [˚A]0 π π π π θ [ r a d ] ≥ W ( θ ) [ k B T ] z m π π π π ˚A θθ x y z z = z m = 7 ˚A (a)(b) Figure 3: The simulation snapshot illustrates the benzene as ligand, a section of the hemisphericalbinding site and part of the water. The binding process is also constrained to one-dimensionaldiffusion along z . The pocket bottom is defined as origin z = 0 such that the pocket mouthis around z m = 7 ˚A. (a) The dependence of the angular potential W ( θ, z ) on z illustrates thepathway upon binding. For z >
10 ˚A a perpendicular orientation to the wall (and towards thewater interface) is energetically favored whereas for z <
10 ˚A it aligns with the pocket bottom andthe x - y -plane. Panel (b) exemplifies sampled data for W ( θ, z ) as gray symbols including blue linedfits for z = 11 . z = 7 ˚A. While the ligand is at z = 6 ˚A the perpendicular orientation for θ = π/ k B T . The example of W ( θ, z ) for z = 6 ˚A is shown inFig. 2 (b) where the gray symbols represent the simulation sampled data, and the dashedblue line is a fit of a shifted hyperbolic tangent. Thus at these ligand separations to thewall benzene partly desolvates if it orients perpendicularly to the x - y plane. Proceedingto smaller z values aligning parallel/lateral to the wall ( θ = 0) is favored because of stericrepulsion with the wall. The example data (gray symbols) and fit (solid blue line) of theangular potential for z = 5 ˚A are again shown in Fig. 2 (b).We now turn to binding to the hydrophobic pocket. In Fig. 3 (a) and (b) we see that thebenzene association to the pocketed binding site is qualitatively similar, but the reorientationfrom perpendicular to lateral occurs on a much broader range in z . While the ligand is10 − −
12 ˚A away from the pocket bottom, it favors the perpendicular orientation. It isactually only slightly favored by little more than half a k B T at z = 11 . z = 7 ˚A benzene favors the lateral orientationby more than half a k B T as shown in panel (b). Hence, benzene favorably aligns with thepocket bottom before it enters the pocket (compare Fig. 3 (a) again).We conclude that binding of benzene to our hydrophobic binding sites involves an en-ergetically favored, perpendicular orientation which is possibly disadvantageous since thefinal bound configuration requires parallel alignment. For slit-like binding sites, this couldbe advantageous if a ligand must enter perpendicularly before binding. We suppose thatthe orientation pathway is steered by water and thus solvation free energy of the ligandbecause the molecule can partly desolvate while orienting out of the water interface. Inthe case of a hemispherically molded binding site, the smeared water interface allows evenearlier desolvation and reorientation upon binding which, on top, energetically weakens theperpendicular orientation.Additional observations on the reorientation pathways cover the association of our remain-ing aromatic compounds benzyl alcohol, phenol, toluene and ethylbenzene to the pocketedbinding site only. All of these ligands comprise an aromatic ring onto which an additionalresidue is attached such as the hydroxyl group to phenol.In Fig. 4 (a) two examples are shown for U ( φ, z ) where ethylbenzene and toluene wereconstrained at z = 11 ˚A. Both exhibit a barrier around φ = 3 π/
8, whereas the barrieris smaller in the case of toluene. This behavior seems to be very significant for thesetwo aromatic compounds which comprise an alkyl residue. Hence, if these ligands are atintermediate positions they partly solvate either the aromatic ring (minimum at φ = π ) orthe alkyl group (minimum at φ = 0). Both of these orientations yield an energetic gain overan unfavored tilted orientation around φ = 3 π/
8. Nevertheless, φ = π is globally favoredbecause the aromatic ring yields higher energetic contributions from the electrostatic energyfrom its partial charges. Fig. 4 (b) and (c) show that the bimodal orientations of tolueneand ethylbenzene range from z = 12 ˚A to 9 ˚A. Overall the orientation pathway funnels alongangles that are larger than π/ z = 4 ˚A,ethylbenzene and toluene are sterically hindered to take orientations other than φ ≈ π/ U ( φ, z = 11 ˚A) for benzyl alcohol and phenol exhibit a favored orientationfor φ = 0. Both of these ligands have polar residues which contain a hydroxyl group. These9 ( φ , z ) [ k B T ] φ [rad] p h e n . b e n z y l a l c . t o l u . e t h y l b e n z . z = 11 ˚A 4 6 8 10 z [˚A]0 π π π π φ [ r a d ] z [˚A] π π π π z m ≥ U ( φ , z ) [ k B T ] z m π π π π φ R zz m (a) (b) (c)(d) (e) (f) Figure 4: The upper sketch schematically represents the pocket and ligand connected by the z -axis. Also the pocket mouth z m = 7 ˚A. The angle φ is taken to be the angle between therespective ligand’s residue R and the z -axis. Plot (a) shows two examples for the angular potentials U ( φ, z = 11 ˚A) from umbrella sampling of ethylbenzene (dashed) and toluene (solid). (The graysymbols plot the sampled data and the blue lines represent smooth interpolation functions.) Thecolor map plots show the angular potentials from all umbrella windows for (b) toluene and (c)ethylbenzene. Panel (d) also shows the examples of U ( φ, z = 11 ˚A) for benzyl alcohol (dashed)and phenol (solid). For (e) phenol and (f) benzyl alcohol we also show the full angular potentials U ( φ, z ) as color map plots. impose two potential hydrogen bonding sites: the oxygen and the associated hydrogen atom.As a consequence, the favorable solvation of the hydroxyl group yields a strong orientationto φ = 0. This angle is even more favored at closer distances z , such that the barrieraround φ = π increases to several k B T . In Fig. 4 (e) and (f) for phenol and benzyl alcohol,respectively, one can see that the energy to take an almost perpendicular orientation exceedsthe 1 . k B T plotted scale. Overall these plots make evident that benzyl alcohol and phenolfavorably solvate their hydroxyl groups while they approach the pocket. Finally both ligandsorient to φ = π/ z = 4 ˚A). In summary, the orientation of all ligandssuggests that their pathways are driven by solvation free energy such that the parts whichhave the highest energetic costs upon hydration are primarily desolvated.10 (a) (b.1) (c.1) π π π π (b.2) π π π π (c.2) h N i z [˚A] benzyl alc.phen.tolu.ethylbenz. h N ih N i φ φ Figure 5: Panel (a) shows the dewetting transition of the pocket plotting the average pocketwater occupancy against ligand distance to the pocket for benzyl alcohol, phenol, toluene, andethylbenzene. The gray vertical lines serve as the legend to the reference positions that are used asline styles in the remaining figure panels. Panel (b.1) and (b.2) plot h N i against ligand orientation φ for toluene and ethylbenzene, respectively, at the indicated positions in panel (a). Panel (c.1)and (c.2) plot the same h N i dependence for phenol and benzyl alcohol, respectively. Note, thatwhile a given ligand’s hydroxyl group orients towards the pocket, the otherwise dewetted bindingsite can be considerably hydrated. In turn, we observe how ligand position and orientation influence the pocket dewetting.Fig. 5 plots the average pocket water occupancy h N i against the two spatial descriptorsposition z and angle φ for ethylbenzene, toluene, phenol, and benzyl alcohol. The transitionfrom wet to dry occurs while a given ligand is at z = 10 ˚A to z = 15 ˚A as shown inFig. 5 (a). An average of six to seven pocket water molecules transition to less than anaverage of one while the ligand is still in front of the pocket at z = 10 ˚A. If the ligand entersthe pocket this minimum pocket water occupancy increases again up to three to four watermolecules which fit into the pocket if the ligand is bound around z = 5 ˚A. Additionally,comparing the dewetting transition for instance of ethylbenzene and phenol, we see thatmore elongated and hydrophobic ligands induce the dewetting farther away. Further, the11ertical gray lines indicate the reference positions for which we observe how h N i dependson the ligand orientation in panels (b.1), (b.2), (c.1), and (c.2) for toluene, ethylbenzene,phenol, and benzyl alcohol, respectively. In panels (b.1) and (b.2), toluene and ethylbenzenepronouncedly induce dewetting if either elongated edge, the aromatic or residue group, reachtowards the pocket and, thus, increase the hydrophobic confinement at φ = π or φ = 0. Incontrast, phenol and benzyl alcohol enhance pocket wetting if their polar residue is orientedtowards the pocket, which essentially pushes the associated solvation layers into the pocket.Hence the hydroxyl groups can considerably hydrate an otherwise dewetted pocket if theyare oriented toward it. B. Dissipative forces and kinetic barriers
Previously we discussed that the pocket water density fluctuations yield additional dissi-pative forces that slow the binding [19, 24, 36]. Firstly, pocket water occupancy fluctuationsand ligand friction were shown to couple [19]. The long time transients of the water fluc-tuations lead to long time transients in the ligand’s force correlations, and for small ligandseparations to the pocket, the water fluctuations increase due to the increased confinement.Secondly, we derived that the hydration fluctuation time scale and the ligand friction di-rectly couple by a proportional relation [36]. In this context, we found that a bimodal natureof hydration fluctuations is sufficient to enhance the ligand friction prior to association tothe pocket. And finally, we demonstrated that results from ligand constraining simulationsmust be corrected by the time transients which occur in unconstrained simulations. Onlythen one can capture the non-Markovian properties, i.e. long-time correlations, for accuratekinetic predictions. We obtain the dissipative forces and time transients (memory) in afriction profile calculated via [37] βξ M ( z ) = ∂ T ( z ) ∂z e − βV ( z ) R z max z d z e − βV ( z ) (2)where the PMF V ( z ) and the MFPT curve T ( z, z f ) are employed. Hence, we can com-bine the results from ligand constraining simulations, i.e., the PMF, and unconstrainedsimulations, i.e., the MFPT. We denote this profile by ξ M ( z ) accounting for the Markovianassumption of Eq. (2). Still, we know from our previous work that this profile non-triviallyincorporates the non-Markovian memory effects by our MFPT input. Rigorously speaking,12e shall not consider ξ M ( z ) as friction profile [19, 24, 36]. We refer to ξ M ( z ) as the kineticprofile which can be incorporated in a rescaled free energy landscape capturing all kineticeffects, whereas we follow the lines of Hinczewski et al. [37].Fig. 6 (a.1) and (a.2) show ξ M ( z ) /ξ ∞ whereas we normalize by the respective bulk frictionconstant to compare our various ligands. The values of the bulk friction values are analyzedand discussed in the SI. For example, the kinetic profile of benzene binding to the wall canbe well assumed to be constant. Moreover, if benzene binds to the pocket, the dominantfeature of the steady-state friction ξ M ( z ) is a Gaussian peaking structure, which we willmodel and fit by ξ ( z ) = ξ ∞ + ∆ ξ e − ( z − z p ) /σ . (3)Thus, the peak height ∆ ξ , position z p and width σ define a peak that adds to the bulk frictionconstant ξ ∞ . This peak roots from pocket hydration fluctuations as we also previouslydiscussed for binding of a spherical ligand to the same pocket [24] which we replot here asgray circles. The key to the additional dissipative forces are the bimodal wet-dry hydrationfluctuations which couple to the ligand. In comparison to the data for the spherical ligand,the dissipative forces for benzene peak wider and shift slightly further into the bulk byroughly half an ˚A which makes their tail reach to z ∼
11 ˚A. This well coincides with theposition where the perpendicular orientation is favored (see Fig. 3 (a)). Hence, benzene isexposed to the increasing friction at larger z -values because it reaches with its extended sidetowards the fluctuating interface.The kinetic profiles for the remaining aromatic compounds ethylbenzene, toluene, phenoland benzyl alcohol are shown in Fig. 6 (a.2) where they are compared to the replotted profileof benzene. The compounds phenol and benzyl alcohol are extended by a hydroxyl groupand a methanol group, respectively, thus offering polar patches. Their size is elongatedcompared to the benzene ring; however, their kinetic profiles match the one of benzene,i.e., their kinetic barriers well coincide. Since the orientation pathways of these two ligandsdominantly expose the aromatic ring to the pocket, the hydration fluctuations yield a similarkinetic profile.Ethylbenzene and toluene are purely hydrophobic compounds made up of a conjugatedcarbon ring that is extended by an ethyl and a methyl group, respectively. Their kineticprofiles are more enhanced and reach further into the bulk. For these two ligands, the peak iseven higher and hints that it contains contributions other than the bimodal pocket hydration.13 M ( z ) / ξ ∞ (a.1) ξ M ( z ) / ξ ∞ (a.1) (a.2) V ( z ) & V ( z ) [ k B T ] z [˚A] (b.1) V ( z ) – symbols V ( z ) – lines (b.1) s l o p e z [˚A] (b.2) V ( Q ) & V ( Q ) [ k B T ] Q [a.u.] V ( Q ) – symbols V ( Q ) – lines (c.1)(c.1) Q [a.u.] (c.2) spherebenz.-pocketbenz.-wall benzyl alc.phenoltolueneethylbenz. -16-12-8-40 4 6 8 10 12 14 4 6 8 10 12 14-16-12-8-404 4 8 12 16 20 24 28 4 8 12 16 20 24 28 Figure 6: (a.1) The kinetic profiles from Eq. (2) exhibit peaks if benzene (green squares) and thespherical ligand (gray circles) bind to the pocketed site. If benzene binds to the wall (red crosses) thefriction can be well assumed constant. The respectively colored lines represent Gaussian functionfits from Eq. (3). (a.2) Comparably peaking profiles can be observed for our remaining aromaticcompounds, whereas those of ethylbenzene and toluene are even more enhanced. Panels (b.1) and(b.2) show the original PMF V ( z ) as colored symbols. Note that the original PMFs of benzeneand the spherical ligand binding to the pocket share the similar attracting slope (blue line) whichsets in around z = 11 ˚A. Panels (c.1) and (c.2) show the rescaled energy landscapes V ( Q ) as lineswhich exhibit the additional kinetic barrier along the rescaled coordinate. We suggest that the additionally bimodally fluctuating orientation adds to the peak of thekinetic profile. In essence, binding of these two ligands involves two degrees of freedom,which bimodally fluctuate and which thus can both add to additional dissipative forces inour one-dimensional description. More importantly the peak positions for ethylbenzene andtoluene shift farther away from the pocket. In comparison to benzene these ligands are evenlonger and are subject to the hydration fluctuations farther outside the pocket.To judge the impact of the kinetic profiles we rescale it into an effective free energy land-scape. We choose a new reaction coordinate Q = Q ( z ), as suggested by Hinczewski et al. [37],such that in the new coordinates the friction is scaled to the constant value 1 k B T ns nm − .14hen the rescaled coordinate is determined by Q = d Q / d z = p ξ M ( z ) /ξ ∞ and the PMFmust be consistently rescaled such that V ( Q ( z )) = V ( z ) + ( β ) − ln( Q ( z ))= V ( z ) + (2 β ) − ln( ξ M ( z ) /ξ ∞ ) (4)In panels (c.1) and (c.2) we plot the rescaled potentials against the new coordinate Q . Therescaled energy landscapes exhibit additional kinetic barriers which naturally origin fromthe kinetic profiles. The rescaled coordinate is calculated as the integral over p ξ M ( z ) /ξ ∞ such that the integration of the Gaussian shaped peak stretches the reaction coordinate. Incomparison to the case for the spherical ligand (gray line), the peak in the kinetic barriersfor the aromatic compounds (colored lines) shift farther away from the pocket, namely toincreasing values of Q . In general, the farther the barrier shifts down the attracting slope,the smaller is its impact because the anyhow attracting slope diminishes the repulsive slopeon the r.h.s. of the barrier. In other words, part of the repulsive slope of the kinetic barrierreaches across the onset of attraction which makes its effect more significant for a slowedassociation. This result is consistent with the MFPT data which we present in the SI, wherebenzene and other aromatic compounds bind slightly slower than the spherical ligand. Thebinding times of ethylbenzene and toluene are even slower than those of benzene. Thebinding speeds of phenol and benzyl alcohol, however, are similar to the binding times ofbenzene.In sum, the size and nature of a ligand can shift and tune the dissipative forces and theresulting kinetic barriers in the ξ M ( z ) profiles. So far we only discussed this qualitativelywith a scientist’s intuition for the shapes of energy landscapes. In the following, we approachour arguments in a quantitative picture by which we explore the full range of the possibleimpact of the kinetic barrier. C. Impact of Steady-State Friction
We already formalized the kinetic profile ξ M ( z ) by its dominant feature the fitted Gaussianpeak from Eq. (3) (see also Fig. 6). The common and dominant feature of the PMF is itssignificantly attracting slope with roughly f = 13 k B T nm − (see blue line in Fig. 6 (b.1)).Thus in the following minimalistic model we use V ( z ) = f ( z − ¯ z ) for z ≤ ¯ z and V ( z ) = 015 k B T k B T k B T k B T f ¯ z (a) . . g o n ( z p , f ) [ n m ] . (b) g o ff ( z p , f ) [ n m ] z p [ ˚A] . f (c) Figure 7: (a) The dominant features of the ligand binding process are the potential V ( z ) = f ( z − ¯ z )(blue), that strongly attracts the ligand given z ≤ ¯ z , and a Gaussian friction peak (gray and green)modeled by Eq. (3). Various other potential slopes are sketched as thin blue lines. (b) The factor g on ( z p , f ) in Eq. (6) strongly depends on the friction peak position. While the negative shift z p − ¯ z decreases, g on ( z p , f ) and thus the impact of the friction peak decreases. Moreover turning to (evenslightly) repulsive slopes f drastically increases the impact of the kinetic barrier (double dotteddashed). (c) The factor g off ( z p , f ) for the unbinding process exponentially increases with frictionpeak position z p and the repulsive slope f . In the example of ethylbenzene (red curve), the potentialslope is the steepest, and the friction peak position lies farthest outside the pocket which is whythe scaling factor is on the order of O (10 ). otherwise, such that ¯ z denotes the inset position of a constant attraction with strength f .For illustration the simplified potential and friction are plotted together in the upper sketchof Fig. 7. In particular, the contribution of the friction peak, i.e., the second summand onthe r.h.s. of Eq. (3), to the binding time in Eq. (1) is given by∆ T = β ∆ ξ Z zz f d z e βV ( z ) − ( z zp )2 σ Z z max z d z e − βV ( z ) (5)16hereas the stepwise definition of V ( z ) has yet to be evaluated.Fixing ¯ z = 11 ˚A, z f = 4 ˚A,and z max = 29 ˚A, leaves the bracket in the integral in Eq. 5dependent on the friction peak position z p , width σ and the force constant f . We lay outthe detailed integration steps in the Appendix A. The result simplifies to∆ T = √ πσ β ∆ ξ · g on ( z p , f ) (6)the product of friction peak height, width and a scaling factor g on ( z p , f ). For the moment wecan neglect the dependence of g on ( z p , f ) on σ because the direct proportionality of ∆ T ∝ σ is the dominating peak width dependence for our values of σ . The factor g on ( z p , f ) quantifiesthe impact of the friction peak on the binding time which is why we will also refer to it asthe scaling factor.If we choose the slope of f = 13 k B T nm − from Fig. 6 (b.1), the scaling factor, shown asblue solid line in Fig. 7 (b), steeply increases with the peak position of the kinetic profile.The broken blue line types indicate how g on ( z p , f ) increases with decreasing potential slope f whereas the thick black line is the case for f = 0. If the force constant even becomesrepulsive the scaling factor increases drastically because the repulsive potential slope and therepulsive kinetic barrier add up. We find that the mean binding time is certainly affected byand thus proportional to the friction peak height ∆ ξ , although, the impact can drasticallydecrease if the peak shifts downward the attracting slope. For repulsive slopes the scalingcan drastically dominate such that for our case the unbinding is dominantly affected.The result of the scaling factor g off for the unbinding process is determined by interchang-ing the integration boundaries z max and z f in Eq. (5). We exemplify the scaling factor for theunbinding in Fig. 7 (c). Generally, the scaling factors exponentially scale with force constant f . Especially for our linearly attractive potential, they scale with g on ∝ exp( βf (¯ z − z f ))for the binding process and g off ∝ exp( βf ( z p − ¯ z )) for the unbinding process. Thus forunbinding it takes values much larger than one, i.e., O (10 ), where for the binding it takesvalues two orders of magnitude smaller, i.e., O (10 − ). In the special case of ethylbenzene,the slope f is the steepest and the friction peak is farthest away from the pocket. Hence,the scaling factor is on the order of O (10 ) for ethylbenzene. See also red curve in Fig. 7 (c).17 . Unbinding In this section we want to compare the impact on average binding times T on = R z max z b d z T ( z , z b )( z max − z b ) (7)where we choose the boundary to be bound as z b = 10 ˚A, and the average unbinding times T off = R z f z ub d z T ( z , z ub )( z max − z ub ) (8)where we choose the boundary to be unbound as z ub = 12 ˚A. The boundary for binding isread from the kinetic barrier peak position. The boundary for the unbinding z ub is chosensuch that the unbinding process can be considered complete by overcoming the kineticbarrier. The respective MFPT curve T ( z, z b / ub ) is calculated from Eq. (1) where the upperintegration boundary is z max for the binding case and z f for the unbinding case. For eachligand we choose two scenarios – one neglecting the kinetic barrier, thus ξ ( z ) = ξ ∞ and oneincorporating the kinetic barrier.The histogram of Fig. 8 (a) plots the ratio T won / T woon of the binding time with and withoutkinetic barrier. The barrier slows the binding time by a factor smaller than two for allligands. In contrast, the average unbinding time is dominantly affected. We estimate thatthe kinetic barrier adds an extra 221 µ s to the unbinding time of ethylbenzene and lessthan 1 µ s to that of phenol. Moreover, the ratio of the average unbinding times with and without kinetic barrier in Fig. 8 (b) yields a factor of more than five for ethylbenzene and isgenerally non-constant for our various ligands. Note, however, that the proper kinetic barrier e t h . - b e n z . t o l u e n e b e n z e n e ph e n o l b e n z . - a l c . T w o n / T w oo n e t h . - b e n z . t o l u e n e b e n z e n e ph e n o l b e n z . - a l c . T w o ff / T w oo ff (a) (b) Figure 8: (a) The ratio of the average binding time with T won and without T woon kinetic barrier exhibitsa constant impact during binding. (b) In contrast, the unbinding times dominantly increase by afactor of five in the ratio of the average unbinding times with T woff and without T wooff kinetic barrier. ξ M ( z )-profiles which we originally extractedfrom the binding process. Hence we neglect possible hysteresis effects. Nevertheless, ourprocedure is most sufficient and efficient for the conclusive interpretation of our estimatesof the unbinding times. Thus, we assume that the conclusions and implications remain thesame since, in particular, the energy landscape or binding affinity mainly steers the scalingof the kinetic barrier, while height modulations of the barriers linearly influence the averageunbinding time (compare Eq. (6) again).On the one hand, the resulting estimates for the unbinding time generally confirm thatmore (extended) hydrophobic ligands reside longer inside the pocket. On the other hand,if we neglect the additional kinetic barrier in front of the pocket, the unbinding estimatescan be wrong by several hundred microseconds. In particular, ethylbenzene would actuallystay for 270 µ s if estimated by Eq. (1) incorporating the kinetic barrier, while it would onlystay for 49 µ s if the kinetic barrier is ignored. In summary, we find a constant impact onthe binding times, while the impact on unbinding times predominantly changes. III. CONCLUSION
The ubiquitous motifs of hydrophobic and hydrophilic groups in active compounds areundoubtedly recognized in biomedical applications and the optimization of the overall effi-cacy of in vitro and in vivo systems. One increasingly appreciated aspect of the optimizationis the kinetics of association and dissociation. In this regard, solvent-mediated interactions,offer novel possibilities for control of synthetic cavitands and drug discovery.In summary, we investigated how binding site hydration influences the binding and un-binding kinetics of aspherical, i.e., aromatic, ligands. Therefore, we compared binding ofbenzene to two different binding sites: our hydrophobic pocket and a planar wall. We foundthat the benzene ring intermediately oriented such that it could maximally desolvate. Thearomatic ring took a perpendicular orientation to the binding site if it was just about to enterthe pocket. In contrast, if it was bound it favored a laterally aligned, i.e., flat, orientationto the binding site. The general rule, however, was apparent from adding the observationsof other aromatic compounds, i.e., ethylbenzene, toluene, phenol, and benzyl alcohol. Theligands would undergo a reorientation process that seemed to be driven by solvation freeenergy. Hydrophobic groups would primarily desolvate by orienting towards the water in-19erface. In the cases of ethylbenzene and toluene, the energetically favored orientation evenbecame bimodal such that two distinct orientations were locally stable. In all cases, the aro-matic ligands underwent a reorientation process in which an intermediate orientation wasorthogonal to the orientation of the final bound state. We found that concerning this or-thogonal reorientation it could be advantageous to bind to a dewetted pocket in comparisonof binding to a wall because the energetic penalties for reorientation were more moderate inthe pocket case. Additionally, the whole reorientation process was stretched over a broaderspatial range also because the pocket was strongly dewetted.These findings complement on previous, more coarse-grained, studies about the solvent-mediated depletion potentials and entropically driven torque from density functional theory(DFT) [33, 34]. These studies on the solvent-mediated association of ideal solutes revealedan association pathway which exhibited an intermediately tilted orientations of the extended,aspherical solutes – neither orthogonal nor laterally aligned. In particular, a ligand wouldthus approach a binding site with a relatively tilted orientation and then lay down into thepocket [34]. In contrast to these DFT studies, we found that for our systems a hydropho-bically driven torque orients the ligand perpendicular to the binding site, which can beconsidered ’disadvantageous’ if the bound state requires the ligand to be parallely alignedto the binding site. The perpendicular orientation could be considered ’advantageous’, ifthe ligand has to enter a narrow, slit-like, elliptical pocket. Moreover, a solvation driventorque should comprise entropic as well as enthalpic contributions in comparison to theaforementioned entropic torque. In sum, the ’solvation’ or ’hydrophobic torque’ originatingfrom solvation free energy is steered by specific chemical groups of the ligands where thespecific behavior of a respective ligand can often be well anticipated such that they implysimple design principles for steered orientation pathways.Additionally, we studied the kinetics by asking how the kinetic coupling impacts thebinding and also the unbinding times. We compared the various ligands to our previousresults of the spherical ligand from Ref. [24]. The aromatic compounds were binding slower,when we compared binding times that were normalized by the bulk friction coefficients,even though some of them exhibited a stronger binding affinity than the sphere. For thediscussion of these findings, we returned to the approach of Eq. (1) to extract the kineticprofiles via Eq. (2). By definition, the kinetic profiles reproduced the correct mean bindingtimes. Nevertheless, we could reinterpret the peaks from dissipative forces as kinetic barriers20hich scaled into new energy landscapes using a methodology depicted in Ref. [37].We rationalized the effect of the kinetic barrier regarding a scaling factor which dom-inantly depended on the barriers position relative to the anyhow attractive slope of therespective PMF. We found that the dissipative forces can have a much higher impact onunbinding times. The binding times were similarly enhanced by various kinetic barriersfor the different ligands; however, the effect on unbinding times scaled from a factor fromless than two to five. In particular, the residence time could be extended by hundreds ofmicroseconds if the theoretical estimates accounted for the kinetic barrier. This slow downwas especially pronounced for the ligands for which the orientation fluctuated bimodally infront of the pocket. In general, the additional degree of freedom of pocket hydration adds todissipative forces that are not captured inside a PMF [36]. In this context, we suggest thatother bimodal degrees of freedom can add to the kinetic barrier such that the bimodallyfluctuating orientation of ethylbenzene and toluene could also increase the effective friction.Nevertheless, the major impact of the kinetic barriers on unbinding times was steered by theslope in the energy landscapes and how far the barriers reached outside the pocket. In thatrespect, extended ligands proved to be appropriate to shift the peak away from the pocketwhile the extended side of elongated compounds orients towards the pocket and inducedsolvation fluctuations while the ligand is even farther outside.As a final notion, we highlight again that one significant model restriction is the one-dimensional treatment along the z -coordinate. In particular, Tiwary et al. [25] criticallyassessed the one-dimensional restraint in a similar MD setup, where they found that theligand least likely enters via a pathway that would include enhanced water fluctuations.This stands in line with our interpretation of the enhanced friction as kinetic barriers.Consequently, a ligand might particularly avoid a route comprising kinetic barriers whichin turn guide the possible pathways in a given system. Our model enables a purified andidealized investigation regarding mechanisms which are certainly not exclusively restrictedto this ligand-pocket setup. Similar results from other model solutes similarly infer the far-reaching consequences [45, 46]. Hence, we focus on this model system to investigate waterfluctuation driven effects which system-dependently influence association kinetics. We leavethe assessment to which extent water fluctuations play a role in a given system to futurestudies which especially deal with realistic association processes. Further, the presenceand relevance of drying transitions in free energy pathways have also been emphasized21n the context of folding and function of proteins. The kinetics in protein folding haspreviously been explored by Hinczewski et al. [37] when they introduced the aforementionedrescaling procedure. One of their main conclusions was that the importance of novel featuresin the rescaled energy landscape especially increased due to explicit water effects whichintroduced new kinetic mechanisms. Our study is oriented along these lines of a fundamentalunderstanding of how solvation impacts kinetic mechanisms while the model nature of oursetup yields the insights for fundamental relationships and controllability. Supporting information
Details on the bulk friction constants ξ ∞ , MFPT, and PMF. Acknowledgments
The authors thank the Deutsche Forschungsgemeinschaft (DFG) for financial support forthis project. P.S. is supported by EMBO IG 3051/2015.
Appendix A: Calculation of scaling factor g on Piecewise evaluation of the inner integral I ( z ) comprises a trivial case, i.e. an integralover unity, while z > ¯ z and an integral over the Boltzmann factor e − βV ( z ) while z ≤ ¯ z . Theresult thus remains piecewise defined, such that I ( z ) = βf (cid:2) e − βf ( z − ¯ z ) − (cid:3) + ( z max − ¯ z ) for z < ¯ zz max − z for z > ¯ z (A1)(Note that the piecewise definition I ( z ) would be lost if V ( z ) would have been discretizedby a Heaviside Step function.) In the next step the outer integral is∆ T ∝ Z zz f d z e − ( z − z p ) /σ e βV ( z ) × I ( z ) (A2a)= Z ¯ zz f d z e − ( z − z p ) /σ (cid:20) − e βf ( z − ¯ z ) βf + ( z max − ¯ z )e βf ( z − ¯ z ) (cid:21) (A2b)+ Z z ¯ z d z e − ( z − z p ) /σ × ( z max − z ) (A2c)where we used the piecewise definitions of V ( z ) and I ( z ) to split the integral from z f to z into an integral from z f to ¯ z and another one from ¯ z to z . Additionally the inverse22oltzmann factor e βV ( z ) in Eq. (A2a) is pulled into the square brackets in Eq. (A2b) ( andis one in Eq. (A2c)). Completing the squares, if necessary, all integrals can be related to aGaussian/Euler-Poisson integral. If we neglect O ( σ ) terms the solution reads∆ T = √ πσ β ∆ ξ n βf (cid:18) erf (cid:20) z p − z f σ (cid:21) − erf (cid:20) z p − ¯ zσ (cid:21)(cid:19) (A3a)+ (cid:18) − βf + ( z max − ¯ z ) (cid:19) e βf ( z p − ¯ z ) erf (cid:20) βf σ z p − z f σ (cid:21) (A3b) − (cid:18) − βf + ( z max − ¯ z ) (cid:19) e βf ( z p − ¯ z ) erf (cid:20) βf σ z p − ¯ zσ (cid:21) (A3c)+ ( z max − z p ) (cid:18) erf (cid:20) z p − ¯ zσ (cid:21) − erf (cid:20) z p − zσ (cid:21)(cid:19) o (A3d) [1] C. U. Kim, W. Lew, M. A. Williams, H. Liu, L. Zhang, S. Swaminathan, N. Bischofberger,M. S. Chen, D. B. Mendel, C. Y. Tai, W. G. Laver, and R. C. Stevens, J. Am. Chem. So. , 681 (1997).[2] I. Creese, D. R. Burt, and S. H. Snyder, Science , 481 (1976).[3] S. Borngraeber, M.-J. Budny, G. Chiellini, S. T. Cunha-Lima, M. Togashi, P. Webb, J. D.Baxter, T. S. Scanlan, and R. J. Fletterick, Proc. Natl. Acad. Sci. USA , 15358 (2003).[4] F. Teixeira-Clerc, B. Julien, P. Grenard, J. T. Van Nhieu, V. Deveaux, L. Li, V. Serriere-Lanneau, C. Ledent, A. Mallat, and S. Lotersztajn, Nat. Med. , 671 (2006).[5] M. B. Hillyer, C. L. D. Gibb, P. Sokkalingam, J. H. Jordan, S. E. Ioup, and B. C. Gibb,Org. Lett. , 4048 (2016).[6] P. P. Wanjari, B. C. Gibb, and H. S. Ashbaugh, Annu. Rev. Phys. Chem. , 617 (2016).[7] M. Xue, Y. Yang, X. Chi, Z. Zhang, and F. Huang, Acc. Chem. Res. , 1294 (2012).[8] D. E. Koshland, Angew. Chem., Int. Ed. Engl. , 2375 (1995).[9] J. Hu, Y. Cheng, Q. Wu, L. Zhao, and T. Xu, J. Phys. Chem. B , 10650 (2009).[10] P. Ball, Chem. Rev. , 74 (2008).[11] T. Young, R. Abel, B. Kim, B. J. Berne, and R. A. Friesner, Proc. Natl. Acad. Sci. (USA) , 808 (2007).[12] T. Young, L. Hua, X. Huang, R. Abel, R. A. Friesner, and B. J. Berne, Proteins , 1856(2010).
13] S. K. Nair, T. L. Calderone, D. M. Christianson, and C. A. Fierke, J. Biol. Chem. , 17320(1991).[14] C. Carey, Y.-K. Cheng, and P. J. Rossky, Chem. Phys. , 415 (2000).[15] P. Setny, Z. Wang, L.-T. Cheng, B. Li, J. A. McCammon, and J. Dzubiella, Phys. Rev. Lett. , 187801 (2009).[16] P. Setny, R. Baron, and J. A. McCammon, J. Chem. Theory Comput. , 2866 (2010).[17] R. Baron, P. Setny, and J. A. McCammon, J. Am. Chem. Soc. , 12091 (2010).[18] A. C. Pan, D. W. Borhani, R. O. Dror, and D. E. Shaw, Drug Dicov. Today , 667 (2013).[19] P. Setny, R. Baron, P. Kekenes-Huskey, J. A. McCammon, and J. Dzubiella,Proc. Natl. Acad. Sci. (USA) , 1197 (2013).[20] P. Setny and M. Geller, J. Chem. Phys. , 14417 (2006).[21] P. Setny, J. Chem. Phys. , 054505 (2007).[22] P. Setny, J. Chem. Phys. , 125105 (2008).[23] J. Mondal, J. A. Morrone, and B. J. , Proc. Natl. Acad. Sci. (USA) , 13277 (2013).[24] R. G. Weiß, P. Setny, and J. Dzubiella, J. Chem. Theory Comput. , 3012 (2017).[25] P. Tiwary, J. Mondal, J. A. Morrone, and B. J. Berne, Proc. Natl. Acad. Sci. (USA) ,12015 (2015).[26] P. Tiwary and B. J. Berne, J. Chem. Phys. , 054113 (2016).[27] C. Ludwig, P. J. A. Michiels, X. Wu, K. L. Kavanagh, E. Pilka, A. Jansson, U. Oppermann,and U. L. G¨unther, J. Med. Chem. , 1 (2007).[28] H.-J. Wittmann and A. Strasser, Naunyn-Schmiedeberg’s Arch. Pharmacol. , 595 (2017).[29] J. B. Bruning, A. A. Parent, G. Gil, M. Zhao, J. Nowak, M. C. Pace, C. L. Smith, P. V.Afonine, P. D. Adams, J. A. Katzenellenbogen, and K. W. Nettles, Nat. Chem. Biol. , 837(2010).[30] Y. Takaoka, M. Ohta, A. Takeuchi, K. Miura, M. Matsuo, T. Sakaeda, A. Sugano, andH. Nishio, J. Biochem. , 25 (2010).[31] L. Benkaidali, F. Andr´e, B. Maouche, P. Siregar, M. Benyettou, F. Maurel, and M. Petitjean,Bioinformatics , 792 (2013).[32] M. Kinoshito, Chem. Phys. Lett. , 47 (2004).[33] R. Roth, R. van Roij, D. Andrienko, K. R. Mecke, and S. Dietrich, Phys. Rev. Lett. (2002).
34] P.-M. K¨onig, R. Roth, and S. Dietrich, Europhys. Lett. (2008).[35] M. D lugosz and J. M. Antosiewicz, J. Phys. Chem. B , 7114 (2016).[36] R. G. Weiß, P. Setny, and J. Dzubiella, J. Phys. Chem. B , 8127 (2016).[37] M. Hinczewski, Y. von Hansen, J. Dzubiella, and R. R. Netz, J. Chem. Phys. , 245103(2010).[38] W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, J. Am. Chem. Soc , 11225 (1996).[39] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem.Phys. , 926 (1983).[40] W. L. Jorgensen and J. D. Madura, Mol. Phys. , 1381 (1985).[41] S. Kumar, D. Bouzida, R. H. Swendsen, P. A. Kollman, and J. M. Rosenberg, J. Com-put. Chem. , 1011 (1992).[42] F. Zhu and G. Hummer, J. Comput. Chem. , 453 (2012).[43] G. H. Weiss, Adv. Chem. Phys. , 1 (1967).[44] A. J. F. Siegert, Phys. Rev. , 617 (1951).[45] J. A. Morrone, J. Li, and B. J. Berne, J. Phys. Chem. B , 378 (2012).[46] J. Li, J. A. Morrone, and B. J. Berne, J. Phys. Chem. B , 11537 (2012). upporting information for: Affinity, kinetics, and pathways ofanisotropic ligands binding to hydrophobic model pockets R. Gregor Weiß,
1, 2, ∗ Richard Chudoba,
1, 2
Piotr Setny, and Joachim Dzubiella
1, 2, ∗ Institut f¨ur Physik, Humboldt-Universit¨at zu Berlin,Newtonstr. 15, D-12489 Berlin, Germany Institut f¨ur Weiche Materie und Funktionale Materialen,Helmholtz-Zentrum Berlin, Hahn-Meitner-Platz 1, D-14109 Berlin, Germany Centre of New Technologies, University of Warsaw, 00-927 Warsaw, Poland ∗ To whom correspondence should be addressed. E-mail: [email protected] [email protected] a r X i v : . [ c ond - m a t . s o f t ] S e p . BULK FRICTION CONSTANTS OF THE AROMATIC LIGANDS Here we briefly summarize the bulk friction constants of aromatic compounds. In generalwe can probe the bulk friction in umbrella windows that are far away from the binding siteby the position autocorrelation [1] βξ ∞ = R ∞ h δz ( t ) δz (0) i d t h δz i (1)where h δz ( t ) δz (0) i is the position auto-correlation function with δz ( t ) = z ( t ) − h z i . For thespherical ligand we know that the bulk friction is around βξ ∞ = 0 . − from our previ-ous work [2]. According to Stokes friction the bulk friction should increase with ligand sizewhich holds for our aromatic compounds. Table S1 lists the bulk friction for ethylbenzene,toluene, benzene, phenol, benzyl alcohol, and the spherical ligand. Throughout the paperwe used these values to normalize the mean binding time and kinetic profiles. This enabledus to compare the effective differences between our various ligands. II. MFPT OF THE AROMATIC LIGANDS
In Fig. S1 we plot the MFPT curves for all aromatic and the spherical ligand. The upperpanels (a) and (b) show the raw data T ( z, z f ) which makes evident that, e.g., ethylbenzenebinds slowest, however, this is a trivial result because the bulk friction constant of ethyl-benzene is the largest. If we rescale the MFPT to the respective bulk friction constants Table S1: The different aromatic compounds have different bulk friction values βξ ∞ . In general,the more elongated/bigger ligands have larger friction constants. βξ ∞ [ns nm − ]ethylbenzene 0 . . . . . . [ n s ] (a.1) T / β ξ ∞ [ n m ] z [ ˚A](b.1) (a.2) z [ ˚A](b.2) V ( z ) [ k B T ] z [ ˚A](c.1) z [ ˚A](c.2)0.00.51.01.52.02.50.00.51.01.52.02.53.0 5 10 15 20 25 5 10 15 20 25 spherebenz.-pocketbenz.-wall -16-12-8-40 4 6 8 10 12 14 benzyl alc.phenoltolueneethylbenz. Figure S1: Panel (a.1) shows the MFPT curves T for benzene binding to the wall and to the pocketas well as binding of the spherical ligand to the pocket. Panel (a.2) compares the MFPT curvesfor benzyl alcohol, phenol, benzene, toluene, and ethylbenzene. In agreement with Table S1 thebigger ligands with larger friction values bind overall slower. In panels (b.1) and (b.2) we plot theMFPT results which are rescaled by the respective bulk friction value such that T /βξ ∞ . Panels(c.1) and (c.2) plot all PMFs for all pairs of ligands and binding sites. the curves for benzene, phenol and benzyl alcohol mostly coincide in panel (d). Hence theoverall rescaled kinetics are in good agreement. This is consistent with the overlapping ki-netic barriers for these ligands in Fig. 5 in the main text. The rescaled kinetics T /βξ ∞ ofethylbenzene and toluene are, however, slightly decelerated in comparison to the benzenebinding times. Thus here the deviations of the kinetic barrier in Fig. 5 qualitatively agreeswith the normalized MFPT curves here. Moreover benzene clearly binds slower than thespherical ligand because its kinetic barrier is shifted further away from the pocket. [1] G. Hummer, New J. Phys. , 34 (2005).[2] R. G. Weiß, P. Setny, and J. Dzubiella, J. Chem. Theory Comput. , 3012 (2017)., 3012 (2017).