Soft-core potentials in thermodynamic integration: comparing one- and two-step transformations

Thermodynamic Decomposition of Solvation Free Energies with Particle Mesh Ewald and Long-Range Lennard-Jones Interactions in Grid Inhomogeneous Solvation Theory

Grid Inhomogeneous Solvation Theory (GIST) maps out solvation thermodynamic properties on a fine meshed grid and provides a statistical mechanical formalism for thermodynamic end-state calculations. However, differences in how long-range nonbonded interactions are calculated in molecular dynamics engines and in the current implementation of GIST have prevented precise comparisons between free energies estimated using GIST and those from other free energy methods such as thermodynamic integration (TI). Here, we address this by presenting PME-GIST, a formalism by which particle mesh Ewald (PME)-based electrostatic energies and long-range Lennard-Jones (LJ) energies are decomposed and assigned to individual atoms and the corresponding voxels they occupy in a manner consistent with the GIST approach. PME-GIST yields potential energy calculations that are precisely consistent with modern simulation engines and performs these calculations at a dramatically faster speed than prior implementations. Here, we apply PME-GIST end-state analyses to 32 small molecules whose solvation free energies are close to evenly distributed from 2 kcal/mol to -17 kcal/mol and obtain solvation energies consistent with TI calculations (R² = 0.99, mean unsigned difference 0.8 kcal/mol). We also estimate the entropy contribution from the second and higher order entropy terms that are truncated in GIST by the differences between entropies calculated in TI and GIST. With a simple correction for the high order entropy terms, PME-GIST obtains solvation free energies that are highly consistent with TI calculations (R² = 0.99, mean unsigned difference = 0.4 kcal/mol) and experimental results (R² = 0.88, mean unsigned difference = 1.4 kcal/mol). The precision of PME-GIST also enables us to show that the solvation free energy of small hydrophobic and hydrophilic molecules can be largely understood based on perturbations of the solvent in a region extending a few solvation shells from the solute. We have integrated PME-GIST into the open-source molecular dynamics analysis software CPPTRAJ.

Parametrization of Trivalent and Tetravalent Metal Ions for the OPC3, OPC, TIP3P-FB, and TIP4P-FB Water Models

Commonly seen in rare-earth chemistry and materials science, highly charged metal ions play key roles in many chemical processes. Computer simulations have become an important tool for scientific research nowadays. Meaningful simulations require reliable parameters. In the present work, we parametrized 18 M(III) and 6 M(IV) metal ions for four new water models (OPC3, OPC, TIP3P-FB, TIP4P-FB) in conjunction with each of the 12-6 and 12-6-4 nonbonded models. Similar to what was observed previously, issues with the 12-6 model can be fixed by using the 12-6-4 model. Moreover, the four new water models showed comparable performance or considerable improvement over the previous water models (TIP3P, SPC/E, and TIP4P_EW) in the same category (3-point or 4-point water models, respectively). Finally, we reported a study of a metalloprotein system demonstrating the capability of the 12-6-4 model to model metalloproteins. The reported parameters will facilitate accurate simulations of highly charged metal ions in aqueous solution.

Variational Method for Networkwide Analysis of Relative Ligand Binding Free Energies with Loop Closure and Experimental Constraints

We describe an efficient method for the simultaneous solution of all free energies within a relative binding free-energy (RBFE) network with cycle closure and experimental/reference constraint conditions using Bennett Acceptance Ratio (BAR) and Multistate BAR (MBAR) analysis. Rather than solving the BAR or MBAR equations for each transformation independently, the simultaneous solution of all transformations are obtained by performing a constrained minimization of a global objective function. The nonlinear optimization of the objective function is subjected to affine linear constraints that couple the free energies between the network edges. The constraints are used to enforce the closure of thermodynamic cycles within the RBFE network, and to enforce an additional set of linear constraint conditions demonstrated here to be subsets of (1 or 2) experimental values. We describe details of the practical implementation of the network BAR/MBAR procedure, including use of generalized coordinates in the minimization of the free-energy objective function, propagation of bootstrap errors from those coordinates, and performance and memory optimization. In some cases it is found that use of restraints in the optimization is more practical than use of generalized coordinates for enforcing constraint conditions. The fast BARnet and MBARnet methods are used to analyze the RBFEs of six prototypical protein-ligand systems, and it is shown that enforcement of cycle closure conditions reduces the error in the predictions only modestly, and further reduction in errors can be achieved when one or two experimental RBFEs are included in the optimization procedure. These methods have been implemented into FE-ToolKit, a new free-energy analysis toolkit. The BARnet/MBARnet framework presented here opens the door to new, more efficient and robust free-energy analysis with enhanced predictive capability for drug discovery applications.

Binding of calcium and magnesium to human cardiac troponin C

Cardiac muscle thin filaments are composed of actin, tropomyosin, and troponin that change conformation in response to Ca²⁺ binding, triggering muscle contraction. Human cardiac troponin C (cTnC) is the Ca²⁺-sensing component of the thin filament. It contains structural sites (III/IV) that bind both Ca²⁺ and Mg²⁺ and a regulatory site (II) that has been thought to bind only Ca²⁺. Binding of Ca²⁺ at this site initiates a series of conformational changes that culminate in force production. However, the mechanisms that underpin the regulation of binding at site II remain unclear. Here, we have quantified the interaction between site II and Ca²⁺/Mg²⁺ through isothermal titration calorimetry and thermodynamic integration simulations. Direct and competitive binding titrations with WT N-terminal cTnC and full-length cTnC indicate that physiologically relevant concentrations of both Ca²⁺/Mg²⁺ interacted with the same locus. Moreover, the D67A/D73A N-terminal cTnC construct in which two coordinating residues within site II were removed was found to have significantly reduced affinity for both cations. In addition, 1 mM Mg²⁺ caused a 1.4-fold lower affinity for Ca²⁺. These experiments strongly suggest that cytosolic-free Mg²⁺ occupies a significant population of the available site II. Interaction of Mg²⁺ with site II of cTnC likely has important functional consequences for the heart both at baseline as well as in diseased states that decrease or increase the availability of Mg²⁺, such as secondary hyperparathyroidism or ischemia, respectively.

