Efficient use of nonequilibrium measurement to estimate free energy differences for molecular systems

Methods for Classical-Mechanical Molecular Simulation in Chemistry: Achievements, Limitations, Perspectives

More than a half century ago it became feasible to simulate, using classical-mechanical equations of motion, the dynamics of molecular systems on a computer. Since then classical-physical molecular simulation has become an integral part of chemical research. It is widely applied in a variety of branches of chemistry and has significantly contributed to the development of chemical knowledge. It offers understanding and interpretation of experimental results, semiquantitative predictions for measurable and nonmeasurable properties of substances, and allows the calculation of properties of molecular systems under conditions that are experimentally inaccessible. Yet, molecular simulation is built on a number of assumptions, approximations, and simplifications which limit its range of applicability and its accuracy. These concern the potential-energy function used, adequate sampling of the vast statistical-mechanical configurational space of a molecular system and the methods used to compute particular properties of chemical systems from statistical-mechanical ensembles. During the past half century various methodological ideas to improve the efficiency and accuracy of classical-physical molecular simulation have been proposed, investigated, evaluated, implemented in general simulation software or were abandoned. The latter because of fundamental flaws or, while being physically sound, computational inefficiency. Some of these methodological ideas are briefly reviewed and the most effective methods are highlighted. Limitations of classical-physical simulation are discussed and perspectives are sketched.

Convergence to the Asymptotic Large Deviation Limit

Large deviation theory offers a powerful and general statistical framework to study the asymptotic dynamical properties of rare events. The application of the formalism to concrete experimental situations is, however, often restricted by finite statistics. Data might not suffice to reach the asymptotic regime or judge whether large deviation estimators converge at all. We here experimentally study the large deviation properties of the stochastic work and heat of a levitated nanoparticle subjected to nonequilibrium feedback control. This setting allows us to determine for each quantity the convergence domain of the large deviation estimators using a criterion that does not require the knowledge of the probability distribution. By extracting both the asymptotic exponential decay and the subexponential prefactors, we demonstrate that singular prefactors significantly restrict the convergence characteristics close to the singularity. Our results provide unique insight into the approach to the asymptotic large deviation limit and underscore the pivotal role of singular prefactors.

Path separation of dissipation-corrected targeted molecular dynamics simulations of protein-ligand unbinding

Protein-ligand (un)binding simulations are a recent focus of biased molecular dynamics simulations. Such binding and unbinding can occur via different pathways in and out of a binding site. Here, we present a theoretical framework on how to compute kinetics along separate paths and on how to combine the path-specific rates into global binding and unbinding rates for comparison with experimental results. Using dissipation-corrected targeted molecular dynamics in combination with temperature-boosted Langevin equation simulations [S. Wolf et al., Nat. Commun. 11, 2918 (2020)] applied to a two-dimensional model and the trypsin-benzamidine complex as test systems, we assess the robustness of the procedure and discuss the aspects of its practical applicability to predict multisecond kinetics of complex biomolecular systems.

Thermodynamic fluctuation theorems govern human sensorimotor learning

The application of thermodynamic reasoning in the study of learning systems has a long tradition. Recently, new tools relating perfect thermodynamic adaptation to the adaptation process have been developed. These results, known as fluctuation theorems, have been tested experimentally in several physical scenarios and, moreover, they have been shown to be valid under broad mathematical conditions. Hence, although not experimentally challenged yet, they are presumed to apply to learning systems as well. Here we address this challenge by testing the applicability of fluctuation theorems in learning systems, more specifically, in human sensorimotor learning. In particular, we relate adaptive movement trajectories in a changing visuomotor rotation task to fully adapted steady-state behavior of individual participants. We find that human adaptive behavior in our task is generally consistent with fluctuation theorem predictions and discuss the merits and limitations of the approach.

A modified Jarzynski free-energy estimator to eliminate non-conservative forces and its application in nanoparticle-membrane interactions

Computational methods to understand interactions in bio-complex systems are however limited to time-scales typically much shorter than in Nature. For example, on the nanoscale level, interactions between nanoparticles (NPs)/molecules/peptides and membranes are central in complex biomolecular processes such as membrane-coated NPs or cellular uptake. This can be remedied by the application of e.g. Jarzynski's equality where thermodynamic properties are extracted from non-equilibrium simulations. Although, the out of equilibrium work leads to non-conservative forces. We here propose a correction Pair Forces method, that removes these forces. Our proposed method is based on the calculation of pulling forces in backward and forward directions for the Jarzynski free-energy estimator using steered molecular dynamics simulation. Our results show that this leads to much improvement for NP-membrane translocation free energies. Although here we have demonstrated the application of the method in molecular dynamics simulation, it could be applied for experimental approaches.

