Chapter 4 Alchemical Free Energy Calculations: Ready for Prime Time?

Calculating Binding Free Energies in Model Host-Guest Systems with Unrestrained Advanced Sampling

Host-guest interactions are important to the design of pharmaceuticals and, more broadly, to soft materials as they can enable targeted, strong, and specific interactions between molecules. The binding process between the host and guest may be classified as a "rare event" when viewing the system at atomic scales, such as those explored in molecular dynamics simulations. To obtain equilibrium binding conformations and dissociation constants from these simulations, it is essential to resolve these rare events. Advanced sampling methods such as the adaptive biasing force (ABF) promote the occurrence of less probable configurations in a system, therefore facilitating the sampling of essential collective variables that characterize the host-guest interactions. Here, we present the application of ABF to a rod-cavitand coarse-grained model of host-guest systems to acquire the potential of mean force. We show that the employment of ABF enables the computation of the configurational and thermodynamic properties of bound and unbound states, including the free energy landscape. Moreover, we identify important dynamic bottlenecks that limit sampling and discuss how these may be addressed in more general systems.

Exploring Free Energies of Specific Protein Conformations Using the Martini Force Field

Coarse-grained (CG) level molecular dynamics simulations are routinely used to study various biomolecular processes. The Martini force field is currently the most widely adopted parameter set for such simulations. The functional form of this and several other CG force fields enforces secondary protein structure support by employing a variety of harmonic potentials or restraints that favor the protein's native conformation. We propose a straightforward method to calculate the energetic consequences of transitions between predefined conformational states in systems in which multiple factors can affect protein conformational equilibria. This method is designed for use within the Martini force field and involves imposing conformational transitions by linking a Martini-inherent elastic network to the coupling parameter λ. We demonstrate the applicability of our method using the example of five biomolecular systems that undergo experimentally characterized conformational transitions between well-defined structures (Staphylococcal nuclease, C-terminal segment of surfactant protein B, LAH4 peptide, and β2-adrenergic receptor) as well as between folded and unfolded states (GCN4 leucine zipper protein). The results show that the relative free energy changes associated with protein conformational transitions, which are affected by various factors, such as pH, mutations, solvent, and lipid membrane composition, are correctly reproduced. The proposed method may be a valuable tool for understanding how different conditions and modifications affect conformational equilibria in proteins.

Expanded ensemble predictions of absolute binding free energies in the SAMPL9 host-guest challenge

As part of the SAMPL9 community-wide blind host-guest challenge, we implemented an expanded ensemble workflow to predict absolute binding free energies for 13 small molecules against pillar[6]arene. Notable features of our protocol include consideration of a variety of protonation and enantiomeric states for both host and guests, optimization of alchemical intermediates, and analysis of free energy estimates and their uncertainty using large numbers of simulation replicates performed using distributed computing. Our predictions of absolute binding free energies resulted in a mean absolute error of 2.29 kcal mol-1 and an R2 of 0.54. Overall, results show that expanded ensemble calculations using all-atom molecular dynamics simulations are a valuable and efficient computational tool in predicting absolute binding free energies.

Quantitative approaches for decoding the specificity of the human T cell repertoire

T cell receptor (TCR)-peptide-major histocompatibility complex (pMHC) interactions play a vital role in initiating immune responses against pathogens, and the specificity of TCRpMHC interactions is crucial for developing optimized therapeutic strategies. The advent of high-throughput immunological and structural evaluation of TCR and pMHC has provided an abundance of data for computational approaches that aim to predict favorable TCR-pMHC interactions. Current models are constructed using information on protein sequence, structures, or a combination of both, and utilize a variety of statistical learning-based approaches for identifying the rules governing specificity. This review examines the current theoretical, computational, and deep learning approaches for identifying TCR-pMHC recognition pairs, placing emphasis on each method's mathematical approach, predictive performance, and limitations.

ACES: Optimized Alchemically Enhanced Sampling

We present an alchemical enhanced sampling (ACES) method implemented in the GPU-accelerated AMBER free energy MD engine. The methods hinges on the creation of an "enhanced sampling state" by reducing or eliminating selected potential energy terms and interactions that lead to kinetic traps and conformational barriers while maintaining those terms that curtail the need to otherwise sample large volumes of phase space. For example, the enhanced sampling state might involve transforming regions of a ligand and/or protein side chain into a noninteracting "dummy state" with internal electrostatics and torsion angle terms turned off. The enhanced sampling state is connected to a real-state end point through a Hamiltonian replica exchange (HREMD) framework that is facilitated by newly developed alchemical transformation pathways and smoothstep softcore potentials. This creates a counterdiffusion of real and enhanced-sampling states along the HREMD network. The effect of a differential response of the environment to the real and enhanced-sampling states is minimized by leveraging the dual topology framework in AMBER to construct a counterbalancing HREMD network in the opposite alchemical direction with the same (or similar) real and enhanced sampling states at inverted end points. The method has been demonstrated in a series of test cases of increasing complexity where traditional MD, and in several cases alternative REST2-like enhanced sampling methods, are shown to fail. The hydration free energy for acetic acid was shown to be independent of the starting conformation, and the values for four additional edge case molecules from the FreeSolv database were shown to have a significantly closer agreement with experiment using ACES. The method was further able to handle different rotamer states in a Cdk2 ligand identified as fractionally occupied in crystal structures. Finally, ACES was applied to T4-lysozyme and demonstrated that the side chain distribution of V111χ1 could be reliably reproduced for the apo state, bound to p-xylene, and in p-xylene→ benzene transformations. In these cases, the ACES method is shown to be highly robust and superior to a REST2-like enhanced sampling implementation alone.

