Search

Large-scale conformational transitions in multi-domain proteins are often essential for their functions. To investigate the transitions, it is necessary to explore multiple potential pathways, which involve different intermediate structures. Here, we present a multi-basin (MB) coarse-grained (CG) structure-based Go̅ model for describing transitions in proteins with more than two moving domains. This model is an extension of our dual-basin Go̅ model in which system-dependent parameters are determined systematically using the multistate Bennett acceptance ratio method. In the MB Go̅ model for multi-domain proteins, we assume that intermediate structures may have partial inter-domain native contacts. This approach allows us to search multiple transition pathways that involve distinct intermediate structures using the CG molecular dynamics (MD) simulations. We apply this scheme to an enzyme, adenylate kinase (AdK), which has three major domains and can move along two different pathways. Using the optimized mixing parameters for each pathway, AdK shows frequent transitions between the Open, Closed, and the intermediate basins and samples a wide variety of conformations within each basin. The explored multiple transition pathways could be compared with experimental data and examined in more detail by atomistic MD simulations.

Free-energy calculations based on atomistic Hamiltonians and sampling are key to a first-principles understanding of biomolecular processes, material properties, and macromolecular chemistry. Here, we generalize the free-energy perturbation method and derive nonlinear Hamiltonian transformation sequences yielding free-energy estimates with minimal mean squared error with respect to the exact values. Our variational approach applies to finite sampling and holds for any finite number of intermediate states. We show that our sequences are also optimal for the Bennett acceptance ratio (BAR) method, thereby generalizing BAR to small sampling sizes and non-Gaussian error distributions.

The binding of a series of small organic molecules, acting as agonists of the cannabinoid receptor CB1, was investigated by means of three methods of computational chemistry. Binding modes were predicted by means of molecular docking, and binding free energy was estimated via docking, molecular-mechanics Poisson-Boltzmann surface area method, and multistate Bennett acceptance ratio. No evident correlation was observed for the molecules between the experimental characteristics of affinity and three computed binding free energy estimates. The reasons for the discrepancy were discussed.

Free energy perturbation (FEP) simulations have been widely applied to obtain predictions of the relative binding free energy for a series of congeneric ligands binding to the same receptor, which is an essential component for the lead optimization process in computer-aided drug discovery. In the case of several congeneric ligands forming a perturbation map involving a closed thermodynamic cycle, the summation of the estimated free energy change along each edge in the cycle using Bennett acceptance ratio (BAR) usually will deviate from zero due to systematic and random errors, which is the hysteresis of cycle closure. In this work, the advanced reweighting techniques binless weighted histogram analysis method (UWHAM) and locally weighted histogram analysis method (LWHAM) are applied to provide statistical estimators of the free energy change along each edge in order to eliminate the hysteresis effect. As an example, we analyze a closed thermodynamic cycle involving four congeneric ligands which bind to HIV-1 integrase, a promising target which has emerged for antiviral therapy. We demonstrate that, compared with FEP and BAR, more accurate and hysteresis-free estimates of free energy differences can be achieved by using UWHAM to find a single estimate of the density of states based on all of the data in the cycle. Furthermore, by comparison of LWHAM results obtained from the inclusion of different numbers of neighboring states with UWHAM estimation involving all the states, we show how to determine the optimal neighborhood size in the LWHAM analysis to balance the trade-offs between computational cost and accuracy of the free energy prediction. Even with the smallest neighborhood, LWHAM can improve the BAR free energy estimates using the same input data as BAR. We introduce an overlapping states matrix that is constructed by using the global jump formula of LWHAM and plot its heat map. The heat map provides a quantitative measure of the overlap between pairs of alchemical/thermodynamic states. We explain how to identify and improve the FEP calculations along the edges that most likely cause large systematic errors by using the heat map of the overlapping states matrix and by comparing the BAR and UWHAM estimates of the free energy change.

In a survey of novel interactions between an IgG1 antibody and different Fcγ receptors (FcγR), molecular dynamics simulations were performed of interactions of monoclonal antibody involved complexes with FcγRs. Free energy simulations were also performed of isolated wild-type and substituted Fc regions bound to FcγRs with the aim of assessing their relative binding affinities. Two different free energy calculation methods, Molecular Mechanical/Generalized Born Molecular Volume (MM/GBMV) and Bennett Acceptance Ratio (BAR), were used to evaluate the known effector substitution G236A that is known to selectively increase antibody dependent cellular phagocytosis. The obtained results for the MM/GBMV binding affinity between different FcγRs are in good agreement with previous experiments, and those obtained using the BAR method for the complete antibody and the Fc-FcγR simulations show increased affinity across all FcγRs when binding to the substituted antibody. The FcγRIIa, a key determinant of antibody agonistic efficacy, shows a 10-fold increase in binding affinity, which is also consistent with the published experimental results. Novel interactions between the Fab region of the antibody and the FcγRs were discovered with this in silico approach, and provide insights into the antibody-FcγR binding mechanism and show promise for future improvements of therapeutic antibodies for preclinical studies of biological drugs.