Reweighted Jarzynski Sampling: Acceleration of Rare Events and Free Energy Calculation with a Bias Potential Learned from Nonequilibrium Work

We introduce a simple enhanced sampling approach for the calculation of free energy differences and barriers along a one-dimensional reaction coordinate. First, a small number of short nonequilibrium simulations are carried out along the reaction coordinate, and the Jarzynski equality is used to learn an approximate free energy surface from the nonequilibrium work distribution. This free energy estimate is represented in a compact form as an artificial neural network and used as an external bias potential to accelerate rare events in a subsequent molecular dynamics simulation. The final free energy estimate is then obtained by reweighting the equilibrium probability distribution of the reaction coordinate sampled under the influence of the external bias. We apply our reweighted Jarzynski sampling recipe to four processes of varying scales and complexities─spanning chemical reaction in the gas phase, pair association in solution, and droplet nucleation in supersaturated vapor. In all cases, we find reweighted Jarzynski sampling to be a very efficient strategy, resulting in rapid convergence of the free energy to high precision.

Single molecule force spectroscopy of a streptomycin-binding RNA aptamer: An out-of-equilibrium molecular dynamics study

Here, we investigate the unfolding behavior of a streptomycin-binding ribonucleic acid (RNA) aptamer under application of force in shear geometry. Using Langevin out-of-equilibrium simulations to emulate the single-molecule force spectroscopy (SMFS) experiment, we were able to understand the hierarchical unfolding process that occurs in the RNA molecule under application of stretching force and the influence of streptomycin modifying this unfolding. Subsequently, the application of the Jarzynski equality to the force profiles obtained in the pulling simulations shows that the free energies for individual systems and the difference of unfolding free energy upon streptomycin binding to the RNA free aptamer are in fair agreement with the experimental values, obtained through SMFS by Nick et al. [J. Phys. Chem. B 120, 6479 (2016)].

Lessons learned about steered molecular dynamics simulations and free energy calculations

The calculation of free energy profiles is central in understanding differential enzymatic activity, for instance, involving chemical reactions that require QM-MM tools, ligand migration, and conformational rearrangements that can be modeled using classical potentials. The use of steered molecular dynamics (sMD) together with the Jarzynski equality is a popular approach in calculating free energy profiles. Here, we first briefly review the application of the Jarzynski equality to sMD simulations, then revisit the so-called stiff-spring approximation and the consequent expectation of Gaussian work distributions and, finally, reiterate the practical utility of the second-order cumulant expansion, as it coincides with the parametric maximum-likelihood estimator in this scenario. We illustrate this procedure using simulations of CO, both in aqueous solution and in a carbon nanotube as a model system for biologically relevant nanoheterogeneous environments. We conclude the use of the second-order cumulant expansion permits the use of faster pulling velocities in sMD simulations, without introducing bias due to large dispersion in the non-equilibrium work distribution.

Organic molecule fluorescence as an experimental test-bed for quantum jumps in thermodynamics

We demonstrate with an experiment how molecules are a natural test bed for probing fundamental quantum thermodynamics. Single-molecule spectroscopy has undergone transformative change in the past decade with the advent of techniques permitting individual molecules to be distinguished and probed. We demonstrate that the quantum Jarzynski equality for heat is satisfied in this set-up by considering the time-resolved emission spectrum of organic molecules as arising from quantum jumps between states. This relates the heat dissipated into the environment to the free energy difference between the initial and final state. We demonstrate also how utilizing the quantum Jarzynski equality allows for the detection of energy shifts within a molecule, beyond the relative shift.

Direction-dependent force-induced dissociation dynamics of an entropic-driven lock-and-key assembly

The unbinding dynamics of a nanosized sphere-and-cavity assembly under the pulling of constant force and constant loading rate is explored by dissipative particle dynamics simulations. The formation of this matched lock-and-key pair in a polymer solution is driven by the depletion attraction. The two-dimensional free energy landscape U(x,z) associated with this assembly is constructed. Our results indicate that the unbinding pathway along the orientation of the assembly is unfavorable due to the relatively high energy barrier compared to that along the tortuous minimum path whose energy barrier is not high. It is also found that the dissociation rate depends on the direction of the external force (θ) with respect to the assembly orientation. The presence of the force component perpendicular to the assembly orientation can reduce the bond lifetime significantly by driving the key particle to approach the minimum path. Moreover, the dissociation dynamics can be facilitated even by a pushing force compared to the spontaneous dissociation (without forces). To elucidate the effective pathway under pulling, the escaping position is analyzed and its mean direction with respect to the assembly orientation rises generally with increasing θ, revealing that the presence of the force component along the minimum pathway is helpful. The importance of the direction of the external pulling has been demonstrated in our simple system. Therefore, this effect should be considered in more complicated unbinding experiments.

