Free-energy landscapes and thermodynamic parameters of complex molecules from nonequilibrium simulation trajectories

Thermodynamic parameters such as free energies and heat capacities are important quantities for understanding processes involving structural transitions in complex molecules such as proteins. Computational investigations provide simulated data that can be used for calculating thermodynamic parameters. However, calculations give accurate results only if the simulations sample all of configuration space with the appropriate temperature-dependent Boltzmann equilibrium probabilities. For many systems, truly comprehensive sampling of configuration space is not computationally feasible. We present an approximation technique for the calculations that will give accurate values for thermodynamic parameters when the data is incomplete. Our work is applicable to systems in which there are two distinct, important regions of configuration space that must be sampled. Importantly, the results are also valid when the system is more complex than two-state systems. Transition pathways that involve intermediate configurations between two stable regions are allowed in this treatment, and therefore the results are valid for multistate systems.

Figure 1

(a) Energy time series for folding and unfolding dynamics of the leucine zipper dimer. This simulation was run at a temperature that allows multiple transitions between low-energy native states and high-energy unfolded states. (b) Magnified view of a section in the MC time series which clearly shows the multistate nature of the transition. (c) Transitions occur between the high-energy disorganized unfolded state, a midenergy transition state, and a low-energy organized native-state dimer.

Figure 2

Native structure of the two-helix bundle showing each helix and the sidechain interactions (dark circles) between the two helices.

Figure 3

Calculated heat capacity as a function of temperature for the two-helix bundle. The open circle (○) τ${}_e$$=$∞ data points are the correct values to be used as reference points. The solid curve is the heat capacity calculated from single-temperature ($T~352\hspace{0.5em}\mathrm{K}$) equilibrium simulations using the histogram technique [54] and helps guide the eye. (a) Uncompensated calculations are shown at $T$$=$337 and 367 K for different τ${}_e$, given in units of M$=$10${}^6$ MC steps. For smaller τ${}_e$, the success rate for simulations decreases (Table ) and the uncompensated $C_v$ deviate further from the reference values, exemplified by the 2M value at 367 K. (b) Heat capacities from the same incomplete sampling simulations used in (a) but calculated using the compensation method of Eq. (11). The compensation method produces results that are very close to the correct value (τ${}_e$$=$∞) for all values of τ${}_e$, including τ${}_e$ that are so short that configuration space is poorly sampled and the success rate for transitions is low.

Figure 4

Free-energy curves $F$($E$) for the two-helix bundle as a function of energy at the same temperatures that are focused on for Fig. 3, 337 K $(<T_c)$ and 367 K $(>T_c)$. Different $F$($E$) curves are calculated using simulations that imposed different cutoff times τ${}_e$. (a) and (c): Uncompensated results showing large deviations from the correct τ${}_e$ $=$ ∞ reference curve. (b) and (d): Compensated calculations produce results that collapse for all τ${}_e$ onto the reference curve τ${}_e$ $=$ ∞.

Figure 5

Free-energy landscapes for the leucine zipper dimer of Fig. 1. Because of the complexity of the system, only nonequilibrium, incomplete simulation trajectories are available. The calculation of $F$ employs the compensation technique of Eq. (11). Red represent low $F$ and blue represent high $F$ (scale shown in kcal mol${}^{-1}$). In both (a) and (b), $F$($q$,$x$) is calculated for different values of the secondary structural parameter, $q$ which represents the fraction of amino acids that are in the native helical state. The other structural parameter $x$ is (a) the fraction of native tertiary contacts $Q$ that are intact. The region with $q$ ∼ 1, $Q$ ∼ 1 corresponds to the native-state basin, whereas $q$ ∼ 0 and $Q$ ∼ 0 correspond to the unfolded, random coil basin. (b) In addition to $q$, the other structural parameter is the center-of-mass separation $d_{\mathrm{cm}}$. The region with $q$ ∼ 1, $d_{\mathrm{cm}}~6$ corresponds to the native-state basin.