The overall charge content and the patterning of charged residues have a profound impact on the conformational ensembles adopted by intrinsically disordered proteins. These parameters can be altered by charge regulation, which refers to the effects of post-translational modifications, pH-dependent changes to charge, and conformational fluctuations that modify the pKa values of ionizable residues. Although atomistic simulations have played a prominent role in uncovering the major sequence-ensemble relationships of IDPs, most simulations assume fixed charge states for ionizable residues. This may lead to erroneous estimates for conformational equilibria if they are linked to charge regulation. Here, we report the development of a new method we term q-canonical Monte Carlo sampling for modeling the linkage between charge regulation and conformational equilibria. The method, which is designed to be interoperable with the ABSINTH implicit solvation model, operates as follows: For a protein sequence with n ionizable residues, we start with all 2n charge microstates and use a criterion based on model compound pKa values to prune down to a subset of thermodynamically relevant charge microstates. This subset is then grouped into mesostates, where all microstates that belong to a mesostate have the same net charge. Conformational distributions, drawn from a canonical ensemble, are generated for each of the charge microstates that make up a mesostate using a method we designate as proton walk sampling. This method combines Metropolis Monte Carlo sampling in conformational space with an auxiliary Markov process that enables interconversions between charge microstates along a mesostate. Proton walk sampling helps identify the most likely charge microstate per mesostate. We then use thermodynamic integration aided by the multistate Bennett acceptance ratio method to estimate the free energies for converting between mesostates. These free energies are then combined with the per-microstate weights along each mesostate to estimate standard state free energies and pH-dependent free energies for all thermodynamically relevant charge microstates. The results provide quantitative estimates of the probabilities and preferred conformations associated with every thermodynamically accessible charge microstate. We showcase the application of q-canonical sampling using two model systems. The results establish the soundness of the method and the importance of charge regulation in systems characterized by conformational heterogeneity.

Accurate methods to estimate free energies play an important role for studying diverse condensed-phase problems in chemistry and biochemistry. The most common methods used in conjunction with molecular dynamics (MD) and Monte Carlo statistical mechanics (MC) simulations are free energy perturbation (FEP) and thermodynamic integration (TI). For common applications featuring the conversion of one molecule to another, simulations are run in stages or multiple "λ-windows" to promote convergence of the results. For computation of absolute free energies of solvation or binding, calculations are needed in which the solute is typically annihilated in the solvent and in the complex. The present work addresses identification of optimal protocols for such calculations, specifically, the creation/annihilation of organic molecules in aqueous solution. As is common practice, decoupling of the perturbations for electrostatic and Lennard-Jones interactions was performed. Consistent with earlier reports, FEP calculations for molecular creations are much more efficient, while annihilations require many more windows and may converge to incorrect values. Strikingly, we find that as few as four windows may be adequate for creation calculations for solutes ranging from argon to ethylbenzene. For a larger druglike molecule, MIF180, which contains 22 non-hydrogen atoms and three rotatable bonds, 10 creation windows are found to be adequate to yield the correct free energy of hydration. Convergence is impeded with procedures that use any sampling in the annihilation direction, and there is no need for postprocessing methods such as the Bennett acceptance ratio (BAR).

Solvation free energies can be essential in the process of evaluating and developing force fields. In addition, they can be used to obtain a diversity of other thermodynamic properties such as Henry's law constants and partition coefficients, for instance. In this paper, we study the solvation of asphaltene-like molecules (polyaromatic hydrocarbons) in both aqueous and organic solvents by using a coarse-grained model known as the SAFT-γ Mie force field. This model relies on a top-down approach in which the force-field parameters are obtained using an equation of state. The use of solvation free energy calculations to evaluate this force field can help in improving the model and increasing the scale in which these simulations can be applied. The results presented here were obtained by carrying out molecular dynamics simulations, with the expanded ensemble method being applied to sample an alchemical path. For this, we employed a softcore variant of the Mie potential. From the output of these simulations, we estimated free-energy differences by means of the Multistate Bennett Acceptance Ratio method. The results with organic solvents exhibited low absolute deviations from experimental data. In turn, hydration free energy calculations required a binary interaction parameter to be estimated from molecular dynamics data, which displays a flaw in the top-down parameterization approach. Fortunately, however, a single binary parameter could be used for all pairs of polyaromatic hydrocarbons with water, proving that the SAFT-γ Mie force field exhibits a suitable transferability property that deserves further investigation.

Previously we proposed the equilibrium and nonequilibrium adaptive alchemical free energy simulation methods Optimum Bennett’s Acceptance Ratio (OBAR) and Optimum Crooks’ Equation (OCE). They are based on the statistically optimal bidirectional reweighting estimator named Bennett’s Acceptance Ratio (BAR) or Crooks’ Equation (CE). They perform initial sampling in the staging alchemical transformation and then determine the importance rank of different states via the time-derivative of the variance (TDV). The method is proven to give speedups compared with the equal time rule. In the current work, we extended the time derivative of variance guided adaptive sampling method to the configurational space, falling in the term of Steered MD (SMD). The SMD approach biasing physically meaningful collective variable (CV) such as one dihedral or one distance to pulling the system from one conformational state to another. By minimizing the variance of the free energy differences along the pathway in an optimized way, a new type of adaptive SMD (ASMD) is introduced. As exhibits in the alchemical case, this adaptive sampling method outperforms the traditional equal-time SMD in nonequilibrium stratification. Also, the method gives much more efficient calculation of potential of mean force than the selection criterion based ASMD scheme, which is proven to be more efficient than traditional SMD. The variance-linearly-dependent minus time derivative of overall variance proposed for OBAR and OCE criterion is extended to determine the importance rank of the nonequilibrium pulling in the configurational space. It is shown that the importance rank given by the standard deviation of the free energy difference is wrong, but by correcting it with the simulation time we obtain the true importance rank in nonequilibrium stratification. The OCE workflow is periodicity-of-CV dependent while ASMD is not. In the non-periodic CV case, the end-state discrimination in the SD rank is eliminated in the TDV rank, while in the periodic CV case the correction introduced in the TDV rank is not that significant. The performance is demonstrated in a dihedral flipping case and two distance pulling cases, accounting for periodic and non-periodic CVs, respectively.

