Estimation of Free-Energy Differences from Computed Work Distributions: An Application of Jarzynski’s Equality

Amodiaquine Nonspecifically Binds Double Stranded and Three-Way Junction DNA Structures

The 4-aminoquinoline class of compounds includes the important antimalarial compounds amodiaquine and chloroquine. Despite their medicinal importance, the mode of action of these compounds is poorly understood. In a previous study we observed these compounds, as well as quinine and mefloquine, tightly bind the DNA cocaine-binding aptamer. Here, we further explore the range of nucleic acid structures bound by these compounds. To gauge a wide range of binding affinities, we used isothermal titration calorimetry to explore high affinity binding (nM to tens of μM) and NMR spectroscopy to assay weak binding biding in the hundreds of micromolar range. We find that amodiaquine tightly binds all double stranded DNA structures explored. Mefloquine binds double stranded DNA duplex molecules tightly and weakly associates with a three-way junction DNA construct. Quinine and chloroquine only weakly bind duplex DNA but do not tightly bind any of the DNA constructs explored. A simulation of the free energy of binding of these ligands to the Dickerson-Drew dodecamer resulted in an excellent agreement between the simulated and experimental free energy. These results provide new insight into the DNA binding of clinically important antimalarial compounds and may play a role in future development of new antimalarials.

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.

Implementation of Telescoping Boxes in Adaptive Steered Molecular Dynamics

Long-time dynamical processes, such as those involving protein unfolding and ligand interactions, can be accelerated and realized through steered molecular dynamics (SMD). The challenge has been the extraction of information from such simulations that generalize for complex nonequilibrium processes. The use of Jarzynski's equality opened the possibility of determining the free energy along the steered coordinate, but sampling over the nonequilibrium trajectories is slow to converge. Adaptive steered molecular dynamics (ASMD) and other related techniques have been introduced to overcome this challenge through the use of stages. Here, we take advantage of these stages to address the numerical cost that arises from the required use of very large solvent boxes. We introduce telescoping box schemes within adaptive steered molecular dynamics (ASMD) in which we adjust the solvent box between stages and thereby vary (and optimize) the required number of solvent molecules. We have benchmarked the method on a relatively long α-helical peptide, Ala₃₀, with respect to the potential of mean force and hydrogen bonds. We show that the use of telescoping boxes introduces little numerical error while significantly reducing the computational cost.

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.

Comparison of umbrella sampling and steered molecular dynamics methods for computing free energy profiles of aromatic substrates through phospholipid bilayers

Understanding the permeation of molecules through lipid membranes is fundamental for predicting the cellular uptake of solutes and drug delivery mechanisms. In molecular simulations, the usual approach is to compute the free energy (FE) profile of a molecule across a model lipid bilayer, which can then be used to estimate the permeability of the molecule. Umbrella Sampling (US), which involves carrying out a series of biased simulations along a defined reaction coordinate (usually the bilayer normal direction), is a popular method for the computation of such FE profiles. However, US can be challenging to implement because the results are dependent on the strength of the biasing potential and the spacing of windows along the reaction coordinate, which, in practice, are usually optimized by an inefficient trial and error approach. The Steered Molecular Dynamics implementation of the Jarzynski Equality (JE-SMD) has been identified as an alternative to equilibrium sampling methods for measuring the FE change across a reaction coordinate. In the JE-SMD approach, equilibrium FE values are evaluated from the average of rapid non-equilibrium trajectories, thus avoiding the practical issues that come with US. Here, we use three different corrections of the JE-SMD method to calculate the FE change for the translocation of two aromatic substrates, phenylalanine and toluene, across a lipid bilayer and compare the accuracy and computational efficiency of these approaches to the results obtained using US. We show evidence that when computing the free energy profile, the JE-SMD approach suffers from insufficient sampling convergence of the bilayer environment and is dependent on the characteristic of the aromatic substrate itself. We deduce that, despite its drawbacks, US remains the more viable approach of the two for computing the FE profile.

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.

Constrained Unfolding of a Helical Peptide: Implicit versus Explicit Solvents

Steered Molecular Dynamics (SMD) has been seen to provide the potential of mean force (PMF) along a peptide unfolding pathway effectively but at significant computational cost, particularly in all-atom solvents. Adaptive steered molecular dynamics (ASMD) has been seen to provide a significant computational advantage by limiting the spread of the trajectories in a staged approach. The contraction of the trajectories at the end of each stage can be performed by taking a structure whose nonequilibrium work is closest to the Jarzynski average (in naive ASMD) or by relaxing the trajectories under a no-work condition (in full-relaxation ASMD--namely, FR-ASMD). Both approaches have been used to determine the energetics and hydrogen-bonding structure along the pathway for unfolding of a benchmark peptide initially constrained as an α-helix in a water environment. The energetics are quite different to those in vacuum, but are found to be similar between implicit and explicit solvents. Surprisingly, the hydrogen-bonding pathways are also similar in the implicit and explicit solvents despite the fact that the solvent contact plays an important role in opening the helix.

High-precision work distributions for extreme nonequilibrium processes in large systems

The distributions of work for strongly nonequilibrium processes are studied using a very general form of a large-deviation approach, which allows one to study distributions down to extremely small probabilities of almost arbitrary quantities of interest for equilibrium, nonequilibrium stationary, and even nonstationary processes. The method is applied to quickly vary the external field in a wide range B = 3 ↔ 0 for a critical (T = 2.269) two-dimensional Ising system of size L × L = 128 × 128. To obtain free-energy differences from the work distributions, they must be studied in ranges where the probabilities are as small as 10^{-240}, which is not possible using direct simulation approaches. By comparison with the exact free energies, which are available for this model for the zero-field case, one sees that the present approach allows one to obtain the free energy with a very high relative precision of 10^{-4}. This works well also for a nonzero field, i.e., for a case where standard umbrella-sampling methods are not efficient to calculate free energies. Furthermore, for the present case it is verified that the resulting distributions of work for forward and backward processes fulfill Crooks theorem with high precision. Finally, the free energy for the Ising magnet as a function of the field strength is obtained.