Predictive Molecular Models for Charged Materials Systems: From Energy Materials to Biomacromolecules

Electrostatic interactions play a dominant role in charged materials systems. Understanding the complex correlation between macroscopic properties with microscopic structures is of critical importance to develop rational design strategies for advanced materials. But the complexity of this challenging task is augmented by interfaces present in the charged materials systems, such as electrode-electrolyte interfaces or biological membranes. Over the last decades, predictive molecular simulations that are founded in fundamental physics and optimized for charged interfacial systems have proven their value in providing molecular understanding of physicochemical properties and functional mechanisms for diverse materials. Novel design strategies utilizing predictive models have been suggested as promising route for the rational design of materials with tailored properties. Here, an overview of recent advances in the understanding of charged interfacial systems aided by predictive molecular simulations is presented. Focusing on three types of charged interfaces found in energy materials and biomacromolecules, how the molecular models characterize ion structure, charge transport, morphology relation to the environment, and the thermodynamics/kinetics of molecular binding at the interfaces is discussed. The critical analysis brings two prominent field of energy materials and biological science under common perspective, to stimulate crossover in both research field that have been largely separated.

A Coarse-Grained Molecular Dynamics Description of Docetaxel-Conjugate Release from PLGA Matrices

Despite the extensive use of poly-lactic-glycolic-acid (PLGA) in biomedical applications, computational research on the mesoscopic characterization of PLGA-based delivery systems is limited. In this study, a computational model for PLGA is proposed, developed, and validated for the reproducibility of transport properties that can influence drug release, the rate of which remains difficult to control. For computational efficiency, coarse-grained (CG) models of the molecular components under consideration were built using the MARTINI force field version 2.2. The translocation free energy barrier ΔG_t^* across the PLGA matrix in the aqueous phase of docetaxel and derivatives of varying sizes and solubilities was predicted via molecular dynamics (MD) simulations and compared with experimental release data. The thermodynamic quantity ΔG_t^* anticipates and can help explain the release kinetics of hydrophobic compounds from the PLGA matrix, albeit within the limit of a drug concentration below a critical aggregation concentration. The proposed computational framework would allow one to predict the pharmacological behavior of polymeric implants loaded with a variety of payloads under different conditions, limiting the experimental workload and associated costs.

Estimating the binding energetics of reversible covalent inhibitors of the SARS-CoV-2 main protease: an in silico study

The main protease (M^pro) of the SARS-CoV-2 virus is an attractive therapeutic target for developing antivirals to combat COVID-19. M^pro is essential for the replication cycle of the SARS-CoV-2 virus, so inhibiting M^pro blocks a vital piece of the cell replication machinery of the virus. A promising strategy to disrupt the viral replication cycle is to design inhibitors that bind to the active site cysteine (Cys145) of the M^pro. Cysteine targeted covalent inhibitors are gaining traction in drug discovery owing to the benefits of improved potency and extended drug-target engagement. An interesting aspect of these inhibitors is that they can be chemically tuned to form a covalent, but reversible bond, with their targets of interest. Several small-molecule cysteine-targeting covalent inhibitors of the M^pro have been discovered-some of which are currently undergoing evaluation in early phase human clinical trials. Understanding the binding energetics of these inhibitors could provide new insights to facilitate the design of potential drug candidates against COVID-19. Motivated by this, we employed rigorous absolute binding free energy calculations and hybrid quantum mechanical/molecular mechanical (QM/MM) calculations to estimate the energetics of binding of some promising reversible covalent inhibitors of the M^pro. We find that the inclusion of enhanced sampling techniques such as replica-exchange algorithm in binding free energy calculations can improve the convergence of predicted non-covalent binding free energy estimates of inhibitors binding to the M^pro target. In addition, our results indicate that binding free energy calculations coupled with multiscale simulations can be a useful approach to employ in ranking covalent inhibitors to their targets. This approach may be valuable in prioritizing and refining covalent inhibitor compounds for lead discovery efforts against COVID-19 and other coronavirus infections.

Fisetin 8-C-glucoside as entry inhibitor in SARS CoV-2 infection: molecular modelling study

Coronaviruses are RNA viruses that infect varied species including humans. TMPRSS2 is gateway for SARS CoV-2 entry into the host cell. It causes proteolytic activation of spike protein and discharge of the peptide into host cell. The TMPRSS2 inhibition could be one of the approaches to stop the viral entry, therefore, interaction pattern and binding energies for Fisetin and TMPRSS2 have been explored in the present study. TMPRSS2 peptide was used for homology modelling and then for further study. Molecular docking score and MMGBSA Binding energy of Fisetin was better than Nafamostat, a known inhibitor of TMPRSS2. Post docking MM-GBSA free energy for Fisetin and Nafamostat was -42.78 and -21.11 kcal/mol, respectively. Fisetin forms H bond with Val 25, His 41, Lys 42, Lys 45, Glu 44, Ser186. Nafamostat formed H bonds with Lys 85, Asp 90, Asp 203. RMSD plots of TMPRSS2, TMPRSS2-Fisetin and TMPRSS2-Nafamostat complex showed stable profile with very small fluctuation during entire simulation of 150 ns. Significant decrease in TMPRSS2-Fisetin and TMPRSS2-Nafamostat complex fluctuation occurred around His 41, Glu 44, Gly 136, Ser 186 in RMSF study. During simulation Fisetin interaction was observed with residues Val 25, His 41, Glu 44, Lys 45, Lys 87, Gly 136, Gln 183, Ser 186 likewise interaction of Nafamostat with Lys 85, Asp 90, Asn 163, Asp 203 and Ser 205. Post simulation MM-GBSA free energy was found to be -51.87 ± 4.3 and -48.23 ± 4.39 kcal/mol for TMPRSS2 with Fisetin and Nafamostat, respectively.Communicated by Ramaswamy H. Sarma.

