Fabio Pietrucci

PLUMED: a portable plugin for free-energy calculations with molecular dynamics

Here we present a program aimed at free-energy calculations in molecular systems. It consists of a series of routines that can be interfaced with the most popular classical molecular dynamics (MD) codes through a simple patching procedure. This leaves the possibility for the user to exploit many different MD engines depending on the system simulated and on the computational resources available. Free-energy calculations can be performed as a function of many collective variables, with a particular focus on biological problems, and using state-of-the-art methods such as metadynamics, umbrella sampling and Jarzynski-equation based steered MD. The present software, written in ANSI-C language, can be easily interfaced with both Fortran and C/C++ codes.

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.

A Collective Variable for the Efficient Exploration of Protein Beta-Sheet Structures: Application to SH3 and GB1

We introduce a new class of collective variables which allow forming efficiently beta-sheet structures in all-atom explicit-solvent simulations of proteins. By this approach we are able to systematically fold a 16-residue beta hairpin using metadynamics on a single replica. Application to the 56-residue SH3 and GB1 proteins show that, starting from extended states, in ∼100 ns tens of structures containing more than 30% beta-sheet are obtained, including parts of the native fold. Using these variables may allow folding moderate size proteins with an accurate explicit solvent description. Moreover, it may allow investigating the presence of misfolded states that are relevant for diseases (e.g., prion and Alzheimer) and studying beta-aggregation (amyloid diseases).

METAGUI. A VMD interface for analyzing metadynamics and molecular dynamics simulations

We present a new computational tool, METAGUI, which extends the VMD program with a graphical user interface that allows constructing a thermodynamic and kinetic model of a given process simulated by large-scale molecular dynamics. The tool is specially designed for analyzing metadynamics based simulations. The huge amount of diverse structures generated during such a simulation is partitioned into a set of microstates (i.e. structures with similar values of the collective variables). Their relative free energies are then computed by a weighted-histogram procedure and the most relevant free energy wells are identified by diagonalization of the rate matrix followed by a commitor analysis. All this procedure leads to a convenient representation of the metastable states and long-time kinetics of the system which can be compared with experimental data. The tool allows to seamlessly switch between a collective variables space representation of microstates and their atomic structure representation, which greatly facilitates the set-up and analysis of molecular dynamics simulations. METAGUI is based on the output format of the PLUMED plugin, making it compatible with a number of different molecular dynamics packages like AMBER, NAMD, GROMACS and several others. The METAGUI source files can be downloaded from the PLUMED web site (http://www.plumed-code.org).

Multidimensional view of amyloid fibril nucleation in atomistic detail.

Starting from a disordered aggregate, we have simulated the formation of ordered amyloid-like beta structures in a system formed by 18 polyvaline chains in explicit solvent, employing molecular dynamics accelerated by bias-exchange metadynamics. We exploited 8 different collective variables to compute the free energy of hundreds of putative aggregate structures, with variable content of parallel and antiparallel β-sheets and different packing among the sheets. This allowed characterizing in detail a possible nucleation pathway for the formation of amyloid fibrils: first the system forms a relatively large ordered nucleus of antiparallel β-sheets, and then a few parallel sheets start appearing. The relevant nucleation process culminates at this point: when a sufficient number of parallel sheets is formed, the free energy starts to decrease toward a new minimum in which this structure is predominant. The complex nucleation pathway we found cannot be described within classical nucleation theory, namely employing a unique simple reaction coordinate like the total content of β-sheets.

Exploring the Universe of Protein Structures beyond the Protein Data Bank

It is currently believed that the atlas of existing protein structures is faithfully represented in the Protein Data Bank. However, whether this atlas covers the full universe of all possible protein structures is still a highly debated issue. By using a sophisticated numerical approach, we performed an exhaustive exploration of the conformational space of a 60 amino acid polypeptide chain described with an accurate all-atom interaction potential. We generated a database of around 30,000 compact folds with at least of secondary structure corresponding to local minima of the potential energy. This ensemble plausibly represents the universe of protein folds of similar length; indeed, all the known folds are represented in the set with good accuracy. However, we discover that the known folds form a rather small subset, which cannot be reproduced by choosing random structures in the database. Rather, natural and possible folds differ by the contact order, on average significantly smaller in the former. This suggests the presence of an evolutionary bias, possibly related to kinetic accessibility, towards structures with shorter loops between contacting residues. Beside their conceptual relevance, the new structures open a range of practical applications such as the development of accurate structure prediction strategies, the optimization of force fields, and the identification and design of novel folds.

The Fate of a Zwitterion in Water from ab Initio Molecular Dynamics: Monoethanolamine (MEA)-CO2

Understanding the fundamental reactions accompanying the capture of carbon dioxide in amine solutions is critical for the design of high-performance solvents and requires an accurate modeling of the solute-solvent interaction. As a first step toward this goal, using ab initio molecular dynamics (Car-Parrinello) simulations, we investigate a zwitterionic carbamate, a species long proposed as intermediate in the formation of a stable carbamate, in a dilute aqueous solution. CO2 release and deprotonation are competitive routes for its dissociation and are both characterized by free-energy barriers of 6-8 kcal/mol. Water molecules play a crucial role in both pathways, resulting in large entropic effects. This is especially true in the case of CO2 release, which is accompanied by a strong reorganization of the solvent beyond the first coordination shell, leading to the formation of a water cage entrapping the solute (hydrophobic effect). Our results contrast with the assumptions of implicit solvent models.

Bridging Static and Dynamical Descriptions of Chemical Reactions: An ab Initio Study of CO2 Interacting with Water Molecules.

Extracting reliable thermochemical parameters from molecular dynamics simulations of chemical reactions, although based on ab initio methods, is generally hampered by difficulties in reproducing the results and controlling the statistical errors. This is a serious drawback with respect to the quantum-chemical description based on potential energy surfaces. This work is an attempt to fill this gap. We apply molecular dynamics, based on density functional theory (DFT) and empowered by path metadynamics (MTD), to simulate the reaction of CO2 with (one, two, and three) water molecules in the gas phase. This study relies on a strategy that ensures a precise control of the accuracy of the reaction coordinates and of the reconstructed free-energy surface on this space, namely, on (i) fully reversible MTD simulations, (ii) a committor probability analysis for the diagnosis of the collective variables, and (iii) a cluster analysis for the characterization of the reconstructed free-energy surfaces. This robust procedure permits a meaningful comparison with more traditional calculations of the potential energy surfaces that we also perform within the same DFT computational scheme. This comparison shows in particular that the reactants and products of systems with only three water molecules can no longer be understood in terms of one structure but must be described as statistical configuration ensembles. Calculations carried out with different prescriptions for the exchange-correlation functionals also allow us to establish their quantitative effect on the activation barriers for the formation and the dissociation of carbonic acid. Their decrease induced by the addition of one water molecule (catalytic effect) is found to be largely independent of the specific functional.

Optimizing the Performance of Bias-Exchange Metadynamics: Folding a 48-Residue LysM Domain Using a Coarse-Grained Model

Computer simulation of complex conformational transitions in biomolecules, such as protein folding, is considered one of the main goals of computational chemistry. A recently developed methodology, bias-exchange metadynamics, was successfully used to reversibly fold some small globular proteins. The objective of this work is to further improve this promising technique. This is accomplished by searching for the optimal set of parameters that enable folding a 48 amino acid protein, 1E0G , in the shortest possible time, using a coarse-grained force field UNRES. It is shown that bias-exchange metadynamics, if appropriately optimized, allows finding the folded state of 1E0G significantly faster than normal replica exchange.