In recent years, the design, development, and evaluation of several inhibitors of the BACE1 enzyme, as part of Alzheimer's treatment, have gathered the scientific community's interest. Here, a linear regression model was built using binding free energy calculations through the Bennett acceptance ratio method for 20 known inhibitors of the BACE1 enzyme, with a Pearson coefficient of R = 0.88 and R2 = 0.78. The validation of this model was verified employing eight additional random inhibitors, which also gave a linear correlation with R = 0.97 and R2 = 0.93. Furthermore, this linear regression model was also used for proposing the structure of four potential BACE1 inhibitors, and the most active of them gave a theoretical Kd = 10 nM. However, these molecules have not been synthesized yet. Our team used a total time of more than 800 ns for the Molecular Dynamics to carry out this study, and all the software used were freely available.

The dual-basin Gō-model is a structural-based coarsegrained model for simulating a conformational transition between two known structures of a protein. Two parameters are required to produce a dual-basin potential mixed using two single-basin potentials, although the determination of mixing parameters is usually not straightforward. Here, we have developed an efficient scheme to determine the mixing parameters using the Multistate Bennett Acceptance Ratio (MBAR) method after short simulations with a set of parameters. In the scheme, MBAR allows us to predict observables at various unsimulated conditions, which are useful to improve the mixing parameters in the next round of iterative simulations. The number of iterations that are necessary for obtaining the converged mixing parameters are significantly reduced in the scheme. We applied the scheme to two proteins, the glutamine binding protein and the ribose binding protein, for showing the effectiveness in the parameter determination. After obtaining the converged parameters, both proteins show frequent conformational transitions between open and closed states, providing the theoretical basis to investigate structure-dynamics-function relationships of the proteins.

The Jarzynski equality is one of the most widely celebrated and scrutinized nonequilibrium work theorems, relating free energy to the external work performed in nonequilibrium transitions. In practice, the required ensemble average of the Boltzmann weights of infinite nonequilibrium transitions is estimated as a finite sample average, resulting in the so-called Jarzynski estimator, (Formula presented.). Alternatively, the second-order approximation of the Jarzynski equality, though seldom invoked, is exact for Gaussian distributions and gives rise to the Fluctuation-Dissipation estimator (Formula presented.). Here we derive the parametric maximum-likelihood estimator (MLE) of the free energy (Formula presented.) considering unidirectional work distributions belonging to Gaussian or Gamma families, and compare this estimator to (Formula presented.). We further consider bidirectional work distributions belonging to the same families, and compare the corresponding bidirectional (Formula presented.) to the Bennett acceptance ratio ((Formula presented.)) estimator. We show that, for Gaussian unidirectional work distributions, (Formula presented.) is in fact the parametric MLE of the free energy, and as such, the most efficient estimator for this statistical family. We observe that (Formula presented.) and (Formula presented.) perform better than (Formula presented.) and (Formula presented.), for unidirectional and bidirectional distributions, respectively. These results illustrate that the characterization of the underlying work distribution permits an optimal use of the Jarzynski equality. Fil: Arrar, Mehrnoosh. Consejo Nacional de Investigaciones Científicas y Técnicas. Oficina de Coordinación Administrativa Ciudad Universitaria. Instituto de Química, Física de los Materiales, Medioambiente y Energía. Universidad de Buenos Aires. Facultad de Ciencias Exactas y Naturales. Instituto de Química, Física de los Materiales, Medioambiente y Energía; Argentina Fil: Boubeta, Fernando Martín. Consejo Nacional de Investigaciones Científicas y Técnicas. Oficina de Coordinación Administrativa Ciudad Universitaria. Instituto de Química, Física de los Materiales, Medioambiente y Energía. Universidad de Buenos Aires. Facultad de Ciencias Exactas y Naturales. Instituto de Química, Física de los Materiales, Medioambiente y Energía; Argentina Fil: Szretter Noste, María Eugenia. Consejo Nacional de Investigaciones Científicas y Técnicas; Argentina. Universidad de Buenos Aires. Facultad de Ciencias Exactas y Naturales. Instituto de Cálculo; Argentina Fil: Sued, Raquel Mariela. Universidad de Buenos Aires. Facultad de Ciencias Exactas y Naturales. Instituto de Cálculo; Argentina. Consejo Nacional de Investigaciones Científicas y Técnicas; Argentina Fil: Boechi, Leonardo. Consejo Nacional de Investigaciones Científicas y Técnicas; Argentina. Universidad de Buenos Aires. Facultad de Ciencias Exactas y Naturales. Instituto de Cálculo; Argentina Fil: Rodriguez, Daniela. Consejo Nacional de Investigaciones Científicas y Técnicas; Argentina. Universidad de Buenos Aires. Facultad de Ciencias Exactas y Naturales. Instituto de Cálculo; Argentina

