Molecular Dynamics and Free Energy Calculations of Dicyclohexano-18-crown-6 Diastereoisomers with Sm2+, Eu2+, Dy2+, Yb2+, Cf2+, and Three Halide Salts in Tetrahydrofuran and Acetonitrile Using the AMOEBA Force Field

With the continual development of lanthanides (Ln) in current technological devices, an efficient separation process is needed that can recover greater amounts of these rare elements. Dicyclohexano-18-crown-6 (DCH18C6) is a crown ether that may be a promising candidate for Ln separation, but additional research is required. As such, molecular dynamics (MD) simulations have been performed on four divalent lanthanide halide salts (Sm2+, Eu2+, Dy2+, and Yb2+) and one divalent actinide halide salt (Cf2+) bound to three diastereoisomers of DCH18C6. Dy2+, Yb2+, Cf2+, DCH18C6, and tetrahydrofuran (THF) solvent were parameterized for the AMOEBA polarizable force field for the first time, whereas existing parameters for Sm2+ and Eu2+ were utilized from our previous efforts. A coordination number (CN) of six for Ln2+/An2+-O solvated in THF indicated that the cations interacted almost entirely with the oxygens of the polyether ring. A CN of one for Ln2+/An2+-N solvated in acetonitrile for systems containing iodide suggested that the N atom of acetonitrile was competitive with I- for cation interactions. Fluctuation between five and six CNs for Dy2+ and Yb2+ suggested that although the cations remained in the polyether ring, the size of the ring may not be an ideal fit as these cations possess comparatively smaller ionic radii. Gibbs binding free energies of Sm2+ in all DCH18C6 diastereoisomers solvated in THF were calculated. The binding free energy of the cis-syn-cis diastereoisomer was the most favorable, followed by cis-anti-cis, and then trans-anti-trans. Finally, two major types of conformation were observed for each diastereoisomer that were related to the electrostatic interactions and charge density of the cations.

Computational methods for calculation of protein-ligand binding affinities in structure-based drug design

Drug design is an expensive and time-consuming process. Any method that allows reducing the time the costs of the drug development project can have great practical value for the pharmaceutical industry. In structure-based drug design, affinity prediction methods are of great importance. The majority of methods used to predict binding free energy in protein-ligand complexes use molecular mechanics methods. However, many limitations of these methods in describing interactions exist. An attempt to go beyond these limits is the application of quantum-mechanical description for all or only part of the analyzed system. However, the extensive use of quantum mechanical (QM) approaches in drug discovery is still a demanding challenge. This chapter briefly reviews selected methods used to calculate protein-ligand binding affinity applied in virtual screening (VS), rescoring of docked poses, and lead optimization stage, including QM methods based on molecular simulations.

Accurate Binding Free Energy Method from End-State MD Simulations

Herein, we introduce a new strategy to estimate binding free energies using end-state molecular dynamics simulation trajectories. The method is adopted from linear interaction energy (LIE) and ANI-2x neural network potentials (machine learning) for the atomic simulation environment (ASE). It predicts the single-point interaction energies between ligand–protein and ligand–solvent pairs at the accuracy of the wb97x/6-31G* level for the conformational space that is sampled by molecular dynamics (MD) simulations. Our results on 54 protein–ligand complexes show that the method can be accurate and have a correlation of R = 0.87–0.88 to the experimental binding free energies, outperforming current end-state methods with reduced computational cost. The method also allows us to compare BFEs of ligands with different scaffolds. The code is available free of charge (documentation and test files) at https://github.com/otayfuroglu/deepQM.

Small-sample limit of the Bennett acceptance ratio method and the variationally derived intermediates

Free energy calculations based on atomistic Hamiltonians provide microscopic insight into the thermodynamic driving forces of biophysical or condensed matter systems. Many approaches use intermediate Hamiltonians interpolating between the two states for which the free energy difference is calculated. The Bennett acceptance ratio (BAR) and variationally derived intermediates (VI) methods are optimal estimator and intermediate states in that the mean-squared error of free energy calculations based on independent sampling is minimized. However, BAR and VI have been derived based on several approximations that do not hold for very few sample points. Analyzing one-dimensional test systems, we show that in such cases BAR and VI are suboptimal and that established uncertainty estimates are inaccurate. Whereas for VI to become optimal, less than seven samples per state suffice in all cases; for BAR the required number increases unboundedly with decreasing configuration space densities overlap of the end states. We show that for BAR, the required number of samples is related to the overlap through an inverse power law. Because this relation seems to hold universally and almost independent of other system properties, these findings can guide the proper choice of estimators for free energy calculations.

Affordable Ab Initio Path Integral for Thermodynamic Properties via Molecular Dynamics Simulations Using Semiempirical Reference Potential

Path integral molecular dynamics (PIMD) is becoming a routinely applied method for incorporating the nuclear quantum effect in computer simulations. However, direct PIMD simulations at an ab initio level of theory are formidably expensive. Using the protonated 1,8-bis(dimethylamino)naphthalene molecule as an example, we show in this work that the computational expense for the intramolecular proton transfer between the two nitrogen atoms can be remarkably reduced by implementing the idea of reference-potential methods. The simulation time can be easily extended to a scale of nanoseconds while maintaining the accuracy on an ab initio level of theory for thermodynamic properties. In addition, postprocessing can be carried out in parallel on massive computer nodes. A 545-fold reduction in the total CPU time can be achieved in this way as compared to a direct PIMD simulation at the same ab initio level of theory.