Modeling the Binding Mechanism of Remdesivir, Favilavir, and Ribavirin to SARS-CoV-2 RNA-Dependent RNA Polymerase

Recent efforts to repurpose drugs to combat COVID-19 have identified Remdesivir as a candidate. It acts on the RNA-dependent, RNA polymerase (RdRp) of the SARS-CoV-2 virus, a protein complex responsible for mediating replication of the virus's genome. However, its exact action mechanism, and that of other nucleotide analogue inhibitors, is not known. In this study, we examine at the molecular level the interaction of this drug and that of similar nucleotide analogue inhibitors, ribavirin and favilavir, by relying on atomistic molecular simulations and advanced sampling. By analyzing the binding free energies of these different drugs, it is found that all of them bind strongly at the active site. Surprisingly, however, ribavirin and favilavir do not bind the nucleotide on the complementary strand as effectively and seem to act by a different mechanism than remdesivir. Remdesivir exhibits similar binding interactions to the natural base adenine. Moreover, by analyzing remdesivir at downstream positions of the RNA, we also find that, consistent with a "delayed" termination mechanism, additional nucleotides can be incorporated after remdesivir is added, and its highly polar 1'-cyano group induces a set of conformational changes that can affect the normal RdRp complex function. By analyzing the fluctuations of residues that are altered by remdesivir binding, and comparing them to those induced by lethal point mutations, we find a possible secondary mechanism in which remdesivir destabilizes the protein complex and its interactions with the RNA strands.

Excess chemical potential of thiophene in [C4MIM] [BF4, Cl, Br, CH3COO] ionic liquids, determined by molecular simulations

The excess chemical potential of thiophene in imidazolium-based ionic liquids [C₄mim][BF₄], [C₄mim][Cl], [C₄mim][Br], and [C₄mim][CH₃COO] were determined by means of molecular dynamics in conjunction with free energy perturbation techniques employing non-polarizable force fields at 300 K and 343.15 K. In addition, energetic and structural analysis were performed such as: interaction energies, averaged noncovalent interactions, radial, and combined distribution functions. The results from this work revealed that the ionic liquids (ILs) presenting the most favorable excess chemical potentials ([C₄mim][BF₄], [C₄mim][CH₃COO]) are associated with the strongest energetic interaction between the thiophene molecule and the ionic liquid anion, and with the weakest energetic interaction between the thiophene molecule and the ionic liquid cation.

Excess chemical potential of thiophene in imidazolium-based ionic liquids [C₄mim][BF₄], [C₄mim][Cl], [C₄mim][Br], and [C₄mim][CH₃COO] determined by molecular simulations.

Computational Prediction of Mutational Effects on SARS-CoV-2 Binding by Relative Free Energy Calculations

The ability of coronaviruses to infect humans is invariably associated with their binding strengths to human receptor proteins. Both SARS-CoV-2, initially named 2019-nCoV, and SARS-CoV were reported to utilize angiotensin-converting enzyme 2 (ACE2) as an entry receptor in human cells. To better understand the interplay between SARS-CoV-2 and ACE2, we performed computational alanine scanning mutagenesis on the "hotspot" residues at protein-protein interfaces using relative free energy calculations. Our data suggest that the mutations in SARS-CoV-2 lead to a greater binding affinity relative to SARS-CoV. In addition, our free energy calculations provide insight into the infectious ability of viruses on a physical basis and also provide useful information for the design of antiviral drugs.

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.

Protein Denaturation, Zero Entropy Temperature, and the Structure of Water around Hydrophobic and Amphiphilic Solutes

The hydrophobic effect plays a key role in many chemical and biological processes, including protein folding. Nonetheless, a comprehensive picture of the effect of temperature on hydrophobic hydration and protein denaturation remains elusive. Here, we study the effect of temperature on the hydration of model hydrophobic and amphiphilic solutes, through molecular dynamics, aiming at getting insight on the singular behavior of water, concerning the zero-entropy temperature, T_S, and entropy convergence, T_S^*, also observed for some proteins, upon denaturation. We show that, similar to hydrocarbons, polar amphiphilic solutes exhibit a T_S, although strongly dependent on solute-water interactions, opposite to hydrocarbons. Further, the temperature dependence of the hydration entropy, normalized by the solvent accessible surface area, is shown to be nearly solute size independent for hydrophobic but not for amphiphilic solutes, for similar reasons. These results are further discussed in the light of information theory (IT) and the structure of water around hydrophobic groups. The latter shows that the tetrahedral enhancement of some water molecules around hydrophobic groups, associated with the reduction of water defects, leads to the strengthening of the weakest hydrogen bonds, relative to bulk water. In addition, a larger tetrahedrality is found in low density water populations, demonstrating that pure water has encoded structural information, similar to that associated with hydrophobic hydration. The reversal of the hydration entropy dependence on the solute size, above T_S^*, is also analyzed and shown to be associated with a greater loss of water molecules exhibiting enhanced orientational order, in the coordination sphere of large solutes. Finally, the source of the differences between Kauzmann's "hydrocarbon model" on protein denaturation and hydrophobic hydration is discussed, with relatively large amphiphilic hydrocarbons seemingly displaying a more similar behavior to some globular proteins than aliphatic hydrocarbons.

