Ernesto Suárez

Entropy Calculations of Single Molecules by Combining the Rigid-Rotor and Harmonic-Oscillator Approximations with Conformational Entropy Estimations from Molecular Dynamics Simulations.

As shown by previous theoretical and computational work, absolute entropies of small molecules that populate different conformers can be predicted accurately on the basis of the partitioning of the intramolecular entropy into vibrational and conformational contributions. Herein, we further elaborate on this idea and propose a protocol for entropy calculations of single molecules that combines the rigid rotor harmonic oscillator (RRHO) entropies with the direct sampling of the molecular conformational space by means of classical molecular dynamics simulations. In this approach, the conformational states are characterized by discretizing the time evolution of internal rotations about single bonds, and subsequently, the mutual information expansion (MIE) is used to approach the full conformational entropy from the converged probability density functions of the individual torsion angles, pairs of torsions, triads, and so on. This RRHO&MIE protocol could have broad applicability, as suggested by our test calculations on systems ranging from hydrocarbon molecules in the gas phase to a polypeptide molecule in aqueous solution. For the hydrocarbon molecules, the ability of the RRHO&MIE protocol to predict absolute entropies is assessed by carefully comparing theoretical and experimental values in the gas phase. For the rest of the test systems, we analyze the advantages and limitations of the RRHO&MIE approach in order to capture high order correlation effects and yield converged conformational entropies within a reasonable simulation time. Altogether, our results suggest that the RRHO&MIE strategy could be useful for estimating absolute and/or relative entropies of single molecules either in the gas phase or in solution.

WESTPA: An Interoperable, Highly Scalable Software Package for Weighted Ensemble Simulation and Analysis

The weighted ensemble (WE) path sampling approach orchestrates an ensemble of parallel calculations with intermittent communication to enhance the sampling of rare events, such as molecular associations or conformational changes in proteins or peptides. Trajectories are replicated and pruned in a way that focuses computational effort on underexplored regions of configuration space while maintaining rigorous kinetics. To enable the simulation of rare events at any scale (e.g., atomistic, cellular), we have developed an open-source, interoperable, and highly scalable software package for the execution and analysis of WE simulations: WESTPA (The Weighted Ensemble Simulation Toolkit with Parallelization and Analysis). WESTPA scales to thousands of CPU cores and includes a suite of analysis tools that have been implemented in a massively parallel fashion. The software has been designed to interface conveniently with any dynamics engine and has already been used with a variety of molecular dynamics (e.g., GROMACS, NAMD, OpenMM, AMBER) and cell-modeling packages (e.g., BioNetGen, MCell). WESTPA has been in production use for over a year, and its utility has been demonstrated for a broad set of problems, ranging from atomically detailed host–guest associations to nonspatial chemical kinetics of cellular signaling networks. The following describes the design and features of WESTPA, including the facilities it provides for running WE simulations and storing and analyzing WE simulation data, as well as examples of input and output.

Simultaneous Computation of Dynamical and Equilibrium Information Using a Weighted Ensemble of Trajectories.

Equilibrium formally can be represented as an ensemble of uncoupled systems undergoing unbiased dynamics in which detailed balance is maintained. Many nonequilibrium processes can be described by suitable subsets of the equilibrium ensemble. Here, we employ the “weighted ensemble” (WE) simulation protocol [Huber and Kim, Biophys. J.1996, 70, 97–110] to generate equilibrium trajectory ensembles and extract nonequilibrium subsets for computing kinetic quantities. States do not need to be chosen in advance. The procedure formally allows estimation of kinetic rates between arbitrary states chosen after the simulation, along with their equilibrium populations. We also describe a related history-dependent matrix procedure for estimating equilibrium and nonequilibrium observables when phase space has been divided into arbitrary non-Markovian regions, whether in WE or ordinary simulation. In this proof-of-principle study, these methods are successfully applied and validated on two molecular systems: explicitly solvated methane association and the implicitly solvated Ala4 peptide. We comment on challenges remaining in WE calculations.