Replica-Exchange Enveloping Distribution Sampling Using Generalized AMBER Force-Field Topologies: Application to Relative Hydration Free-Energy Calculations for Large Sets of Molecules

Free-energy differences between pairs of end-states can be estimated based on molecular dynamics (MD) simulations using standard pathway-dependent methods such as thermodynamic integration (TI), free-energy perturbation, or Bennett’s acceptance ratio. Replica-exchange enveloping distribution sampling (RE-EDS), on the other hand, allows for the sampling of multiple end-states in a single simulation without the specification of any pathways. In this work, we use the RE-EDS method as implemented in GROMOS together with generalized AMBER force-field (GAFF) topologies, converted to a GROMOS-compatible format with a newly developed GROMOS++ program amber2gromos, to compute relative hydration free energies for a series of benzene derivatives. The results obtained with RE-EDS are compared to the experimental data as well as calculated values from the literature. In addition, the estimated free-energy differences in water and in vacuum are compared to values from TI calculations carried out with GROMACS. The hydration free energies obtained using RE-EDS for multiple molecules are found to be in good agreement with both the experimental data and the results calculated using other free-energy methods. While all considered free-energy methods delivered accurate results, the RE-EDS calculations required the least amount of total simulation time. This work serves as a validation for the use of GAFF topologies with the GROMOS simulation package and the RE-EDS approach. Furthermore, the performance of RE-EDS for a large set of 28 end-states is assessed with promising results.

First Principles Calculation of Protein-Protein Dimer Affinities of ALS-Associated SOD1 Mutants

Cu,Zn superoxide dismutase (SOD1) is a 32 kDa homodimer that converts toxic oxygen radicals in neurons to less harmful species. The dimerization of SOD1 is essential to the stability of the protein. Monomerization increases the likelihood of SOD1 misfolding into conformations associated with aggregation, cellular toxicity, and neuronal death in familial amyotrophic lateral sclerosis (fALS). The ubiquity of disease-associated mutations throughout the primary sequence of SOD1 suggests an important role of physicochemical processes, including monomerization of SOD1, in the pathology of the disease. Herein, we use a first-principles statistical mechanics method to systematically calculate the free energy of dimer binding for SOD1 using molecular dynamics, which involves sequentially computing conformational, orientational, and separation distance contributions to the binding free energy. We consider the effects of two ALS-associated mutations in SOD1 protein on dimer stability, A4V and D101N, as well as the role of metal binding and disulfide bond formation. We find that the penalty for dimer formation arising from the conformational entropy of disordered loops in SOD1 is significantly larger than that for other protein–protein interactions previously considered. In the case of the disulfide-reduced protein, this leads to a bound complex whose formation is energetically disfavored. Somewhat surprisingly, the loop free energy penalty upon dimerization is still significant for the holoprotein, despite the increased structural order induced by the bound metal cations. This resulted in a surprisingly modest increase in dimer binding free energy of only about 1.5 kcal/mol upon metalation of the protein, suggesting that the most significant stabilizing effects of metalation are on folding stability rather than dimer binding stability. The mutant A4V has an unstable dimer due to weakened monomer-monomer interactions, which are manifested in the calculation by a separation free energy surface with a lower barrier. The mutant D101N has a stable dimer partially due to an unusually rigid β-barrel in the free monomer. D101N also exhibits anticooperativity in loop folding upon dimerization. These computational calculations are, to our knowledge, the most quantitatively accurate calculations of dimer binding stability in SOD1 to date.

PIGNet: a physics-informed deep learning model toward generalized drug-target interaction predictions

Recently, deep neural network (DNN)-based drug–target interaction (DTI) models were highlighted for their high accuracy with affordable computational costs. Yet, the models' insufficient generalization remains a challenging problem in the practice of in silico drug discovery. We propose two key strategies to enhance generalization in the DTI model. The first is to predict the atom–atom pairwise interactions via physics-informed equations parameterized with neural networks and provides the total binding affinity of a protein–ligand complex as their sum. We further improved the model generalization by augmenting a broader range of binding poses and ligands to training data. We validated our model, PIGNet, in the comparative assessment of scoring functions (CASF) 2016, demonstrating the outperforming docking and screening powers than previous methods. Our physics-informing strategy also enables the interpretation of predicted affinities by visualizing the contribution of ligand substructures, providing insights for further ligand optimization.

PIGNet, a deep neural network-based drug–target interaction model guided by physics and extensive data augmentation, shows significantly improved generalization ability and model performance.

Perturbation Free-Energy Toolkit: An Automated Alchemical Topology Builder