Efficient Alchemical Intermediate States in Free Energy Calculations Using λ-Enveloping Distribution Sampling

Alchemical free energy calculations generally require intermediate states along a coupling parameter λ to establish sufficient phase space overlap for obtaining converged results. Such intermediate states can also be engineered to lower the energy barriers and, consequently, reduce the required sampling time. The recently introduced λ-enveloping distribution sampling (λ-EDS) scheme combines the properties of the minimum variance pathway and the EDS methods to improve sampling and allow for larger steps along the alchemical pathway compared to conventional approaches. This scheme also eliminates the need for soft-core potentials and retains the behavior of conventional λ-intermediate states as a limiting case. In this study, an automated procedure is developed to select the parameters of λ-EDS for optimal performance. The underlying theory is illustrated based on simulations of simple test systems (bond length changes in harmonic oscillators, mutations of dihedral angles, and charge creation in water), as well as on the calculation of the absolute hydration free energies of 12 small organic molecules.

Accelerated Computation of Free Energy Profile at Ab Initio Quantum Mechanical/Molecular Mechanics Accuracy via a Semiempirical Reference Potential. 4. Adaptive QM/MM

Although quantum mechanical/molecular mechanics (QM/MM) methods are now routinely applied to the studies of chemical reactions in condensed phases and enzymatic reactions, they may experience technical difficulties when the reactive region is varying over time. For instance, when the solvent molecules are directly participating in the reaction, the exchange of water molecules between the QM and MM regions may occur on a time scale comparable to the reaction time. To cope with this situation, several adaptive QM/MM schemes have been proposed. However, these methods either add significantly to the computational cost or introduce artificial restraints to the system. In this work, we developed a novel adaptive QM/MM scheme and applied it to the study of a nucleophilic addition reaction. In this scheme, the configuration sampling was performed with a small QM region (without solvent molecules), and the thermodynamic properties under another potential energy function with a larger QM region (with a certain number of solvent molecules and/or different levels of QM theory) are computed via extrapolation using the reference-potential method. Our simulation results show that this adaptive QM/MM scheme is numerically stable, at least for the case studied in this work. Furthermore, this method also offers an inexpensive way to examine the convergence of the QM/MM calculation with respect to the size of the QM region.

Connecting dynamic reweighting Algorithms: Derivation of the dynamic reweighting family tree

Thermally driven processes of molecular systems include transitions of energy barriers on the microsecond timescales and higher. Sufficient sampling of such processes with molecular dynamics simulations is challenging and often requires accelerating slow transitions using external biasing potentials. Different dynamic reweighting algorithms have been proposed in the past few years to recover the unbiased kinetics from biased systems. However, it remains an open question if and how these dynamic reweighting approaches are connected. In this work, we establish the link between the two main reweighting types, i.e., path-based and energy-based reweighting. We derive a path-based correction factor for the energy-based dynamic histogram analysis method, thus connecting the previously separate reweighting types. We show that the correction factor can be used to combine the advantages of path-based and energy-based reweighting algorithms: it is integrator independent, more robust, and at the same time able to reweight time-dependent biases. We can furthermore demonstrate the relationship between two independently derived path-based reweighting algorithms. Our theoretical findings are verified on a one-dimensional four-well system. By connecting different dynamic reweighting algorithms, this work helps to clarify the strengths and limitations of the different methods and enables a more robust usage of the combined types.

Alchemical Binding Free Energy Calculations in AMBER20: Advances and Best Practices for Drug Discovery

Predicting protein-ligand binding affinities and the associated thermodynamics of biomolecular recognition is a primary objective of structure-based drug design. Alchemical free energy simulations offer a highly accurate and computationally efficient route to achieving this goal. While the AMBER molecular dynamics package has successfully been used for alchemical free energy simulations in academic research groups for decades, widespread impact in industrial drug discovery settings has been minimal because of the previous limitations within the AMBER alchemical code, coupled with challenges in system setup and postprocessing workflows. Through a close academia-industry collaboration we have addressed many of the previous limitations with an aim to improve accuracy, efficiency, and robustness of alchemical binding free energy simulations in industrial drug discovery applications. Here, we highlight some of the recent advances in AMBER20 with a focus on alchemical binding free energy (BFE) calculations, which are less computationally intensive than alternative binding free energy methods where full binding/unbinding paths are explored. In addition to scientific and technical advances in AMBER20, we also describe the essential practical aspects associated with running relative alchemical BFE calculations, along with recommendations for best practices, highlighting the importance not only of the alchemical simulation code but also the auxiliary functionalities and expertise required to obtain accurate and reliable results. This work is intended to provide a contemporary overview of the scientific, technical, and practical issues associated with running relative BFE simulations in AMBER20, with a focus on real-world drug discovery applications.

Accelerated Computation of Free Energy Profile at Ab Initio Quantum Mechanical/Molecular Mechanics Accuracy via a Semiempirical Reference Potential. 3. Gaussian Smoothing on Density-of-States