Opposing Intermolecular Tuning of Ca2+ Affinity for Calmodulin by Neurogranin and CaMKII Peptides

We investigated the impact of bound calmodulin (CaM)-target compound structure on the affinity of calcium (Ca²⁺) by integrating coarse-grained models and all-atomistic simulations with nonequilibrium physics. We focused on binding between CaM and two specific targets, Ca²⁺/CaM-dependent protein kinase II (CaMKII) and neurogranin (Ng), as they both regulate CaM-dependent Ca²⁺ signaling pathways in neurons. It was shown experimentally that Ca²⁺/CaM (holoCaM) binds to the CaMKII peptide with overwhelmingly higher affinity than Ca²⁺-free CaM (apoCaM); the binding of CaMKII peptide to CaM in return increases the Ca²⁺ affinity for CaM. However, this reciprocal relation was not observed in the Ng peptide (Ng_13–49), which binds to apoCaM or holoCaM with binding affinities of the same order of magnitude. Unlike the holoCaM-CaMKII peptide, whose structure can be determined by crystallography, the structural description of the apoCaM-Ng_13–49 is unknown due to low binding affinity, therefore we computationally generated an ensemble of apoCaM-Ng_13–49 structures by matching the changes in the chemical shifts of CaM upon Ng_13–49 binding from nuclear magnetic resonance experiments. Next, we computed the changes in Ca²⁺ affinity for CaM with and without binding targets in atomistic models using Jarzynski’s equality. We discovered the molecular underpinnings of lowered affinity of Ca²⁺ for CaM in the presence of Ng_13–49 by showing that the N-terminal acidic region of Ng peptide pries open the β-sheet structure between the Ca²⁺ binding loops particularly at C-domain of CaM, enabling Ca²⁺ release. In contrast, CaMKII peptide increases Ca²⁺ affinity for the C-domain of CaM by stabilizing the two Ca²⁺ binding loops. We speculate that the distinctive structural difference in the bound complexes of apoCaM-Ng_13–49 and holoCaM-CaMKII delineates the importance of CaM’s progressive mechanism of target binding on its Ca²⁺ binding affinities.

The role of intramolecular nonbonded interaction and angle sampling in single-step free energy perturbation

Single-step free energy perturbation (sFEP) has often been proposed as an efficient tool for a quick free energy scan due to its straightforward protocol and the ability to recycle an existing molecular dynamics trajectory for free energy calculations. Although sFEP is expected to fail when the sampling of a system is inefficient, it is often expected to hold for an alchemical transformation between ligands with a moderate difference in their sizes, e.g., transforming a benzene into an ethylbenzene. Yet, exceptions were observed in calculations for anisole and methylaniline, which have similar physical sizes as ethylbenzene. In this study, we show that such exceptions arise from the sampling inefficiency on an unexpected rigid degree of freedom, namely, the bond angle θ. The distributions of θ differ dramatically between two end states of a sFEP calculation, i.e., the conformation of the ligand changes significantly during the alchemical transformation process. Our investigation also reveals the interrelation between the ligand conformation and the intramolecular nonbonded interactions. This knowledge suggests a best combination of the ghost ligand potential and the dual topology setting, which improves the accuracy in a single reference sFEP calculation by bringing down its error from around 5k_BT to k_BT.

Molecular Dynamics Simulations of the Human Glucose Transporter GLUT1

Glucose transporters (GLUTs) provide a pathway for glucose transport across membranes. Human GLUTs are implicated in devastating diseases such as heart disease, hyper- and hypo-glycemia, type 2 diabetes and cancer. The human GLUT1 has been recently crystalized in the inward-facing open conformation. However, there is no other structural information for other conformations. The X-ray structures of E. coli Xylose permease (XylE), a glucose transporter homolog, are available in multiple conformations with and without the substrates D-xylose and D-glucose. XylE has high sequence homology to human GLUT1 and key residues in the sugar-binding pocket are conserved. Here we construct a homology model for human GLUT1 based on the available XylE crystal structure in the partially occluded outward-facing conformation. A long unbiased all atom molecular dynamics simulation starting from the model can capture a new fully opened outward-facing conformation. Our investigation of molecular interactions at the interface between the transmembrane (TM) domains and the intracellular helices (ICH) domain in the outward- and inward-facing conformation supports that the ICH domain likely stabilizes the outward-facing conformation in GLUT1. Furthermore, inducing a conformational transition, our simulations manifest a global asymmetric rocker switch motion and detailed molecular interactions between the substrate and residues through the water-filled selective pore along a pathway from the extracellular to the intracellular side. The results presented here are consistent with previously published biochemical, mutagenesis and functional studies. Together, this study shed light on the structure and functional relationships of GLUT1 in multiple conformational states.