Free-energy calculations play an important role in the application of computational chemistry to a range of fields, including protein biochemistry, rational drug design, or materials science. Importantly, the free-energy difference is directly related to experimentally measurable quantities such as partition and adsorption coefficients, water activity, and binding affinities. Among several techniques aimed at predicting free-energy differences, perturbation approaches, involving the alchemical transformation of one molecule into another through intermediate states, stand out as rigorous methods based on statistical mechanics. However, despite the importance of free-energy calculations, the applicability of the perturbation approaches is still largely impeded by a number of challenges, including the definition of the perturbation path, i.e., alchemical changes leading to the transformation of one molecule to the other. To address this, an automatic perturbation topology builder based on a graph-matching algorithm is developed, which can identify the maximum common substructure (MCS) of two or multiple molecules and provide the perturbation topologies suitable for free-energy calculations using the GROMOS and the GROMACS simulation packages. Various MCS search options are presented leading to alternative definitions of the perturbation pathway. Moreover, perturbation topologies generated using the default multistate MCS search are used to calculate the changes in free energy between lysine and its two post-translational modifications, 3-methyllysine and acetyllysine. The pairwise free-energy calculations performed on this test system led to a cycle closure of 0.5 ± 0.3 and 0.2 ± 0.2 kJ mol-1, with GROMOS and GROMACS simulation packages, respectively. The same relative free energies between the three states are obtained by employing the enveloping distribution sampling (EDS) approach when compared to the pairwise perturbations. Importantly, this toolkit is made available online as an open-source Python package (https://github.com/drazen-petrov/SMArt).

Molecular Mechanism Exploration of Potent Fluorinated PI3K Inhibitors with a Triazine Scaffold: Unveiling the Unusual Synergistic Effect of Pyridine-to-Pyrimidine Ring Interconversion and CF3 Defluorination

The phosphatidylinostitol-3-kinase (PI3K)/AKT/mammalian target of rapamycin signaling pathway is a vital regulator of cell proliferation, growth, and survival, which is frequently overactivated in many human cancers. To this effect, PI3K, which is an important mediator of this pathway, has been pinpointed as a crucial target in cancer therapy and hence the importance of PI3K inhibitors. It was recently reported that defluorination and pyridine-to-pyrimidine ring interconversion increase the potency of specific small-molecule inhibitors of PI3K. Compound 4, an inhibitor with the difluorinated pyrimidine motif, was found to be eight times more potent against PI3K than compound 1, an inhibitor with the trifluorinated pyridine motif. This observation presents the need to rationally resolve the differential inhibitory mechanisms exhibited by both compounds. In this present work, we employed multiple computational approaches to investigate and distinguish the binding modes of 1 and 4 in addition to the effects they mediate on the secondary structure of PI3K. Likewise, we evaluated two other derivatives, compounds 2 with the difluorinated pyridine motif and 3 with the trifluorinated pyrimidine motif, to investigate the cooperativity effect between the defluorination of CF₃ and pyridine-to-pyrimidine ring interconversion. Findings revealed that PI3K, upon interaction with 4, exhibited a series of structural changes that favored the binding of the inhibitor at the active-site region. Furthermore, a positive (synergistic) cooperativity effect was observed between CF₃ defluorination and pyridine-to-pyrimidine ring interconversion. Moreover, there was a good correlation between the binding free energy estimated and the biological activity reported experimentally. Energy decomposition analysis revealed that the major contributing force to binding affinity variations between 1 and 4 is the electrostatic energy. Per-residue energy-based hierarchical clustering analysis further identified four hot-spot residues ASP841, TYR867, ASP964, and LYS833 and four warm-spot residues ASP836, SER806, ASP837, and LYS808, which essentially mediate the optimal and higher-affinity binding of compound 4 to PI3K relative to 1. This study therefore provides rational insights into the mechanisms by which 4 exhibited superior PI3K-inhibitory activities over 1, which is vital for future structure-based drug discovery efforts in PI3K targeting.

Predicting solvation free energies in non-polar solvents using classical density functional theory based on the PC-SAFT equation of state

We propose a predictive Density Functional Theory (DFT) for the calculation of solvation free energies. Our approach is based on a Helmholtz free-energy functional that is consistent with the Perturbed-Chain Statistical Associating Fluid Theory (PC-SAFT) equation of state. This allows for a coarse-grained description of the solvent based on an inhomogeneous density of PC-SAFT segments. The solute, on the other hand, is described in full detail by atomistic Lennard-Jones interaction sites. The approach is entirely predictive as it only takes the PC-SAFT parameters of the solvent and the force-field parameters of the solute as input. No adjustable parameters or empirical corrections are involved. The framework is applied to study self-solvation of n-alkanes and to the calculation of residual chemical potentials in binary solvent mixtures. Our DFT approach accurately predicts solvation free energies of small molecular solutes in three different non-polar solvents, namely n-hexane, cyclohexane, and benzene. Additionally, we show that the calculated solvation free energies agree well with those obtained by molecular dynamics simulations and with the residual chemical potential calculated by the bulk PC-SAFT equation of state. We observe higher deviations for the solvation free energy of systems with significant solute-solvent Coulomb interactions.

Structure-based in silico approaches for drug discovery against Mycobacterium tuberculosis

Mycobacterium tuberculosis is the causative agent of TB and was estimated to cause 1.4 million death in 2019, alongside 10 million new infections. Drug resistance is a growing issue, with multi-drug resistant infections representing 3.3% of all new infections, hence novel antimycobacterial drugs are urgently required to combat this growing health emergency. Alongside this, increased knowledge of gene essentiality in the pathogenic organism and larger compound databases can aid in the discovery of new drug compounds. The number of protein structures, X-ray based and modelled, is increasing and now accounts for greater than > 80% of all predicted M. tuberculosis proteins; allowing novel targets to be investigated. This review will focus on structure-based in silico approaches for drug discovery, covering a range of complexities and computational demands, with associated antimycobacterial examples. This includes molecular docking, molecular dynamic simulations, ensemble docking and free energy calculations. Applications of machine learning onto each of these approaches will be discussed. The need for experimental validation of computational hits is an essential component, which is unfortunately missing from many current studies. The future outlooks of these approaches will also be discussed.

Accurate absolute free energies for ligand-protein binding based on non-equilibrium approaches

The accurate calculation of the binding free energy for arbitrary ligand-protein pairs is a considerable challenge in computer-aided drug discovery. Recently, it has been demonstrated that current state-of-the-art molecular dynamics (MD) based methods are capable of making highly accurate predictions. Conventional MD-based approaches rely on the first principles of statistical mechanics and assume equilibrium sampling of the phase space. In the current work we demonstrate that accurate absolute binding free energies (ABFE) can also be obtained via theoretically rigorous non-equilibrium approaches. Our investigation of ligands binding to bromodomains and T4 lysozyme reveals that both equilibrium and non-equilibrium approaches converge to the same results. The non-equilibrium approach achieves the same level of accuracy and convergence as an equilibrium free energy perturbation (FEP) method enhanced by Hamiltonian replica exchange. We also compare uni- and bi-directional non-equilibrium approaches and demonstrate that considering the work distributions from both forward and reverse directions provides substantial accuracy gains. In summary, non-equilibrium ABFE calculations are shown to yield reliable and well-converged estimates of protein-ligand binding affinity.

Challenges Encountered Applying Equilibrium and Nonequilibrium Binding Free Energy Calculations

Binding free energy calculations have become increasingly valuable to drive decision making in drug discovery projects. However, among other issues, inadequate sampling can reduce accuracy, limiting the value of the technique. In this paper, we apply absolute binding free energy calculations to ligands binding to T4 lysozyme L99A and HSP90 using equilibrium and nonequilibrium approaches. We highlight sampling problems encountered in these systems, such as slow side chain rearrangements and slow changes of water placement upon ligand binding. These same types of challenges are also likely to show up in other protein-ligand systems, and we propose some strategies to diagnose and test for such problems in alchemical free energy calculations. We also explore similarities and differences in how the equilibrium and the nonequilibrium approaches handle these problems. Our results show the large amount of work still to be done to make free energy calculations robust and reliable and provide insight for future research in this area.

Estimating the Roles of Protonation and Electronic Polarization in Absolute Binding Affinity Simulations

Accurate prediction of binding free energies is critical to streamlining the drug development and protein design process. With the advent of GPU acceleration, absolute alchemical methods, which simulate the removal of ligand electrostatics and van der Waals interactions with the protein, have become routinely accessible and provide a physically rigorous approach that enables full consideration of flexibility and solvent interaction. However, standard explicit solvent simulations are unable to model protonation or electronic polarization changes upon ligand transfer from water to the protein interior, leading to inaccurate prediction of binding affinities for charged molecules. Here, we perform extensive simulation totaling ∼540 μs to benchmark the impact of modeling conditions on predictive accuracy for absolute alchemical simulations. Binding to urokinase plasminogen activator (UPA), a protein frequently overexpressed in metastatic tumors, is evaluated for a set of 10 inhibitors with extended flexibility, highly charged character, and titratable properties. We demonstrate that the alchemical simulations can be adapted to utilize the MBAR/PBSA method to improve the accuracy upon incorporating electronic polarization, highlighting the importance of polarization in alchemical simulations of binding affinities. Comparison of binding energy prediction at various protonation states indicates that proper electrostatic setup is also crucial in binding affinity prediction of charged systems, prompting us to propose an alternative binding mode with protonated ligand phenol and Hid-46 at the binding site, a testable hypothesis for future experimental validation.

Large-scale binding affinity calculations on commodity compute clouds

In recent years, it has become possible to calculate binding affinities of compounds bound to proteins via rapid, accurate, precise and reproducible free energy calculations. This is imperative in drug discovery as well as personalized medicine. This approach is based on molecular dynamics (MD) simulations and draws on sequence and structural information of the protein and compound concerned. Free energies are determined by ensemble averages of many MD replicas, each of which requires hundreds of cores and/or GPU accelerators, which are now available on commodity cloud computing platforms; there are also requirements for initial model building and subsequent data analysis stages. To automate the process, we have developed a workflow known as the binding affinity calculator. In this paper, we focus on the software infrastructure and interfaces that we have developed to automate the overall workflow and execute it on commodity cloud platforms, in order to reliably predict their binding affinities on time scales relevant to the domains of application, and illustrate its application to two free energy methods.

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.

Computing Absolute Free Energy with Deep Generative Models

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 intermediates. Here, 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.

Molecular dynamics simulation and free energy calculation studies of Coagulin L as dipeptidyl peptidase-4 inhibitor

Plant derived product can be used as other alternatives to currently used drugs for controlling chronic diseases like Diabetes mellitus. The potential of Coagulin L (a constituent of Withania coagulans) as dipeptidyl peptidase-4 (DPP-4) inhibitor was evaluated by molecular modelling study. It was observed that amino acid residues such as Glu205, Glu206, Tyr 547, His 740, and Try662 interacts with Coagulin L and Saxagliptin (a known DPP-4 inhibitor). Other nonbonded interactions of Coagulin L and Saxagliptin with DPP-4 binding residues were also found similar. The docking energy of Coagulin L was found to be -7.69 Kcal/mol whereas -8.44 kcal/mol was recorded for Saxagliptin. MD simulation study revealed stable binding throughout 100 ns simulation. RMSD plot of the complex was stabilized in 43 ns and remained stable during entire simulation(100 ns). RMSF plot of DPP-4 Coagulin L interaction showed major fluctuations at residue 246 and 766, however, Arg 125, Glu 205, Ser 209 and His 740 showed no major perturbations. Principal Component Analysis showed that important dynamics of the protein remain unchanged during entire simulation since the non-polar, van der waals, ionic interaction and solvation energy, altogether play important role in the complex stability. The molecular modelling study of DPP-4 with Coagulin L was an effort to establish correlation with traditional practices of Withania coagulans as antidiabetic agent in Indian subcontinent.Communicated by Ramaswamy H. Sarma.

Solubility Advantage of Amorphous Ketoprofen. Thermodynamic and Kinetic Aspects by Molecular Dynamics and Free Energy Approaches

Thermodynamic and kinetic aspects of crystalline (c-KTP) and amorphous (a-KTP) ketoprofen dissolution in water have been investigated by molecular dynamics simulation focusing on free energy properties. Absolute free energies of all relevant species and phases have been determined by thermodynamic integration on a novel path, first connecting the harmonic to the anharmonic system Hamiltonian at low T and then extending the result to the temperature of interest. The free energy required to transfer one ketoprofen molecule from the crystal to the solution is in fair agreement with the experimental value. The absolute free energy of the amorphous form is 19.58 kJ/mol higher than for the crystal, greatly enhancing the ketoprofen concentration in water, although as a metastable species in supersaturated solution. The kinetics of the dissolution process has been analyzed by computing the free energy profile along a reaction coordinate bringing one ketoprofen molecule from the crystal or amorphous phase to the solvated state. This computation confirms that, compared to the crystal form, the dissolution rate is nearly 7 orders of magnitude faster for the amorphous form, providing one further advantage to the latter in terms of bioavailability. The problem of drug solubility, of great practical importance, is used here as a test bed for a refined method to compute absolute free energies, which could be of great interest in biophysics and drug discovery in particular.

In silico design of peptides as potential ligands to resistin

Resistin is a hormone of biological interest due to its connection with several diseases of worldwide concern. This work aims to design a series of cyclic peptides as "lead compounds" to identify potential ligands to resistin. To this end, we propose an approach based on a peptide design algorithm plus a two-stage selection which accounts for selectivity, one of the most forgotten steps in the design of ligands. Following this approach, we have been able to identify several peptides as strong candidates for the design of elements of bio-recognition. Those peptides present low scoring binding energy to albumin, good water solubility, stability in water at 300 K, and high scoring binding energy to resistin. Among those peptides, two were chosen, to perform a more rigorous calculation of binding free energy based on the Alchemical Absolute Binding Free Energy method. We were able to establish a methodological route for the development of strong candidates for the design of ligands to resistin. Graphical Abstract Combined MD + MC + AABFE approach to design and screening of high-affinity binders to resistin.

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

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

Implementation of adaptive integration method for free energy calculations in molecular systems

Estimating free energy differences by computer simulation is useful for a wide variety of applications such as virtual screening for drug design and for understanding how amino acid mutations modify protein interactions. However, calculating free energy differences remains challenging and often requires extensive trial and error and very long simulation times in order to achieve converged results. Here, we present an implementation of the adaptive integration method (AIM). We tested our implementation on two molecular systems and compared results from AIM to those from a suite of other methods. The model systems tested here include calculating the solvation free energy of methane, and the free energy of mutating the peptide GAG to GVG. We show that AIM is more efficient than other tested methods for these systems, that is, AIM results converge to a higher level of accuracy and precision for a given simulation time.

Multiple Free Energy Calculations from Single State Point Continuous Fractional Component Monte Carlo Simulation Using Umbrella Sampling

We introduce an alternative method to perform free energy calculations for mixtures at multiple temperatures and pressures from a single simulation, by combining umbrella sampling and the continuous fractional component Monte Carlo method. One can perform a simulation of a mixture at a certain pressure and temperature and accurately compute the chemical potential at other pressures and temperatures close to the simulation conditions. This method has the following advantages: (1) Accurate estimates of the chemical potential as a function of pressure and temperature are obtained from a single state simulation without additional postprocessing. This can potentially reduce the number of simulations of a system for free energy calculations for a specific temperature and/or pressure range. (2) Partial molar volumes and enthalpies are obtained directly from the estimated chemical potentials. We tested our method for a Lennard-Jones system, aqueous mixtures of methanol at T = 298 K and P = 1 bar, and a mixture of ammonia, nitrogen, and hydrogen at T = 573 K and P = 800 bar. For pure methanol (N = 410 molecules), we observed that the estimated chemical potentials from umbrella sampling are in excellent agreement with the reference values obtained from independent simulations, for ΔT = ±15 K and ΔP = 100 bar (with respect to the simulated system). For larger systems, this range becomes smaller because the relative fluctuations of energy and volume become smaller. Without sufficient overlap, the performance of the method will become poor especially for nonlinear variations of the chemical potential.

Critical temperature, surface tension and vapor density estimations of methanol using the scaled model

In the present work, we applied scaled model to calculate surface tension, vapor densities and the critical temperatures of four different models of methanol: namely, H1, J1, J2 and L1 models. The scaled model is based on calculating the free energy of the system. Free energy calculations were performed by applying the Bennet acceptance ratio (BAR) using Monte-Carlo simulations at low temperature range of 220K-280K. The BAR is based on calculating the free energy difference of n-molecules and (n-1)-molecules plus a free probe on methanol. Estimations of vapor densities are based on extrapolating the intercept of the scaled free energy linear line as number of molecules approaches infinity, which requires a pre-known values for liquid densities. To accomplish this, a series of molecular dynamic simulations were performed at low temperature range of 200K-300K with steps of 10K. All the estimated properties were in excellent agreement with experimental published data.

Molecular simulation as a computational pharmaceutics tool to predict drug solubility, solubilization processes and partitioning

In this review we will discuss how computational methods, and in particular classical molecular dynamics simulations, can be used to calculate solubility of pharmaceutically relevant molecules and systems. To the extent possible, we focus on the non-technical details of these calculations, and try to show also the added value of a more thorough and detailed understanding of the solubilization process obtained by using computational simulations. Although the main focus is on classical molecular dynamics simulations, we also provide the reader with some insights into other computational techniques, such as the COSMO-method, and also discuss Flory-Huggins theory and solubility parameters. We hope that this review will serve as a valuable starting point for any pharmaceutical researcher, who has not yet fully explored the possibilities offered by computational approaches to solubility calculations.

Alchemical Free-Energy Calculations by Multiple-Replica λ-Dynamics: The Conveyor Belt Thermodynamic Integration Scheme

A new method is proposed to calculate alchemical free-energy differences based on molecular dynamics (MD) simulations, called the conveyor belt thermodynamic integration (CBTI) scheme. As in thermodynamic integration (TI), K replicas of the system are simulated at different values of the alchemical coupling parameter λ. The number K is taken to be even, and the replicas are equally spaced on a forward-turn-backward-turn path, akin to a conveyor belt (CB) between the two physical end-states; and as in λ-dynamics (λD), the λ-values associated with the individual systems evolve in time along the simulation. However, they do so in a concerted fashion, determined by the evolution of a single dynamical variable Λ of period 2π controlling the advance of the entire CB. Thus, a change of Λ is always associated with K/2 equispaced replicas moving forward and K/2 equispaced replicas moving backward along λ. As a result, the effective free-energy profile of the replica system along Λ is periodic of period 2 πK^-1, and the magnitude of its variations decreases rapidly upon increasing K, at least as K^-1 in the limit of large K. When a sufficient number of replicas is used, these variations become small, which enables a complete and quasi-homogeneous coverage of the λ-range by the replica system, without application of any biasing potential. If desired, a memory-based biasing potential can still be added to further homogenize the sampling, the preoptimization of which is computationally inexpensive. The final free-energy profile along λ is calculated similarly to TI, by binning of the Hamiltonian λ-derivative as a function of λ considering all replicas simultaneously, followed by quadrature integration. The associated quadrature error can be kept very low owing to the continuous and quasi-homogeneous λ-sampling. The CBTI scheme can be viewed as a continuous/deterministic/dynamical analog of the Hamiltonian replica-exchange/permutation (HRE/HRP) schemes or as a correlated multiple-replica analog of the λD or λ-local elevation umbrella sampling (λ-LEUS) schemes. Compared to TI, it shares the advantage of the latter schemes in terms of enhanced orthogonal sampling, i.e. the availability of variable-λ paths to circumvent conformational barriers present at specific λ-values. Compared to HRE/HRP, it permits a deterministic and continuous sampling of the λ-range, is expected to be less sensitive to possible artifacts of the thermo- and barostating schemes, and bypasses the need to carefully preselect a λ-ladder and a swapping-attempt frequency. Compared to λ-LEUS, it eliminates (or drastically reduces) the dead time associated with the preoptimization of a biasing potential. The goal of this article is to provide the mathematical/physical formulation of the proposed CBTI scheme, along with an initial application of the method to the calculation of the hydration free energy of methanol.

Challenges in the use of atomistic simulations to predict solubilities of drug-like molecules

Background: Solubility is a physical property of high importance to the pharmaceutical industry, the prediction of which for potential drugs has so far been a hard task. We attempted to predict the solubility of acetylsalicylic acid (ASA) by estimating the absolute chemical potentials of its most stable polymorph and of solutions with different concentrations of the drug molecule. Methods: Chemical potentials were estimated from all-atom molecular dynamics simulations.  We used the Einstein molecule method (EMM) to predict the absolute chemical potential of the solid and solvation free energy calculations to predict the excess chemical potentials of the liquid-phase systems. Results: Reliable estimations of the chemical potentials for the solid and for a single ASA molecule using the EMM required an extremely large number of intermediate states for the free energy calculations, meaning that the calculations were extremely demanding computationally. Despite the computational cost, however, the computed value did not agree well with the experimental value, potentially due to limitations with the underlying energy model. Perhaps better values could be obtained with a better energy model; however, it seems likely computational cost may remain a limiting factor for use of this particular approach to solubility estimation.    Conclusions: Solubility prediction of drug-like solids remains computationally challenging, and it appears that both the underlying energy model and the computational approach applied may need improvement before the approach is suitable for routine use.

Multiscale Modeling and Simulation Approaches to Lipid-Protein Interactions

Lipid membranes play a crucial role in living systems by compartmentalizing biological processes and forming a barrier between these processes and the environment. Naturally, a large apparatus of biomolecules is responsible for construction, maintenance, transport, and degradation of these lipid barriers. Additional classes of biomolecules are tasked with transport of specific substances or transduction of signals from the environment across lipid membranes. In this article, we intend to describe a set of techniques that enable one to build accurate models of lipid systems and their associated proteins, and to simulate their dynamics over a variety of time and length scales. We discuss the methods and challenges that allow us to derive structural, mechanistic, and thermodynamic information from these models. We also show how these models have recently been applied in research to study some of the most complex lipid-protein systems to date, including bacterial and viral envelopes, neuronal membranes, and mammalian signaling systems.

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.

Assimilating Radial Distribution Functions To Build Water Models with Improved Structural Properties

The structural properties of three- and four-site water models are improved by extending the ForceBalance parametrization code to include a new methodology allowing for the targeting of any radial distribution function (RDF) during the parametrization of a force field. The mean squared difference (MSD) between the experimental and simulated RDFs contributes to an objective function, allowing for the systematic optimization of force field parameters to reach closer overall agreement with experiment. RDF fitting is applied to develop modified versions of the TIP3P and TIP4P/2005 water models in which the Lennard-Jones potential is replaced by a Buckingham potential. The optimized TIP3P-Buckingham and TIP4P-Buckingham potentials feature 93 and 98% lower MSDs in the OO RDF compared to the TIP3P and TIP4P/2005 models respectively, with marked decreases in the height of the first peak. Additionally, these Buckingham models predict the entropy of water more accurately, reducing the error in the entropy of TIP3P from 11 to 3% and the error in the entropy of TIP4P/2005 from 11 to 2%. These new Buckingham models have improved predictive power for many nonfitted properties particularly in the case of TIP3P. Our work directly demonstrates how the Buckingham potential can improve the description of water's structural properties beyond the Lennard-Jones potential. Moreover, adding a Buckingham potential is a favorable alternative to adding interaction sites in terms of computational speed on modern GPU hardware.

Estimating Atomic Contributions to Hydration and Binding Using Free Energy Perturbation

We present a general method called atom-wise free energy perturbation (AFEP), which extends a conventional molecular dynamics free energy perturbation (FEP) simulation to give the contribution to a free energy change from each atom. AFEP is derived from an expansion of the Zwanzig equation used in the exponential averaging method by defining that the system total energy can be partitioned into contributions from each atom. A partitioning method is assumed and used to group terms in the expansion to correspond to individual atoms. AFEP is applied to six example free energy changes to demonstrate the method. Firstly, the hydration free energies of methane, methanol, methylamine, methanethiol, and caffeine in water. AFEP highlights the atoms in the molecules that interact favorably or unfavorably with water. Finally AFEP is applied to the binding free energy of human immunodeficiency virus type 1 protease to lopinavir, and AFEP reveals the contribution of each atom to the binding free energy, indicating candidate areas of the molecule to improve to produce a more strongly binding inhibitor. FEP gives a single value for the free energy change and is already a very useful method. AFEP gives a free energy change for each "part" of the system being simulated, where part can mean individual atoms, chemical groups, amino acids, or larger partitions depending on what the user is trying to measure. This method should have various applications in molecular dynamics studies of physical, chemical, or biochemical phenomena, specifically in the field of computational drug discovery.

The Role of the Active Site Flap in Streptavidin/Biotin Complex Formation

Obtaining a detailed description of how active site flap motion affects substrate or ligand binding will advance structure-based drug design (SBDD) efforts on systems including the kinases, HSP90, HIV protease, ureases, etc. Through this understanding, we will be able to design better inhibitors and better proteins that have desired functions. Herein we address this issue by generating the relevant configurational states of a protein flap on the molecular energy landscape using an approach we call MT_Flex-b and then following this with a procedure to estimate the free energy associated with the motion of the flap region. To illustrate our overall workflow, we explored the free energy changes in the streptavidin/biotin system upon introducing conformational flexibility in loop_3-4 in the biotin unbound ( apo) and bound ( holo) state. The free energy surfaces were created using the Movable Type free energy method, and for further validation, we compared them to potential of mean force (PMF) generated free energy surfaces using MD simulations employing the FF99SBILDN and FF14SB force fields. We also estimated the free energy thermodynamic cycle using an ensemble of closed-like and open-like end states for the ligand unbound and bound states and estimated the binding free energy to be approximately -16.2 kcal/mol (experimental -18.3 kcal/mol). The good agreement between MT_Flex-b in combination with the MT method with experiment and MD simulations supports the effectiveness of our strategy in obtaining unique insights into the motions in proteins that can then be used in a range of biological and biomedical applications.