Calculations of the free energy profile, also known as potential of mean force (PMF), along a chosen collective variable (CV) are now routinely applied in the studies of chemical processes, such as enzymatic reactions and chemical reactions in condensed phases. However, if the ab initio quantum mechanical/molecular mechanics (QM/MM) level of accuracy is required for the PMF, it can be formidably demanding even with the most advanced enhanced sampling methods, such as umbrella sampling. To ameliorate this difficulty, we developed a novel method for the computation of the free energy profile based on the reference-potential method recently, in which a low-level reference Hamiltonian is employed for phase space sampling and the free energy profile can be corrected to the level of interest (the target Hamiltonian) by energy reweighting in a nonparametric way. However, when the reference Hamiltonian is very different from the target Hamiltonian, the calculated ensemble averages, including the PMF, often suffer from numerical instability, which mainly comes from the overestimation of the density-of-states (DoS) in the low-energy region. Stochastic samplings of these low-energy configurations are rare events, and some low-energy conformations may get oversampled in simulations of a finite length. In this work, an assumption of Gaussian distribution is applied to the DoS in each CV bin, and the weight of each configuration is rescaled according to the accumulated DoS. The results show that this smoothing process can remarkably reduce the ruggedness of the PMF and increase the reliability of the reference-potential method.

An Alternative to Conventional λ-Intermediate States in Alchemical Free Energy Calculations: λ-Enveloping Distribution Sampling

Alchemical free energy calculations typically rely on intermediate states to bridge between the relevant phase spaces of the two end states. These intermediate states are usually created by mixing the energies or parameters of the end states according to a coupling parameter λ. The choice of the procedure has a strong impact on the efficiency of the calculation, as it affects both the encountered energy barriers and the phase space overlap between the states. The present work builds on the connection between the minimum variance pathway (MVP) and enveloping distribution sampling (EDS). It is shown that both methods can be regarded as special cases of a common scheme referred to as λ-EDS, which can also reproduce the behavior of conventional λ-intermediate states. A particularly attractive feature of λ-EDS is its ability to emulate the use of soft core potentials (SCP) while avoiding the associated computational overhead when applying efficient free energy estimators such as the multistate Bennett's acceptance ratio (MBAR). The method is illustrated for both relative and absolute free energy calculations considering five benchmark systems. The first two systems (charge inversion and cavity creation in a dipolar solvent) demonstrate the use of λ-EDS as an alternative coupling scheme in the context of thermodynamic integration (TI). The three other systems (change of bond length, change of dihedral angles, and cavity creation in water) investigate the efficiency and optimal choice of parameters in the context of free energy perturbation (FEP) and Bennett's acceptance ratio (BAR). It is shown that λ-EDS allows larger steps along the alchemical pathway than conventional intermediate states.

Prediction of aqueous free energies of solvation using coupled QM and MM explicit solvent simulations

A method based on molecular dynamics simulations which employ two distinct levels of theory is proposed and tested for the prediction of Gibbs free energies of solvation for non-ionic solutes in water. The method consists of two additive contributions: (i) an evaluation of the free energy of solvation predicted by a computationally efficient molecular mechanics (MM) method; and (ii) an evaluation of the free energy difference between the potential energy surface of the MM method and that of a more computationally intensive first-principles quantum-mechanical (QM) method. The latter is computed by a thermodynamic integration method based on a series of shorter molecular dynamics simulations that employ weighted averages of the QM and MM force evaluations. The combined computational approach is tested against the experimental free energies of aqueous solvation for four solutes. For solute-solvent interactions that are found to be described qualitatively well by the MM method, the QM correction makes a modest improvement in the predicted free energy of aqueous solvation. However, for solutes that are found to not be adequately described by the MM method, the QM correction does not improve agreement with experiment. These preliminary results provide valuable insights into the novel concept of implementing thermodynamic integration between two model chemistries, suggesting that it is possible to use QM methods to improve upon the MM predictions of free energies of aqueous solvation.

BAR-based multi-dimensional nonequilibrium pulling for indirect construction of QM/MM free energy landscapes: from semi-empirical to ab initio

The indirect method for the construction of quantum mechanics (QM)/molecular mechanics (MM) free energy landscapes provides a cheaper alternative for free energy simulations at the QM level. The indirect method features a direct calculation of the free energy profile with a computationally efficient but less accurate Hamiltonian (i.e. low-level Hamiltonian) and a low-level-to-high-level correction. In the thermodynamic cycle, the direct low-level calculation along the physically meaningful reaction coordinate is corrected via the alchemical method, which is often achieved with perturbation-based techniques. In our previous work, a multi-dimensional nonequilibrium pulling framework is proposed for the indirect construction of QM/MM free energy landscapes. Previously, we focused on obtaining semi-empirical QM (SQM) results indirectly from direct MM simulations and MM to SQM corrections. In this work, we apply this method to obtain results under ab initio QM Hamiltonians by combining direct SQM results and SQM to QM corrections. A series of SQM and QM Hamiltonians are benchmarked. It is observed that PM6 achieves the best performance among the low-level Hamiltonians. Therefore, we recommend using PM6 as the low-level theory in the indirect free energy simulation. Considering its higher similarity to the high-level Hamiltonians, PM6 corrected with the bond charge correction could be more accurate than the existing AM1-BCC model. Another central result in the current work is a basic protocol of choosing the strength of restraints and an appropriate time step in nonequilibrium free energy simulation at the stiff spring limit. We provide theoretical derivations to emphasize the importance of using a sufficiently large force constant and choosing an appropriate time step. It is worth noting that a general rule of thumb for choosing the time step, according to our derivation, is that a time step of 1 fs or smaller should be used, as long as the stiff spring approximation is employed, even in simulations with constraints on bonds involving hydrogen atoms.