A fast and high-quality charge model for the next generation general AMBER force field

The General AMBER Force Field (GAFF) has been broadly used by researchers all over the world to perform in silico simulations and modelings on diverse scientific topics, especially in the field of computer-aided drug design whose primary task is to accurately predict the affinity and selectivity of receptor-ligand binding. The atomic partial charges in GAFF and the second generation of GAFF (GAFF2) were originally developed with the quantum mechanics derived restrained electrostatic potential charge, but in practice, users usually adopt an efficient charge method, Austin Model 1-bond charge corrections (AM1-BCC), based on which, without expensive ab initio calculations, the atomic charges could be efficiently and conveniently obtained with the ANTECHAMBER module implemented in the AMBER software package. In this work, we developed a new set of BCC parameters specifically for GAFF2 using 442 neutral organic solutes covering diverse functional groups in aqueous solution. Compared to the original BCC parameter set, the new parameter set significantly reduced the mean unsigned error (MUE) of hydration free energies from 1.03 kcal/mol to 0.37 kcal/mol. More excitingly, this new AM1-BCC model also showed excellent performance in the solvation free energy (SFE) calculation on diverse solutes in various organic solvents across a range of different dielectric constants. In this large-scale test with totally 895 neutral organic solvent-solute systems, the new parameter set led to accurate SFE predictions with the MUE and the root-mean-square-error of 0.51 kcal/mol and 0.65 kcal/mol, respectively. This newly developed charge model, ABCG2, paved a promising path for the next generation GAFF development.

Improved Alchemical Free Energy Calculations with Optimized Smoothstep Softcore Potentials

Progress in the development of GPU-accelerated free energy simulation software has enabled practical applications on complex biological systems and fueled efforts to develop more accurate and robust predictive methods. In particular, this work re-examines concerted (a.k.a., one-step or unified) alchemical transformations commonly used in the prediction of hydration and relative binding free energies (RBFEs). We first classify several known challenges in these calculations into three categories: endpoint catastrophes, particle collapse, and large gradient-jumps. While endpoint catastrophes have long been addressed using softcore potentials, the remaining two problems occur much more sporadically and can result in either numerical instability (i.e., complete failure of a simulation) or inconsistent estimation (i.e., stochastic convergence to an incorrect result). The particle collapse problem stems from an imbalance in short-range electrostatic and repulsive interactions and can, in principle, be solved by appropriately balancing the respective softcore parameters. However, the large gradient-jump problem itself arises from the sensitivity of the free energy to large values of the softcore parameters, as might be used in trying to solve the particle collapse issue. Often, no satisfactory compromise exists with the existing softcore potential form. As a framework for solving these problems, we developed a new family of smoothstep softcore (SSC) potentials motivated by an analysis of the derivatives along the alchemical path. The smoothstep polynomials generalize the monomial functions that are used in most implementations and provide an additional path-dependent smoothing parameter. The effectiveness of this approach is demonstrated on simple yet pathological cases that illustrate the three problems outlined. With appropriate parameter selection, we find that a second-order SSC(2) potential does at least as well as the conventional approach and provides vast improvement in terms of consistency across all cases. Last, we compare the concerted SSC(2) approach against the gold-standard stepwise (a.k.a., decoupled or multistep) scheme over a large set of RBFE calculations as might be encountered in drug discovery.

Molecular characterization of ebselen binding activity to SARS-CoV-2 main protease

There is an urgent need to repurpose drugs against severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). Recent computational-experimental screenings have identified several existing drugs that could serve as effective inhibitors of the virus' main protease, M^pro, which is involved in gene expression and replication. Among these, ebselen (2-phenyl-1,2-benzoselenazol-3-one) appears to be particularly promising. Here, we examine, at a molecular level, the potential of ebselen to decrease M^pro activity. We find that it exhibits a distinct affinity for the catalytic region. Our results reveal a higher-affinity, previously unknown binding site localized between the II and III domains of the protein. A detailed strain analysis indicates that, on such a site, ebselen exerts a pronounced allosteric effect that regulates catalytic site access through surface-loop interactions, thereby inducing a reconfiguration of water hotspots. Together, these findings highlight the promise of ebselen as a repurposed drug against SARS-CoV-2.

Validation of Free Energy Methods in AMBER

Herein we provide high-precision validation tests of the latest GPU-accelerated free energy code in AMBER. We demonstrate that consistent free energy results are obtained in both the gas phase and in solution. We first show, in the context of thermodynamic integration (TI), that the results are invariant with respect to "split" (e.g., stepwise decharge-vdW-recharge) versus "unified" protocols. This brought to light a subtle inconsistency in previous versions of AMBER that was traced to the improper treatment of 1-4 vdW and electrostatic interactions involving atoms across the softcore boundary. We illustrate that under the assumption that the ensembles produced by different legs of the alchemical transformation between molecules A and B in the gas phase and aqueous phase are very small, the inconsistency in the relative hydration free energy ΔΔG_hydr[A → B] = ΔG_aq[A → B] - ΔG_gas[A → B] is minimal. However, for general cases where the ensembles are shown to be substantially different, as expected in ligand-protein binding applications, these errors can be large. Finally, we demonstrate that results for relative hydration free energy simulations are independent of TI or multistate Bennett's acceptance ratio (MBAR) analysis, invariant to the specific choice of the softcore region, and in agreement with results derived from absolute hydration free energy values.