Thermochemical Fragment Energy Method for Biomolecules: Application to a Collagen Model Peptide.

Herein, we first review different methodologies that have been proposed for computing the quantum mechanical (QM) energy and other molecular properties of large systems through a linear combination of subsystem (fragment) energies, which can be computed using conventional QM packages. Particularly, we emphasize the similarities among the different methods that can be considered as variants of the multibody expansion technique. Nevertheless, on the basis of thermochemical arguments, we propose yet another variant of the fragment energy methods, which could be useful for, and readily applicable to, biomolecules using either QM or hybrid quantum mechanical/molecular mechanics methods. The proposed computational scheme is applied to investigate the stability of a triple-helical collagen model peptide. To better address the actual applicability of the fragment QM method and to properly compare with experimental data, we compute average energies by carrying out single-point fragment QM calculations on structures generated by a classical molecular dynamics simulation. The QM calculations are done using a density functional level of theory combined with an implicit solvent model. Other free-energy terms such as attractive dispersion interactions or thermal contributions are included using molecular mechanics. The importance of correcting both the intermolecular and intramolecular basis set superposition error (BSSE) in the QM calculations is also discussed in detail. On the basis of the favorable comparison of our fragment-based energies with experimental data and former theoretical results, we conclude that the fragment QM energy strategy could be an interesting addition to the multimethod toolbox for biomolecular simulations in order to investigate those situations (e.g., interactions with metal clusters) that are beyond the range of applicability of common molecular mechanics methods.

CENCALC: a computational tool for conformational entropy calculations from molecular simulations.

We present the CENCALC software that has been designed to estimate the conformational entropy of single molecules from extended Molecular Dynamics (MD) simulations in the gas‐phase or in solution. CENCALC uses both trajectory coordinates and topology information in order to characterize the conformational states of the molecule of interest by discretizing the time evolution of internal rotations. The implemented entropy methods are based on the mutual information expansion, which is built upon the converged probability density functions of the individual torsion angles, pairs of torsions, triads, and so on. Particularly, the correlation‐corrected multibody local approximation selects an optimum cutoff in order to retrieve the maximum amount of genuine correlation from a given MD trajectory. We illustrate these capabilities by carrying out conformational entropy calculations for a decapeptide molecule either in its unbound form or in complex with a metalloprotease enzyme. CENCALC is distributed under the GNU public license at http://sourceforge.net/projects/cencalc/.

Entropic Control of the Relative Stability of Triple-helical Collagen Peptide Models

Herein, we show that current methodologies in atomistic simulations can yield reliable standard free energy values in aqueous solution for the transition from the dissociated monomeric form to the triple-helix state of collagen model peptides. The calculations are performed on a prototypical highly stable triple-helical peptide, [(Pro-Hyp-Gly)10]3 (POG10), and on the so-called T3-785 triple-helix mimicking a fragment from the type III human collagen, which is more thermally labile. On the basis of extensive MD simulations in explicit solvent followed by molecular-mechanical and electrostatic Poisson-Boltzmann calculations complemented with an accurate estimation of the nonpolar contributions to solvation, the computed free energy change for the aggregation processes of the POG10 and T3-785 peptides leading to their triple-helices is -6.6 and -6.1 kcal/mol, respectively. For POG10, this value is in agreement with differential scanning calorimetric data. However, it is shown that conformational entropy, which is estimated by means of an expansion of mutual information functions, preferentially destabilizes the triple-helical state of T3-785 by around 4.6 kcal/mol, thus explaining its lower thermal stability. Altogether, our computational results allow us to ascertain, for the first time, the actual thermodynamic forces controlling the absolute and relative stability of collagen model peptides.

Multibody local approximation: Application to conformational entropy calculations on biomolecules