An Albumin Sandwich Enhances in Vivo Circulation and Stability of Metabolically Labile Peptides

The effectiveness of numerous molecular drugs is hampered by their poor pharmacokinetics. Different from previous approaches with limited effectiveness, most recently, emerging high-affinity albumin binding moieties (ABMs) for in vivo hitchhiking of endogenous albumin opens up an avenue to chaperone small molecules for long-acting therapeutics. Although several FDA-approved fatty acids have shown prolonged residence and therapeutic effect, an easily synthesized, water-soluble, and high-efficiency ABM with versatile drug loading ability is urgently needed to improve the therapeutic efficacy of short-lived constructs. We herein identified an ideal bivalent Evans blue derivative, denoted as N(tEB)₂, as a smart ABM-delivery platform to chaperone short-lived molecules, through both computational modeling screening and efficient synthetic schemes. The optimal N(tEB)₂ could reversibly link two molecules of albumin through its two binding heads with a preferable spacer, resulting in significantly extended circulation half-life of a preloaded cargo and water-soluble. Notably, this in situ dimerization of albumin was able to sandwich peptide therapeutics to protect them from proteolysis. As an application, we conjugated N(tEB)₂ with exendin-4 for long-acting glucose control in a diabetic mouse model, and it was superior to both previously tested NtEB-exendin-4 (Abextide) and the newly FDA-approved semaglutide, which has been arguably the best commercial weekly formula so far. Hence, this novel albumin binder has excellent clinical potential for next-generation biomimetic drug delivery systems.

Exploring Gatekeeper Mutations in EGFR through Computer Simulations

The emergence of resistance against drugs that inhibit a particular protein is a major problem in targeted therapy. There is a clear need for rigorous methods to predict the likelihood of specific drug-resistance mutations arising in response to the binding of a drug. In this work we attempt to develop a robust computational protocol for predicting drug resistant mutations at the gatekeeper position (T790) in EGFR. We explore how mutations at this site affects interactions with ATP and three drugs that are currently used in clinics. We found, as expected, that certain mutations are not tolerated structurally, while some other mutations interfere with the natural substrate and hence are unlikely to be selected for. However, we found five possible mutations that are well tolerated structurally and energetically. Two of these mutations were predicted to have increased affinity for the drugs over ATP, as has been reported earlier. By reproducing the trends in the experimental binding affinities of the data, the methods chosen here are able to correctly predict the effects of these mutations on the binding affinities of the drugs. However, the increased affinity does not always translate into increased efficacy, because the efficacy is affected by several other factors such as binding kinetics, competition with ATP, and residence times. The computational methods used in the current study are able to reproduce or predict the effects of mutations on the binding affinities. However, a different set of methods is required to predict the kinetics of drug binding.

BAR-based multi-dimensional nonequilibrium pulling for indirect construction of a QM/MM free energy landscape

Construction of free energy landscapes at the quantum mechanics (QM) level is computationally demanding. As shown in previous studies, by employing an indirect scheme (i.e. constructing a thermodynamic cycle connecting QM states via an alchemical pathway), simulations are converged with much less computational burden. The indirect scheme makes QM/molecular mechanics (MM) free energy simulation orders of magnitude faster than the direct QM/MM schemes. However, the indirect QM/MM simulations were mostly equilibrium sampling based and the nonequilibrium methods were merely exploited in one-dimensional alchemical QM/MM end-state correction at two end states. In this work, we represent a multi-dimensional nonequilibrium pulling scheme for indirect QM/MM free energy simulations, where the whole free energy simulation is performed only with nonequilibrium methods. The collective variable (CV) space we explore is a combination of one alchemical CV and one physically meaningful CV. The current nonequilibrium indirect QM/MM simulation method can be seen as the generalization of equilibrium perturbation based indirect QM/MM methods. The test systems include one backbone dihedral case and one distance case. The two cases are significantly different in size, enabling us to investigate the dependence of the speedup of the indirect scheme on the size of the system. It is shown that the speedup becomes larger when the size of the system becomes larger, which is consistent with the scaling behavior of QM Hamiltonians.

Solvation Free Energy Calculation Using a Fixed-Charge Model: Implicit and Explicit Treatments of the Polarization Effect

In this work, IPolQ-Mod charges and the reference potential scheme are used to calculate the solvation free energies of a set of organic molecules. Both methods could capture the phase transfer of a solute with accompanying polarization cost utilizing a fixed-charge model. The IPolQ-Mod charges, which are the average of two charge sets fitted in a vacuum state and a condensed phase, take account of the polarization effect implicitly. For the reference potential method, the quantum mechanics polarization corrections are calculated explicitly by thermodynamic perturbation. The polarization effect captured by the IPolQ-Mod charges is an approximation to that of the reference potential method theoretically. In the present study, the reference potential method shows a slight improvement over the classical restrained electrostatic potential (RESP) charges, which perform pretty well in predicting the solvation free energy. However, IPolQ-Mod(MP2) shows a poor agreement with the experimental data. Compared with IPolQ-Mod(MP2), IPolQ-Mod(M06-2X) or IPolQ-Mod(ωB97X) is found to give more appropriate prediction of the molecule's dipole and the solvation free energies calculated by IPolQ-Mod(M06-2X) or IPolQ-Mod(ωB97X) are more compatible with those of the RESP charges. If the other force field parameters remain unchanged, M06-2X or ωB97X is recommended to derive the IPolQ-Mod charges.