Repulsive Soft-Core Potentials for Efficient Alchemical Free Energy Calculations

In alchemical free energy (FE) simulations, annihilation and creation of atoms are generally achieved with the soft-core potential that shifts the interparticle separations. While this soft-core potential eliminates the numerical instability occurring near the two end states of the transformation, it makes the hybrid Hamiltonian vary nonlinearly with respect to the parameter λ, which interpolates between the Hamiltonians representing the two end states. This complicates FE estimation by Bennett acceptance ratio (BAR), free energy perturbation (FEP), and thermodynamic integration (TI) methods, thus reducing their calculation efficiency. In this work, we develop a new type of repulsive soft-core potential, called Gaussian soft-core (GSC) potential, with two parameters controlling its maximum and width. The main advantage of this potential is the linearity of the hybrid Hamiltonian with respect to λ, thus permitting the direct application of BAR, FEP, TI, and other variant FE methods. The accuracy and efficiency of the GSC potential are demonstrated by comparing the free energies of annihilation determined for 13 small molecules and an alchemical mutation of a protein side chain. In addition, in combination with a TI integrand (∂H/∂λ) estimation strategy, we show that GSC can considerably reduce the number of λ simulations compared to the commonly used separation-shifted soft-core potential.

Mechanism of Fully Reversible, pH-Sensitive Inhibition of Human Glutamine Synthetase by Tyrosine Nitration

Glutamine synthetase (GS) catalyzes an ATP-dependent condensation of glutamate and ammonia to form glutamine. This reaction-and therefore GS-are indispensable for the hepatic nitrogen metabolism. Nitration of tyrosine 336 (Y336) inhibits human GS activity. GS nitration and the consequent loss of GS function are associated with a broad range of neurological diseases. The mechanism by which Y336 nitration inhibits GS, however, is not understood. Here, we show by means of unbiased MD simulations, binding, and configurational free energy computations that Y336 nitration hampers ATP binding but only in the deprotonated and negatively charged state of residue 336. By contrast, for the protonated and neutral state, our computations indicate an increased binding affinity for ATP. pK_a computations of nitrated Y336 within GS predict a pK_a of ∼5.3. Thus, at physiological pH, nitrated Y336 exists almost exclusively in the deprotonated and negatively charged state. In vitro experiments confirm these predictions, in that, the catalytic activity of nitrated GS is decreased at pH 7 and 6 but not at pH 4. These results indicate a novel, fully reversible, pH-sensitive mechanism for the regulation of GS activity by tyrosine nitration.

A Polarization-Consistent Model for Alcohols to Predict Solvation Free Energies

Classical nonpolarizable models, normally based on a combination of Lennard-Jones sites and point charges, are extensively used to model thermodynamic properties of fluids, including solvation. An important shortcoming of these models is that they do not explicitly account for polarization effects, i.e., a description of how the electron density responds to changes in the molecular environment. Instead, polarization is implicitly included, in a mean-field sense, into the parameters of the model, usually by fitting to pure liquid properties (e.g., density). This causes problems when trying to describe thermodynamic properties that involve a change of phase (e.g., enthalpy of vaporization), that directly depend on the electronic response of the medium (e.g., dielectric constant), and that require mixing or solvation in different media (e.g., solvation free energies). Fully polarizable models present a natural route for addressing these limitations but at the price of a much higher computational cost. In this work, we combine the best of those two approaches by running fast simulations using nonpolarizable models and applying post facto corrections to the computed properties in order to account for the effects of polarization. By applying this new paradigm, a new united-atom force field for alcohols is developed that is able to predict both pure liquid properties, including dielectric constant, and solvation free energies in different solvents with a high degree of accuracy. This paves the way for the development of a generic classical nonpolarizable force field that can predict solvation of drug-like molecules in a variety of solvents.

Best Practices for Alchemical Free Energy Calculations [Article v1.0]

Alchemical free energy calculations are a useful tool for predicting free energy differences associated with the transfer of molecules from one environment to another. The hallmark of these methods is the use of "bridging" potential energy functions representing alchemical intermediate states that cannot exist as real chemical species. The data collected from these bridging alchemical thermodynamic states allows the efficient computation of transfer free energies (or differences in transfer free energies) with orders of magnitude less simulation time than simulating the transfer process directly. While these methods are highly flexible, care must be taken in avoiding common pitfalls to ensure that computed free energy differences can be robust and reproducible for the chosen force field, and that appropriate corrections are included to permit direct comparison with experimental data. In this paper, we review current best practices for several popular application domains of alchemical free energy calculations performed with equilibrium simulations, in particular relative and absolute small molecule binding free energy calculations to biomolecular targets.

Atomic insights into the effects of pathological mutants through the disruption of hydrophobic core in the prion protein

Destabilization of prion protein induces a conformational change from normal prion protein (PrPC) to abnormal prion protein (PrPSC). Hydrophobic interaction is the main driving force for protein folding, and critically affects the stability and solvability. To examine the importance of the hydrophobic core in the PrP, we chose six amino acids (V176, V180, T183, V210, I215, and Y218) that make up the hydrophobic core at the middle of the H2-H3 bundle. A few pathological mutants of these amino acids have been reported, such as V176G, V180I, T183A, V210I, I215V, and Y218N. We focused on how these pathologic mutations affect the hydrophobic core and thermostability of PrP. For this, we ran a temperature-based replica-exchange molecular dynamics (T-REMD) simulation, with a cumulative simulation time of 28 μs, for extensive ensemble sampling. From the T-REMD ensemble, we calculated the protein folding free energy difference between wild-type and mutant PrP using the thermodynamic integration (TI) method. Our results showed that pathological mutants V176G, T183A, I215V, and Y218N decrease the PrP stability. At the atomic level, we examined the change in pair-wise hydrophobic interactions from valine-valine to valine-isoleucine (and vice versa), which is induced by mutation V180I, V210I (I215V) at the 180th-210th (176th-215th) pair. Finally, we investigated the importance of the π-stacking between Y218 and F175.