Computer simulations of molecular models are powerful technique that have improved the under- standing of many biochemical phenomena. The method is frequently applied to study the motions of biological macromolecules such as protein and nucleic acids, which can be useful for interpreting the results of certain biophysical experiments. In this work, we have estimated the solvation free energy for light alkane (methane, ethane, propane and n-butane) dissolved in water and methanol respectively over a broad range of temperatures, from 275 K to 375 K, using molecular dynamics simulations. The alkane (methane, ethane, propane and n-butane), and methanol molecules are described by the OPLS-AA (Optimized Potentials for Liquid Simulations-All Atom) potential, while water is modeled by TIP3P (Transferable Intermolecular Potential with 3-Points) model. We have used the free energy perturbation method (Bennett Acceptance Ratio (BAR) method) for the calculation of free energy of solvation. The estimated values of solvation free energy of alkane in the corresponding solvents agree well with the available experimental data.BIBECHANA 16 (2019) 91-104

We redevelop the variational free energy profile (vFEP) method using a cardinal B-spline basis to extend the method for analyzing free energy surfaces (FESs) involving 3-or-more reaction coordinates. We also implemented software for evaluating high-dimensional profiles based on the multistate Bennett acceptance ratio (MBAR) method which constructs an unbiased probability density from global reweighting of the observed samples. The MBAR method takes advantage of a fast algorithm for solving the unbinned weighted histogram (UWHAM)/MBAR equations which replaces the solution of simultaneous equations with a nonlinear optimization of a convex function. We make use of cardinal B-splines and multiquadric radial basis functions to obtain smooth, differentiable MBAR profiles in arbitrary high dimensions. The cardinal B-spline vFEP and MBAR methods are compared using three example systems that examine 1-, 2-, and 3-dimensional profiles. Both methods are found to be useful and produce nearly indistinguishable results. The vFEP method is found to be 150 times faster than MBAR when applied to periodic 2-dimensional profiles, but the MBAR method is 4.5 times faster than vFEP when evaluating unbounded 3-dimensional profiles. In agreement with previous comparisons, we find the vFEP method produces superior FESs when the overlap between umbrella window simulations decreases. Finally, the associative reaction mechanism of hammerhead ribozyme is characterized using 3-, 4-, and 6-dimensional profiles, and the higher-dimensional profiles are found to have smaller reaction barriers by as much as 1.5 kcal/mol. The methods presented here have been implemented into the FE-ToolKit software package along with new methods for network-wide free energy analysis in drug discovery.

The fast and accurate calculation of standard binding free energy has many important applications. Existing methodologies struggle at balancing accuracy and efficiency. We introduce a new method to compute binding free energy using deep generative models and the Bennett acceptance ratio method (DeepBAR). Compared to the rigorous potential of mean force (PMF) approach that requires sampling from intermediate states, DeepBAR is an order-of-magnitude more efficient as demonstrated in a series of host-guest systems. Notably, DeepBAR is exact and does not suffer from approximations for entropic contributions used in methods such as the molecular mechanics energy combined with the generalized Born and surface area continuum solvation (MM/GBSA). We anticipate DeepBAR to be a valuable tool for computing standard binding free energy used in drug design.

Free energy calculations based on atomistic Hamiltonians provide microscopic insight into the thermodynamicdriving forces of biophysical or condensed matter systems. Many approaches use intermediate Hamiltoniansinterpolating between the two states for which the free energy difference is calculated. The Bennett acceptanceratio (BAR) and variationally derived intermediates (VI) methods are optimal estimator and intermediate states inthat 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 thatestablished uncertainty estimates are inaccurate. Whereas for VI to become optimal, less than seven samplesper state suffice in all cases; for BAR the required number increases unboundedly with decreasing configurationspace densities overlap of the end states. We show that for BAR, the required number of samples is relatedto the overlap through an inverse power law. Because this relation seems to hold universally and almostindependent of other system properties, these findings can guide the proper choice of estimators for free energycalculations.

All-atom molecular dynamics simulations with stratified alchemical free energy calculations were used to predict the octanol-water partition coefficient [Formula: see text] of eleven small molecules as part of the SAMPL6-[Formula: see text] blind prediction challenge using four different force field parametrizations: standard OPLS-AA with transferable charges, OPLS-AA with non-transferable CM1A charges, AMBER/GAFF, and CHARMM/CGenFF. Octanol parameters for OPLS-AA, GAFF and CHARMM were validated by comparing the density as a function of temperature, the chemical potential, and the hydration free energy to experimental values. The partition coefficients were calculated from the solvation free energy for the compounds in water and pure ("dry") octanol or "wet" octanol with 27 mol% water dissolved. Absolute solvation free energies were computed by thermodynamic integration (TI) and the multistate Bennett acceptance ratio with uncorrelated samples from data generated by an established protocol using 5-ns windowed alchemical free energy perturbation (FEP) calculations with the Gromacs molecular dynamics package. Equilibration of sets of FEP simulations was quantified by a new measure of convergence based on the analysis of forward and time-reversed trajectories. The accuracy of the [Formula: see text] predictions was assessed by descriptive statistical measures such as the root mean square error (RMSE) of the data set compared to the experimental values. Discarding the first 1 ns of each 5-ns window as an equilibration phase had a large effect on the GAFF data, where it improved the RMSE by up to 0.8 log units, while the effect for other data sets was smaller or marginally worsened the agreement. Overall, CGenFF gave the best prediction with RMSE 1.2 log units, although for only eight molecules because the current CGenFF workflow for Gromacs does not generate files for certain halogen-containing compounds. Over all eleven compounds, GAFF gave an RMSE of 1.5. The effect of using a mixed water/octanol solvent slightly decreased the accuracy for CGenFF and GAFF and slightly increased it for OPLS-AA. The GAFF and OPLS-AA results displayed a systematic error where molecules were too hydrophobic whereas CGenFF appeared to be more balanced, at least on this small data set.

