w-REXAMD: A Hamiltonian Replica Exchange Approach to Improve Free Energy Calculations for Systems with Kinetically Trapped Conformations.

Free energy governs the equilibrium extent of many biological processes. High barriers separating free energy minima often limit the sampling in molecular dynamics (MD) simulations, leading to inaccurate free energies. Here, we demonstrate enhanced sampling and improved free energy calculations, relative to conventional MD, using windowed accelerated MD within a Hamiltonian replica exchange framework (w-REXAMD). We show that for a case in which multiple conformations are separated by large free energy barriers, w-REXAMD is a useful enhanced sampling technique, without any necessary reweighting.

Introduction

The change in the Gibbs free energy is the thermodynamic property that drives equilibrium aspects of complex biological processes like the folding of proteins, drug binding, and molecular recognition of proteins in cell signaling pathways.[1−3] A great deal of progress has been made in the area of free energy calculations in recent years.[4−8]

One of the obstacles in using computer simulation to calculate free energy changes, for example in the solvation of small molecules like ibuprofen, or in comparing protein–drug binding interactions, is the limit of conformational sampling.[9,10] A transformation of interest, perhaps changing a substitutent on a drug-like molecule, or a mutation in a protein active site, may involve a conformation change that is inaccessible on the time scale of the simulation. Without considering the mutation in both possible conformations, the calculated free energy change for the process is inaccurate. This sampling problem can be particularly difficult to solve if the two (or more) conformations are not known beforehand.

Enhanced sampling techniques are often leveraged to mitigate the sampling problem.[11−13] Accelerated molecular dynamics (aMD)[14] is a Hamiltonian-modifying technique that has been used extensively in our group to simulate micro- to millisecond time scale conformational changes in biomolecules.[15−23] The original form of aMD adds a bias to potential energy minima, making barriers easier to cross, without defining a reaction coordinate a priori. While this method has been shown to be a useful tool in exploring the conformational space of large proteins, the simulation must be reweighted to recover the Boltzmann distribution. In this exponential reweighting, the conformations at energy minima make significantly larger contributions than all other conformations to the average value, and the effective number of data points is reduced, introducing statistical noise in ensemble averages. By using the original aMD equation in Hamiltonian replica exchange (REXAMD), Fajer et al. showed that the exponential reweighting problem could be ameliorated by simulating replicas at various states of acceleration, including an unaccelerated “ground” state, and collecting data only from the unaccelerated state.[24]

Recently a new “windowed” aMD (w-aMD) method has been introduced in which the potential energy barriers are scaled down while preserving the original features of the landscape (eqs 1 and 2).[25] The advantage of this implementation relative to the original aMD implementation is that the structures already at energy minima remain on the unbiased potential energy surface, and energy barriers are scaled without affecting the overall shape of the energy landscape. While this method was shown to be an improvement on the previous barrier-lowering aMD method,[26] one limitation of the windowed aMD method is the introduction of additional parameters that determine the degree of acceleration: E1 and E2 define the window in which the dynamics are accelerated, while α1 and α2 are tuning parameters (see Supporting Information for details). Here, we implement the windowed aMD equation in a Hamiltonian replica exchange framework (w-REXAMD), in which multiple sets of acceleration parameters distinguish the Hamiltonians of exchanging replicas.where

In our discussion of w-REXAMD simulations, a “state” refers to a particular set of aMD parameters introduced in eq 2. A “replica” refers to an evolving configuration which may exchange its state with that of its neighbor’s, according to the Metropolis criterion. In an analogous temperature replica exchange simulation, the “state” would be defined by the temperature, while a “replica” may visit different temperatures.

We anticipate that the additional conformational sampling from multiple replicas at different acceleration states will enrich the sampling in the unaccelerated state ensemble. As with any replica exchange method though, we must acknowledge the increase in computational cost associated with the multiple replicas. To account for both the advantage in sampling simply from multiple independent simulations and the limitation of increased computational cost, we compare w-REXAMD to conventional MD in the same replica exchange framework, which we refer to as REGREX. In a REGREX simulation, multiple independent conventional MD replicas exchange their configurations at the same frequency at which exchanges are attempted in a w-REXAMD simulation. The states in a REGREX simulation are equivalent to each other, and each state contains equivalent portions of independent MD trajectories. One REGREX state, then, may also benefit from conformational space sampled in multiple replica trajectories but will still contain the same number of structures as one w-REXAMD state.