Multibody type expansions like mutual information expansions are widely used for computing or analyzing properties of large composite systems. The power of such expansions stems from their generality. Their weaknesses, however, are the large computational cost of including high order terms due to the combinatorial explosion and the fact that truncation errors do not decrease strictly with the expansion order. Herein, we take advantage of the redundancy of multibody expansions in order to derive an efficient reformulation that captures implicitly all-order correlation effects within a given cutoff, avoiding the combinatory explosion. This approach, which is cutoff dependent rather than order dependent, keeps the generality of the original expansions and simultaneously mitigates their limitations provided that a reasonable cutoff can be used. An application of particular interest can be the computation of the conformational entropy of flexible peptide molecules from molecular dynamics trajectories. By combining the multibody local estimations of conformational entropy with average values of the rigid-rotor and harmonic-oscillator entropic contributions, we obtain by far a tighter upper bound of the absolute entropy than the one obtained by the broadly used quasi-harmonic method.

Accurate Estimation of Protein Folding and Unfolding Times: Beyond Markov State Models.

Because standard molecular dynamics (MD) simulations are unable to access time scales of interest in complex biomolecular systems, it is common to “stitch together” information from multiple shorter trajectories using approximate Markov state model (MSM) analysis. However, MSMs may require significant tuning and can yield biased results. Here, by analyzing some of the longest protein MD data sets available (>100 μs per protein), we show that estimators constructed based on exact non-Markovian (NM) principles can yield significantly improved mean first-passage times (MFPTs) for protein folding and unfolding. In some cases, MSM bias of more than an order of magnitude can be corrected when identical trajectory data are reanalyzed by non-Markovian approaches. The NM analysis includes “history” information, higher order time correlations compared to MSMs, that is available in every MD trajectory. The NM strategy is insensitive to fine details of the states used and works well when a fine time-discretization (i.e., small “lag time”) is used.

Estimating First Passage Time Distributions from Weighted Ensemble Simulations and non‐Markovian Analyses

First‐passage times (FPTs) are widely used to characterize stochastic processes such as chemical reactions, protein folding, diffusion processes or triggering a stock option. In previous work (Suarez et al., JCTC 2014;10:2658‐2667), we demonstrated a non‐Markovian analysis approach that, with a sufficient subset of history information, yields unbiased mean first‐passage times from weighted‐ensemble (WE) simulations. The estimation of the distribution of the first‐passage times is, however, a more ambitious goal since it cannot be obtained by direct observation in WE trajectories. Likewise, a large number of events would be required to make a good estimation of the distribution from a regular “brute force” simulation. Here, we show how the previously developed non‐Markovian analysis can generate approximate, but highly accurate, FPT distributions from WE data. The analysis can also be applied to any other unbiased trajectories, such as from standard molecular dynamics simulations. The present study employs a range of systems with independent verification of the distributions to demonstrate the success and limitations of the approach. By comparison to a standard Markov analysis, the non‐Markovian approach is less sensitive to the user‐defined discretization of configuration space.

The role of CuCl on the mechanism of dibenzo-p-dioxin formation from poly-chlorophenol precursors: A computational study.

A computational study is performed for the elucidation of the role played by CuCl in the condensation of two polychlorophenol molecules to yield PCDDs. The mechanism found consists of six sequential steps, which allow the final recuperation of the CuCl molecule, and applies for phenol molecules with an ortho chlorine. In the temperature range of 453-473 K (previously reported as adequate to diminish PCDDs formation in the post-combustion area), CuCl is able to softly retain chlorophenol molecules, mainly those less chlorinated. After a first HCl release, Cu(I) remains bonded to phenol oxygen atom, thus avoiding the formation of phenoxy radicals and the subsequent radical processes. A temperature raise up to 1200 K destabilizes the initial CuCl-chlorophenol complexes and causes that the rate limiting step change from the formation of the first oxygen bridge to HCl elimination. It has been checked that tetra and penta-chlorophenols undergo essentially the same reaction process of 2-chlorophenol. In view of our results and trying to arrive at a practical way to diminish the rate of formation of PCDDs, we propose that an extra addition of powdered CuCl to the post-combustion zone, cooled down to temperatures lower than 473 K, could act as an inhibitor in the formation of these pollutants.