Alchemical free energy calculations using conventional molecular dynamics and thermodynamic integration rely on simulations performed at fixed values of the coupling parameter λ. When multiple conformers in equilibrium are separated by high barriers in the space orthogonal to λ, proper convergence may require extremely long simulations. Four main strategies can be employed to address this orthogonal-sampling problem: (a) λ-variations, where λ can change along the simulations to circumvent barriers; (b) λ-extrapolations, where statistical information is transferred between λ-points; (c) specific biasing, where orthogonal barriers are reduced using a biasing potential designed specifically for the system; and (d) generic biasing, where orthogonal barriers are reduced using a generic approach. Here, we investigate the relative merits of the first three strategies considering two benchmark systems. The KXK system involves a mutation of the central residue in a tripeptide to a glycine and the XTP system involves a hydrogen-to-bromine mutation in the base of a nucleotide. Three sampling methods are compared, the latter two involving λ-variations: molecular dynamics simulations with fixed λ-points, Hamiltonian replica exchange, and the recently introduced conveyor belt method. Two free energy estimators are applied, the second one involving λ-extrapolations: thermodynamic integration with Simpson quadrature and the multistate Bennett acceptance ratio. Finally, three different seeding schemes are considered for the generation of the initial configurations. For both benchmark systems, λ-extrapolations are found to provide little gain, whereas λ-variations can significantly enhance the convergence. They are sufficient on their own if the orthogonal barriers are low in at least one state (e.g., the glycine state in KXK). However, if the orthogonal barriers are high over the entire λ-range (e.g., the XTP system), λ-variations are only effective when applied together with a specific biasing for introducing such a low-barrier state.

We calculate the absolute binding free energies of tetra-methylated octa-acids host-guest systems as a part of the SAMPL6 blind challenge (receipt ID vq30p). We employed two different free energy simulation methods, i.e., the umbrella sampling (US) and double decoupling method (DDM). The US method was used with the weighted histogram analysis method (WHAM) (US-WHAM scheme). In the DDM scheme, Hamiltonian replica-exchange method (HREM) was combined with the Bennett acceptance ratio (BAR) (HREM-BAR scheme). We obtained initial binding poses via molecular docking using GalaxyDock-HG program, which is developed for the SAMPL challenge. The root mean square deviation (RMSD) and the mean absolute deviations (MAD) using US-WHAM scheme were 1.33 and 1.02 kcal/mol, respectively. The MAD was the top among all submissions, however the correlation with respect to experiment was unexceptional. While the RMSD and MAD via HREM-BAR scheme were greater than US-WHAM scheme, (i.e., 2.09 and 1.76 kcal/mol), their correlations were slightly better than US-WHAM. The correlation between the two methods was high. Further discussion on the DDM method can be found in a companion paper by Han et al. (receipt ID 3z83m) in the same issue.

We investigate the ability of several free-energy calculation methods to combine two alchemical changes. We use Bennett acceptance ratio (BAR), thermodynamic integration (TI), extended TI (X-TI), and enveloping distribution sampling (EDS) to perturb a water molecule, which is restrained to an amino acid that is also being perturbed. In addition to these pairwise methods, we present two two-dimensional approaches, EDS-TI and two-dimensional TI (2D-TI). We compare feasibility, efficiency and usability of these methods in regard to our simple model system, which mimics the displacement of a water molecule in the active site of a protein on residue mutation. The correct treatment of structural water has been shown to greatly aid binding affinity calculations in some cases that remained elusive otherwise. This is of broad interest in, for example, drug design, and we conclude that thus far, the pairwise method BAR and also the newer X-TI remain the most suitable methods to treat this problem as long as few end states are involved. © 2018 Wiley Periodicals, Inc.

This study reports the results of binding free energy calculations for CB[8] host-guest systems in the SAMPL6 blind challenge (receipt ID 3z83m). Force-field parameters were developed specific for each of host and guest molecules to improve configurational sampling. We used quantum mechanical (QM) implicit solvent calculations and QM force matching to determine non-bonded (partial atomic charges) and bonded terms, respectively. Free energy calculations were carried out using the double-decoupling method (DDM) combined with Hamiltonian replica exchange method (HREM) and Bennett acceptance ratio (BAR) method. The root mean square error (RMSE) of the predicted values using DDM with respect to the experimental results was 4.32 kcal/mol. The coefficient of determination (R(2)) and Kendall rank coefficient (τ) were 0.49 and 0.52, respectively, highest of all submissions. In addition, these were compared to the results obtained by umbrella sampling (US) and weighted histogram analysis method (WHAM). Overall, DDM achieved a higher prediction accuracy than the US method. Results are discussed in terms of parameterization and free energy simulations.