Prior-predictive value from fast-growth simulations: error analysis and bias estimation

Variants of fluctuation theorems recently discovered in the statistical mechanics of nonequilibrium processes may be used for the efficient determination of high-dimensional integrals as typically occurring in Bayesian data analysis. In particular for multimodal distributions, Monte Carlo procedures not relying on perfect equilibration are advantageous. We provide a comprehensive statistical error analysis for the determination of the prior-predictive value (the evidence) in a Bayes problem, building on a variant of the Jarzynski equation. Special care is devoted to the characterization of the bias intrinsic to the method and statistical errors arising from exponential averages. We also discuss the determination of averages over multimodal posterior distributions with the help of a consequence of the Crooks relation. All our findings are verified by extensive numerical simulations of two model systems with bimodal likelihoods.

Correcting for the free energy costs of bond or angle constraints in molecular dynamics simulations

BACKGROUND

Free energy simulations are an important tool in the arsenal of computational biophysics, allowing the calculation of thermodynamic properties of binding or enzymatic reactions. This paper introduces methods to increase the accuracy and precision of free energy calculations by calculating the free energy costs of constraints during post-processing. The primary purpose of employing constraints for these free energy methods is to increase the phase space overlap between ensembles, which is required for accuracy and convergence.

METHODS

The free energy costs of applying or removing constraints are calculated as additional explicit steps in the free energy cycle. The new techniques focus on hard degrees of freedom and use both gradients and Hessian estimation. Enthalpy, vibrational entropy, and Jacobian free energy terms are considered.

RESULTS

We demonstrate the utility of this method with simple classical systems involving harmonic and anharmonic oscillators, four-atomic benchmark systems, an alchemical mutation of ethane to methanol, and free energy simulations between alanine and serine. The errors for the analytical test cases are all below 0.0007kcal/mol, and the accuracy of the free energy results of ethane to methanol is improved from 0.15 to 0.04kcal/mol. For the alanine to serine case, the phase space overlaps of the unconstrained simulations range between 0.15 and 0.9%. The introduction of constraints increases the overlap up to 2.05%. On average, the overlap increases by 94% relative to the unconstrained value and precision is doubled.

CONCLUSIONS

The approach reduces errors arising from constraints by about an order of magnitude. Free energy simulations benefit from the use of constraints through enhanced convergence and higher precision.

GENERAL SIGNIFICANCE

The primary utility of this approach is to calculate free energies for systems with disparate energy surfaces and bonded terms, especially in multi-scale molecular mechanics/quantum mechanics simulations. This article is part of a Special Issue entitled Recent developments of molecular dynamics.

Proton-coupled sugar transport in the prototypical major facilitator superfamily protein XylE

The major facilitator superfamily (MFS) is the largest collection of structurally related membrane proteins that transport a wide array of substrates. The proton-coupled sugar transporter XylE is the first member of the MFS that has been structurally characterized in multiple transporting conformations, including both the outward and inward-facing states. Here we report the crystal structure of XylE in a new inward-facing open conformation, allowing us to visualize the rocker-switch movement of the N-domain against the C-domain during the transport cycle. Using molecular dynamics simulation, and functional transport assays, we describe the movement of XylE that facilitates sugar translocation across a lipid membrane and identify the likely candidate proton-coupling residues as the conserved Asp27 and Arg133. This study addresses the structural basis for proton-coupled substrate transport and release mechanism for the sugar porter family of proteins.

Glucose transporters are a medically important class of membrane proteins often deregulated in diseases such as Type 2 diabetes. Here, Wisedchaisri et al. report the crystal structure of XylE in an inward-facing open conformation to provide a general mechanism of substrate transport for the sugar porter family of proteins.

Practical Aspects of Free-Energy Calculations: A Review

Free-energy calculations in the framework of classical molecular dynamics simulations are nowadays used in a wide range of research areas including solvation thermodynamics, molecular recognition, and protein folding. The basic components of a free-energy calculation, that is, a suitable model Hamiltonian, a sampling protocol, and an estimator for the free energy, are independent of the specific application. However, the attention that one has to pay to these components depends considerably on the specific application. Here, we review six different areas of application and discuss the relative importance of the three main components to provide the reader with an organigram and to make nonexperts aware of the many pitfalls present in free energy calculations.

