Path-breaking schemes for nonequilibrium free energy calculations

Seeding the multi-dimensional nonequilibrium pulling for Hamiltonian variation: indirect nonequilibrium free energy simulations at QM levels

The combination of free energy simulations in the alchemical and configurational spaces provides a feasible route to access the thermodynamic profiles under a computationally demanding target Hamiltonian. Normally, due to the significant differences between the computational cost of ab initio quantum mechanics (QM) calculations and those of semi-empirical quantum mechanics (SQM) and molecular mechanics (MM), this indirect method could be used to obtain the QM thermodynamics by combining the SQM or MM results and the SQM-to-QM or MM-to-QM corrections. In our previous work, a multi-dimensional nonequilibrium pulling framework for Hamiltonian variations was introduced based on bidirectional pulling and bidirectional reweighting. The method performs nonequilibrium free energy simulations in the configurational space to obtain the thermodynamic profile along the conformational change pathway under a selected computationally efficient Hamiltonian, and uses the nonequilibrium alchemical method to correct or perturb the thermodynamic profile to that under the target Hamiltonian. The BAR-based method is designed to achieve the best generality and transferability and thus leads to modest (∼20 fold) speedup. In this work, we explore the possibility of further accelerating the nonequilibrium free energy simulation by employing unidirectional pulling and using the selection criterion to obtain the initial configurations used to initiate nonequilibrium trajectories following the idea of adaptive steered molecular dynamics (ASMD). A single initial condition is used to seed the whole multi-dimensional nonequilibrium free energy simulation and the sampling is performed fully in the nonequilibrium ensemble. Introducing very short ps-length equilibrium sampling to grab more initial seeds could also be helpful. The ASMD scheme estimates the free energy difference with the unidirectional exponential average (EXP), but it does not follow exactly the requirements of the EXP estimator. Another deficiency of the seeding simulation is the inherently sequential or serial pulling due to the inter-segment dependency, which triggers some problems in the parallelizability of the simulation. Numerical tests are performed to grasp some insights and guidelines for using this selection-criterion-based ASMD scheme. The presented selection-criterion-based multi-dimensional ASMD scheme follows the same perturbation network of the BAR-based method, and thus could be used in various Hamiltonian-variation cases.

Conformational Fluctuations in GTP-Bound K-Ras: A Metadynamics Perspective with Harmonic Linear Discriminant Analysis

Biomacromolecules often undergo significant conformational rearrangements during function. In proteins, these motions typically consist in nontrivial, concerted rearrangement of multiple flexible regions. Mechanistic, thermodynamics, and kinetic predictions can be obtained via molecular dynamics simulations, provided that the simulation time is at least comparable to the relevant time scale of the process of interest. Because of the substantial computational cost, however, plain MD simulations often have difficulty in obtaining sufficient statistics for converged estimates, requiring the use of more-advanced techniques. Central in many enhanced sampling methods is the definition of a small set of relevant degrees of freedom (collective variables) that are able to describe the transitions between different metastable states of the system. The harmonic linear discriminant analysis (HLDA) has been shown to be useful for constructing low-dimensional collective variables in various complex systems. Here, we apply HLDA to study the free-energy landscape of a monomeric protein around its native state. More precisely, we study the K-Ras protein bound to GTP, focusing on two flexible loops and on the region associated with oncogenic mutations. We perform microsecond-long biased simulations on the wild type and on G12C, G12D, G12 V mutants, describe the resulting free-energy landscapes, and compare our predictions with previous experimental and computational studies. The fast interconversion between open and closed macroscopic states and their similar thermodynamic stabilities are observed. The mutation-induced effects include the alternations of the relative stabilities of different conformational states and the introduction of many microscopic metastable states. Together, our results demonstrate the applicability of the HLDA-based protocol for the conformational sampling of multiple flexible regions in folded proteins.

On the Consistency of the Exfoliation Free Energy of Graphenes by Molecular Simulations

Monolayer graphene is now produced at significant yields, by liquid phase exfoliation of graphites in solvents. This has increased the interest in molecular simulation studies to give new insights in the field. We use decoupling simulations to compute the exfoliation free energy of graphenes in a liquid environment. Starting from a bilayer graphene configuration, we decouple the Van der Waals interactions of a graphene monolayer in the presence of saline water. Then, we introduce the monolayer back into water by coupling its interactions with water molecules and ions. A different approach to compute the graphene exfoliation free energy is to use umbrella sampling. We apply umbrella sampling after pulling the graphene monolayer on the shear direction up to a distance from a bilayer. We show that the decoupling and umbrella methods give highly consistent free energy results for three bilayer graphene samples with different size. This strongly suggests that the systems in both methods remain closely in equilibrium as we move between the states before and after the exfoliation. Therefore, the amount of nonequilibrium work needed to peel the two layers apart is minimized efficiently.