A Comparison of QM/MM Simulations with and without the Drude Oscillator Model Based on Hydration Free Energies of Simple Solutes

Maintaining a proper balance between specific intermolecular interactions and non-specific solvent interactions is of critical importance in molecular simulations, especially when predicting binding affinities or reaction rates in the condensed phase. The most rigorous metric for characterizing solvent affinity are solvation free energies, which correspond to a transfer from the gas phase into solution. Due to the drastic change of the electrostatic environment during this process, it is also a stringent test of polarization response in the model. Here, we employ both the CHARMM fixed charge and polarizable force fields to predict hydration free energies of twelve simple solutes. The resulting classical ensembles are then reweighted to obtain QM/MM hydration free energies using a variety of QM methods, including MP2, Hartree⁻Fock, density functional methods (BLYP, B3LYP, M06-2X) and semi-empirical methods (OM2 and AM1 ). Our simulations test the compatibility of quantum-mechanical methods with molecular-mechanical water models and solute Lennard⁻Jones parameters. In all cases, the resulting QM/MM hydration free energies were inferior to purely classical results, with the QM/MM Drude force field predictions being only marginally better than the QM/MM fixed charge results. In addition, the QM/MM results for different quantum methods are highly divergent, with almost inverted trends for polarizable and fixed charge water models. While this does not necessarily imply deficiencies in the QM models themselves, it underscores the need to develop consistent and balanced QM/MM interactions. Both the QM and the MM component of a QM/MM simulation have to match, in order to avoid artifacts due to biased solute⁻solvent interactions. Finally, we discuss strategies to improve the convergence and efficiency of multi-scale free energy simulations by automatically adapting the molecular-mechanics force field to the target quantum method.

GPU-Accelerated Molecular Dynamics and Free Energy Methods in Amber18: Performance Enhancements and New Features

We report progress in graphics processing unit (GPU)-accelerated molecular dynamics and free energy methods in Amber18. Of particular interest is the development of alchemical free energy algorithms, including free energy perturbation and thermodynamic integration methods with support for nonlinear soft-core potential and parameter interpolation transformation pathways. These methods can be used in conjunction with enhanced sampling techniques such as replica exchange, constant-pH molecular dynamics, and new 12-6-4 potentials for metal ions. Additional performance enhancements have been made that enable appreciable speed-up on GPUs relative to the previous software release.

On the convergence of multi-scale free energy simulations

In this work we employ simple model systems to evaluate the relative performance of two of the most important free energy methods: The Zwanzig equation (also known as "Free energy perturbation") and Bennett's acceptance ratio method (BAR). Although our examples should be transferable to other kinds of free energy simulations, we focus on applications of multi-scale free energy simulations. Such calculations are especially complex, since they connect two different levels of theory with very different requirements in terms of speed, accuracy, sampling and parallelizability. We try to reconcile all those different factors by developing some simple criteria to guide the early stages of the development of a free energy protocol. This is accomplished by quantifying how many λ intermediate steps and how many potential energy evaluations are necessary in order to reach a certain level of convergence.

Relative Ligand-Binding Free Energies Calculated from Multiple Short QM/MM MD Simulations

We have devised a new efficient approach to compute combined quantum mechanical (QM) and molecular mechanical (MM, i.e. QM/MM) ligand-binding relative free energies. Our method employs the reference-potential approach with free-energy perturbation both at the MM level (between the two ligands) and from MM to QM/MM (for each ligand). To ensure that converged results are obtained for the MM → QM/MM perturbations, explicit QM/MM molecular dynamics (MD) simulations are performed with two intermediate mixed states. To speed up the calculations, we utilize the fact that the phase space can be extensively sampled at the MM level. Therefore, we run many short QM/MM MD simulations started from snapshots of the MM simulations, instead of a single long simulation. As a test case, we study the binding of nine cyclic carboxylate ligands to the octa-acid deep cavitand. Only the ligand is in the QM system, treated with the semiempirical PM6-DH+ method. We show that for eight of the ligands, we obtain well converged results with short MD simulations (1-15 ps). However, in one case, the convergence is slower (∼50 ps) owing to a mismatch between the conformational preferences of the MM and QM/MM potentials. We test the effect of initial minimization, the need of equilibration, and how many independent simulations are needed to reach a certain precision. The results show that the present approach is about four times faster than using standard MM → QM/MM free-energy perturbations with the same accuracy and precision.

Toward Fast and Accurate Binding Affinity Prediction with pmemdGTI: An Efficient Implementation of GPU-Accelerated Thermodynamic Integration

We report the implementation of the thermodynamic integration method on the pmemd module of the AMBER 16 package on GPUs (pmemdGTI). The pmemdGTI code typically delivers over 2 orders of magnitude of speed-up relative to a single CPU core for the calculation of ligand-protein binding affinities with no statistically significant numerical differences and thus provides a powerful new tool for drug discovery applications.

An Estimation of Hybrid Quantum Mechanical Molecular Mechanical Polarization Energies for Small Molecules Using Polarizable Force-Field Approaches

