A kinetic model of trp-cage folding from multiple biased molecular dynamics simulations

Abstract

Trp-cage is a designed 20-residue polypeptide that, in spite of its size, shares several features with larger globular proteins.Although the system has been intensively investigated experimentally and theoretically, its folding mechanism is not yet fully understood. Indeed, some experiments suggest a two-state behavior, while others point to the presence of intermediates. In this work we show that the results of a bias-exchange metadynamics simulation can be used for constructing a detailed thermodynamic and kinetic model of the system. The model, although constructed from a biased simulation, has a quality similar to those extracted from the analysis of long unbiased molecular dynamics trajectories. This is demonstrated by a careful benchmark of the approach on a smaller system, the solvated Ace-Ala3-Nme peptide. For theTrp-cage folding, the model predicts that the relaxation time of 3100 ns observed experimentally is due to the presence of a compact molten globule-like conformation. This state has an occupancy of only 3% at 300 K, but acts as a kinetic trap.Instead, non-compact structures relax to the folded state on the sub-microsecond timescale. The model also predicts the presence of a state at Calpha-RMSD of 4.4 A from the NMR structure in which the Trp strongly interacts with Pro12. This state can explain the abnormal temperature dependence of the Pro12-delta3 and Gly11-alpha3 chemical shifts. The structures of the two most stable misfolded intermediates are in agreement with NMR experiments on the unfolded protein. Our work shows that, using biased molecular dynamics trajectories, it is possible to construct a model describing in detail the Trp-cage folding kinetics and thermodynamics in agreement with experimental data.

Figure 1. Bins free energies of Ala₃ from BE and from MD.

Correlation between the bins free energies calculated using Eq. 6 applied on BE simulations data and using the standard thermodynamics relation on MD results. A bin size of 30° has been used. In the inset it is shown the distribution of the deviations between the bins free energies calculated from BE and from MD, divided by the estimated error on the MD free energy. A Gaussian fit of the distribution is also shown.

Figure 2. Mean first passage times between the free energy basins of Ala₃.

Panel A: correlation between the MFPT among the four regions in , , , and , and among the eight attractors (see text and Fig. S1), obtained by MD simulations and by KMC using the kinetic model. The MFPT are calculated as the average time to go from one region to another, without passing through different regions. The error bars due to the statistical error in the MD simulations are also displayed. Large bins have a cubic side of 36°, while when not specified a cubic side of 30° is used. Panel B: distribution of FPTs from to for MD and the kinetic model. Panel C: distribution of FPTs from to for MD and the kinetic model. For panel B and C a cubic side of 30° and a time lag of 16 ps was used for calculating the diffusion matrix .

Figure 3. Dependence of the diffusion coefficient of Ala₃ on the time lag and the trajectory length.

Dependence of the slope of the line fitting the pairs of mean first passage times (see text and Fig. 2A) from the parameters used in the fit of the diffusion matrix : the length of the MD run and the time lag . For converges to the optimal value 1 (dashed line). A cubic side of 30° was used.

Figure 4. Metastable kinetic clusters of Trp-cage.

Panel A: metastable sets (clusters) detected by MCL method using . The colored spheres correspond to the lowest free energy bins of each cluster. The corresponding structures are shown with the same color code. Panel B: occupancy as a function of temperature of cluster 1, 2, and 5.

Figure 5. Trp6 interactions in the clusters reference structures of Trp-cage.

Hydrophobic contacts within 3.9 Å and hydrogen bonds(Å) are displayed. The distances(Å) between Leu7, Pro12, Arg16 and Trp6 selected protons are shown for the 3 most populated clusters. The corresponding values can be compared with the unfolded state NOE contact distances reported in Ref. . The nearest hyperpolarized Trp6 protons in the NMR experiment are selected for measuring distances. Short Ile4-Trp6 proton distances (4–5 Å) are not reported in the figure since they are found mostly in open random-coil like structures and in some more compact cluster with population <1%. This figure was generated using the program LIGPLOT .

Figure 6. Simulated NMR chemical shift deviations and ring current shifts in Trp-cage.

Panel A: correlation between experimental and calculated protons CSD for the cluster 1 (black circles), the lowest free energy bin (empty circles), and the ensemble obtained from a simulation started from the NMR structure at 282 K (black squares) and 300 K (empty squares). The continuous and dashed lines are obtained from a linear regression on the black circles and the squares, respectively. The thin dashed line corresponds to a proportionality factor of 1 between experiment and theory. Panel B: correlation between protons ring current shift temperature derivative and the corresponding ring current shift value evaluated at 298 K. Results are shown for protons (empty circles) and side chain protons (black circles). Ring current shift temperature derivative is calculated as a finite difference between 298 and 303 K using the chemical shift temperature extrapolation obtained using Eq. 9 and 10.

Figure 7. Schematic representation of the Trp-cage folding dynamics.

Times (inverse of rates) for the transitions between the relevant clusters are shown on the arrows. The uncertainty on each transition time due to both the error on the free energies and the position-dependence of is at most 40%. Only the clusters whose population is higher than 1% are shown. Continuous arrows correspond to direct transitions between clusters that occur on a time smaller than 1 µs. Dashed arrows correspond instead to transition that occur on a time larger than 1 µs or taking place through other intermediate low-populated clusters, not represented in the Figure.