Four thermodynamic properties of pure N2 and O2 fluids were calculated over a continuous temperature range in the near-critical region employing molecular dynamics combined with the multistate Bennett acceptance ratio (MBAR) technique. The non-asymptotic critical anomalies of the isothermal compressibility, the thermal expansion coefficient and the heat capacities at constant pressure and constant volume were determined in the 0.80Tc–1.20Tc temperature range along the 1.00Pc and 1.15Pc isobars. A single-site united-atom model and a two-site rigid rotor force fields were used in the simulations. The results show that effective force fields developed to describe accurately the behavior of liquids and dense gases in general can be used to obtain near-critical anomalies of the thermodynamic properties. Furthermore, the computed values are in very good agreement with NIST data when available. The results fill existing gaps of the NIST database in the near-critical vicinity of the critical temperature along the critical isobar, especially for the heat capacities. The MBAR technique was crucial to describe the continuous and sharp peak of each property as functions of the temperature. The results indicated that even with a simple force field such as the united-atom one-site approach, very good results could be obtained for the non-asymptotic critical anomalies of thermodynamic properties without using renormalization group techniques. The computation of accurate results was possible due to the availability of precise critical constants especially for the simpler one-site force field. The inclusion of the anisotropy in the force field does not necessarily improve the computed thermodynamic properties in the near critical region.

Following the previously proposed equilibrate-state sampling based adaptive sampling regime Optimum Bennett Acceptance Ratio (OBAR), we introduce its nonequilibrium extension, the Optimum Crooks' Equation (OCE) in the current work. The efficiency of the NonEquilibrium Work (NEW) stratification is improved by adaptively manipulating the significance of each nonequilibrium realization followed by importance sampling. As is exhibited in the equilibrium case, the nonequilibrium extension outperforms the simple equal time rule used in nonequilibrium stratification in the sense of minimizing the total variance of the free energy estimate. The speedup of this non-equal time rule is more than 1-fold. The Time Derivative of total Variance (TDV) proposed for the OBAR criterion is extended to determine the importance of each nonequilibrium transformation, which is linearly dependent on the variance. The TDV in the nonequilibrium case gives a totally different importance rank from the standard errors of the free energy differences and OBAR TDV due to the duration of nonequilibrium pulling being added into the OCE equation. The performance of the OCE workflow is demonstrated in the solvation of several small molecules with a series of lambda increments and relaxation times between successive perturbations. To the best of our knowledge, such a nonequilibrium adaptive sampling regime in alchemical transformation is unprecedented.

Understanding the solubility of flavonoids is of relevance in view of the broad application such as antioxidant and antitumor agents. The solubility of flavonoid rutin has been experimentally investigated in a homologous series of alcohols, acetone, n-heptane and water solvents. In this work theoretical evaluation of Gibbs free energy of solvation was performed at the Density Functional Theory (DFT) and Molecular Dynamics (MD) levels, with the ω B97x-D functional and multistate Bennett acceptance ratio (MBAR) method, respectively. The degree of solubility of a given solute among various solvents is dictated by solute-solvent interactions, which may be represented by the solvation energy. Our calculated relative energy values due to solvent effects, using quantum chemical methods (PCM and SMD solvent models) do not follow strictly the observed solubility trend. However, the PCM model including solute-solvent dispersion interaction, solute-solvent repulsion interaction and solute cavitation energy contribution to the total energy (named DIS,REP,CAV) yielded an overall agreement with experimental solubility profile, with exception of ethanol solvent. A similar satisfactory accordance with experiment was found for MBresults, with a discrepancy also found for 1-butanol, which deserve further investigation. Analysis of spatial and radial distribution functions (SDF and RDF) for pure solvents and the analysis of the contributions of vdW and electrostatic terms for the total free energy reveal that a balance between dispersive and electrostatic interactions between the rutin and the solvent is necessary to solubilize the flavonoid, which corroborates with the values of these energy contributions to the solvation-free energy calculated by the MBAR method. The use of explicit solvent molecules in DFT calculations is precluded due to the size, flexibility, and many sites for solvent interactions of rutin molecule. Our results indicate that in general DFT-PCM (DIS,REP,CAV) and MBAR levels of calculations can satisfactorily predict the solvation free energy and solubility trend for flavonoid rutin along a series of solvents, and also provide a reasonable estimate of relative solubilities of other flavonoids, information that can be relevant for further molecular modeling studies.

The indirect construction of the free energy landscape at Quantum mechanics (QM)/ molecular mechanics (MM) levels provides a feasible alternative to the direct QM/MM free energy simulations. The main idea under the indirect method is constructing a thermodynamic cycle, exploring the configurational space at a computationally efficient but less accurate low-level Hamiltonian, and performing an alchemical correction to obtain the thermodynamics at an accurate but computationally demanding high-level Hamiltonian. In our previous works, we developed a multi-dimensional nonequilibrium free energy simulation framework to obtain QM/MM free energy landscapes indirectly. Specifically, we considered obtaining semi-empirical QM (SQM) results by combining the MM results and the MM-to-SQM correction and obtaining the QM results by combining the SQM results and the SQM-to-QM correction. In this work, we explore the possibility of changing the region for electronic structure calculations in the multi-scale QM/MM treatment, which could also be considered as a change of the level of theory. More generally, the multi-dimensional nonequilibrium Hamiltonian-variation/perturbation framework could be used to obtain transformations between different Hamiltonians of interest, such as changing the QM theory, the size of the QM region, and the basis set simultaneously.