In this work, we report two polarizable molecular mechanics (polMM) force field models for estimating the polarization energy in hybrid quantum mechanical molecular mechanical (QM/MM) calculations. These two models, named the potential of atomic charges (PAC) and potential of atomic dipoles (PAD), are formulated from the ab initio quantum mechanical (QM) response kernels for the prediction of the QM density response to an external molecular mechanical (MM) environment (as described by external point charges). The PAC model is similar to fluctuating charge (FQ) models because the energy depends on external electrostatic potential values at QM atomic sites; the PAD energy depends on external electrostatic field values at QM atomic sites, resembling induced dipole (ID) models. To demonstrate their uses, we apply the PAC and PAD models to 12 small molecules, which are solvated by TIP3P water. The PAC model reproduces the QM/MM polarization energy with a R² value of 0.71 for aniline (in 10,000 TIP3P water configurations) and 0.87 or higher for other 11 solute molecules, while the PAD model has a much better performance with R² values of 0.98 or higher. The PAC model reproduces reference QM/MM hydration free energies for 12 solute molecules with a RMSD of 0.59 kcal/mol. The PAD model is even more accurate, with a much smaller RMSD of 0.12 kcal/mol, with respect to the reference. This suggests that polarization effects, including both local charge distortion and intramolecular charge transfer, can be well captured by induced dipole type models with proper parametrization.

Binding free energies in the SAMPL5 octa-acid host-guest challenge calculated with DFT-D3 and CCSD(T)

We have tried to calculate the free energy for the binding of six small ligands to two variants of the octa-acid deep cavitand host in the SAMPL5 blind challenge. We employed structures minimised with dispersion-corrected density-functional theory with small basis sets and energies were calculated using large basis sets. Solvation energies were calculated with continuum methods and thermostatistical corrections were obtained from frequencies calculated at the HF-3c level. Care was taken to minimise the effects of the flexibility of the host by keeping the complexes as symmetric and similar as possible. In some calculations, the large net charge of the host was reduced by removing the propionate and benzoate groups. In addition, the effect of a restricted molecular dynamics sampling of structures was tested. Finally, we tried to improve the energies by using the DLPNO-CCSD(T) approach. Unfortunately, results of quite poor quality were obtained, with no correlation to the experimental data, systematically too positive affinities (by ~50 kJ/mol) and a mean absolute error (after removal of the systematic error) of 11-16 kJ/mol. DLPNO-CCSD(T) did not improve the results, so the accuracy is not limited by the energy function. Instead, four likely sources of errors were identified: first, the minimised structures were often incorrect, owing to the omission of explicit solvent. They could be partly improved by performing the minimisations in a continuum solvent with four water molecules around the charged groups of the ligands. Second, some ligands could bind in several different conformations, requiring sampling of reasonable structures. Third, there is an indication the continuum-solvation model has problems to accurately describe the binding of both the negatively and positively charged guest molecules. Fourth, different methods to calculate the thermostatistical corrections gave results that differed by up to 30 kJ/mol and there is an indication that HF-3c overestimates the entropy term. In conclusion, it is a challenge to calculate binding affinities for this octa-acid system with quantum-mechanical methods.

Absolute binding free energy calculations of CBClip host-guest systems in the SAMPL5 blind challenge

Herein, we report the absolute binding free energy calculations of CBClip complexes in the SAMPL5 blind challenge. Initial conformations of CBClip complexes were obtained using docking and molecular dynamics simulations. Free energy calculations were performed using thermodynamic integration (TI) with soft-core potentials and Bennett's acceptance ratio (BAR) method based on a serial insertion scheme. We compared the results obtained with TI simulations with soft-core potentials and Hamiltonian replica exchange simulations with the serial insertion method combined with the BAR method. The results show that the difference between the two methods can be mainly attributed to the van der Waals free energies, suggesting that either the simulations used for TI or the simulations used for BAR, or both are not fully converged and the two sets of simulations may have sampled difference phase space regions. The penalty scores of force field parameters of the 10 guest molecules provided by CHARMM Generalized Force Field can be an indicator of the accuracy of binding free energy calculations. Among our submissions, the combination of docking and TI performed best, which yielded the root mean square deviation of 2.94 kcal/mol and an average unsigned error of 3.41 kcal/mol for the ten guest molecules. These values were best overall among all participants. However, our submissions had little correlation with experiments.

Blind prediction of distribution in the SAMPL5 challenge with QM based protomer and pK a corrections

The computation of distribution coefficients between polar and apolar phases requires both an accurate characterization of transfer free energies between phases and proper accounting of ionization and protomerization. We present a protocol for accurately predicting partition coefficients between two immiscible phases, and then apply it to 53 drug-like molecules in the SAMPL5 blind prediction challenge. Our results combine implicit solvent QM calculations with classical MD simulations using the non-Boltzmann Bennett free energy estimator. The OLYP/DZP/SMD method yields predictions that have a small deviation from experiment (RMSD = 2.3 [Formula: see text] D units), relative to other participants in the challenge. Our free energy corrections based on QM protomer and [Formula: see text] calculations increase the correlation between predicted and experimental distribution coefficients, for all methods used. Unfortunately, these corrections are overly hydrophilic, and fail to account for additional effects such as aggregation, water dragging and the presence of polar impurities in the apolar phase. We show that, although expensive, QM-NBB free energy calculations offer an accurate and robust method that is superior to standard MM and QM techniques alone.