Molecular dynamics simulations suggest ligand's binding to nicotinamidase/pyrazinamidase

The research on the binding process of ligand to pyrazinamidase (PncA) is crucial for elucidating the inherent relationship between resistance of Mycobacterium tuberculosis and PncA’s activity. In the present study, molecular dynamics (MD) simulation methods were performed to investigate the unbinding process of nicotinamide (NAM) from two PncA enzymes, which is the reverse of the corresponding binding process. The calculated potential of mean force (PMF) based on the steered molecular dynamics (SMD) simulations sheds light on an optimal binding/unbinding pathway of the ligand. The comparative analyses between two PncAs clearly exhibit the consistency of the binding/unbinding pathway in the two enzymes, implying the universality of the pathway in all kinds of PncAs. Several important residues dominating the pathway were also determined by the calculation of interaction energies. The structural change of the proteins induced by NAM’s unbinding or binding shows the great extent interior motion in some homologous region adjacent to the active sites of the two PncAs. The structure comparison substantiates that this region should be very important for the ligand’s binding in all PncAs. Additionally, MD simulations also show that the coordination position of the ligand is displaced by one water molecule in the unliganded enzymes. These results could provide the more penetrating understanding of drug resistance of M. tuberculosis and be helpful for the development of new antituberculosis drugs.

Delineation of the unbinding pathway of α-conotoxin ImI from the α7 nicotinic acetylcholine receptor

α-Conotoxins potently and specifically inhibit isoforms of nicotinic acetylcholine receptors (nAChRs) and are used as molecular probes and as drugs or drug leads. Interactions occurring during binding and unbinding events are linked to binding kinetics, and knowledge of these interactions could help in the development of α-conotoxins as drugs. Here, the unbinding process for the prototypical α-conotoxin ImI/α7-nAChR system was investigated theoretically, and three exit routes were identified using random accelerated molecular dynamics simulations. The route involving the smallest conformation perturbation was further divided into three subpathways, which were studied using steered molecular dynamics simulations. Of the three subpathways, two had better experimental support and lower potential of mean force, indicating that they might be sampled more frequently. Additionally, these subpathways were supported by previous experimental studies. Several pairwise interactions, including a cation-π interaction and charge and hydrogen bond interactions, were identified as potentially playing important roles in the unbinding event.

Investigation of indazole unbinding pathways in CYP2E1 by molecular dynamics simulations

Human microsomal cytochrome P450 2E1 (CYP2E1) can oxidize not only low molecular weight xenobiotic compounds such as ethanol, but also many endogenous fatty acids. The crystal structure of CYP2E1 in complex with indazole reveals that the active site is deeply buried into the protein center. Thus, the unbinding pathways and associated unbinding mechanisms remain elusive. In this study, random acceleration molecular dynamics simulations combined with steered molecular dynamics and potential of mean force calculations were performed to identify the possible unbinding pathways in CYP2E1. The results show that channel 2c and 2a are most likely the unbinding channels of CYP2E1. The former channel is located between helices G and I and the B-C loop, and the latter resides between the region formed by the F-G loop, the B-C loop and the β1 sheet. Phe298 and Phe478 act as the gate keeper during indazole unbinding along channel 2c and 2a, respectively. Previous site-directed mutagenesis experiments also supported these findings.

Probing carbohydrate product expulsion from a processive cellulase with multiple absolute binding free energy methods

Understanding the enzymatic mechanism that cellulases employ to degrade cellulose is critical to efforts to efficiently utilize plant biomass as a sustainable energy resource. A key component of cellulase action on cellulose is product inhibition from monosaccharide and disaccharides in the product site of cellulase tunnel. The absolute binding free energy of cellobiose and glucose to the product site of the catalytic tunnel of the Family 7 cellobiohydrolase (Cel7A) of Trichoderma reesei (Hypocrea jecorina) was calculated using two different approaches: steered molecular dynamics (SMD) simulations and alchemical free energy perturbation molecular dynamics (FEP/MD) simulations. For the SMD approach, three methods based on Jarzynski's equality were used to construct the potential of mean force from multiple pulling trajectories. The calculated binding free energies, -14.4 kcal/mol using SMD and -11.2 kcal/mol using FEP/MD, are in good qualitative agreement. Analysis of the SMD pulling trajectories suggests that several protein residues (Arg-251, Asp-259, Asp-262, Trp-376, and Tyr-381) play key roles in cellobiose and glucose binding to the catalytic tunnel. Five mutations (R251A, D259A, D262A, W376A, and Y381A) were made computationally to measure the changes in free energy during the product expulsion process. The absolute binding free energies of cellobiose to the catalytic tunnel of these five mutants are -13.1, -6.0, -11.5, -7.5, and -8.8 kcal/mol, respectively. The results demonstrated that all of the mutants tested can lower the binding free energy of cellobiose, which provides potential applications in engineering the enzyme to accelerate the product expulsion process and improve the efficiency of biomass conversion.

