pubs.acs.org/JPCL
Direct Measurement of the Kinetics and Thermodynamics of Association of Hydrophobic Molecules from Molecular Dynamics Simulations Andrew S. Thomas and Adrian H. Elcock* Department of Biochemistry, University of Iowa, 51 Newton Road, Iowa City, Iowa 52242, United States
ABSTRACT A detailed study of the kinetics and thermodynamics of associations of model hydrophobic molecules is likely to be valuable for understanding the fundamental driving forces for processes such as protein folding and proteinprotein association. To this end, we present results from a series of 500 ns long molecular dynamics (MD) simulations examining associations of 13 types of different alkane pairs in explicit water. In addition to providing accurate measurements of the association thermodynamics, the unbiased nature of the configurational sampling in the MD simulations allows the association and dissociation kinetics to be directly quantified. We show that by choosing a suitable reaction coordinate, the computed free energies of all of the alkane-alkane complexes can be linearly related to their buried molecular surface areas, that their dissociation kinetics can be reliably estimated from the height of the barrier on the computed free energy surfaces, and that their association kinetics are effectively diffusionlimited. SECTION Biophysical Chemistry
required. Zhang and McCammon15 provided a framework for extracting both the thermodynamics and kinetics from unbiased molecular dynamics (MD) simulations and illustrated the approach by application to the methane-methane association in water. Here we show that a thorough sampling of more complex molecular associations with unbiased MD simulations is also now achievable, enabling a detailed examination of the connection between the thermodynamics and kinetics of hydrophobic associations and their dependence on geometric features of the interacting solutes. We have studied the bimolecular associations of the following pairs of alkanes in explicit solvent using unbiased 500 ns MD simulations: methane-methane (M-M), methane-ethane (M-E), methane-propane (M-P), methane-butane (M-B), methane-pentane (M-PNT), ethane-ethane (E-E), ethanepropane (E-P), ethane-butane (E-B), propane-propane (P-P), propane-butane (P-B), butane-butane (B-B), neopentane-neopentane (N-N), and n-pentane-n-pentane (PNT-PNT). The computed interaction free energies for all 13 pairs of alkanes are shown in Figure 1a,b, plotted versus two different definitions of reaction coordinate. The data shown in Figure 1a are plotted versus the distance between the centers of mass of the two associating molecules; this is perhaps the most natural and certainly the most commonly used reaction
M
olecular simulations have been used for many years to study the association thermodynamics of nonpolar molecules in water and to gain thereby an understanding of hydrophobic interactions. Early studies focused on hydrophobic solutes that could be considered essentially spherical (e.g., methane) and whose association could therefore be adequately described by a 1-D interaction free energy surface or potential of mean force (PMF).1-6 In more recent years, the association thermodynamics of larger nonspherical hydrophobic molecules (e.g., butane and amino acid analogs) have been computed, with the work of the Scheraga group being especially broad in scope.7-11 A common feature in many of these previous studies has been the use of biasing potentials (e.g., those used in umbrellasampling)12 that restrict sampling so that the thermodynamics of association are measured only along selected pathways such as approaches parallel or perpendicular to predefined molecular geometries.8-11 Whereas the use of biasing potentials has the advantage that the computationally expensive problem of sampling the entire configurational space open to the associating molecules can be largely avoided, a drawback is that it complicates the extraction of kinetic information, such as dissociation rate constants, from the simulations; a number of laboratories have, therefore, explored ways to derive indirectly kinetic information from simulations,13,14 but for many situations, current computer power is now sufficient that kinetic information can often be obtained directly from unbiased (“brute force”) sampling, with a minimum of assumptions being
r 2010 American Chemical Society
Received Date: November 2, 2010 Accepted Date: December 7, 2010 Published on Web Date: December 13, 2010
19
DOI: 10.1021/jz1014899 |J. Phys. Chem. Lett. 2011, 2, 19–24
pubs.acs.org/JPCL
Figure 1. Association thermodynamics of the 13 alkane pairs. Free energy curves for the association of the alkane pairs are plotted using the same simulation data but with different reaction coordinates: (a) the center-of-mass distance and (b) the closest carbon-carbon distance. The inset of panel b highlights the region of the contact minimum. (c) Correlation between the free energy of the contact minimum (from panel b) and the average amount of solvent-accessible surface area buried upon association.
coordinate in the literature (e.g., refs 8-11). The data shown in Figure 1b are instead plotted versus the shortest carboncarbon distance measurable between the two associating molecules; for two molecules that contain N and M carbon atoms, respectively, there are N M such distances to measure for every sampled configuration of the system. (See the Computational Details.) Both reaction coordinates yield free energy curves that exhibit the usual features of molecular associations in aqueous solution: (1) a contact minimum, (2) a desolvation barrier, and (3) a solvent-separated minimum. Comparison of the two panels, however, immediately suggests that the use of the closest carbon-carbon distance as the reaction coordinate leads to a more readily interpreted collection of free energy curves. As described below, this idea receives further support when attempts are made to interpret the computed thermodynamics and kinetics in terms of geometric properties of the associating molecules. The free energies of the contact minima shown in Figure 1b (all of which occur at a minimum carbon-carbon distance of ∼3.9 Å) show the expected size dependence: the smallest alkane pair, methane-methane, shows the weakest interaction at contact (ΔG = -0.74 kcal/mol), whereas the two largest alkane pairs studied, those of n-pentane (PNT-PNT) and neopentane (N-N), show the strongest interactions (ΔG = -1.31 and -1.25 kcal/mol, respectively). This strengthening of the association with increasing alkane size is consistent with the Scheraga group's measurements of PMFs for the association of alkanes in specified relative orientations.8,9 In conjunction with the strengthening of the contact minimum with increasing alkane size, however, there is also a clear increase in the height of the desolvation barrier (which occurs at a closest carbon-carbon distance of ∼5.9 Å). Notably, neither of these trends can be readily observed when the center-of-mass distance is used as the reaction coordinate because the distances at which these features appear on the free energy surfaces are so variable (Figure 1a). One very commonly used method to compute rapidly the free energy contribution made by hydrophobic interactions to biomolecular associations is to measure the hydrophobic surface area (SA) that is buried upon association. Such a method relies on the assumption that the free energy
r 2010 American Chemical Society
contribution is linearly dependent on the buried surface area, that is, that ΔG = γ ΔSA, where γ is a proportionality constant.16-18 The results compiled in Figure 1b provide us with an important new opportunity to test directly this assumption in a well-defined molecular system. Figure 1c shows the free energies of the contact minima shown in Figure 1b plotted versus the average amount of solventaccessible surface area buried, ΔSASA, upon association of each of the 13 alkane-pairs studied. A high quality linear fit is obtained (r2 = 0.975) with a slope, γ=19.3 cal/mol/Å2; interestingly, this regressed value is in good agreement both with γ values derived from surface-area based studies of association thermodynamics19 and with γ values derived from studies of hydrophobic hydration. (See ref 20 and refs therein.) Also, notable from Figure 1c is the fact that even with the largest alkane pairs studied here we observe no significant deviation from a linear dependence on the ΔSASA. Others have shown convincingly that a fundamental shift in the nature of the hydrophobic hydration occurs once a nonpolar solute's radius exceeds a threshold value,21 with estimates of this threshold ranging from a radius of ∼2.5 Å (roughly that of neopentane) to 4 Å;21-23 our data do not argue against this concept but do suggest that under the present simulation conditions the pentane isomers are not large enough to have “crossed over” into this large-size regime. It is encouraging that such a simple surface-area-based model can so accurately reproduce the interaction free energies of the alkane pairs studied here, and it argues in favor of the use of this approach in situations where rapidly estimating the free energy contribution of hydrophobic burial is paramount. It should be noted, however, as has been shown by others previously,24,25 that a ΔSASA-based approach cannot reproduce all features of the free energy curves shown in Figure 1b because it does not predict the thermodynamically unfavorable desolvation barrier (data not shown). In addition to yielding thermodynamic information, the unbiased MD simulations reported here allow new insight into the kinetics of hydrophobic interactions in solution to be obtained. Because the dissociation of alkane pairs can be considered to be effectively a unimolecular process, we expect the “survival function” of the associated complex15
20
DOI: 10.1021/jz1014899 |J. Phys. Chem. Lett. 2011, 2, 19–24
pubs.acs.org/JPCL
a linear relationship should be obtained when the directly computed dissociation rate constant is plotted versus exp(-ΔG‡dissociation/RT). Figure 2b shows the results of such a comparison when ΔG‡dissociation is calculated from the interaction free energy profile shown in Figure 1a, that is, using the center-of-mass distance as a reaction coordinate. Figure 2c shows the plot obtained when ΔG‡dissociation is instead calculated from the interaction free energy profile shown in Figure 1b, that is, using the closest carbon-carbon distance as a reaction coordinate; in this case, ΔG‡dissociation is calculated as the difference between the free energy at the desolvation barrier (∼5.9 Å) and the free energy at the contact minimum (∼3.9 Å). Clearly, only the latter plot is linear (r2 = 0.991). This indicates that the relative dissociation kinetics of the alkane complexes can be nicely described by the Arrhenius equation and that they can therefore be predicted adequately from purely thermodynamic calculations but only when the more appropriate closest-carbon reaction coordinate is used. It is worth noting, however, that, as suggested by a reviewer, the center-of-mass reaction coordinate can nevertheless be used to predict the dissociation kinetics using an Arrhenius-like expression if instead of using the free energy difference between the desolvation barrier and the contact minimum one uses the free energy difference between the associated and dissociated states obtained by numerically integrating the free energy curves15,26,27 shown in Figure 1a. (See Figure S1 in the Supporting Information.) Given the preceding paragraph and given that the association thermodynamics were previously shown to be simply related to the ΔSASA of the alkane-alkane complexes, it is perhaps not surprising to find that the directly computed dissociation rate constants are also straightforwardly related to the ΔSASA of the complexes (Figure 2d). In this case, the relationship is an exponential, which is as expected given the Arrhenius-like behavior demonstrated above and the linear relationship previously noted between ΔG and ΔSASA. Finally we turn to the association kinetics. As noted by Zhang and McCammon, the calculation of association rate constants from MD simulation data is less straightforward, and a number of approaches might be used.15 Here we have opted to define the unbound and bound states as having closest carbon-carbon distances of >10 and