Binding Thermodynamics and Interaction Patterns of Inhibitor-Major Urinary Protein-I Binding from Extensive Free-Energy Calculations: Benchmarking AMBER Force Fields

Mouse major urinary protein (MUP) plays a key role in the pheromone communication system. The one-end-closed β-barrel of MUP-I forms a small, deep, and hydrophobic central cavity, which could accommodate structurally diverse ligands. Previous computational studies employed old protein force fields and short simulation times to determine the binding thermodynamics or investigated only a small number of structurally similar ligands, which resulted in sampled regions far from the experimental structure, nonconverged sampling outcomes, and limited understanding of the possible interaction patterns that the cavity could produce. In this work, extensive end-point and alchemical free-energy calculations with advanced protein force fields were performed to determine the binding thermodynamics of a series of MUP-inhibitor systems and investigate the inter- and intramolecular interaction patterns. Three series of inhibitors with a total of 14 ligands were simulated. We independently simulated the MUP-inhibitor complexes under two advanced AMBER force fields. Our benchmark test showed that the advanced AMBER force fields including AMBER19SB and AMBER14SB provided better descriptions of the system, and the backbone root-mean-square deviation (RMSD) was significantly lowered compared with previous computational studies with old protein force fields. Surprisingly, although the latest AMBER force field AMBER19SB provided better descriptions of various observables, it neither improved the binding thermodynamics nor lowered the backbone RMSD compared with the previously proposed and widely used AMBER14SB. The older but widely used AMBER14SB actually achieved better performance in the prediction of binding affinities from the alchemical and end-point free-energy calculations. We further analyzed the protein-ligand interaction networks to identify important residues stabilizing the bound structure. Six residues including PHE38, LEU40, PHE90, ALA103, LEU105, and TYR120 were found to contribute the most significant part of protein-ligand interactions, and 10 residues were found to provide favorable interactions stabilizing the bound state. The two AMBER force fields gave extremely similar interaction networks, and the secondary structures also showed similar behavior. Thus, the intra- and intermolecular interaction networks described with the two AMBER force fields are similar. Therefore, AMBER14SB could still be the default option in free-energy calculations to achieve highly accurate binding thermodynamics and interaction patterns.

Sulfur-substitution-induced base flipping in the DNA duplex

Base flipping is widely observed in a number of important biological processes. The genetic codes deposited inside the DNA duplex become accessible to external agents upon base flipping. The sulfur substitution of guanine leads to thioguanine, which alters the thermodynamic stability of the GC base pairs and the GT mismatches. Experimental studies conclude that the sulfur substitution decreases the lifetime of the GC base pair. In this work, under three AMBER force fields for nucleotide systems, we firstly performed equilibrium and nonequilibrium free energy simulations to investigate the variation of the thermodynamic profiles in base flipping upon sulfur substitution. It is found that the bsc0 modification, the bsc1 modification and the OL15 modification of AMBER force fields are able to qualitatively describe the sulfur-substitution dependent behavior of the thermodynamics. However, only the two last-generation AMBER force fields are able to provide quantitatively correct predictions. The second computational study on the sulfur substitutions focused on the relative stability of the S6G-C base pair and the S6G-T mismatch. Two conflicting experimental observations were reported by the same authors. One suggested that the S6G-C base pair was more stable, while the other concludes that the S6G-T mismatch was more stable. We answered this question by constructing the free energy profiles along the base flipping pathway computationally.

BAR-based optimum adaptive steered MD for configurational sampling

The equilibrium and nonequilibrium adaptive alchemical free energy simulation methods optimum Bennett's acceptance ratio and optimum crooks' equation (OCE), based on the statistically optimal bidirectional reweighting estimator named Bennett's Acceptance Ratio or Crooks' equation, perform initial sampling in the staging alchemical transformation and then determine the importance rank of different states via the time-derivative of the variance. The method is proven to give speedups compared with the equal time rule. In the current work, we extend the time derivative of variance guided adaptive sampling method to the configurational space, falling in the term of steered MD (SMD). The SMD approach biasing physically meaningful collective variable (CV) such as one dihedral or one distance to pulling the system from one conformational state to another. By minimizing the variance of the free energy differences along the pathway in an optimized way, a new type of adaptive SMD (ASMD) is introduced. As exhibits in the alchemical case, this adaptive sampling method outperforms the traditional equal-time SMD in nonequilibrium stratification. Also, the method gives much more efficient calculation of potential of mean force than the selection criterion-based ASMD scheme, which is proven to be more efficient than traditional SMD. The OCE workflow is periodicity-of-CV dependent while ASMD is not. The performance is demonstrated in a dihedral flipping case and two distance pulling cases, accounting for periodic and nonperiodic CVs, respectively. © 2019 Wiley Periodicals, Inc.