Molecular dynamics study of a heteroditopic-calix[4]diquinone-assisted transfer of KCl and dopamine through a water-chloroform liquid-liquid interface

The ability of two heteroditopic calix[4]diquinone receptors to transport a KCl ion-pair and a dopamine zwitterion through a water-chloroform interface was investigated via molecular dynamics (MD) simulations. Gas-phase conformational analysis has been carried on KCl and dopamine receptor binding associations and the lowest energy structures found in both cases show that the recognition of KCl and dopamine zwitterion occurs through multiple and cooperative N-H...anion and O...cation bonding interactions, with the receptor adopting equivalent folded conformations stabilized by pi-stacking interactions. The unconstrained MD simulations performed on KCl and dopamine complexes inserted in either the chloroform or water phase revealed that receptors are preferentially located at the interface with the hydrophobic tert-butyl groups of the calix[4]diquinone moiety immersed in the chloroform bulk while the polar anion binding cavity is directed toward the water phase. When the KCl complex is placed in chloroform, the release of the ion-pair occurs only after the first contact with the water interface, being a nonsimultaneous event, with the chloride anion leaving the receptor before the potassium cation. The dopamine, via the -NH(3)(+) binding entity, remains bound to the receptor during the entire time of the MD simulation (10 ns). In contrast, when both complexes were inserted in the water bulk, the full release of KCl and dopamine are fast events. The potentials of mean force (PMFs), associated with the migration of the complexes from chloroform to water through the interface, were calculated from steered molecular dynamics (SMD) simulations. The PMFs for the free KCl and zwitterionic dopamine migrations were also obtained for comparison purposes. The transport of KCl from water to chloroform (the reverse path) mediated by the receptor has a free energy barrier estimated in 6.50 kcal mol(-1), which is 3.0 kcal mol(-1) smaller than that found for the free KCl. The transport of dopamine complex along the reverse path is characterized by downhill energy profile, with a small free energy barrier of 6.56 kcal mol(-1).

Optimum protocol for fast-switching free-energy calculations

Free-energy differences computed from fast-switching simulations or measurements according to the Jarzynski equation are independent of the particular protocol specifying how the control parameter is changed in time. In contrast, the average work carried out on the system as well the accuracy of the resulting free energy strongly depend on the protocol. Recently, Schmiedl and Seifert [Phys. Rev. Lett. 98, 108301 (2007)] found that protocols that minimize the average work for a given duration of the switching process have discrete steps at the beginning and the end. Here we determine numerically the protocols that minimize the statistical error in the free energy estimate and find that such minimum error protocols have similar discrete jumps. Our analysis shows that the reduction in computational effort achieved by the use of steplike protocols can be considerable. Such large savings of computing time, however, typically occur for parameter ranges in which an application of the Jarzynski equation is impractical due to large statistical errors arising from the exponential work average.

Computational insights into the mechanism of ligand unbinding and selectivity of estrogen receptors

Estrogen receptors (ER) belong to the nuclear receptor superfamily, and two subtypes, ERalpha and ERbeta, have been identified to date. The differentiated functions and receptor expressions of ERalpha and ERbeta made it attracted to discover subtype-specified ligands with high selectivity. However, these two subtypes are highly homologous and only two residues differ in the ligand binding pocket. Therefore, the mechanism of ligand selectivity has become an important issue in searching selective ligands of ER subtypes. In this study, steered molecular dynamics simulations were carried out to investigate the unbinding pathways of two selective ERbeta ligands from the binding pocket of both ERalpha and ERbeta, which demonstrated that the pathway between the H11 helix and the H7 approximately H8 loop was the most probable for ligand escaping. Then potentials of mean force for ligands unbinding along this pathway were calculated in order to gain insights into the molecular basis for energetics of ligand unbinding and find clues of ligand selectivity. The results indicated that His524/475 in ERalpha/ERbeta acted as a "gatekeeper" during the ligand unbinding. Especially, the H7 approximately H8 loop of ERbeta acted as a polar "transmitter" that controlled the ligand unbinding from the binding site and contributed to the ligand selectivity. Finally, the mechanism of ligand selectivity of ER subtypes was discussed from a kinetic perspective and suggestions for improving the ligand selectivity of ERbeta were also presented. These findings could be helpful for rational design of highly selective ERbeta ligands.

