Microsecond simulations of the folding/unfolding thermodynamics of the Trp-cage miniprotein

Abstract

We study the unbiased folding/unfolding thermodynamics of the Trp-cage miniprotein using detailed molecular dynamics simulations of an all-atom model of the protein in explicit solvent using the Amberff99SB force field. Replica-exchange molecular dynamics simulations are used to sample the protein ensembles over a broad range of temperatures covering the folded and unfolded states at two densities. The obtained ensembles are shown to reach equilibrium in the 1 mus/replica timescale. The total simulation time used in the calculations exceeds 100 mus. Ensemble averages of the fraction folded, pressure, and energy differences between the folded and unfolded states as a function of temperature are used to model the free energy of the folding transition, DeltaG(P, T), over the whole region of temperatures and pressures sampled in the simulations. The DeltaG(P, T) diagram describes an ellipse over the range of temperatures and pressures sampled, predicting that the system can undergo pressure-induced unfolding and cold denaturation at low temperatures and high pressures, and unfolding at low pressures and high temperatures. The calculated free energy function exhibits remarkably good agreement with the experimental folding transition temperature (T(f) = 321 K), free energy, and specific heat changes. However, changes in enthalpy and entropy are significantly different than the experimental values. We speculate that these differences may be due to the simplicity of the semiempirical force field used in the simulations and that more elaborate force fields may be required to describe appropriately the thermodynamics of proteins.

Figure 1. Time evolution of the REMD trajectories

A) Fraction of the time that each replica (1-40) spends in the folded state (rmsd<0.22 nm) over 1 ns time block averages. All replicas were started from an unfolded configuration. B) First folding time for replicas along the folding trajectory (red line). The profile shows an exponential growth with exponent of the order on 100 ns. The green line shows the first time a replica has folded, unfolded past 0.6 nm, and refolded. The curve is very similar to the first folding time, expect for a time delay of tens of ns. C) Fraction of replicas that are folded at any time during the simulation as a function of time (red line(UIC) and green line(FIC)). All replicas, at different temperatures are averaged together. We also show the fraction of replicas that are folded for a simulation in which all replicas were started in the folded state. The two ensembles show similar averages after 200 ns. The blue curve in the inset of c) is the time correlation function of the fraction folded time-history shows. Estimates of the correlation time for the fraction for replicas folded are 50 ns. We use simulation segments if twice this time, 100 ns, to estimate uncertainties in the averages.

Figure 2

Temperature dependence of distributions of a) rmsd, b) Trp Nε to Asp backbone carbonyl hydrogen bond, and c) Arg- Asp ion pair distances. All distances are in nm.

Figure 3. Superposition of representative structures belonging to the folded basins. Four representative structures are obtained by clustering the structures sampled at the lowest temperature (280K), during the last 10 ns of the REMD simulation

The top most populated clusters correspond to structures at 0.075 nm (48%), 0.69 nm (7%), 0.15 nm (6%), 0.17 nm (5%), and 0.21 nm (4%) from structure 1 of the NMR ensemble. Cluster 2 is assigned as an unfolded state. The plots show the superposition of the structures in the folded basin (rmsd < 0.22 nm). We used the Daura clustering method implemented in Gromacs. A) Tube representation of the backbone. B) All non hydrogen atoms. The structures representative of the top four clusters in the folded basins are labeled in cyan first cluster), blue (third cluster), red (fourth cluster), and orange (fifth cluster). The main differences in the structures are in the 3-10 turn (Pro 12- Ser13-Ser-14-Gly15) region and the N and C termini.

Figure 4. Contour maps of the free energy of the Trp-cage as a function of the rmsd and the radius of gyration (R_g) at 280K, 320K and 350K

The contour shows that the system occupies two main basins- - one corresponding to the folded state (low rmsd and compact states) and another at larger rmsd > 0.4 nm and less compact, corresponding to the unfolded state. At even higher temperatures the system becomes extended and only the large rmsd and R_g basin is occupied. The folded basin for low rmsd and R_g has two basins (substates) at low temperature. The legends use a logarithmic scale, with -10 corresponding to the highest probability and 0 to the lowest.

Figure 5. Fraction of replicas folded as a function of temperature

A) Comparison of the fraction of Trp-cage protein folded as a function of temperature for two different force fields—ff99SB (red and green curves) and ff94 (blue curve). The folding profiles for the ff99SB are practically indistinguishable for REMD simulation starting from completely unfolded states (red) and folded states (green). The ff99SB profiles cross the 50% fraction folded at T=321 K, in close agreement with experimental data (T_exp = 317 K). However, the calculated fraction of folded states is lower than measured values at low T. B) Fraction folded as a function of other order parameters that are representative of the folded Trp-cage folded structure. These parameters are rmsd (red curve--folded is rmsd< 0.22 nm), fraction of a Trp-Nε – Asp backbone carbonyl hydrogen bond (blue curve-- d < .4 nm), Arg-Asp ion pair (magenta curve -- d<0 .6 nm), and number of alpha helical amino acids (divided by five to fit the same scale in the plot, green curve). An alpha helical amino acid is an amino acid with three consecutive phi-psi angles in the alpha helical region ( ϕ= −60 +/− 30 Deg. and φ = −47+/− 30 Deg) . The maximum number of alpha helical segments observed is ten.

Figure 6. Comparison of the computed (solid line) and measured (dotted line)

a) Free energy, b) fraction folded, c) entropy and d) enthalpy, as a function of temperature.

Figure 7. Stability of the Trp-cage protein as function of temperature and pressure

A) Average energy difference between unfolded and folded state ensembles calculated at various temperatures for a system at low pressure. B) Average energy difference between unfolded and folded state ensembles calculated at various temperatures for a system at high pressure. The uncertainties are calculated over five, 100 ns long, time segments. C) Fraction of the ensemble folded as a function of temperature. The two curves show the results obtained from simulations at low and high average pressures. The solid lines are the best χ² fit of the thermodynamics of the system, following a fitting procedure described in the text. The parameters describing the changes in Gibbs free energy are shown in table 1. D) Pressure-Temperature stability diagram for the Trp-cage protein. The contours for ΔG(P,T)=0 (contour closer to P=0 and T=320K), −5, −10,−15 kJ/mol are elliptical. The free energy profile shows small pressure dependence and the protein will require high pressures to unfold at low temperatures. The dotted lines represent the states sampled during the REMD simulations at low (blue) and high (red) densities. The black solid line represents the ΔV=0 isochore. The red solid line represents ΔS=0. All states above this curve decrease stability upon cooling, and have higher entropy in the folded state than in the unfolded state.