Correspondence between light-absorption spectrum and nonequilibrium work distribution as a mean to access free energy differences between electronic states

The problem of recovering the free energy difference between two electronic states has been investigated by Frezzato [Chem. Phys. Lett. 533, 106 (2012)], exploring the equivalence between light-absorption spectra and work distribution, hence opening to the application of a spectroscopic version of the Jarzynski equality (JE) [Phys. Rev. Lett. 78, 2690 (1997)]. Here, assuming the validity of the time-dependent perturbation theory, we demonstrate that such equivalence does not lead to the known form of the JE. This is ascribed to the fact that light-absorption processes cannot be described as stochastic processes. To emphasize such an aspect, we devise a stochastic model for the UV-vis (ultraviolet and visible) absorption, suitable for determining the free energy difference between two generic quantum manifolds in a JE-like fashion. However, the model would require explicit knowledge of the transition dipole moments, which are in general not available. Nonetheless, we derive a spectroscopic version of the JE that allows us to recover the free energy difference between the ground and an excited electronic state when the latter state is the only one observed in the spectrum.

Extensive Assessment of Various Computational Methods for Aspartate's pKa Shift

A series of computational methods for pK_a shift prediction are extensively tested on a set of benchmark protein systems, aiming at identifying pitfalls and evaluating their performance on high variants. Including 19 ASP residues in 10 protein systems, the benchmark set consists of both residues with highly shifted pK_a values as well as those varying little from the reference value, with an experimental RMS free energy differences of 2.49 kcal/mol with respect to blocked amino acid, namely the RMS pK_a shift being 1.82 pK_a units. The constant pH molecular dynamics (MD), alchemical methods, PROPKA3.1, and multiconformation continuum electrostatics give RMSDs of 1.52, 2.58, 1.37, and 3.52 pK_a units, respectively, on the benchmark set. The empirical scoring method is the most accurate one with extremely low computational cost, and the pH-dependent model is also able to provide accurate results, while the accuracy of MD sampling incorporating alchemical free energy simulation is prohibited by convergence achievement and the performance of conformational search incorporating multiconformation continuum electrostatics is bad. Former research works did not define statistical uncertainty with care and yielded the questionable conclusion that alchemical methods perform well in most benchmarks. In this work the traditional alchemical methods are thoroughly tested for high variants. We also performed the first application of nonequilibrium alchemical methods to the pK_a cases.

Evaluation of protein-ligand affinity prediction using steered molecular dynamics simulations

In computational drug design, ranking a series of compound analogs in a manner that is consistent with experimental affinities remains a challenge. In this study, we evaluated the prediction of protein-ligand binding affinities using steered molecular dynamics simulations. First, we investigated the appropriate conditions for accurate predictions in these simulations. A conic harmonic restraint was applied to the system for efficient sampling of work values on the ligand unbinding pathway. We found that pulling velocity significantly influenced affinity predictions, but that the number of collectable trajectories was less influential. We identified the appropriate pulling velocity and collectable trajectories for binding affinity predictions as 1.25 Å/ns and 100, respectively, and these parameters were used to evaluate three target proteins (FK506 binding protein, trypsin, and cyclin-dependent kinase 2). For these proteins using our parameters, the accuracy of affinity prediction was higher and more stable when Jarzynski's equality was employed compared with the second-order cumulant expansion equation of Jarzynski's equality. Our results showed that steered molecular dynamics simulations are effective for predicting the rank order of ligands; thus, they are a potential tool for compound selection in hit-to-lead and lead optimization processes.

Jarzynski matrix equality: Calculating the free-energy difference by nonequilibrium simulations with an arbitrary initial distribution

The Jarzynski equality (JE) method, which relates the work of a nonequilibrium process to the free-energy difference between its initial and final states, provides an efficient way to calculate free energies of thermodynamic systems in simulations or experiments. However, more extensive applications of the JE are hindered by the requirement that the initial state must be in equilibrium. In this work we extend the JE method to be the Jarzynski matrix equality (JME) method, which relates the work of trajectories connecting metastable conformational regions to their local free energies, and thus we can estimate the free energy from the nonequilibrium trajectories starting from an almost arbitrary initial distribution. We then apply the JME to toy models, Lennard-Jones fluids, and polymer chain models, demonstrating its efficiency in free-energy calculations with satisfactory accuracy. The JME extends the applicability of the nonequilibrium methods to complex systems whose initial equilibrium states are difficult to reach.

Elastic Barrier Dynamical Freezing in Free Energy Calculations: A Way To Speed Up Nonequilibrium Molecular Dynamics Simulations by Orders of Magnitude