The original reaction move for the reaction ensemble Monte Carlo (RxMC) method is adapted to align both the position and orientation of inserted product molecules and deleted reactant molecules. The accuracy and efficiency of this move is demonstrated for xylene isomerization in vapor, liquid, and supercritical phases. Classical RxMC requires the ideal gas free energy of reaction ΔGrxnideal as an input. We compare three methods for computing ΔGrxnideal: using tabulated enthalpies and entropies of formation, using the harmonic oscillator and rigid rotor approximations and using QM/MM alchemical transformation combined with multistate Bennett acceptance ratio. We find that the tabulated free energies of reaction give the best agreement with experimental equilibrium compositions in bulk fluids. RxMC simulations in a carbon nanotube with an inner diameter of approximately 6 A show that p-xylene becomes the dominant isomer under confinement, an effect consistent with the production of p-xylene in the zeolite Z...

Calcium and magnesium ions play important roles in many physicochemical processes. To facilitate the investigation of phenomena related to these ions that occur over large length and time scales, a coarse-grained force field (CGFF) is developed for MgCl2 and CaCl2 aqueous solutions. The ions are modeled by CG beads with characteristics of hydration shells. To accurately describe the nonideal behavior of the solutions, osmotic coefficients in a wide range of concentrations were used as guidance for parametrization. The osmotic coefficients were obtained from the chemical potential increments of water calculated using the Bennett acceptance ratio (BAR) method. The result CGFF accurately reproduces experimental osmotic coefficients, densities, surface tensions, and cation–anion separations of calcium chloride and magnesium chloride solutions at molalities up to 3.0 mol/kg. As a preliminary application, the force field is applied to simulate aggregations of sodium dodecyl sulfate (SDS) in calcium chloride sol...

A method that can reliably predict protein–ligand binding free energies is essential for rational drug design. Much effort has been devoted to this field, but it remains challenging especially for flexible ligands. In this work, both a molecular mechanical (MM) method and a hybrid quantum mechanical/molecular mechanical (QM/MM) method have been applied in the study of the binding affinities of methyl-α-L-fucoside to Ralstonia solanacearum lectins. The free energy at the MM level was calculated using the double-decoupling method (DDM) and the free energy change at each step was calculated via a series of intermediate states using the Bennett acceptance ratio (BAR). The binding free energy agrees well with the experimental measurement, no matter whether the general AMBER force field or GLYCAM06j was applied to the ligand. Nonetheless, slow convergence for some intermediate states has been observed, which requires substantially longer simulations than were used in many other studies. The QM/MM free energy was calculated by thermodynamic perturbation (TP) from the MM states. This strategy has been shown to yield minimal variance for the calculated free energy without direct sampling at the QM/MM level in a previous study. However, after this MM-to-QM/MM correction, the agreement with the experimental value decreased. This study serves as an implication of the demand for substantially longer simulations for the alchemical process than those that were used in many other studies and for further improvement of QM/MM methods, especially the description of interactions between the QM and MM regions.

We have developed an enhanced conformational sampling method combining replica-exchange umbrella sampling (REUS) with Gaussian accelerated molecular dynamics (GaMD). REUS enhances the sampling along predefined reaction coordinates, while GaMD accelerates the conformational dynamics by adding a boost potential to the system energy. The method, which we call GaREUS (Gaussian accelerated replica-exchange umbrella sampling), enhances the sampling more efficiently than REUS or GaMD, while the computational resource for GaREUS is the same as that required for REUS. The two-step reweighting procedure using the multistate Bennett acceptance ratio method and the cumulant expansion for the exponential average is applied to the simulation trajectories for obtaining the unbiased free-energy landscapes. We apply GaREUS to the calculations of free-energy landscapes for three different cases: conformational equilibria of N-glycan, folding of chignolin, and conformational change of adenyl kinase. We show that GaREUS speeds up the convergences of free-energy calculations using the same amount of computational resources as REUS. The free-energy landscapes reweighted from the trajectories of GaREUS agree with previously reported ones. GaREUS is applicable to free-energy calculations of various biomolecular dynamics and functions with reasonable computational costs.

We introduce a means to implement the interface potential approach for computing wetting properties within a molecular dynamics framework. The general approach provides a means to determine the contact angle of a liquid droplet on a solid substrate in a mother vapor. We present a framework for implementing "spreading" and "drying" versions of the method within an isothermal-isobaric ensemble. Two free energy methods are considered: cumulative integration of average force profile and multistate Bennett acceptance ratio. An umbrella sampling strategy is used to restrain volume fluctuations and to ensure adequate sampling of a broad volume range. We explore implementation of the approach with the GROningen MAchine for Chemical Simulations and the Large-scale Atomic/Molecular Massively Parallel Simulator. We test the accuracy and efficiency of the method with models consisting of a monoatomic Lennard-Jones fluid in the vicinity of a structureless or atomistically detailed substrate. Our results show that one can successfully generate the drying potential within the framework pursued here. The efficiency of the method is strongly dependent upon how one handles the dynamics of the two confining walls. These decisions impact the rate of volume fluctuations, and therefore, the quality of the volume distributions collected. Our efforts to implement the spreading method with molecular dynamics alone proved unsuccessful. The rate at which the configuration space of the vapor phase evolves is insufficient. We show how one can overcome this challenge by implementing a coupled molecular dynamics/Monte Carlo approach. Finally, we show how one can determine the variation in interfacial properties with temperature and substrate strength.

