L44
Biophysical Journal Volume 98 May 2010 L44–L46
Linking Well-Tempered Metadynamics Simulations with Experiments Alessandro Barducci,* Massimiliano Bonomi, and Michele Parrinello Computational Science, Department of Chemistry and Applied Biosciences, ETH Zurich, USI Campus, Lugano, Switzerland
ABSTRACT Linking experiments with the atomistic resolution provided by molecular dynamics simulations can shed light on the structure and dynamics of protein-disordered states. The sampling limitations of classical molecular dynamics can be overcome using metadynamics, which is based on the introduction of a history-dependent bias on a small number of suitably chosen collective variables. Even if such bias distorts the probability distribution of the other degrees of freedom, the equilibrium Boltzmann distribution can be reconstructed using a recently developed reweighting algorithm. Quantitative comparison with experimental data is thus possible. Here we show the potential of this combined approach by characterizing the conformational ensemble explored by a 13-residue helix-forming peptide by means of a well-tempered metadynamics/parallel tempering approach and comparing the reconstructed nuclear magnetic resonance scalar couplings with experimental data.
Received for publication 10 December 2009 and in final form 28 January 2010. *Correspondence:
[email protected]
Molecular dynamics (MD) simulations are a valuable tool for studying complex biomolecular systems because they can provide an atomistic description of their structure and dynamics. This is particularly important for protein conformationally disordered states (1,2). Characterizing the structure and the dynamics of these states with standard experimental tools is challenging, as the heterogeneous nature of these states can only be represented by ensembles of conformations. Much is to be gained by combining MD simulations with experiments. This can be achieved either by restraining the MD simulations to reproduce the experimental data (3) or by calculating the experimental observables from extensive MD sampling. Unfortunately, MD predictive power is severely limited due to the complexity of the biomolecule free energy landscape, which prevents exhaustive sampling by means of standard MD. Many methods have been proposed to overcome the sampling limitations; here we focus on metadynamics (4,5), which has been successfully applied to a variety of complex biomolecular processes (6–8). Metadynamics relies on the introduction of a history-dependent potential acting on a selected number of slow degrees of freedom, the collective variables (CVs). If properly applied, metadynamics can greatly speedup sampling and reconstruct the free energy surface (FES) associated to the CVs through the relation Vðs; t/NÞ ¼
T þ DT FðsÞ; DT
(1)
where s are the CVs, V(s, t) is the bias potential, T is the system temperature, and DT is an input parameter representing an effective sampling temperature of the CV space (5). The power of metadynamics can be further enhanced by combining it with parallel tempering (PT) (9): the combined method (PTMetaD) has been demonstrated to greatly improve the performance of both PT and single-replica meta-
dynamics (10,11). The introduction of metadynamics bias potential leads to a correct (see Eq. 1) probability distribution along the CVs but distorts that of the other degrees of freedom. Because CVs often do not correspond to experimental observables, this feature has severely limited the possibility for quantitatively comparing metadynamics simulations to experiments. Very recently, we have been able to solve this problem devising a method to recover, on-thefly, the correct Boltzmann distribution from metadynamics trajectories (12). Here, we show that the conformational ensemble spanned by a C-peptide, corresponding to the first 13 N-terminal residues of RNase A, can be easily characterized by means of PTMetaD. We validate our results via comparison of the reconstructed nuclear magnetic resonance (NMR) scalar couplings with experiments. This peptide has attracted considerable attention due to its propensity to form a-helical structures, which is unusual for such a short chain. Circular dichroism (CD) shows that the average helicity of C-peptide depends on both temperature and pH and exhibits a maximum at pH 5.25 and T ¼ 276 K (13). Even under these optimal conditions, NMR data suggested the presence of multiple conformations of different helical content (14). We choose a mutant sequence of C-peptide whose 3JHNHa scalar coupling constants and CD spectra were determined at T ¼ 277 K and pH ¼ 4.5 (15). The peptide was solvated in 1192 TIP3P (16) water molecules and the AMBER99SB (17) force field was adopted. The system was simulated with a PTMetaD approach for 20 ns using GROMACS4 (18) and PLUMED (19). We used 64 replicas spanning the
Editor: Gregory A. Voth. Ó 2010 by the Biophysical Society doi: 10.1016/j.bpj.2010.01.033
Biophysical Letters
L45
temperature interval between 274 and 650 K. The stochastic thermostat of Bussi et al. (20) was used for this. Additional technical details can be found in the Supporting Material. The metadynamics bias was applied to a single CV that measures the population of a-helical backbone H-bonds, defined as Sa ¼
6 9 X 1 ðdðOi ; HNi þ 4 Þ=3:0Þ
8;
i¼1
1 ðdðOi ; HNi þ 4 Þ=3:0Þ
(2)
where d(Oi, HNiþ4) is the distance between the backbone oxygen atom of residue i and the backbone amide hydrogen of residue i þ 4. It must be stressed that Sa does not represent the only slow mode of our system. However, combining metadynamics with PT greatly accelerates the sampling of degrees of freedom that are transverse to the CVs (10). To assess the reliability of the simulation protocol, we performed two PTMetaD simulations, starting from different initial configurations. A first run (i.e., Open) started from a random coil structure, whereas the second one (i.e., Closed) started from an a-helical configuration corresponding to ribonuclease A NMR structure (see Supporting Material). In both of the simulations, the entire CV range (from Sa x 0 to Sa x 9) was rapidly explored and the estimate of F(Sa) converged in ~10 ns. We can arbitrarily subdivide the Sa space into three regions: an almost random coil state (AC, Sa < 4); a state of an a-helicity roughly similar to that of the ribonuclease NMR structure (PH, 4 < Sa < 7); and an almost fully a-helical state (FH, Sa > 7). In Fig. 1 the estimate of the free-energy differences among those states as a function of the simulation time is reported for replica with T ¼ 278 K. The free-energy differences indicate that AC is the most stable state. However, configurations with partially formed a-helices are significantly populated (DGPH-AC x 0.5 kcal/mol), whereas fully folded a-helical structures are less likely (DGFH-AC x 2 kcal/mol). Such findings are in qualitative agreement with the average helicity measured with CD spectroscopy (13). NMR experiments performed on wild-type C-peptide suggested that in similar
FIGURE 1 Estimates of free-energy difference among the three Sa regions (see text) as a function of the simulation time for two independent PTMetaD simulations (open and closed).
conditions three sets of conformations were populated: a set of largely extended conformations; a set of largely helical forms; and a set of partially a-helical structures with a saltbridge between the Glu2 and Arg10 side chains (14). Such interaction is believed to induce a kink in the backbone preventing the N-terminal residues from adopting an a-helical structure. To get a deeper insight into the interplay between a-helix formation and Glu2-Arg10 side-chain interaction, we used the reweighting algorithm for reconstructing the FES as a function of both Sa and the salt-bridge coordinate SSB. The resulting FES is shown in Fig. 2. The landscape reconstructed from the simulations provides a detailed picture of the conformational ensemble which remarkably matches the hypothesis of Osterhout et al. (14). Indeed, the suggested sets of configurations correspond to local free-energy minima. The free-energy minimum at (Sa x 6, SSB x 0.8) corresponds to conformations with a formed salt-bridge and whose helicity is slightly less than NMR structure. Fully a-helical configurations require breaking the salt-bridge and are less stable. To make quantitative comparison with experimental data, we focus here on NMR scalar couplings. These observables are related to dihedral angles by the classical Karplus relations (21). In particular, 3JHNHa values can be expressed as a function of the F backbone angle J(F) ¼ Acos2(F þ D) þ Bcos(F þ D) þ C, where A, B, C, and D are empirical parameters. Different parameters have been proposed, but here we focus on two empirical sets: the parameters calibrated on ubiquitin (22) (UBIQ) and another set that was calibrated on flavodoxin (23) (FLAV). As done in other studies (24– 26), we measured the agreement between experiments and simulations using NX ¼ 13 2 c2 Jsim ; J exp ¼ N 1 hJi isim Ji;exp =s2 ;
(3)
i¼1
where s ¼ 0.91 is the estimated systematic error (25,26). The results for the two PTMetaD simulations are reported in Table 1. Although the agreement is in general good, its
FIGURE 2 FES as a function of Sa and SSB from the open simulation. SSB ¼ [1 – (d/5.0)6]/[1 – (d/5.0)12], where d is the distance in ˚ ngstroms between the Cd of Glu2 and the central carbon atom in A the guanidine group of Arg10. Contour lines are plotted every kBT. Biophysical Journal 98(9) L44–L46
L46 TABLE 1
Biophysical Letters c2 values
Open Closed MD
UBIQ
FLAV
1.71 1.88 4.56
0.88 0.92 3.29
c2 values using UBIQ (22) and FLAV (23) parameter set for the PTMetaD simulations and for a standard 20-ns MD simulation averaged over eight independent runs (see Supporting Material).
quality depends on the parameter choice. However, the scalar couplings obtained in the two simulations are in good agreement within the same parameter set (c2(Jclosed, Jopen) < 0.09). This makes us confident that sampling has been adequate and the remaining error is due to force-field inaccuracies and/or J(F) parameterization. In this letter, we presented the characterization of the disordered structural ensemble explored by a mutant of C-peptide by means of PTMetaD simulation. Using a reweighting algorithm, we have been able to compare our results to NMR data. Such comparison allowed us both to validate our simulations and to gain a deeper insight into the conformationally averaged information provided by the experiments. We believe that such a combined approach may provide a much better understanding of the protein conformationally heterogeneous states. Future developments will point toward the reconstruction of the kinetics relative to the conformational transitions by means of stochastic dynamics simulations on the reconstructed FESs. SUPPORTING MATERIAL One figure and one table are available at http://www.biophysj.org/biophysj/ supplemental/S0006-3495(10)00198-0.
ACKNOWLEDGMENTS The authors thank ‘‘Centro Enrico Fermi’’ for providing the computational resources.
REFERENCES and FOOTNOTES 1. Dyson, H. J., and P. E. Wright. 2005. Intrinsically unstructured proteins and their functions. Nat. Rev. Mol. Cell Biol. 6:197–208. 2. Chiti, F., and C. M. Dobson. 2006. Protein misfolding, functional amyloid, and human disease. Annu. Rev. Biochem. 75:333–366. 3. Vendruscolo, M. 2007. Determination of conformationally heterogeneous states of proteins. Curr. Opin. Struct. Biol. 17:15–20. 4. Laio, A., and M. Parrinello. 2002. Escaping free-energy minima. Proc. Natl. Acad. Sci. USA. 99:12562–12566. 5. Barducci, A., G. Bussi, and M. Parrinello. 2008. Well-tempered metadynamics: a smoothly converging and tunable free-energy method. Phys. Rev. Lett. 100:020603.
Biophysical Journal 98(9) L44–L46
6. Parrinello, M. 2008. Physical Biology: From Atoms to Medicine A. H. Zewail, editor. Imperial College Press, London, UK. 7. Laio, A., and F. L. Gervasio. 2008. Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science. Rep. Prog. Phys. 71:126601. 8. Pfaendtner, J., D. Branduardi, ., G. A. Voth. 2009. Nucleotide-dependent conformational states of actin. Proc. Natl. Acad. Sci. USA. 106:12723–12728. 9. Hansmann, U. H. E. 1997. Parallel tempering algorithm for conformational studies of biological molecules. Chem. Phys. Lett. 281:140–150. 10. Bussi, G., F. L. Gervasio, ., M. Parrinello. 2006. Free-energy landscape for b-hairpin folding from combined parallel tempering and metadynamics. J. Am. Chem. Soc. 128:13435–13441. 11. Bonomi, M., F. L. Gervasio, ., M. Parrinello. 2007. Insight into the folding inhibition of the HIV-1 protease by a small peptide. Biophys. J. 93:2813–2821. 12. Bonomi, M., A. Barducci, and M. Parrinello. 2009. Reconstructing the equilibrium Boltzmann distribution from well-tempered metadynamics. J. Comput. Chem. 30:1615–1621. 13. Shoemaker, K. R., P. S. Kim, ., R. L. Baldwin. 1987. Tests of the helix dipole model for stabilization of a-helices. Nature. 326:563–567. 14. Osterhout, Jr., J. J., R. L. Baldwin, ., P. E. Wright. 1989. 1H NMR studies of the solution conformations of an analogue of the C-peptide of ribonuclease A. Biochemistry. 28:7059–7064. 15. Caballero-Herrera, A., K. Nordstrand, and K. Berndt. 2005. Effect of urea on peptide conformation in water: molecular dynamics and experimental characterization. Biophys. J. 89:842–857. 16. Jorgensen, W. L., J. Chandrasekhar, ., M. L. Klein. 1983. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79:926–935. 17. Hornak, V., R. Abel, ., C. Simmerling. 2006. Comparison of multiple AMBER force fields and development of improved protein backbone parameters. Proteins. 65:712–725. 18. Hess, B., C. Kutzner, ., E. Lindahl. 2008. GROMACS 4: algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 4:435–447. 19. Bonomi, M., D. Branduardi, ., M. Parrinello. 2009. PLUMED: a portable plugin for free-energy calculations with molecular dynamics. Comput. Phys. Commun. 180:1961–1972. 20. Bussi, G., D. Donadio, and M. Parrinello. 2007. Canonical sampling through velocity rescaling. J. Chem. Phys. 126:014101. 21. Karplus, M. 1959. Contact electron-spin coupling of nuclear magnetic moments. J. Chem. Phys. 30:11–15. 22. Hu, J. S., and A. Bax. 1997. Determination of f and c1 angles in proteins from C13-C13 three-bond J couplings measured by threedimensional heteronuclear NMR. how planar is the peptide bond? J. Am. Chem. Soc. 119:6360–6368. 23. Schmidt, J. M., M. Blumel, ., H. Ruterjans. 1999. Self-consistent 3J coupling analysis for the joint calibration of Karplus coefficients and evaluation of torsion angles. J. Biomol. NMR. 14:1–12. 24. Best, R. B., N. V. Buchete, and G. Hummer. 2008. Are current molecular dynamics force fields too helical? Biophys. J. 95:L07–L09. 25. Best, R. B., and G. Hummer. 2009. Optimized molecular dynamics force fields applied to the helix-coil transition of polypeptides. J. Phys. Chem. B. 113:9004–9015. 26. Wickstrom, L., A. Okur, and C. Simmerling. 2009. Evaluating the performance of the ff99SB force field based on NMR scalar coupling data. Biophys. J. 97:853–856.