Many important processes - like chemical reactions, transformations of nanostructures, or protein folding - are characterized by different metastable states interconverting very slowly. This is due to the presence of barriers on the free-energy landscape. Molecular dynamics simulations of such processes, being limited in terms of time scale, require the use of enhanced sampling techniques to explore in an efficient way.
We develop new approaches to enhance the sampling thus reconstructing free energies and kinetic rates, and we implement them into freely available software, e.g., the plugin plumed, compatible with many MD codes, the vmd-interface metagui, or the tool piv_clustering to analyze liquids as well as ordered/disordered solids.
Bias exchange metadynamics
A powerful enhanced sampling approach is bias-exchange metadynamics, developed by A. Laio and S. Piana [J.Phys.Chem.B 2007]. In short, a complex, high-dimensional free energy landscape is efficiently explored by performing MD simulations of several replicas of the system, each one biased with metadynamics along a different collective variable, and periodically attempting replica exchanges with a Metropolis rule.
The low-dimensional bias acting on each replica leads the system to quickly overcome free energy barriers, whereas the exchanges between replicas guarantee a global exploration of the landscape in all directions. When all bias potentials start growing in a stationary way, they achieve convergence to the corresponding equilibrium free energy projections.
We extended this method to the reconstruction of high-dimensional free energy landscapes (see below), we studied its efficiency and convergence as a function of the simulation parameters and boundary conditions [J.Phys.Chem.B 2010] [Curr.Phys.Chem. 2012] [Comp.Phys.Comm. 2012], and we applied it to a number of challenging problems in biophysics, including protein folding, protein-ligand interaction, and amyloid aggregation.
Reconstructing the thermodynamics and kinetics
We developed a method which allows exploiting a collection of biased MD trajectories to compute accurate thermodynamic and kinetic properties [PLoS Comput.Biol. 2009]. A high-dimensional free energy surface in the space of all reaction coordinates is obtained by exploiting weighted histogram analysis. On top of the landscape, an effective Markov state kinetic model (master equation) is constructed, easily allowing to calculate the transition rates of the system, e.g. the folding rate of a protein or the binding rate of an enzyme. This makes it possible to direct compare with experiments over time scales exceeding the hours. The method follows 4 steps:
1) bias exchange trajectories are generated until convergence, i.e., until each bias profile shows self-parallel growth [Phys.Rev.E 2010]. 2) The space of reaction coordinates is partitioned into a set of bins (microstates). 3) The free energy (equilibrium probability) of each bin is computed using a weighted histogram technique. 4) For each bin, the transition rates towards neighbours are estimated assuming a Markovian diffusive behaviour. As a result, hundreds of very long kinetic monte carlo trajectories can be easily generated.
We implemented this approach in the user-friendly interface METAGUI for VMD. The method has been validated on the Ala3 test system, obtaining excellent agreement with ergodic MD trajectories, and it has been successfully applied to a number of problems, e.g., the folding of small proteins (trp-cage, insulin, huntingtin, etc.), the binding mechanism of HIV-1 protease [J.Am.Chem.Soc. 2009] and of the prion protein, and the mechanism of amyloid nucleation. All these works were based on accurate all-atom explicit-solvent force fields.
New collective variables to track secondary structure changes
A crucial ingredient of enhanced sampling approaches based on a bias potential - from umbrella sampling to steered MD to metadynamics etc. - is the identification of collective variables acting as order parameters, discriminating metastable states of the system. When such variables include all relevant slow degrees of freedom, they are good reaction coordinates which allow to obtain reliable transition pathways and free energy barriers [Annu.Rev.Phys.Chem. 2002].
I introduced a new class of reaction coordinates to tackle the formation/disruption of protein secondary structure [J.Chem.Theory.Comput. 2009], like alpha helices and beta sheets, the latter in particular being extremely challenging to sample with MD simulations. The new approach is extremely efficient and allowed us and other groups to obtain some unprecedented results, including the following:o reconstructing the atlas of experimental protein folds up to 70 amino acids [PLoS Comput.Biol. 2010]
o mapping the accurate conformational ensemble of globular proteins like SH3 and GB1 [J.Phys.Chem.Lett. 2013] [PNAS 2013]
o exploring alpha to beta transitions in human amylin monomer, in excellent agreement with transition path sampling simulations [J.Chem.Phys. 2013] [Biophys.J. 2013]
o understanding folding upon binding in intrinsically disordered proteins [PNAS 2013]
o exploring the complex nucleation mechanism of amyloid aggregates [J.Am.Chem.Soc. 2012] [Phys.Rev.Lett. 2013]
for details of publications see the section Papers