Correct calculation of the variation of free energy upon base flipping is crucial in understanding the dynamics of DNA systems. The free-energy landscape along the flipping pathway gives the thermodynamic stability and the flexibility of base-paired states. Although numerous free-energy simulations are performed in the base flipping cases, no theoretically rigorous nonequilibrium techniques are devised and employed to investigate the thermodynamics of base flipping. In the current work, we report a general nonequilibrium stratification scheme for the efficient calculation of the free-energy landscape of base flipping in DNA duplex. We carefully monitor the convergence behavior of the equilibrium sampling based free-energy simulation and the nonequilibrium stratification and determine the empirical length of time blocks required for converged sampling. Comparison between the performances of the equilibrium umbrella sampling and the nonequilibrium stratification is given. The results show that nonequilibrium free-energy simulation achieves similar accuracy and efficiency compared with the equilibrium enhanced sampling technique in the base flipping cases. We further test a convergence criterion we previously proposed and it comes out that the convergence determined by this criterion agrees with those given by the time-invariant behavior of PMF and the nonlinear dependence of standard deviation on the sample size.

In this work, we present a methodology to determine phase coexistence lines for atomic and rigid molecular systems with an emphasis on solid-fluid and on solid-solid equilibria. Phase coexistence points are found by computing the absolute free energy for each candidate phase separately. For solid phases, a combination of the extended Einstein crystal and the Einstein molecule method is presented which constitutes a convenient way to compute the absolute free energy with fixed center of mass. We compare results from equilibrium methods-thermodynamic integration and reweighting using the multistate Bennett acceptance ratio estimator (MBAR)-with simulations using a nonequilibrium method and discuss their advantages and disadvantages. Once absolute free energies of different phases are available, they are combined with simulations performed in the isothermal isobaric ensemble and MBAR, which enables efficient, iterative tracing of coexistence lines. The method is applicable to both liquid-solid as well as solid-solid transitions and is comparably simple and convenient to apply since the same method (MBAR) is used to compute free energies and to trace the coexistence line. Furthermore, statistical uncertainties can readily be computed in a transparent manner. We apply the method to an atomic solid (fcc argon) as well as small molecular systems (methanol and water) using the LAMMPS simulation package. Our study shows that all methods can be used to reliably compute the absolute free energy of solid phases, while MBAR is the most flexible method with high statistical efficiency. We find the nonequilibrium method is an attractive choice since it is simple to set up and to postprocess and is, hence, less prone to errors. The presented workflow provides a flexible, efficient, and robust way to compute phase diagrams using openly available software.

The multistate Bennett acceptance ratio method (MBAR) and unbinned weighted histogram analysis method (UWHAM) are widely employed approaches to calculate relative free energies of multiple thermodynamic states that gain statistical precision by employing free energy contributions from configurations sampled at each of the simulated λ states. With the increasing availability of high throughput computing resources, a large number of configurations can be sampled from hundreds or even thousands of states. Combining sampled configurations from all states to calculate relative free energies requires the iterative solution of large scale MBAR/UWHAM equations. In the current work, we describe the development of a fast solver to iteratively solve these large scale MBAR/UWHAM equations utilizing our previous findings that the MBAR/UWHAM equations can be derived as a Rao-Blackwell estimator. The solver is implemented and distributed as a Python module called FastMBAR. Our benchmark results show that FastMBAR is more than 2 times faster than the currently, and widely used solver, pymbar, when it runs on a central processing unit (CPU) and more than 100 times faster than pymbar when it runs on a graphical processing unit (GPU). The significant speedup achieved by FastMBAR running on a GPU is useful not only for solving large scale MBAR/UWHAM equations but also for estimating uncertainty of calculated free energies using bootstrapping where the MBAR/UWHAM equations need to be solved multiple times.

Construction of free energy landscapes at 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.

Fast and accurate evaluation of free energy has broad applications from drug design to material engineering. Computing the absolute free energy is of particular interest since it allows the assessment of the relative stability between states without the use of intermediates. In this letter, we introduce a general framework for calculating the absolute free energy of a state. A key step of the calculation is the definition of a reference state with tractable deep generative models using locally sampled configurations. The absolute free energy of this reference state is zero by design. The free energy for the state of interest can then be determined as the difference from the reference. We applied this approach to both discrete and continuous systems and demonstrated its effectiveness. It was found that the Bennett acceptance ratio method provides more accurate and efficient free energy estimations than approximate expressions based on work. We anticipate the method presented here to be a valuable strategy for computing free energy differences.

Free energy of solvation (i.e. residual chemical potential) and Henry's law solubility constants were estimated for mono ethylene glycol (MEG), diethylene glycol (DEG) and triethylene glycol (TEG) with water and methane as solvents at both 298 K and 1 atm and at 283 K and 80 atm using molecular dynamics simulations. Three methods of calculating free energies of solvation were compared: free energy perturbation (FEP), thermodynamic integration (TI) and Bennett acceptance ratio (BAR) method. The solubility constants were estimated based on the obtained free energies and were compared to the previously reported values. We found that free energy differences obtained using the three applied methods agreed very well. For both MEG and DEG the estimated solubility constants were within the same order of magnitude as previously reported solubility constants and therefore in good agreement. The estimated solubility constants of TEG were found to be about an order of magnitude smaller than the lowest reported value, and would still yield a reasonable estimate of residual chemical potential due to the logarithmic dependence.