Calculating distribution coefficients based on multi-scale free energy simulations: an evaluation of MM and QM/MM explicit solvent simulations of water-cyclohexane transfer in the SAMPL5 challenge

One of the central aspects of biomolecular recognition is the hydrophobic effect, which is experimentally evaluated by measuring the distribution coefficients of compounds between polar and apolar phases. We use our predictions of the distribution coefficients between water and cyclohexane from the SAMPL5 challenge to estimate the hydrophobicity of different explicit solvent simulation techniques. Based on molecular dynamics trajectories with the CHARMM General Force Field, we compare pure molecular mechanics (MM) with quantum-mechanical (QM) calculations based on QM/MM schemes that treat the solvent at the MM level. We perform QM/MM with both density functional theory (BLYP) and semi-empirical methods (OM1, OM2, OM3, PM3). The calculations also serve to test the sensitivity of partition coefficients to solute polarizability as well as the interplay of the quantum-mechanical region with the fixed-charge molecular mechanics environment. Our results indicate that QM/MM with both BLYP and OM2 outperforms pure MM. However, this observation is limited to a subset of cases where convergence of the free energy can be achieved.

Extended Thermodynamic Integration: Efficient Prediction of Lambda Derivatives at Nonsimulated Points

Thermodynamic integration (TI) is one of the most commonly used free-energy calculation methods. The derivative of the Hamiltonian with respect to lambda, ⟨∂H/∂λ⟩, is determined at multiple λ-points. Because a numerical integration step is necessary, high curvature regions require simulations at densely spaced λ-points. Here, the principle of extended TI is introduced, where ⟨∂H/∂λ⟩ values are predicted at nonsimulated λ-points. On the basis of three model systems, it is shown that extended TI requires significantly fewer λ-points than regular TI to obtain similar accuracy.

An efficient protocol for obtaining accurate hydration free energies using quantum chemistry and reweighting from molecular dynamics simulations

The non-Boltzmann Bennett (NBB) free energy estimator method is applied to 21 molecules from the blind subset of the SAMPL4 challenge. When NBB is applied with the SMD implicit solvent model, and the OLYP/DZP level of quantum chemistry, highly accurate hydration free energy calculations are obtained with respect to experiment (RMSD=0.89kcal·mol^-1). Other quantum chemical methods are also tested, and the effects of solvent model, density functional, basis set are explored in this benchmarking study, providing a framework for improvements in calculating hydration free energies. We provide a practical guide for using the best QM-NBB protocols that are consistently more accurate than either pure QM or pure MM alone. In situations where high accuracy hydration free energy predictions are needed, the QM-NBB method with SMD implicit solvent should be the first choice of quantum chemists.

Methods for Efficiently and Accurately Computing Quantum Mechanical Free Energies for Enzyme Catalysis

Enzyme activity is inherently linked to free energies of transition states, ligand binding, protonation/deprotonation, etc.; these free energies, and thus enzyme function, can be affected by residue mutations, allosterically induced conformational changes, and much more. Therefore, being able to predict free energies associated with enzymatic processes is critical to understanding and predicting their function. Free energy simulation (FES) has historically been a computational challenge as it requires both the accurate description of inter- and intramolecular interactions and adequate sampling of all relevant conformational degrees of freedom. The hybrid quantum mechanical molecular mechanical (QM/MM) framework is the current tool of choice when accurate computations of macromolecular systems are essential. Unfortunately, robust and efficient approaches that employ the high levels of computational theory needed to accurately describe many reactive processes (ie, ab initio, DFT), while also including explicit solvation effects and accounting for extensive conformational sampling are essentially nonexistent. In this chapter, we will give a brief overview of two recently developed methods that mitigate several major challenges associated with QM/MM FES: the QM non-Boltzmann Bennett's acceptance ratio method and the QM nonequilibrium work method. We will also describe usage of these methods to calculate free energies associated with (1) relative properties and (2) along reaction paths, using simple test cases with relevance to enzymes examples.

Converging ligand-binding free energies obtained with free-energy perturbations at the quantum mechanical level

In this article, the convergence of quantum mechanical (QM) free-energy simulations based on molecular dynamics simulations at the molecular mechanics (MM) level has been investigated. We have estimated relative free energies for the binding of nine cyclic carboxylate ligands to the octa-acid deep-cavity host, including the host, the ligand, and all water molecules within 4.5 Å of the ligand in the QM calculations (158-224 atoms). We use single-step exponential averaging (ssEA) and the non-Boltzmann Bennett acceptance ratio (NBB) methods to estimate QM/MM free energy with the semi-empirical PM6-DH2X method, both based on interaction energies. We show that ssEA with cumulant expansion gives a better convergence and uses half as many QM calculations as NBB, although the two methods give consistent results. With 720,000 QM calculations per transformation, QM/MM free-energy estimates with a precision of 1 kJ/mol can be obtained for all eight relative energies with ssEA, showing that this approach can be used to calculate converged QM/MM binding free energies for realistic systems and large QM partitions. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.

Comparison of Methods To Reweight from Classical Molecular Simulations to QM/MM Potentials