Kinase Activation by Small Conformational Changes

Protein kinases (PKs) are allosteric enzymes that play an essential role in signal transduction by regulating a variety of key cellular processes. Most PKs suffer conformational rearrangements upon phosphorylation that strongly enhance the catalytic activity. Generally, it involves the movement of the phosphorylated loop toward the active site and the rotation of the whole C-terminal lobe. However, not all kinases undergo such a large configurational change: The MAPK extracellular signal-regulated protein kinases ERK1 and ERK2 achieve a 50 000 fold increase in kinase activity with only a small motion of the C-terminal region. In the present work, we used a combination of molecular simulation tools to characterize the conformational landscape of ERK2 in the active (phosphorylated) and inactive (unphosphorylated) states in solution in agreement with NMR experiments. We show that the chemical reaction barrier is strongly dependent on ATP conformation and that the "active" low-barrier configuration is subtly regulated by phosphorylation, which stabilizes a key salt bridge between the conserved Lys52 and Glu69 belonging to helix-C and promotes binding of a second Mg ion. Our study highlights that the on-off switch embedded in the kinase fold can be regulated by small, medium, and large conformational changes.

Improved Alignment of PEDOT:PSS Induced by in-situ Crystallization of "Green" Dimethylsulfone Molecules to Enhance the Polymer Thermoelectric Performance

Dimethylsulfone (DMSO2), a small organic molecule, was observed to induce the alignment of poly(3,4-ethylenedioxythiophene): poly(4-styrenesulfonate) (PEDOT:PSS) via in-situ crystallization in PEDOT:PSS mixture, which was verified by field emission scanning electron microscopy (FESEM), X-ray diffraction (XRD) and atomic force microscopy (AFM). A chemically stable dopant, DMSO2, remarkably raised the electrical conductivity of the PEDOT:PSS film, which was fabricated from pre-mixed solution of PEDOT:PSS and DMSO2, up to 1080 S/cm, and more importantly, such a PEDOT:PSS film showed a long-term humidity stability and it retained near 90% electric conductivity after 60 days, suggesting DMSO2 is promising for an eco-friendly alternative to replace dimethyl sulfoxide (DMSO), ethylene glycol (EG) and various acids dopants that have been widely employed to dope and post-treat PEDOT:PSS. Pairwise interaction energies and free energy of solvation between PEDOT:PSS and DMSO2 were calculated by first-principles and molecular mechanics, respectively, revealing the mechanism of DMSO2 in enhancing the electrical conductivity.

Exploring the Effectiveness of Binding Free Energy Calculations

Increasing the accuracy of the evaluation of ligand-binding energies is one of the most important tasks of current computational biology. Here we explore the accuracy of free energy perturbation (FEP) approaches by comparing the performance of a "regular" FEP method to the one using replica exchange to enhance the sampling on a well-defined benchmark. The examination was limited to the so-called alchemical perturbations which are restricted to a fragment of the drug, and therefore, the calculation is a relative one rather than the absolute binding energy of the drug. Overall, our calculations reach the 1 kcal/mol accuracy limit. It is also shown that the accurate prediction of the position of water molecules around the binding pocket is important for FEP calculations. Interestingly, the replica exchange method does not significantly improve the accuracy of binding energies, suggesting that we reach the limit where the force field quality is a critical factor for accurate calculations.

Blinded prediction of protein-ligand binding affinity using Amber thermodynamic integration for the 2018 D3R grand challenge 4