Absolute free energies and equilibrium ensembles of dense fluids computed from a nondynamic growth method

We demonstrate a nondynamical Monte Carlo method to compute free energies and generate equilibrium ensembles of dense fluids. In this method, based on step-by-step polymer growth algorithms, an ensemble of n+1 particles is obtained from an ensemble of n particles by generating configurations of the n+1st particle. A statistically rigorous resampling scheme is utilized to remove configurations with low weights and to avoid a combinatorial explosion; the free energy is obtained from the sum of the weights. In addition to the free energy, the method generates an equilibrium ensemble of the full system. We consider two different system sizes for a Lennard-Jones fluid and compare the results with conventional Monte Carlo methods.

Assessment of nonequilibrium free energy methods

One of the factors preventing the general application of free energy methods in rational drug design remains the lack of sufficient computational resources. Many nonequilibrium (NE) free energy methods, however, are easily made embarrassingly parallel in comparison to equilibrium methods and may be conveniently run on desktop computers using distributed computing software. In recent years, there has been a proliferation of NE methods, but the general applicability of these approaches has not been determined. In this study, a subset including only those NE methods which are easily parallelised were considered for examination, with a view to their application to the prediction of protein-ligand binding affinities. A number of test systems were examined, including harmonic oscillator (HO) systems and the calculation of relative free energies of hydration of water-methane. The latter system uses identical potentials to the protein ligand case and is therefore an appropriate model system on which methods may be tested. As well as investigating existing protocols, a replica exchange NE approach was developed, which was found to offer advantages over conventional methods. It was found that Rosenbluth-based approaches to optimizing the NE work values used in NE free energy estimates were not consistent in the improvements in accuracy achieved and that, given their computational cost, the simple approach of taking each work value in an unbiased way is to be preferred. Of the two free energy estimators examined, Bennett's acceptance ratio was the most consistent and is, therefore, to be preferred over the Jarzynski estimator. The recommended protocols may be run very efficiently within a distributed computing environment and are of similar accuracy and precision to equilibrium free energy methods.

Bayesian estimates of free energies from nonequilibrium work data in the presence of instrument noise

The Jarzynski equality and the fluctuation theorem relate equilibrium free energy differences to nonequilibrium measurements of the work. These relations extend to single-molecule experiments that have probed the finite-time thermodynamics of proteins and nucleic acids. The effects of experimental error and instrument noise have not been considered previously. Here, we present a Bayesian formalism for estimating free energy changes from nonequilibrium work measurements that compensates for instrument noise and combines data from multiple driving protocols. We reanalyze a recent set of experiments in which a single RNA hairpin is unfolded and refolded using optical tweezers at three different rates. Interestingly, the fastest and farthest-from-equilibrium measurements contain the least instrumental noise and, therefore, provide a more accurate estimate of the free energies than a few slow, more noisy, near-equilibrium measurements. The methods we propose here will extend the scope of single-molecule experiments; they can be used in the analysis of data from measurements with atomic force microscopy, optical, and magnetic tweezers.

Protein-protein interaction investigated by steered molecular dynamics: the TCR-pMHC complex

We present a novel steered molecular dynamics scheme to induce the dissociation of large protein-protein complexes. We apply this scheme to study the interaction of a T cell receptor (TCR) with a major histocompatibility complex (MHC) presenting a peptide (p). Two TCR-pMHC complexes are considered, which only differ by the mutation of a single amino acid on the peptide; one is a strong agonist that produces T cell activation in vivo, while the other is an antagonist. We investigate the interaction mechanism from a large number of unbinding trajectories by analyzing van der Waals and electrostatic interactions and by computing energy changes in proteins and solvent. In addition, dissociation potentials of mean force are calculated with the Jarzynski identity, using an averaging method developed for our steering scheme. We analyze the convergence of the Jarzynski exponential average, which is hampered by the large amount of dissipative work involved and the complexity of the system. The resulting dissociation free energies largely underestimate experimental values, but the simulations are able to clearly differentiate between wild-type and mutated TCR-pMHC and give insights into the dissociation mechanism.

Potential of mean force calculations of ligand binding to ion channels from Jarzynski's equality and umbrella sampling