An important issue concerning computer simulations addressed to free energy estimates via nonequilibrium work theorems, such as the Jarzynski equality [Phys. Rev. Lett. 1997, 78, 2690], is the computational effort required to achieve results with acceptable accuracy. In this respect, the dynamical freezing approach [Phys. Rev. E 2009, 80, 041124] has been shown to improve the efficiency of this kind of simulations, by blocking the dynamics of particles located outside an established mobility region. In this report, we show that dynamical freezing produces a systematic spurious decrease of the particle density inside the mobility region. As a consequence, the requirements to apply nonequilibrium work theorems are only approximately met. Starting from these considerations, we have developed a simulation scheme, called "elastic barrier dynamical freezing", according to which a stiff potential-energy barrier is enforced at the boundaries of the mobility region, preventing the particles from leaving this region of space during the nonequilibrium trajectories. The method, tested on the calculation of the distance-dependent free energy of a dimer immersed into a Lennard-Jones fluid, provides an accuracy comparable to the conventional steered molecular dynamics, with a computational speedup exceeding a few orders of magnitude.

Nonequilibrium Candidate Monte Carlo Simulations with Configurational Freezing Schemes

Nonequilibrium Candidate Monte Carlo simulation [Nilmeier et al., Proc. Natl. Acad. Sci. U.S.A. 2011, 108, E1009-E1018] is a tool devised to design Monte Carlo moves with high acceptance probabilities that connect uncorrelated configurations. Such moves are generated through nonequilibrium driven dynamics, producing candidate configurations accepted with a Monte Carlo-like criterion that preserves the equilibrium distribution. The probability of accepting a candidate configuration as the next sample in the Markov chain basically depends on the work performed on the system during the nonequilibrium trajectory and increases with decreasing such a work. It is thus strategically relevant to find ways of producing nonequilibrium moves with low work, namely moves where dissipation is as low as possible. This is the goal of our methodology, in which we combine Nonequilibrium Candidate Monte Carlo with Configurational Freezing schemes developed by Nicolini et al. (J. Chem. Theory Comput. 2011, 7, 582-593). The idea is to limit the configurational sampling to particles of a well-established region of the simulation sample, namely the region where dissipation occurs, while leaving fixed the other particles. This allows to make the system relaxation faster around the region perturbed by the finite-time switching move and hence to reduce the dissipated work, eventually enhancing the probability of accepting the generated move. Our combined approach enhances significantly configurational sampling, as shown by the case of a bistable dimer immersed in a dense fluid.

A Self-Adaptive Steered Molecular Dynamics Method Based on Minimization of Stretching Force Reveals the Binding Affinity of Protein-Ligand Complexes

Binding affinity prediction of protein–ligand complexes has attracted widespread interest. In this study, a self-adaptive steered molecular dynamics (SMD) method is proposed to reveal the binding affinity of protein–ligand complexes. The SMD method is executed through adjusting pulling direction to find an optimum trajectory of ligand dissociation, which is realized by minimizing the stretching force automatically. The SMD method is then used to simulate the dissociations of 19 common protein–ligand complexes which are derived from two homology families, and the binding free energy values are gained through experimental techniques. Results show that the proposed SMD method follows a different dissociation pathway with lower a rupture force and energy barrier when compared with the conventional SMD method, and further analysis indicates the rupture forces of the complexes in the same protein family correlate well with their binding free energy, which reveals the possibility of using the proposed SMD method to identify the active ligand.

Towards bulk thermodynamics via non-equilibrium methods: gaseous methane as a case study

The equation of state of bulk materials is achieved via thermodynamic derivatives of the free energy yielded by nonequilibrium transformations and Jarzynski equality.

Thermodynamics of Deca-alanine Folding in Water

The determination of the folding dynamics of polypeptides and proteins is critical in characterizing their functions in biological systems. Numerous computational models and methods have been developed for studying structure formation at the atomic level. Due to its small size and simple structure, deca-alanine is used as a model system in molecular dynamics (MD) simulations. The free energy of unfolding in vacuum has been studied extensively using the end-to-end distance of the peptide as the reaction coordinate. However, few studies have been conducted in the presence of explicit solvent. Previous results show a significant decrease in the free energy of extended conformations in water, but the α-helical state is still notably favored over the extended state. Although sufficient in vacuum, we show that end-to-end distance is incapable of capturing the full complexity of deca-alanine folding in water. Using α-helical content as a second reaction coordinate, we deduce a more descriptive free-energy landscape, one which reveals a second energy minimum in the extended conformations that is of comparable free energy to the α-helical state. Equilibrium simulations demonstrate the relative stability of the extended and α-helical states in water as well as the transition between the two states. This work reveals both the necessity and challenge of determining a proper reaction coordinate to fully characterize a given process.