In the framework of the 2018 Drug Design Data Resource grand challenge 4, blinded predictions on relative binding free energy were performed for a set of 39 ligands of the Cathepsin S protein. We leveraged the GPU-accelerated thermodynamic integration of Amber 18 to advance our computational prediction. When our entry was compared to experimental results, a good correlation was observed (Kendall's τ: 0.62, Spearman's ρ: 0.80 and Pearson's R: 0.82). We designed a parallelized transformation map that placed ligands into several groups based on common alchemical substructures; TI transformations were carried out for each ligand to the relevant substructure, and between substructures. Our calculations were all conducted using the linear potential scaling scheme in Amber TI because we believe the softcore potential/dual-topology approach as implemented in current Amber TI is highly fault-prone for some transformations. The issue is illustrated by using two examples in which typical preparation for the dual-topology approach of Amber TI fails. Overall, the high accuracy of our prediction is a result of recent advances in force fields (ff14SB and GAFF), as well as rapid calculation of ensemble averages enabled by the GPU implementation of Amber. The success shown here in a blinded prediction strongly suggests that alchemical free energy calculation in Amber is a promising tool for future commercial drug design.

Using AMBER18 for Relative Free Energy Calculations

With renewed interest in free energy methods in contemporary structure-based drug design, there is a pressing need to validate against multiple targets and force fields to assess the overall ability of these methods to accurately predict relative binding free energies. We computed relative binding free energies using graphics processing unit accelerated thermodynamic integration (GPU-TI) on a data set originally assembled by Schrödinger, Inc. Using their GPU free energy code (FEP+) and the OPLS2.1 force field combined with the REST2 enhanced sampling approach, these authors obtained an overall MUE of 0.9 kcal/mol and an overall RMSD of 1.14 kcal/mol. In our study using GPU-TI from AMBER with the AMBER14SB/GAFF1.8 force field but without enhanced sampling, we obtained an overall MUE of 1.17 kcal/mol and an overall RMSD of 1.50 kcal/mol for the 330 perturbations contained in this data set. A more detailed analysis of our results suggested that the observed differences between the two studies arise from differences in sampling protocols along with differences in the force fields employed. Future work should address the problem of establishing benchmark quality results with robust statistical error bars obtained through multiple independent runs and enhanced sampling, which is possible with the GPU-accelerated features in AMBER.

Energetic, Structural and Dynamic Properties of Nucleobase-Urea Interactions that Aid in Urea Assisted RNA Unfolding

Understanding the structure-function relationships of RNA has become increasingly important given the realization of its functional role in various cellular processes. Chemical denaturation of RNA by urea has been shown to be beneficial in investigating RNA stability and folding. Elucidation of the mechanism of unfolding of RNA by urea is important for understanding the folding pathways. In addition to studying denaturation of RNA in aqueous urea, it is important to understand the nature and strength of interactions of the building blocks of RNA. In this study, a systematic examination of the structural features and energetic factors involving interactions between nucleobases and urea is presented. Results from molecular dynamics (MD) simulations on each of the five DNA/RNA bases in water and eight different concentrations of aqueous urea, and free energy calculations using the thermodynamic integration method are presented. The interaction energies between all the nucleobases with the solvent environment and the transfer free energies become more favorable with respect to increase in the concentration of urea. Preferential interactions of urea versus water molecules with all model systems determined using Kirkwood-Buff integrals and two-domain models indicate preference of urea by nucleobases in comparison to water. The modes of interaction between urea and the nucleobases were analyzed in detail. In addition to the previously identified hydrogen bonding and stacking interactions between urea and nucleobases that stabilize the unfolded states of RNA in aqueous solution, NH-π interactions are proposed to be important. Dynamic properties of each of these three modes of interactions have been presented. The study provides fundamental insights into the nature of interaction of urea molecules with nucleobases and how it disrupts nucleic acids.

Understanding Aldose Reductase-Inhibitors interactions with free energy simulation

Aldose Reductase (AR) reduces a variety of substrates, such as aldehydes, aldoses and corticosteroids. It is the first and rate-limiting enzyme of the polyol pathway and is an important target enzyme for diabetes-associated complications, including retinopathy, neuropathy, and nephropathy. Inhibitors targeting this enzyme are structurally different and some of them have side effects. In existing publications, computational techniques are applied to investigate the binding affinities of existing inhibitors and predicting the affinities of newly designed ligands. However, these calculations only employ coarse and approximated methods such as docking and MM/PBSA. Brute force simulations are employed to study the dynamics of the system but no converged statistics is obtained. As a result, these computations provide results not consistent with experimental values and large discrepancies exist. In the current work, we employ the enhanced sampling technique of alchemical free energy simulation to calculate the binding affinities of several ligands targeting AR. The statistical error is defined with care and the correlation in the time-series data is fully considered. The statistically optimal estimators are used to extract the free energy estimates and the predicted results are in agreement with the experimental values. Less computationally demanding end-point free energy methods are also performed to compare their efficiency with the alchemical methods. As is expected, the end-point methods are of less accuracy and reliability compared with the alchemical free energy methods. The decomposition of the free energy difference in each alchemical transformation into the enthalpic and entropic components gives further insights on the thermodynamics. The enthalpy-entropy compensation is observed in this case. As the structural data obtained from experiments are only snapshots and more details are needed to understand the dynamics of the protein-ligand system, the conformational ensemble is analyzed. We identify important residues involved in the protein-ligand binding case and short-lived interactions formed due to fluctuations in the conformational ensemble. The current work shed light on the atomic detailed understanding of the dynamics of AR-inhibitors interactions.

Molecular insight into chymotrypsin inhibitor 2 resisting proteolytic degradation

Chymotrypsin inhibitor 2 (CI2) is a special serine protease inhibitor which can resist hydrolysis for several days with a rapid equilibrium between the Michaelis complex and acyl-enzyme intermediate. The energies and conformational changes for subtilisin-catalyzed proteolysis of CI2 were examined in this paper for the first time by employing pseudo bond ab initio QM/MM MD simulations. In the acylation reaction, a low-barrier hydrogen bond between His64 and Asp32 in the transition state together with the lack of covalent backbone constraints makes the peptide bonds of CI2 break more easily than in other serine protease inhibitors. After acyl-enzyme formation, molecular dynamics simulations showed that the access of hydrolytic water to the active site requires partial dissociation of the leaving group. However, retention of the leaving group mainly by the intra- and inter-molecular H-bonding networks hinders the access of water and retards the deacylation reaction. Instead of the dissociation constant of inhibitors, we suggest employing the free energy at the acyl-enzyme state to predict the relative hydrolysis rates of CI2 mutants, which are testified by the experimental relative hydrolysis rates.

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.

Rapid and accurate structure-based therapeutic peptide design using GPU accelerated thermodynamic integration

Peptide-based therapeutics are an alternative to small molecule drugs as they offer superior specificity, lower toxicity, and easy synthesis. Here we present an approach that leverages the dramatic performance increase afforded by the recent arrival of GPU accelerated thermodynamic integration (TI). GPU TI facilitates very fast, highly accurate binding affinity optimization of peptides against therapeutic targets. We benchmarked TI predictions using published peptide binding optimization studies. Prediction of mutations involving charged side-chains was found to be less accurate than for non-charged, and use of a more complex 3-step TI protocol was found to boost accuracy in these cases. Using the 3-step protocol for non-charged side-chains either had no effect or was detrimental. We use the benchmarked pipeline to optimize a peptide binding to our recently discovered cancer target: EME1. TI calculations predict beneficial mutations using both canonical and non-canonical amino acids. We validate these predictions using fluorescence polarization and confirm that binding affinity is increased. We further demonstrate that this increase translates to a significant reduction in pancreatic cancer cell viability.

Reproducibility of Free Energy Calculations across Different Molecular Simulation Software Packages

Alchemical free energy calculations are an increasingly important modern simulation technique to calculate free energy changes on binding or solvation. Contemporary molecular simulation software such as AMBER, CHARMM, GROMACS, and SOMD include support for the method. Implementation details vary among those codes, but users expect reliability and reproducibility, i.e., for a given molecular model and set of force field parameters, comparable free energy differences should be obtained within statistical bounds regardless of the code used. Relative alchemical free energy (RAFE) simulation is increasingly used to support molecule discovery projects, yet the reproducibility of the methodology has been less well tested than its absolute counterpart. Here we present RAFE calculations of hydration free energies for a set of small organic molecules and demonstrate that free energies can be reproduced to within about 0.2 kcal/mol with the aforementioned codes. Absolute alchemical free energy simulations have been carried out as a reference. Achieving this level of reproducibility requires considerable attention to detail and package-specific simulation protocols, and no universally applicable protocol emerges. The benchmarks and protocols reported here should be useful for the community to validate new and future versions of software for free energy calculations.

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.

Binding free energies in the SAMPL6 octa-acid host-guest challenge calculated with MM and QM methods

We have estimated free energies for the binding of eight carboxylate ligands to two variants of the octa-acid deep-cavity host in the SAMPL6 blind-test challenge (with or without endo methyl groups on the four upper-rim benzoate groups, OAM and OAH, respectively). We employed free-energy perturbation (FEP) for relative binding energies at the molecular mechanics (MM) and the combined quantum mechanical (QM) and MM (QM/MM) levels, the latter obtained with the reference-potential approach with QM/MM sampling for the MM → QM/MM FEP. The semiempirical QM method PM6-DH+ was employed for the ligand in the latter calculations. Moreover, binding free energies were also estimated from QM/MM optimised structures, combined with COSMO-RS estimates of the solvation energy and thermostatistical corrections from MM frequencies. They were performed at the PM6-DH+ level of theory with the full host and guest molecule in the QM system (and also four water molecules in the geometry optimisations) for 10-20 snapshots from molecular dynamics simulations of the complex. Finally, the structure with the lowest free energy was recalculated using the dispersion-corrected density-functional theory method TPSS-D3, for both the structure and the energy. The two FEP approaches gave similar results (PM6-DH+/MM slightly better for OAM), which were among the five submissions with the best performance in the challenge and gave the best results without any fit to data from the SAMPL5 challenge, with mean absolute deviations (MAD) of 2.4-5.2 kJ/mol and a correlation coefficient (R²) of 0.77-0.93. This is the first time QM/MM approaches give binding free energies that are competitive to those obtained with MM for the octa-acid host. The QM/MM-optimised structures gave somewhat worse performance (MAD = 3-8 kJ/mol and R² = 0.1-0.9), but the results were improved compared to previous studies of this system with similar methods.

Charge transfer through a fragment of the respiratory complex I and its regulation: an atomistic simulation approach

We simulate electron transfer within a fragment of the NADH:ubiquinone oxidoreductase (respiratory complex I) of the hyperthermophilic bacterium Aquifex aeolicus. We apply molecular dynamics simulations, thermodynamic integration, and a thermodynamic network least squares analysis to compute two key parameters of Marcus' theory of charge transfer, the thermodynamic driving force and the reorganization energy. Intramolecular contributions to the Gibbs free energy differences of electron and hydrogen transfer processes, ΔG, are accessed by calibrating against experimental redox titration data. This approach permits the computation of the interactions between the species NAD+, FMNH2, N1a-, and N3-, and the construction of a free energy surface for the flow of electrons within the fragment. We find NAD+ to be a strong candidate for the regulation of charge transfer.

Assessing the stability of free-energy perturbation calculations by performing variations in the method

We have calculated relative binding affinities for eight tetrafluorophenyl-triazole-thiogalactoside inhibitors of galectin-3 with the alchemical free-energy perturbation approach. We obtain a mean absolute deviation from experimental estimates of only 2-3 kJ/mol and a correlation coefficient (R²) of 0.5-0.8 for seven relative affinities spanning a range of up to 11 kJ/mol. We also studied the effect of using different methods to calculate the charges of the inhibitor and different sizes of the perturbed group (the atoms that are described by soft-core potentials and are allowed to have differing coordinates). However, the various approaches gave rather similar results and it is not possible to point out one approach as consistently and significantly better than the others. Instead, we suggest that such small and reasonable variations in the computational method can be used to check how stable the calculated results are and to obtain a more accurate estimate of the uncertainty than if performing only one calculation with a single computational setup.

A GPU-Accelerated Parameter Interpolation Thermodynamic Integration Free Energy Method

There has been a resurgence of interest in free energy methods motivated by the performance enhancements offered by molecular dynamics (MD) software written for specialized hardware, such as graphics processing units (GPUs). In this work, we exploit the properties of a parameter-interpolated thermodynamic integration (PI-TI) method to connect states by their molecular mechanical (MM) parameter values. This pathway is shown to be better behaved for Mg2+ → Ca2+ transformations than traditional linear alchemical pathways (with and without soft-core potentials). The PI-TI method has the practical advantage that no modification of the MD code is required to propagate the dynamics, and unlike with linear alchemical mixing, only one electrostatic evaluation is needed (e.g., single call to particle-mesh Ewald) leading to better performance. In the case of AMBER, this enables all the performance benefits of GPU-acceleration to be realized, in addition to unlocking the full spectrum of features available within the MD software, such as Hamiltonian replica exchange (HREM). The TI derivative evaluation can be accomplished efficiently in a post-processing step by reanalyzing the statistically independent trajectory frames in parallel for high throughput. We also show how one can evaluate the particle mesh Ewald contribution to the TI derivative evaluation without needing to perform two reciprocal space calculations. We apply the PI-TI method with HREM on GPUs in AMBER to predict p Ka values in double stranded RNA molecules and make comparison with experiments. Convergence to under 0.25 units for these systems required 100 ns or more of sampling per window and coupling of windows with HREM. We find that MM charges derived from ab initio QM/MM fragment calculations improve the agreement between calculation and experimental results.

Saturation Mutagenesis by Efficient Free-Energy Calculation

Single-point mutations in proteins can greatly influence protein stability, binding affinity, protein function or its expression per se. Here, we present accurate and efficient predictions of the free energy of mutation of amino acids. We divided the complete mutational free energy into an uncharging step, which we approximate by a third-power fitting (TPF) approach, and an annihilation step, which we approximate using the one-step perturbation (OSP) method. As a diverse set of test systems, we computed the solvation free energy of all amino acid side chain analogues and obtained an excellent agreement with thermodynamic integration (TI) data. Moreover, we calculated mutational free energies in model tripeptides and established an efficient protocol involving a single reference state. Again, the approximate methods agreed excellently with the TI references, with a root-mean-square error of only 3.6 kJ/mol over 17 mutations. Our combined TPF+OSP approach does show not only a very good agreement but also a 2-fold higher efficiency than full blown TI calculations.

Towards a coarse-grained model of the peptoid backbone: the case of N,N-dimethylacetamide

Coarse-grained model of DMA, containing the basic motif of the peptoid backbone, based on short ranged many-body ranged interactions.

Absolute Alchemical Free Energy Calculations for Ligand Binding: A Beginner's Guide

Many thermodynamic quantities can be extracted from computer simulations that generate an ensemble of microstates according to the principles of statistical mechanics. Among these quantities is the free energy of binding of a small molecule to a macromolecule, such as a protein. Here, we present an introductory overview of a protocol that allows for the estimation of ligand binding free energies via molecular dynamics simulations. While we focus on the binding of organic molecules to proteins, the approach is in principle transferable to any pair of molecules.

Binding affinities of the farnesoid X receptor in the D3R Grand Challenge 2 estimated by free-energy perturbation and docking

We have studied the binding of 102 ligands to the farnesoid X receptor within the D3R Grand Challenge 2016 blind-prediction competition. First, we employed docking with five different docking software and scoring functions. The selected docked poses gave an average root-mean-squared deviation of 4.2 Å. Consensus scoring gave decent results with a Kendall's τ of 0.26 ± 0.06 and a Spearman's ρ of 0.41 ± 0.08. For a subset of 33 ligands, we calculated relative binding free energies with free-energy perturbation. Five transformations between the ligands involved a change of the net charge and we implemented and benchmarked a semi-analytic correction (Rocklin et al., J Chem Phys 139:184103, 2013) for artifacts caused by the periodic boundary conditions and Ewald summation. The results gave a mean absolute deviation of 7.5 kJ/mol compared to the experimental estimates and a correlation coefficient of R 2 = 0.1. These results were among the four best in this competition out of 22 submissions. The charge corrections were significant (7-8 kJ/mol) and always improved the results. By employing 23 intermediate states in the free-energy perturbation, there was a proper overlap between all states and the precision was 0.1-0.7 kJ/mol. However, thermodynamic cycles indicate that the sampling was insufficient in some of the perturbations.

Molecular Structure Inhibiting Synergism in Charged Surfactant Mixtures: An Atomistic Molecular Dynamics Simulation Study

Synergistic and nonsynergistic surfactant-water mixtures of sodium dodecyl sulfate (SDS), lauryl betaine (C12B), and cocoamidopropyl betaine (CAPB) systems are studied using molecular simulation to understand the role of interactions among headgroups, tailgroups, and water on structural and thermodynamic properties at the air-water interface. SDS is an anionic surfactant, while C12B and CAPB are zwitterionic; CAPB differs from C12B by an amide group in the tail. While the lowest surface tensions at high surface concentrations in the SDS-C12B synergistic system could not be reproduced by simulation, estimated partitioning between surface and bulk shows trends consistent with synergism. Structural analysis shows the influence of the SDS headgroup pulling C12B to the surface, resulting in closely packed structures compared to their respective homomolecular-surfactant systems. The SDS-CAPB system, on the other hand, is nonsynergistic when the surfactants are mixed on account of the tilted structure of the CAPB tail. The translational excess entropy due to the tailgroup interactions discriminates between the synergistic and nonsynergistic systems. The implications of such interactions on surfactant effects in complex, multicomponent atmospheric aerosols are discussed.