Potential of mean force (PMF) calculations provide a reliable method for determination of the absolute binding free energies for protein-ligand systems. The common method used for this purpose -- umbrella sampling with weighted histogram analysis -- is computationally very laborious, which limits its applications. Recently, a much simpler alternative for PMF calculations has become available, namely, using Jarzynski's equality in steered molecular dynamics simulations. So far, there have been a few comparisons of the two methods and mostly in simple systems that do not reflect the complexities of protein-ligand systems. Here, we use both methods to calculate the PMF for ion permeation and ligand binding to ion channels. Comparison of results indicate that Jarzynski's method suffers from relaxation problems in complex systems and would require much longer simulation times to yield reliable PMFs for protein-ligand systems.

Comparison of free energy methods for molecular systems

We present a detailed comparison of computational efficiency and precision for several free energy difference (DeltaF) methods. The analysis includes both equilibrium and nonequilibrium approaches, and distinguishes between unidirectional and bidirectional methodologies. We are primarily interested in comparing two recently proposed approaches, adaptive integration, and single-ensemble path sampling to more established methodologies. As test cases, we study relative solvation free energies of large changes to the size or charge of a Lennard-Jones particle in explicit water. The results show that, for the systems used in this study, both adaptive integration and path sampling offer unique advantages over the more traditional approaches. Specifically, adaptive integration is found to provide very precise long-simulation DeltaF estimates as compared to other methods used in this report, while also offering rapid estimation of DeltaF. The results demonstrate that the adaptive integration approach is the best overall method for the systems studied here. The single-ensemble path sampling approach is found to be superior to ordinary Jarzynski averaging for the unidirectional, "fast-growth" nonequilibrium case. Closer examination of the path sampling approach on a two-dimensional system suggests it may be the overall method of choice when conformational sampling barriers are high. However, it appears that the free energy landscapes for the systems used in this study have rather modest configurational sampling barriers.

Gibbs free-energy estimates from direct path-sampling computations

We have implemented a path-sampling scheme enabling a direct estimation of Gibbs free energy. This scheme consists of a Monte Carlo sampling of constant-pressure Langevin paths, followed by an ensemble averaging carried out over the Markov chain of paths. In practice, we sample an umbrella path ensemble, which requires to rigorously define a statistical weight for the paths, equivalent of the Boltzmann weight. This statistical weight is a function of an effective work related to the path. The umbrella ensemble is chosen so that its work histogram overlaps with the histograms corresponding to the ensembles of forward and backward paths. We have finally investigated the relations between numerical efficiency and overlapping properties of the various work histograms. This analysis yields a built-in criterion for diagnosing the convergence during a single-run simulation.

Peptide conformational equilibria computed via a single-stage shifting protocol

We study the conformational equilibria of two peptides using a novel statistical mechanics approach designed for calculating free energy differences between highly dissimilar conformational states. Our results elucidate the contrasting roles of entropy in implicitly solvated leucine dipeptide and decaglycine. The method extends earlier work by Voter and overcomes the notorious "overlap" problem in free energy computations by constructing a mathematically equivalent calculation with high conformational similarity. The approach requires only equilibrium simulations of the two states of interest, without the need for sampling transition states. We discuss possible extensions and optimizations of the approach.

Phase-space overlap measures. I. Fail-safe bias detection in free energies calculated by molecular simulation

We consider ways to quantify the overlap of the parts of phase space important to two systems, labeled A and B. Of interest is how much of the A-important phase space lies in that important to B, and how much of B lies in A. Two measures are proposed. The first considers four total-energy distributions, formed from all combinations made by tabulating either the A-system or the B-system energy when sampling either the A or B system. Measures for A in B and B in A are given by two overlap integrals defined on pairs of these distributions. The second measure is based on information theory, and defines two relative entropies which are conveniently expressed in terms of the dissipated work for free-energy perturbation (FEP) calculations in the A-->B and B-->A directions, respectively. Phase-space overlap is an important consideration in the performance of free-energy calculations. To demonstrate this connection, we examine bias in FEP calculations applied to a system of independent particles in a harmonic potential. Systems are selected to represent a range of overlap situations, including extreme subset, subset, partial overlap, and nonoverlap. The magnitude and symmetry of the bias (A-->B vs B-->A) are shown to correlate well with the overlap, and consequently with the overlap measures. The relative entropies are used to scale the amount of sampling to obtain a universal bias curve. This result leads to develop a simple heuristic that can be applied to determine whether a work-based free-energy measurement is free of bias. The heuristic is based in part on the measured free energy, but we argue that it is fail-safe inasmuch as any bias in the measurement will not promote a false indication of accuracy.