We first compare the extent of conformational sampling in w-REXAMD and REGREX simulations of two systems: butane and cyclohexane, both in explicit water solvent. For the case of butane, conventional MD, and thus REGREX simulations, can exhaustively sample the easily characterizable conformational space of the system (Figure 1A) and provide a reference for a converged result. In the case of cyclohexane, however, a large (≈10 kcal/mol) barrier separates the two stable chair conformations (Figure 2A), making the interconversion between the two conformations, so-called “chair-flipping,” a rare event that occurs on the microsecond time scale.[27] Since cyclohexane is a model system in which the time scales are no longer accessible for reference REGREX simulations, we instead use umbrella sampling calculations as an appropriate reference to which we compare ensembles generated in w-REXAMD simulations.

Figure 1

Conformational sampling of butane. (A) A schematic of butane conformational landscape according to the C1–C2–C3–C4 dihedral angle (the cis and eclipsed (+) structures are shifted slightly so that the opposed groups are visible). Fractions of state ensembles in the trans (black series), gauche(−) (red series), gauche(+) (blue series), and eclipsed (green series) are shown for (B) REGREX, (C) w-REXAMD, and (D) w-aMD simulations. Total number of frames corresponds to 500 ps of simulation. Error bars show standard deviation from three independent simulations.

Figure 2

Conformational sampling of cyclohexane. (A) Umbrella sampling along the ε coordinate shows a high barrier between chair conformations. (B) Observed probability of chair and boat-like conformations from w-REXAMD (black series), windowed-aMD (red series), and REGREX (magenta series) simulations. Error bars show standard deviation from three independent simulations. Green bars mark the expected probability from umbrella sampling calculation in A as a reference. (C) ε conformation coordinate during windowed-aMD (red series), w-REXAMD (black series), and REGREX (magenta series) shows interconversion of cyclohexane conformations.

We then present an application of w-REXAMD to free energy calculations using an alchemical identity transformation. As mentioned earlier, multiple kinetically trapped conformations can be a major obstacle in free energy calculations. Here, we discuss a simple example of such asymmetry in the transformation of bromocyclohexane to itself. This identity transformation, in which the bromine is displaced from one carbon to another (shown in Figure 3), should result in a free energy change of zero and serves as a force-field independent benchmark. With the bromine substituent, the two chair conformations of the cyclohexane are no longer equivalent. The bromine occupies either an equatorial or axial position, and interconversion of these two conformers occurs on the microsecond time scale. In our alchemical self-transformation, the bromine in equatorial-substituted bromocyclohexane is displaced to the axial position at another carbon. The axial and equatorial conformers exist in equilibrium, so constructing the transformation in this way requires sampling of both conformers according to their thermodynamic probabilities in order to reach the analytical result of zero. On the typical nanosecond time scale of MD simulations, the axial and equatorial configurations are kinetically trapped, and the identity transformation is a good test case for the convergence problem encountered when multiple conformations are separated by high barriers. We show that the enhanced sampling in w-REXAMD simulations is useful in mitigating this convergence issue in the free energy calculation.

Figure 3

Thermodynamic cycle of alchemical identity transformation of bromocyclohexane. In alchemical transformations 1 and 2, the bromine substituent is displaced from C1 to C4, either from an equatorial to an axial position (transformation 1: “C1,eq-to-C4,ax”) or vice versa (transformation 2: “C1,ax-to-C4,eq”). Transformations 3 and 4 represent the slow conformational change isolating the two conformers in conventional MD. ΔG1 + ΔG4 = ΔG3 + ΔG2 = 0 kcal/mol.

Results and Discussion

Conformational Sampling

We distinguish between state trajectories (corresponding to the ensemble of structures sampled with a certain set of acceleration parameters) and replica trajectories (corresponding to the trajectory of a replica as it visited different states).

We use the dihedral angle (C1–C2–C3–C4) of butane during the simulations to monitor the interconversion between gauche(∓) and trans conformations of butane (Figure 1A). In order to probe convergence in conformational sampling in butane, we clustered multiple independent unaccelerated state trajectories into gauche(+), gauche(−), or trans conformations and plot the cumulative fraction for each conformation during the time scale of the w-REXAMD, REGREX, and w-aMD simulations. In the w-REXAMD simulations, all conformation populations converge to fractional occupancies that lie within the error intervals reached in reference REGREX simulations (Figure 1B and C). The two equivalent gauche conformations quickly converge in the w-REXAMD simulations in considerably fewer steps than in the REGREX simulations. For the nonexchanging w-aMD simulations, the equivalent gauche conformations are sampled exhaustively (Figure 1D), but the less favorable eclipsed conformation is sampled to a greater extent, indicating the effectiveness of the windowed aMD equation in scaling down the barriers but illustrating the need for Boltzmann reweighting in order to attain information about the equilibrium ensemble.

