Alessandro Laio

Escaping free-energy minima

We introduce a powerful method for exploring the properties of the multidimensional free energy surfaces (FESs) of complex many-body systems by means of coarse-grained non-Markovian dynamics in the space defined by a few collective coordinates. A characteristic feature of these dynamics is the presence of a history-dependent potential term that, in time, fills the minima in the FES, allowing the efficient exploration and accurate determination of the FES as a function of the collective coordinates. We demonstrate the usefulness of this approach in the case of the dissociation of a NaCl molecule in water and in the study of the conformational changes of a dialanine in solution.

Clustering by fast search and find of density peaks

Discerning clusters of data points Cluster analysis is used in many disciplines to group objects according to a defined measure of distance. Numerous algorithms exist, some based on the analysis of the local density of data points, and others on predefined probability distributions. Rodriguez and Laio devised a method in which the cluster centers are recognized as local density maxima that are far away from any points of higher density. The algorithm depends only on the relative densities rather than their absolute values. The authors tested the method on a series of data sets, and its performance compared favorably to that of established techniques. Science, this issue p. 1492 Local density of points is ranked and analyzed to categorize data. Cluster analysis is aimed at classifying elements into categories on the basis of their similarity. Its applications range from astronomy to bioinformatics, bibliometrics, and pattern recognition. We propose an approach based on the idea that cluster centers are characterized by a higher density than their neighbors and by a relatively large distance from points with higher densities. This idea forms the basis of a clustering procedure in which the number of clusters arises intuitively, outliers are automatically spotted and excluded from the analysis, and clusters are recognized regardless of their shape and of the dimensionality of the space in which they are embedded. We demonstrate the power of the algorithm on several test cases.

Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science

Metadynamics is a powerful algorithm that can be used both for reconstructing the free energy and for accelerating rare events in systems described by complex Hamiltonians, at the classical or at the quantum level. In the algorithm the normal evolution of the system is biased by a history-dependent potential constructed as a sum of Gaussians centered along the trajectory followed by a suitably chosen set of collective variables. The sum of Gaussians is exploited for reconstructing iteratively an estimator of the free energy and forcing the system to escape from local minima. This review is intended to provide a comprehensive description of the algorithm, with a focus on the practical aspects that need to be addressed when one attempts to apply metadynamics to a new system: (i) the choice of the appropriate set of collective variables; (ii) the optimal choice of the metadynamics parameters and (iii) how to control the error and ensure convergence of the algorithm.

A Hamiltonian electrostatic coupling scheme for hybrid Car-Parrinello molecular dynamics simulations

We present a fully Hamiltonian and computationally efficient scheme to include the electrostatic effects due to the classical environment in a Car–Parrinello mixed quantum Mechanics/molecular mechanics (QM/MM) method. The polarization due to the MM atoms close to the quantum system is described by a Coulombic potential modified at short range. We show that the functional form of this potential has to be chosen carefully in order to obtain the correct interaction properties and to prevent an unphysical escape of the electronic density to the MM atoms (the so-called spill-out effect). The interaction between the QM system and the more distant MM atoms is modeled by a Hamiltonian term explicitly coupling the multipole moments of the quantum charge distribution with the classical point charges. Our approach remedies some of the well known deficiencies of current electrostatic coupling schemes in QM/MM methods, allowing molecular dynamics simulations of mixed systems within a fully consistent and energy conser...

Predicting crystal structures: the Parrinello-Rahman method revisited.

By suitably adapting a recent approach [A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U.S.A. 99, 12 562 (2002)]] we develop a powerful molecular dynamics method for the study of pressure-induced structural transformations. We use the edges of the simulation cell as collective variables and define a metadynamics that drives the system away from the local minimum towards a new crystal structure. In contrast to the Parrinello-Rahman method, our approach shows no hysteresis, and crystal structure transformations can occur at the equilibrium pressure. We illustrate the power of the method by studying the pressure-induced diamond to simple hexagonal phase transition in a model of silicon.

Equilibrium Free Energies from Nonequilibrium Metadynamics

In this Letter we propose a new formalism to map history-dependent metadynamics in a Markovian process. We apply this formalism to model Langevin dynamics and determine the equilibrium distribution of a collection of simulations. We demonstrate that the reconstructed free energy is an unbiased estimate of the underlying free energy and analytically derive an expression for the error. The present results can be applied to other history-dependent stochastic processes, such as Wang-Landau sampling.

Anisotropy of Earth's D'' layer and stacking faults in the MgSiO3 post-perovskite phase.

The post-perovskite phase of (Mg,Fe)SiO3 is believed to be the main mineral phase of the Earths lowermost mantle (the D″ layer). Its properties explain numerous geophysical observations associated with this layer—for example, the D″ discontinuity, its topography and seismic anisotropy within the layer. Here we use a novel simulation technique, first-principles metadynamics, to identify a family of low-energy polytypic stacking-fault structures intermediate between the perovskite and post-perovskite phases. Metadynamics trajectories identify plane sliding involving the formation of stacking faults as the most favourable pathway for the phase transition, and as a likely mechanism for plastic deformation of perovskite and post-perovskite. In particular, the predicted slip planes are {010} for perovskite (consistent with experiment) and {110} for post-perovskite (in contrast to the previously expected {010} slip planes). Dominant slip planes define the lattice preferred orientation and elastic anisotropy of the texture. The {110} slip planes in post-perovskite require a much smaller degree of lattice preferred orientation to explain geophysical observations of shear-wave anisotropy in the D″ layer.

A Kinetic Model of Trp-Cage Folding from Multiple Biased Molecular Dynamics Simulations

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 the Trp-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 of 4.4 Å from the NMR structure in which the Trp strongly interacts with Pro12. This state can explain the abnormal temperature dependence of the and 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.

Substrate Binding Mechanism of HIV-1 Protease from Explicit-Solvent Atomistic Simulations

The binding mechanism of a peptide substrate (Thr-Ile-Met-Met-Gln-Arg, cleavage site p2-NC of the viral polyprotein) to wild-type HIV-1 protease has been investigated by 1.6 micros biased all-atom molecular dynamics simulations in explicit water. The configuration space has been explored biasing seven reaction coordinates by the bias-exchange metadynamics technique. The structure of the Michaelis complex is obtained starting from the substrate outside the enzyme within a backbone rmsd of 0.9 A. The calculated free energy of binding is -6 kcal/mol, and the kinetic constants for association and dissociation are 1.3 x 10(6) M(-1) s(-1) and 57 s(-1), respectively, consistent with experiments. In the main binding pathway, the flaps of the protease do not open sizably. The substrate slides inside the enzyme cavity from the tight lateral channel. This may contrast with the natural polyprotein substrate which is expected to bind by opening the flaps. Thus, mutations might influence differently the binding kinetics of peptidomimetic ligands and of the natural substrate.

Reconstructing the density of states by history-dependent metadynamics.

We present a novel method for the calculation of the energy density of states D(E) for systems described by classical statistical mechanics. The method builds on an extension of a recently proposed strategy that allows the free-energy profile of a canonical system to be recovered within a preassigned accuracy [Proc. Natl. Acad. Sci. U.S.A. 99, 12562 (2002)]]. The method allows a good control over the error on the recovered system entropy. This fact is exploited to obtain D(E) more efficiently by combining measurements at different temperatures. The accuracy and efficiency of the method are tested for the two-dimensional Ising model (up to size 50 x 50) by comparison with both exact results and previous studies. This method is a general one and should be applicable to more realistic model systems.