We examine methods to reweight classical molecular mechanics solvation calculations to more expensive QM/MM energy functions. We first consider the solvation free energy difference between ethane and methanol in a QM/MM Hamiltonian from configurations generated in a cheaper MM potential. The solute molecules in the QM/MM Hamiltonian are treated with B3LYP/6-31G*, and the solvent water molecules are treated classically. The free energy difference in the QM/MM Hamiltonian is estimated using Boltzmann reweighting with both the non-Boltzmann Bennett method (NBB) and the multistate Bennett acceptance ratio (MBAR), and the variance of each method is directly compared for an identical data set. For this system, MBAR-derived methods are found to produce smaller overall uncertainties than NBB-based methods. Additionally, we show that to reduce the variance in the overall free energy difference estimate in this system for a fixed amount of QM/MM calculations, the energy re-evaluations in the Boltzmann reweighting step should be concentrated on the physical MM states with the highest overlap to the QM/MM states, rather than allocated equally over all sampled MM states. We also show that reallocating the QM/MM re-evaluations can be used to diagnose poor overlap between the sampled and target state. The solvation free energies for molecules in the SAMPL4 solvation data set are also calculated in the QM/MM Hamiltonian with NBB and MBAR, and the variances are marginally smaller for MBAR. Overall, NBB and MBAR produce similar variances for systems with poor sampling efficiency, and MBAR provides smaller variances than NBB in systems with high sampling efficiency. Both NBB and MBAR converge to identical solvation free energy estimates in the QM/MM Hamiltonian, and the RMSD to experimental values for molecules in the SAMPL4 solvation data set decreases by approximately 28% when switching from the MM Hamiltonian to the QM/MM Hamiltonian.

Computation of Hydration Free Energies Using the Multiple Environment Single System Quantum Mechanical/Molecular Mechanical Method

A recently developed MESS-E-QM/MM method (multiple-environment single-system quantum mechanical molecular/mechanical calculations with a Roothaan-step extrapolation) is applied to the computation of hydration free energies for the blind SAMPL4 test set and for 12 small molecules. First, free energy simulations are performed with a classical molecular mechanics force field using fixed-geometry solute molecules and explicit TIP3P solvent, and then the non-Boltzmann-Bennett method is employed to compute the QM/MM correction (QM/MM-NBB) to the molecular mechanical hydration free energies. For the SAMPL4 set, MESS-E-QM/MM-NBB corrections to the hydration free energy can be obtained 2 or 3 orders of magnitude faster than fully converged QM/MM-NBB corrections, and, on average, the hydration free energies predicted with MESS-E-QM/MM-NBB fall within 0.10-0.20 kcal/mol of full-converged QM/MM-NBB results. Out of five density functionals (BLYP, B3LYP, PBE0, M06-2X, and ωB97X-D), the BLYP functional is found to be most compatible with the TIP3P solvent model and yields the most accurate hydration free energies against experimental values for solute molecules included in this study.

Computational scheme for pH-dependent binding free energy calculation with explicit solvent

We present a computational scheme to compute the pH‐dependence of binding free energy with explicit solvent. Despite the importance of pH, the effect of pH has been generally neglected in binding free energy calculations because of a lack of accurate methods to model it. To address this limitation, we use a constant‐pH methodology to obtain a true ensemble of multiple protonation states of a titratable system at a given pH and analyze the ensemble using the Bennett acceptance ratio (BAR) method. The constant pH method is based on the combination of enveloping distribution sampling (EDS) with the Hamiltonian replica exchange method (HREM), which yields an accurate semi‐grand canonical ensemble of a titratable system. By considering the free energy change of constraining multiple protonation states to a single state or releasing a single protonation state to multiple states, the pH dependent binding free energy profile can be obtained. We perform benchmark simulations of a host‐guest system: cucurbit[7]uril (CB[7]) and benzimidazole (BZ). BZ experiences a large pK_a shift upon complex formation. The pH‐dependent binding free energy profiles of the benchmark system are obtained with three different long‐range interaction calculation schemes: a cutoff, the particle mesh Ewald (PME), and the isotropic periodic sum (IPS) method. Our scheme captures the pH‐dependent behavior of binding free energy successfully. Absolute binding free energy values obtained with the PME and IPS methods are consistent, while cutoff method results are off by 2 kcal mol⁻¹. We also discuss the characteristics of three long‐range interaction calculation methods for constant‐pH simulations.

Binding affinities by alchemical perturbation using QM/MM with a large QM system and polarizable MM model

The most general way to improve the accuracy of binding-affinity calculations for protein-ligand systems is to use quantum-mechanical (QM) methods together with rigorous alchemical-perturbation (AP) methods. We explore this approach by calculating the relative binding free energy of two synthetic disaccharides binding to galectin-3 at a reasonably high QM level (dispersion-corrected density functional theory with a triple-zeta basis set) and with a sufficiently large QM system to include all short-range interactions with the ligand (744-748 atoms). The rest of the protein is treated as a collection of atomic multipoles (up to quadrupoles) and polarizabilities. Several methods for evaluating the binding free energy from the 3600 QM calculations are investigated in terms of stability and accuracy. In particular, methods using QM calculations only at the endpoints of the transformation are compared with the recently proposed non-Boltzmann Bennett acceptance ratio (NBB) method that uses QM calculations at several stages of the transformation. Unfortunately, none of the rigorous approaches give sufficient statistical precision. However, a novel approximate method, involving the direct use of QM energies in the Bennett acceptance ratio method, gives similar results as NBB but with better precision, ∼3 kJ/mol. The statistical error can be further reduced by performing a greater number of QM calculations.

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.