Unlike those of butane, the two stable conformations of cyclohexane are separated by a 10 kcal/mol free energy barrier, making interconversion between the two conformations a rare event.[27] Using a conformation coordinate ε (eq 3), which is a function of the six correlated cyclohexane dihedral angles (τ), we analyzed the sampled conformations of cyclohexane (chair1, boat-like, and chair2).[17]

By clustering the ε coordinate over multiple trajectories, we were able to observe the fractional occupancies of the two chair conformations (ε ≈ ± 50) of cyclohexane. Chair-flipping cannot be observed on reasonable time scales of cMD (0.5 μs, data not shown), and REGREX simulations similarly only sample the initial chair conformation (Figure 2B, magenta series). The probability of each chair conformation for the unaccelerated state from w-REXAMD simulations agrees well with expected probability from reference umbrella calculations (Figure 2A and B, black series). The use of the w-aMD equation allows for enhanced sampling of cyclohexane conformations (Figure 2C, red series). Without the replica exchange framework, w-aMD simulations must be reweighted to obtain the same equilibrium populations (Figure 2B, red series). Although this reweighting is quite successful in recovering equilibrium populations for small systems such as these,[25] it has been shown to introduce statistical noise in ensemble averages.[28] This issue may be more prominent in more complex systems.

Alchemical Free Energy Calculations

In the previous section, we showed that w-REXAMD is a reasonable method to rapidly explore conformational space, particularly in the presence of high barriers separating relevant conformations. We expect, then, that w-REXAMD should also improve the convergence of free energy calculations. We consider the transformation of bromocyclohexane to itself, which has a free energy change of zero. Similar to cyclohexane, bromocyclohexane has a considerable barrier separating its two chair conformations, making sampling of the conformations a microsecond time scale event, presently inaccessible to routine alchemical free energy calculations. The presence of the bromine substituent creates a difference between the two chair conformations; one chair conformation is the axial conformer, while the other is the equatorial conformer. The equatorial conformer is thermodynamically favored, as it results in less steric hindrance between the bulky bromine and neighboring hydrogens, but bromocyclohexane exists in equilibrium between both conformers. To highlight the sampling problem, we let the two conformers be the initial and final states and considered the bromocyclohexane identity transformation, where the bromine is displaced from one carbon to another, in both directions: C1,eq-to-C4,ax, and C1,ax-to-C4,eq. In the C1,eq-to-C4,ax transformation, for example, the equatorial bromine at carbon 1 and an axial hydrogen at carbon 4 are decoupled from the rest of the system, while a bromine atom appears at the axial position on carbon 4 and a hydrogen atom appears at the equatorial position initially occupied by the bromine, at carbon 1 (Figure 3).

Using the REGREX simulations at each of nine lambda intermediates along the alchemical pathway, we found that the observed free energy change for the transformation did not result in zero, and that it depended upon the directionality of the transformation (see Figure 4A). The REGREX C1,ax-to-C4,eq bromocyclohexane transformations converged to a value of approximately −0.2 kcal/mol (Figure 4A, red series), while the C1,eq-to-C4,ax self-transformations converged to approximately +0.3 kcal/mol (Figure 4A, blue series). Because of the presence of the high free energy barrier separating the two conformers, bromocyclohexane did not flip into its equatorial conformation during the alchemical transformation and instead remained exclusively in the axial conformation (see Figure 4B, blue series). In this case, where the relevant conformations are kinetically trapped, we find that the alchemical identity transformation results in a nonzero free energy change. Without any enhanced sampling, the alchemical self-transformations only complete either transformation 1 or 2 of the thermodynamic cycle in Figure 3, with ΔG1 = −ΔG2, within statistical uncertainty. In order for the analytical result to be reached with conventional MD, the conformational changes depicted by transformations 3 and 4 of the thermodynamic cycle in Figure 3 must also be completed. On the other hand, with the w-REXAMD method the initial and final states are able to sample both axial and equatorial conformations, and as a result the identity transformation quickly converges to the analytical result of 0 kcal/mol (Figure 4A, black series). It is important to note that for this small system, reweighting of the most accelerated w-aMD state trajectory also converges to the analytical result (Figure 4A, maroon series).

Figure 4

Free energy of identity transformation for bromocyclohexane. (A) REGREX transformations were done in both the C1,eq-to-C4,ax (blue series) and C1,ax-to-C4,eq (red series) directions. w-REXAMD transformations were performed for just the C1,eq-to-C4,ax direction. The free energy from the unaccelerated state of the w-REXAMD simulations is shown in the black series, and the reweighted most accelerated state is shown in maroon. Each series shows the average of three independent transformations, weighted by their statistical uncertainties. (B) Representative values for ε at λ = 0.7 for the transformations show that only one chair conformation is sampled in REGREX C1,eq-to-C4,ax (blue series) and C1,ax-to-C4,eq (red series) transformations, whereas both chair conformations are sampled in the w-REXAMD unaccelerated state of the C1,eq-to-C4,ax transformation (black series).

Conclusions

We presented an enhanced sampling method, w-REXAMD, which is capable of sampling conformations that would normally interchange on microsecond time scales. We showed that w-REXAMD could be used to improve alchemical free energy calculations, particularly when the conformation of the initial or final state is kinetically trapped. Similar to the use of accelerated molecular dynamics, an advantage of the w-REXAMD method is the gain in sampling without having to define any reaction coordinate a priori. The replica exchange framework offers the additional advantage of (1) avoiding any reweighting of data in obtaining equilibrium ensemble averages relevant for free energy calculations and (2) using multiple sets of acceleration parameters. To optimize the w-REXAMD method, future studies should focus on the distribution of the sets of acceleration parameters in order to both maximize sampling and retain efficient diffusion rates among replicas.

The replica exchange framework does come with the need for additional parallel computing resources, and this need will be amplified as system size increases. In systems with multiple, slow-exchanging conformational states though, a common protocol is to calculate the change in free energy for some transformation, e.g., a mutation in an enzyme active site, for each constrained conformation separately, and then combine the transformations in a thermodynamic cycle to estimate the actual free energy of the mutation. We believe that in such cases, and particularly in cases where the conformational landscapes are not well characterized, w-REXAMD will be a useful enhanced sampling method.

Computational Methods

Hamiltonian Replica Exchange

Each thermodynamic state in the w-REXAMD simulation corresponds to a set of four parameters (E1, E2, α1, and α2) that determine the degree of acceleration of a replica at that state. Four replicas were used in each simulation, and the least accelerated state was always the unmodified Hamiltonian. Every 300 MD steps (0.3 ps), exchanges between adjacent replicas were attempted. A description of how parameters were selected can be found in the Supporting Information. Only torsional potential energy terms were accelerated in these test cases. The Amber ff99SB force field[29] was used for all simulations, with the TIP3P water model.[30]

Thermodynamic Integration

All atoms being coupled or decoupled according to the parameter λ were treated with soft-core potentials. Nine evenly spaced λ points were chosen between 0 and 1 for each transformation, and 1.5 ns of REGREX and w-REXAMD was run at each λ point (NVT) after an equilibration period of 150 ps (NPT). ⟨dV/dλ⟩λ from one state (the unaccelerated state in w-REXAMD simulations) was used for the free energy calculations. Free energies were calculated by numerically integrating the ⟨dV/dλ⟩λ values from 0 to 1 using the Trapezoid rule. The uncertainty associated with each of the nine ⟨dV/dλ⟩λ values was calculated aswhere N is the number of dV/dλ values used to compute the average and g is the statistical inefficiency, as defined by Chodera et al.,[31] which contains the integrated autocorrelation time τ (g = 1 + 2τ). The variance of the free energy from one transformation was estimated by propagating the uncertainty at each of the nine lambda values. Three independent C1,ax-to-C4,eq and C1,eq-to-C4,ax transformations were performed using REGREX, and three C1,eq-to-C4,ax transformations were done using w-REXAMD. The final free energies reported are averages weighted by each transformation’s uncertainty:where w = 1/δ2.

Umbrella Sampling

The umbrella sampling simulations were performed in NAMD[32] using the Amber ff99SB[29] force field. At each of 63 evenly spaced ε centers, ranging from ε = −62° to ε = 62°, a 2 ns MD simulation was performed with a harmonic potential centered at the value of ε. A force constant of 1.0 kcal/mol/deg2 was used for all simulations, and only the last 1.5 ns were used to construct the free energy profile. Initial structures were taken from pre-equilibrated structures from cyclohexane NVT simulations. The free energy profile was constructed using the python implementation of MBAR.[33]