Basic ingredients of free energy calculations: A review

Convergence criteria for single-step free-energy calculations: the relation between the Π bias measure and the sample variance

Free energy calculations play a crucial role in simulating chemical processes, enzymatic reactions, and drug design. However, assessing the reliability and convergence of these calculations remains a challenge. This study focuses on single-step free-energy calculations using thermodynamic perturbation. It explores how the sample distributions influence the estimated results and evaluates the reliability of various convergence criteria, including Kofke's bias measure Π and the standard deviation of the energy difference ΔU, σ ΔU . The findings reveal that for Gaussian distributions, there is a straightforward relationship between Π and σ ΔU , free energies can be accurately approximated using a second-order cumulant expansion, and reliable results are attainable for σ ΔU up to 25 kcal mol-1. However, interpreting non-Gaussian distributions is more complex. If the distribution is skewed towards more positive values than a Gaussian, converging the free energy becomes easier, rendering standard convergence criteria overly stringent. Conversely, distributions that are skewed towards more negative values than a Gaussian present greater challenges in achieving convergence, making standard criteria unreliable. We propose a practical approach to assess the convergence of estimated free energies.

OGRe: Optimal Grid Refinement Protocol for Accurate Free Energy Surfaces and Its Application in Proton Hopping in Zeolites and 2D COF Stacking

While free energy surfaces are the crux of our understanding of many chemical and biological processes, their accuracy is generally unknown. Moreover, many developments to improve their accuracy are often complicated, limiting their general use. Luckily, several tools and guidelines are already in place to identify these shortcomings, but they are typically lacking in flexibility or fail to systematically determine how to improve the accuracy of the free energy calculation. To overcome these limitations, this work introduces OGRe, a Python package for optimal grid refinement in an arbitrary number of dimensions. OGRe is based on three metrics that gauge the confinement, consistency, and overlap of each simulation in a series of umbrella sampling (US) simulations, an enhanced sampling technique ubiquitously adopted to construct free energy surfaces for hindered processes. As these three metrics are fundamentally linked to the accuracy of the weighted histogram analysis method adopted to generate free energy surfaces from US simulations, they facilitate the systematic construction of accurate free energy profiles, where each metric is driven by a specific umbrella parameter. This allows for the derivation of a consistent and optimal collection of umbrellas for each simulation, largely independent of the initial values, thereby dramatically increasing the ease-of-use toward accurate free energy surfaces. As such, OGRe is particularly suited to determine complex free energy surfaces with large activation barriers and shallow minima, which underpin many physical and chemical transformations and hence to further our fundamental understanding of these processes.

An interoperable implementation of collective-variable based enhanced sampling methods in extended phase space within the OpenMM package

Collective variable (CV)-based enhanced sampling techniques are widely used today for accelerating barrier-crossing events in molecular simulations. A class of these methods, which includes temperature accelerated molecular dynamics (TAMD)/driven-adiabatic free energy dynamics (d-AFED), unified free energy dynamics (UFED), and temperature accelerated sliced sampling (TASS), uses an extended variable formalism to achieve quick exploration of conformational space. These techniques are powerful, as they enhance the sampling of a large number of CVs simultaneously compared to other techniques. Extended variables are kept at a much higher temperature than the physical temperature by ensuring adiabatic separation between the extended and physical subsystems and employing rigorous thermostatting. In this work, we present a computational platform to perform extended phase space enhanced sampling simulations using the open-source molecular dynamics engine OpenMM. The implementation allows users to have interoperability of sampling techniques, as well as employ state-of-the-art thermostats and multiple time-stepping. This work also presents protocols for determining the critical parameters and procedures for reconstructing high-dimensional free energy surfaces. As a demonstration, we present simulation results on the high dimensional conformational landscapes of the alanine tripeptide in vacuo, tetra-N-methylglycine (tetra-sarcosine) peptoid in implicit solvent, and the Trp-cage mini protein in explicit water.

Advances in Computational Approaches for Estimating Passive Permeability in Drug Discovery

Passive permeation of cellular membranes is a key feature of many therapeutics. The relevance of passive permeability spans all biological systems as they all employ biomembranes for compartmentalization. A variety of computational techniques are currently utilized and under active development to facilitate the characterization of passive permeability. These methods include lipophilicity relations, molecular dynamics simulations, and machine learning, which vary in accuracy, complexity, and computational cost. This review briefly introduces the underlying theories, such as the prominent inhomogeneous solubility diffusion model, and covers a number of recent applications. Various machine-learning applications, which have demonstrated good potential for high-volume, data-driven permeability predictions, are also discussed. Due to the confluence of novel computational methods and next-generation exascale computers, we anticipate an exciting future for computationally driven permeability predictions.

Temperature Accelerated Sliced Sampling to Probe Ligand Dissociation from Protein

Modeling ligand unbinding in proteins to estimate the free energy of binding and probing the mechanism presents several challenges. They primarily pertain to the entropic bottlenecks resulting from protein and solvent conformations. While exploring the unbinding processes using enhanced sampling techniques, very long simulations are required to sample all of the conformational states as the system gets trapped in local free energy minima along transverse coordinates. Here, we demonstrate that temperature accelerated sliced sampling (TASS) is an ideal approach to overcome some of the difficulties faced by conventional sampling methods in studying ligand unbinding. Using TASS, we study the unbinding of avibactam inhibitor molecules from the Class C β-lactamase (CBL) active site. Extracting CBL-avibactam unbinding free energetics, unbinding pathways, and identifying critical interactions from the TASS simulations are demonstrated.

Revisiting MMPBSA by Adoption of MC-Based Surface Area/Volume, ANI-ML Potentials, and Two-Valued Interior Dielectric Constant

Here, we report the accuracy improvements of molecular mechanics Poisson-Boltzmann surface area (MMPBSA) calculations by adoption of ANI-ML potentials in replacement of MM terms, the use of solvent-accessible surface area (SASA) and volume (SAV) values from the Monte Carlo sampling of the probe, and introducing two different interior dielectric constants for electrostatic interactions of protein-ligand (P-L) and polar solvation term in the MMPBSA calculations. Our results show that the Pearson correlation coefficients of MMPBSA-calculated values with respect to experimental binding free energies can be drastically improved from 0.48 to 0.90 by adoption of ANI-ML potentials in replacement of MM energy terms in the equation, referred to as ANI-PBSA. Moreover, we show that the SASA/SAV-combined equation in the scaled particle theory (SPT) can be a better choice to model nonpolar solvation term, reaching nearly the same accuracy by ANI-PBSA calculations. Finally, we introduce two different values of interior dielectric constants, which could be an alternative strategy between the single and variable constant definitions.

Quantum Mechanical-Cluster Approach to Solve the Bioisosteric Replacement Problem in Drug Design

Bioisosteres are molecules that differ in substituents but still have very similar shapes. Bioisosteric replacements are ubiquitous in modern drug design, where they are used to alter metabolism, change bioavailability, or modify activity of the lead compound. Prediction of relative affinities of bioisosteres with computational methods is a long-standing task; however, the very shape closeness makes bioisosteric substitutions almost intractable for computational methods, which use standard force fields. Here, we design a quantum mechanical (QM)-cluster approach based on the GFN2-xTB semi-empirical quantum-chemical method and apply it to a set of H → F bioisosteric replacements. The proposed methodology enables advanced prediction of biological activity change upon bioisosteric substitution of -H with -F, with the standard deviation of 0.60 kcal/mol, surpassing the ChemPLP scoring function (0.83 kcal/mol), and making QM-based ΔΔG estimation comparable to ∼0.42 kcal/mol standard deviation of in vitro experiment. The speed of the method and lack of tunable parameters makes it affordable in current drug research.

Structural elucidation of ETHR-A and ETHR-B from Plutella xylostella and insight into non-conservative mutations in transmembrane helix-6

The development of Diamondback moth (DBM) depends on the ecdysis triggering hormone receptor (ETHR); a neuronal membrane G-protein coupled receptor (GPCR) connected to the metamorphosis cascade. Lepidopteran insect DBM is an infamous pest of cruciferous plants. This study examined the full-length coding sequences (CDS) of PxETHR-A and PxETHR-B from the DBM genome. The three-dimensional (3 D) models of both receptors were generated in an inactive state. The behaviour and stability of receptors were examined using molecular dynamics simulations in a lipid membrane system for 300 ns and established a GPCR family-based view. Secondary interactions within receptors were studied to know more about factors contributing to their stability. Multiple sequence alignment revealed conserved features of insect ETHRs those compared to the GPCR family proteins. These features were helpful during the evaluation of the molecular models of both receptors. Side-chain orientation of conserved residues, non-conserved and conserved hydrogen-bond networks (HBN) and hydrophobic clusters were examined in the structures of both receptors. The non-conserved residues L6.35, T6.39, C/S6.43, and L6.48, are present in a conserved position on the transmembrane helix-6 (TM6) of ETHRs. In TM6, PxETHR-A and PxETHR-B differ at positions C/S6.43 and Y/F6.51, both being part of the HBN.Communicated by Ramaswamy H. Sarma.

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.

Accelerated Enveloping Distribution Sampling (AEDS) Allows for Efficient Sampling of Orthogonal Degrees of Freedom

One of the most challenging aspects in the molecular simulation of proteins is the study of slowly relaxing processes, such as loop rearrangements or ligands that adopt different conformations in the binding site. State-of-the-art methods used to calculate binding free energies rely on performing several short simulations (lambda steps), in which the ligand is slowly transformed into the endstates of interest. This makes capturing the slowly relaxing processes even more difficult, as they would need to be observed in all of the lambda steps. One attractive alternative is the use of a reference state capable of sampling all of the endstates of interest in a single simulation. However, the energy barriers between the states can be too high to overcome, thus hindering the sampling of all of the relevant conformations. Accelerated enveloping distribution sampling (AEDS) is a recently developed reference state technique that circumvents the high-energy-barrier challenge by adding a boosting potential that flattens the energy landscape without distorting the energy minima. In the present work, we apply AEDS to the well-studied benchmark system T4 lysozyme L99A. The T4 lysozyme L99A mutant contains a hydrophobic pocket in which there is a valine (valine 111), whose conformation influences the binding efficiencies of the different ligands. Incorrectly sampling the dihedral angle can lead to errors in predicted binding free energies of up to 16 kJ mol^-1. This protein constitutes an ideal scenario to showcase the advantages and challenges when using AEDS in the presence of slow relaxing processes. We show that AEDS is able to successfully sample the relevant degrees of freedom, providing accurate binding free energies, without the need of previous information of the system in the form of collective variables. Additionally, we showcase the capabilities of AEDS to efficiently screen a set of ligands. These results represent a promising first step toward the development of free-energy methods that can respond to more intricate molecular events.

Multireference Generalization of the Weighted Thermodynamic Perturbation Method

We describe the generalized weighted thermodynamic perturbation (gwTP) method for estimating the free energy surface of an expensive "high-level" potential energy function from the umbrella sampling performed with multiple inexpensive "low-level" reference potentials. The gwTP method is a generalization of the weighted thermodynamic perturbation (wTP) method developed by Li and co-workers [J. Chem. Theory Comput. 2018, 14, 5583-5596] that uses a single "low-level" reference potential. The gwTP method offers new possibilities in model design whereby the sampling generated from several low-level potentials may be combined (e.g., specific reaction parameter models that might have variable accuracy at different stages of a multistep reaction). The gwTP method is especially well suited for use with machine learning potentials (MLPs) that are trained against computationally expensive ab initio quantum mechanical/molecular mechanical (QM/MM) energies and forces using active learning procedures that naturally produce multiple distinct neural network potentials. Simulations can be performed with greater sampling using the fast MLPs and then corrected to the ab initio level using gwTP. The capabilities of the gwTP method are demonstrated by creating reference potentials based on the MNDO/d and DFTB2/MIO semiempirical models supplemented with the "range-corrected deep potential" (DPRc). The DPRc parameters are trained to ab initio QM/MM data, and the potentials are used to calculate the free energy surface of stepwise mechanisms for nonenzymatic RNA 2'-O-transesterification model reactions. The extended sampling made possible by the reference potentials allows one to identify unequilibrated portions of the simulations that are not always evident from the short time scale commonly used with ab initio QM/MM potentials. We show that the reference potential approach can yield more accurate ab initio free energy predictions than the wTP method or what can be reasonably afforded from explicit ab initio QM/MM sampling.

Assessing the effect of forcefield parameter sets on the accuracy of relative binding free energy calculations

Software for accurate prediction of protein-ligand binding affinity can be a key enabling tool for small molecule drug discovery. Free energy perturbation (FEP) is a computational technique that can be used to compute binding affinity differences between molecules in a congeneric series. It has shown promise in reliably generating accurate predictions and is now widely used in the pharmaceutical industry. However, the high computational cost and use of commercial software, together with the technical challenges to setup, run, and analyze the simulations, limits the usage of FEP. Here, we use an automated FEP workflow which uses the open-source OpenMM package. To enable effective application of FEP, we compared the performance of different water models, partial charge assignments, and AMBER protein forcefields in eight benchmark test cases previously assembled for FEP validation studies.

From free-energy profiles to activation free energies

Given a chemical reaction going from reactant (R) to the product (P) on a potential energy surface (PES) and a collective variable (CV) discriminating between R and P, we define the free-energy profile (FEP) as the logarithm of the marginal Boltzmann distribution of the CV. This FEP is not a true free energy. Nevertheless, it is common to treat the FEP as the “free-energy” analog of the minimum potential energy path and to take the activation free energy, [Formula: see text], as the difference between the maximum at the transition state and the minimum at R. We show that this approximation can result in large errors. The FEP depends on the CV and is, therefore, not unique. For the same reaction, different discriminating CVs can yield different [Formula: see text]. We derive an exact expression for the activation free energy that avoids this ambiguity. We find [Formula: see text] to be a combination of the probability of the system being in the reactant state, the probability density on the dividing surface, and the thermal de Broglie wavelength associated with the transition. We apply our formalism to simple analytic models and realistic chemical systems and show that the FEP-based approximation applies only at low temperatures for CVs with a small effective mass. Most chemical reactions occur on complex, high-dimensional PES that cannot be treated analytically and pose the added challenge of choosing a good CV. We study the influence of that choice and find that, while the reaction free energy is largely unaffected, [Formula: see text] is quite sensitive.

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.

Discovery of Histone Deacetylase Inhibitor Using Molecular Modeling and Free Energy Calculations

The histone acetylation-deacetylation at lysine regulates the functions of many cellular proteins. An increased expression of HDAC6 can cause an increased amount of deacetylated histones, which leads to an inhibition of gene expression and has been associated with cancer cell proliferation. The present study screened the ZINC database to find novel HDAC6 inhibitors using virtual high-throughput screening techniques. The docking score, free energy, and binding pattern of the complexes were used to select a best ligand for further study. Molecular dynamic simulations, binding interactions, and the stability of docked conformations were investigated. Several parameters that determine protein-ligand interactions, such as root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), radius of gyration (Rg), and binding pattern, were observed. Hydrogen bonds were observed at His 573 and Gly 582 after a 150 ns simulation with identified compound ZINC000002845205, and they were similar to known inhibitor Panobinostat. The molecular mechanics with generalised Born and surface area solvation (MM/GBSA) free energy was comparable to known inhibitor Panobinostat. ZINC000002845205 qualifies drug-likeness according to Lipinski's rule-of-five, rule-of-three, and the World Drug Index (WDI)-like rule, but there is one violation in the lead-like rule.

Exploration of high dimensional free energy landscapes by a combination of temperature-accelerated sliced sampling and parallel biasing

Temperature-accelerated sliced sampling (TASS) is an enhanced sampling method for achieving accelerated and controlled exploration of high-dimensional free energy landscapes in molecular dynamics simulations. With the aid of umbrella bias potentials, the TASS method realizes a controlled exploration and divide-and-conquer strategy for computing high-dimensional free energy surfaces. In TASS, diffusion of the system in the collective variable (CV) space is enhanced with the help of metadynamics bias and elevated-temperature of the auxiliary degrees of freedom (DOF) that are coupled to the CVs. Usually, a low-dimensional metadynamics bias is applied in TASS. In order to further improve the performance of TASS, we propose here to use a highdimensional metadynamics bias, in the same form as in a parallel bias metadynamics scheme. Here, a modified reweighting scheme, in combination with artificial neural network is used for computing unbiased probability distribution of CVs and projections of high-dimensional free energy surfaces. We first validate the accuracy and efficiency of our method in computing the four-dimensional free energy landscape for alanine tripeptide in vacuo. Subsequently, we employ the approach to calculate the eight-dimensional free energy landscape of alanine pentapeptide in vacuo. Finally, the method is applied to a more realistic problem wherein we compute the broad four-dimensional free energy surface corresponding to the deacylation of a drug molecule which is covalently complexed with a β-lactamase enzyme. We demonstrate that using parallel bias in TASS improves the efficiency of exploration of high-dimensional free energy landscapes.

Multiple Molecular Dynamics Simulations and Energy Analysis Unravel the Dynamic Properties and Binding Mechanism of Mutants HIV-1 Protease with DRV and CA-p2

PRS17, a variant of human immunodeficiency virus type I protease (HIV-1 PR), has 17 mutated residues showing high levels of multidrug resistance. To describe the effects of these mutated residues on the dynamic properties and the binding mechanism of PR with substrate and inhibitor, focused on six systems (two complexes of WT PR and PRS17 with inhibitor Darunavir (DRV), two complexes of WT PR and PRS17 with substrate analogue CA-p2, two unligand WT PR and PRS17), we performed multiple molecular dynamics (MD) simulations combined with MM-PBSA and solvated interaction energy (SIE) methods. For both the unligand PRs and ligand-PR complexes, the results from simulations revealed 17 mutated residues alter the flap-flap distance, the distance from flap regions to catalytic sites, and the curling degree of the flap tips. These mutated residues changed the flexibility of the flap region in PR, and thus affected its binding energy with DRV and CA-p2, resulting in differences in sensitivity. Hydrophobic cavity makes an important contribution to the binding of PR and ligands. And most noticeable of all, the binding of the guanidine group in CA-p2 and Arg8' of PRS17 is useful for increasing their binding ability. These results have important guidance for the further design of drugs against multidrug resistant PR. IMPORTANCE Developing effective anti-HIV inhibitors is the current requirement to cope with the emergence of the resistance of mutants. Compared with the experiments, MD simulations along with energy calculations help reduce the time and cost of designing new inhibitors. Based on our simulation results, we propose two factors that may help design effective inhibitors against HIV-1 PR: (i) importance of hydrophobic cavity, and (ii) introduction of polar groups similar to the guanidine group.

Optimizing the Calculation of Free Energy Differences in Nonequilibrium Work SQM/MM Switching Simulations

A key step during indirect alchemical free energy simulations using quantum mechanical/molecular mechanical (QM/MM) hybrid potential energy functions is the calculation of the free energy difference ΔAlow→high between the low level (e.g., pure MM) and the high level of theory (QM/MM). A reliable approach uses nonequilibrium work (NEW) switching simulations in combination with Jarzynski's equation; however, it is computationally expensive. In this study, we investigate whether it is more efficient to use more shorter switches or fewer but longer switches. We compare results obtained with various protocols to reference free energy differences calculated with Crooks' equation. The central finding is that fewer longer switches give better converged results. As few as 200 sufficiently long switches lead to ΔAlow→high values in good agreement with the reference results. This optimized protocol reduces the computational cost by a factor of 40 compared to earlier work. We also describe two tools/ways of analyzing the raw data to detect sources of poor convergence. Specifically, we find it helpful to analyze the raw data (work values from the NEW switching simulations) in a quasi-time series-like manner. Principal component analysis helps to detect cases where one or more conformational degrees of freedom are different at the low and high level of theory.

Multiscale Modeling of Bio-Nano Interactions of Zero-Valent Silver Nanoparticles

Understanding the specifics of interaction between the protein and nanomaterial is crucial for designing efficient, safe, and selective nanoplatforms, such as biosensor or nanocarrier systems. Routing experimental screening for the most suitable complementary pair of biomolecule and nanomaterial used in such nanoplatforms might be a resource-intensive task. While a range of computational tools are available for prescreening libraries of proteins for their interactions with small molecular ligands, choices for high-throughput screening of protein libraries for binding affinities to new and existing nanomaterials are very limited. In the current work, we present the results of the systematic computational study of interaction of various biomolecules with pristine zero-valent noble metal nanoparticles, namely, AgNPs, by using the UnitedAtom multiscale approach. A set of blood plasma and dietary proteins for which the interaction with AgNPs was described experimentally were examined computationally to evaluate the performance of the UnitedAtom method. A set of interfacial descriptors (log P^NM, adsorption affinities, and adsorption affinity ranking), which can characterize the relative hydrophobicity/hydrophilicity/lipophilicity of the nanosized silver and its ability to form bio(eco)corona, was evaluated for future use in nano-QSAR/QSPR studies.

Boosting the conformational sampling by combining replica exchange with solute tempering and well-sliced metadynamics

Methods that combine collective variable (CV) based enhanced sampling and global tempering approaches are used in speeding-up the conformational sampling and free energy calculation of large and soft systems with a plethora of energy minima. In this paper, a new method of this kind is proposed in which the well-sliced metadynamics approach (WSMTD) is united with replica exchange with solute tempering (REST2) method. WSMTD employs a divide-and-conquer strategy wherein high-dimensional slices of a free energy surface are independently sampled and combined. The method enables one to accomplish a controlled exploration of the CV-space with a restraining bias as in umbrella sampling, and enhance-sampling of one or more orthogonal CVs using a metadynamics like bias. The new hybrid method proposed here enables boosting the sampling of more slow degrees of freedom in WSMTD simulations, without the need to specify associated CVs, through a replica exchange scheme within the framework of REST2. The high-dimensional slices of the probability distributions of CVs computed from the united WSMTD and REST2 simulations are subsequently combined using the weighted histogram analysis method to obtain the free energy surface. We show that the new method proposed here is accurate, improves the conformational sampling, and achieves quick convergence in free energy estimates. We demonstrate this by computing the conformational free energy landscapes of solvated alanine tripeptide and Trp-cage mini protein in explicit water.

The automated optimisation of a coarse-grained force field using free energy data

Atomistic models provide a detailed representation of molecular systems, but are sometimes inadequate for simulations of large systems over long timescales. Coarse-grained models enable accelerated simulations by reducing the number of degrees of freedom, at the cost of reduced accuracy. New optimisation processes to parameterise these models could improve their quality and range of applicability. We present an automated approach for the optimisation of coarse-grained force fields, by reproducing free energy data derived from atomistic molecular simulations. To illustrate the approach, we implemented hydration free energy gradients as a new target for force field optimisation in ForceBalance and applied it successfully to optimise the un-charged side-chains and the protein backbone in the SIRAH protein coarse-grain force field. The optimised parameters closely reproduced hydration free energies of atomistic models and gave improved agreement with experiment.

Machine Learning Applied to the Modeling of Pharmacological and ADMET Endpoints

The well-known concept of quantitative structure-activity relationships (QSAR) has been gaining significant interest in the recent years. Data, descriptors, and algorithms are the main pillars to build useful models that support more efficient drug discovery processes with in silico methods. Significant advances in all three areas are the reason for the regained interest in these models. In this book chapter we review various machine learning (ML) approaches that make use of measured in vitro/in vivo data of many compounds. We put these in context with other digital drug discovery methods and present some application examples.

New perspectives in cancer drug development: computational advances with an eye to design

Computational chemistry has come of age in drug discovery. Indeed, most pharmaceutical development programs rely on computer-based data and results at some point. Herein, we discuss recent applications of advanced simulation techniques to difficult challenges in drug discovery. These entail the characterization of allosteric mechanisms and the identification of allosteric sites or cryptic pockets determined by protein motions, which are not immediately evident in the experimental structure of the target; the study of ligand binding mechanisms and their kinetic profiles; and the evaluation of drug-target affinities. We analyze different approaches to tackle challenging and emerging biological targets. Finally, we discuss the possible perspectives of future application of computation in drug discovery.

Trends in application of advancing computational approaches in GPCR ligand discovery

G protein-coupled receptors (GPCRs) comprise the most important superfamily of protein targets in current ligand discovery and drug development. GPCRs are integral membrane proteins that play key roles in various cellular signaling processes. Therefore, GPCR signaling pathways are closely associated with numerous diseases, including cancer and several neurological, immunological, and hematological disorders. Computer-aided drug design (CADD) can expedite the process of GPCR drug discovery and potentially reduce the actual cost of research and development. Increasing knowledge of biological structures, as well as improvements on computer power and algorithms, have led to unprecedented use of CADD for the discovery of novel GPCR modulators. Similarly, machine learning approaches are now widely applied in various fields of drug target research. This review briefly summarizes the application of rising CADD methodologies, as well as novel machine learning techniques, in GPCR structural studies and bioligand discovery in the past few years. Recent novel computational strategies and feasible workflows are updated, and representative cases addressing challenging issues on olfactory receptors, biased agonism, and drug-induced cardiotoxic effects are highlighted to provide insights into future GPCR drug discovery.

Machine learning-accelerated quantum mechanics-based atomistic simulations for industrial applications

Atomistic simulations have become an invaluable tool for industrial applications ranging from the optimization of protein-ligand interactions for drug discovery to the design of new materials for energy applications. Here we review recent advances in the use of machine learning (ML) methods for accelerated simulations based on a quantum mechanical (QM) description of the system. We show how recent progress in ML methods has dramatically extended the applicability range of conventional QM-based simulations, allowing to calculate industrially relevant properties with enhanced accuracy, at reduced computational cost, and for length and time scales that would have otherwise not been accessible. We illustrate the benefits of ML-accelerated atomistic simulations for industrial R&D processes by showcasing relevant applications from two very different areas, drug discovery (pharmaceuticals) and energy materials. Writing from the perspective of both a molecular and a materials modeling scientist, this review aims to provide a unified picture of the impact of ML-accelerated atomistic simulations on the pharmaceutical, chemical, and materials industries and gives an outlook on the exciting opportunities that could emerge in the future.

Recent progress in simulating microscopic ion transport mechanisms at liquid-liquid interfaces

Transport of ions through liquid-liquid interfaces is of fundamental importance to a wide variety of applications. However, since it is quite challenging for experimentalists to directly and selectively observe molecules at the interfaces, microscopic mechanisms of ion transport have been largely presumed from kinetic information. This Perspective illustrates recent examples that molecular dynamics simulations with proper free energy surfaces clarified mechanistic pictures of ion transport. The key is a proper choice of coordinates and defining/calculating free energy surfaces in multidimensional space. Once the free energy surfaces for realistic systems are available, they naturally provide new insight into the ion transport in unprecedented details, including water finger, transient ion pairing, and electron transfer.

Ensembler: A Simple Package for Fast Prototyping and Teaching Molecular Simulations

Ensembler is a Python package that enables method prototyping using 1D and 2D model systems and allows deepening of the understanding of different molecular dynamics (MD) methods, starting from basic techniques to enhanced sampling and free-energy approaches. The ease of installing and using the package increases shareability, comparability, and reproducibility of scientific code developments. Here, we describe the implementation and usage of the package and provide an application example for free-energy calculation. The code of Ensembler is freely available on GitHub at https://github.com/rinikerlab/Ensembler.

Ligand Affinities within the Open-Boundary Molecular Mechanics/Coarse-Grained Framework (I): Alchemical Transformations within the Hamiltonian Adaptive Resolution Scheme

Our recently developed Open-Boundary Molecular Mechanics/Coarse Grained (OB-MM/CG) framework predicts ligand poses in important pharmaceutical targets, such as G-protein Coupled Receptors, even when experimental structural information is lacking. The approach, which is based on GROMOS and AMBER force fields, allows for grand-canonical simulations of protein-ligand complexes by using the Hamiltonian Adaptive Resolution Scheme (H-AdResS) for the solvent. Here, we present a key step toward the estimation of ligand binding affinities for their targets within this approach. This is the implementation of the H-AdResS in the GROMACS code. The accuracy of our implementation is established by calculating hydration free energies of several molecules in water by means of alchemical transformations. The deviations of the GROMOS- and AMBER-based H-AdResS results from the reference fully atomistic simulations are smaller than the accuracy of the force field and/or they are in the range of the published results. Importantly, our predictions are in good agreement with experimental data. The current implementation paves the way to the use of the OB-MM/CG framework for the study of large biological systems.

Automation of absolute protein-ligand binding free energy calculations for docking refinement and compound evaluation

Absolute binding free energy calculations with explicit solvent molecular simulations can provide estimates of protein-ligand affinities, and thus reduce the time and costs needed to find new drug candidates. However, these calculations can be complex to implement and perform. Here, we introduce the software BAT.py, a Python tool that invokes the AMBER simulation package to automate the calculation of binding free energies for a protein with a series of ligands. The software supports the attach-pull-release (APR) and double decoupling (DD) binding free energy methods, as well as the simultaneous decoupling-recoupling (SDR) method, a variant of double decoupling that avoids numerical artifacts associated with charged ligands. We report encouraging initial test applications of this software both to re-rank docked poses and to estimate overall binding free energies. We also show that it is practical to carry out these calculations cheaply by using graphical processing units in common machines that can be built for this purpose. The combination of automation and low cost positions this procedure to be applied in a relatively high-throughput mode and thus stands to enable new applications in early-stage drug discovery.

Application and Comprehensive Analysis of Neighbor Approximated Information Theoretic Configurational Entropy Methods to Protein-Ligand Binding Cases

The binding entropy is an important thermodynamic quantity which has numerous applications in studies of the biophysical process, and configurational entropy is often one of the major contributors in it. Therefore, its accurate estimation is important, though it is challenging mostly due to sampling limitations, anharmonicity, and multimodality of atomic fluctuations. The present work reports a Neighbor Approximated Maximum Information Spanning Tree (A-MIST) method for conformational entropy and presents its performance and computational advantage over conventional Mutual Information Expansion (MIE) and Maximum Information Spanning Tree (MIST) for two protein-ligand binding cases: indirubin-5-sulfonate to Plasmodium falciparum Protein Kinase 5 (PfPK5) and P. falciparum RON2-peptide to P. falciparum Apical Membrane Antigen 1 (PfAMA1). Important structural regions considering binding configurational entropy are identified, and physical origins for such are discussed. A thorough performance evaluation is done of a set of four entropy estimators (Maximum Likelihood (ML), Miller-Madow (MM), Chao-Shen (CS), and James and Stein shrinkage (JS)) with known varying degrees of sensitivity of the entropy estimate on the extent of sampling, each with two schemes for discretization of fluctuation data of Degrees of Freedom (DFs) to estimate Probability Density Functions (PDFs). Our comprehensive evaluation of influences of variations of parameters shows Neighbor Approximated MIE (A-MIE) outperforms MIE in terms of convergence and computational efficiency. In the case of A-MIE/MIE, results are sensitive to the choice of root atoms, graph search algorithm used for the Bond-Angle-Torsion (BAT) conversion, and entropy estimator, while A-MIST/MIST are not. A-MIST yields binding entropy within 0.5 kcal/mol of MIST with only 20-30% computation. Moreover, all these methods have been implemented in an OpenMP/MPI hybrid parallel C++11 code, and also a python package for data preprocessing and entropy contribution analysis is developed and made available. A comparative analysis of features of current implementation and existing tools is also presented.

Variationally derived intermediates for correlated free-energy estimates between intermediate states

Free-energy difference calculations based on atomistic simulations generally improve in accuracy when sampling from a sequence of intermediate equilibrium thermodynamic states that bridge the configuration space between two states of interest. For reasons of efficiency, usually the same samples are used to calculate the stepwise difference of such an intermediate to both adjacent intermediates. However, this procedure violates the assumption of uncorrelated estimates that is necessary to derive both the optimal sequence of intermediate states and the widely used Bennett acceptance ratio estimator. In this work, via a variational approach, we derive the sequence of intermediate states and the corresponding estimator with minimal mean-squared error that account for these correlations and assess its accuracy.

On the faithfulness of molecular mechanics representations of proteins towards quantum-mechanical energy surfaces

Force fields based on molecular mechanics (MM) are the main computational tool to study the relationship between protein structure and function at the molecular level. To validate the quality of such force fields, high-level quantum-mechanical (QM) data are employed to test their capability to reproduce the features of all major conformational substates of a series of blocked amino acids. The phase-space overlap between MM and QM is quantified in terms of the average structural reorganization energies over all energy minima. Here, the structural reorganization energy is the MM potential-energy difference between the structure of the respective QM energy minimum and the structure of the closest MM energy minimum. Thus, it serves as a measure for the relative probability of visiting the QM minimum during an MM simulation. We evaluate variants of the AMBER, CHARMM, GROMOS and OPLS biomolecular force fields. In addition, the two blocked amino acids alanine and serine are used to demonstrate the dependence of the measured agreement on the QM method, the phase, and the conformational preferences. Blocked serine serves as an example to discuss possible improvements of the force fields, such as including polarization with Drude particles, or using tailored force fields. The results show that none of the evaluated force fields satisfactorily reproduces all energy minima. By decomposing the average structural reorganization energies in terms of individual energy terms, we can further assess the individual weaknesses of the parametrization strategies of each force field. The dominant problem for most force fields appears to be the van der Waals parameters, followed to a lesser degree by dihedral and bonded terms. Our results show that performing a simple QM energy optimization from an MM-optimized structure can be a first test of the validity of a force field for a particular target molecule.

Enhanced Jarzynski free energy calculations using weighted ensemble

The free energy of transitions between stable states is the key thermodynamic quantity that governs the relative probabilities of the forward and reverse reactions and the ratio of state probabilities at equilibrium. The binding free energy of a drug and its receptor is of particular interest, as it serves as an optimization function for drug design. Over the years, many computational methods have been developed to calculate binding free energies, and while many of these methods have a long history, issues such as convergence of free energy estimates and the projection of a binding process onto order parameters remain. Over 20 years ago, the Jarzynski equality was derived with the promise to calculate equilibrium free energies by measuring the work applied to short nonequilibrium trajectories. However, these calculations were found to be dominated by trajectories with low applied work that occur with extremely low probability. Here, we examine the combination of weighted ensemble algorithms with the Jarzynski equality. In this combined method, an ensemble of nonequilibrium trajectories are run in parallel, and cloning and merging operations are used to preferentially sample low-work trajectories that dominate the free energy calculations. Two additional methods are also examined: (i) a novel weighted ensemble resampler that samples trajectories directly according to their importance to the work of work and (ii) the diffusion Monte Carlo method using the applied work as the selection potential. We thoroughly examine both the accuracy and efficiency of unbinding free energy calculations for a series of model Lennard-Jones atom pairs with interaction strengths ranging from 2 kcal/mol to 20 kcal/mol. We find that weighted ensemble calculations can more efficiently determine accurate binding free energies, especially for deeper Lennard-Jones well depths.

Chirality Effects in Biomolecular Systems: Calculation of the Relative Free Energies by Molecular Dynamics Simulations

Chirality plays an essential role in chemical and biological sciences. At the molecular level, the effects associated with this phenomenon can be studied by using the well-established technique of molecular dynamics simulations. In this work, we present several approaches suited for the molecular dynamics-based free energy calculation in chiral systems. In particular, we have proposed and tested the following strategies relying on the application of general, enhanced sampling methods: (i) biased sampling in the two-dimensional space, along the coordinates defined by the values of the selected torsional angles; (ii) biased sampling in the one- or two-dimensional space, along the path-based coordinate(s); (iii) rational alteration of the system's Hamiltonian in order to enable the interconversion between stereoisomers and reweighting the biased distribution of configurations; (iv) using the free energy landscape generated within approaches (i) or (ii) as time-independent bias in order to further improve sampling efficiency and simultaneously account for multiple chiral centers. All approaches have been tested on a set of model compounds (fenoterol, fructofuranose, and bromochlorofluoromethane), demonstrating the good performance but also some differences in the range of their applicabilities.

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.

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

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

BRIDGE: An Open Platform for Reproducible High-Throughput Free Energy Simulations

Biomolecular Reaction and Interaction Dynamics Global Environment (BRIDGE) is an open-source web platform developed with the aim to provide an environment for the design of reliable methods to conduct reproducible biomolecular simulations. It is built on the well-known Galaxy bioinformatics platform. Through this, BRIDGE hosts computational chemistry tools on public web servers for internet use and provides machine- and operating-system-independent portability using the Docker container platform for local use. This construction improves the accessibility, shareability, and reproducibility of computational methods for molecular simulations. Here we present integrated free energy tools (or apps) to calculate absolute binding free energies (ABFEs) and relative binding free energies (RBFEs), as illustrated through use cases. We present free energy perturbation (FEP) methods contained in various software packages such as GROMACS and YANK that are independent of hardware configuration, software libraries, or operating systems when ported in the Galaxy-BRIDGE Docker container platform. By performing cyclin-dependent kinase 2 (CDK2) FEP calculations on geographically dispersed web servers (in Africa and Europe), we illustrate that large-scale computations can be performed using the exact same codes and methodology by collaborating groups through publicly shared protocols and workflows. The ease of public sharing and independent reproduction of simulations via BRIDGE makes possible an open review of methods and complete simulation protocols. This makes the discovery of inhibitors for drug targets accessible to nonexperts and the computer experiments that are used to arrive at leads verifiable by experts and reviewers. We illustrate this on β-galactoside α-2,3-sialyltransferase I (ST3Gal-I), a breast cancer drug target, where a combination of RBFE and ABFE methods are used to compute the binding free energies of three inhibitors.

Evaluating the Performance of a Non-Bonded Cu2+ Model Including Jahn-Teller Effect into the Binding of Tyrosinase Inhibitors

Tyrosinase (TYR) is a metalloenzyme classified as a type-3 copper protein, which is involved in the synthesis of melanin through a catalytic process beginning with the conversion of the amino acid l-Tyrosine (l-Tyr) to l-3,4-dihydroxyphenylalanine (l-DOPA). It plays an important role in the mechanism of melanogenesis in various organisms including mammals, plants, and fungi. Herein, we used a combination of computational molecular modeling techniques including molecular dynamic (MD) simulations and the linear interaction energy (LIE) model to evaluate the binding free energy of a set of analogs of kojic acid (KA) in complex with TYR. For the MD simulations, we used a dummy model including the description of the Jahn–Teller effect for Cu²⁺ ions in the active site of this enzyme. Our results show that the LIE model predicts the TYR binding affinities of the inhibitor in close agreement to experimental results. Overall, we demonstrate that the classical model provides a suitable description of the main interactions between analogs of KA and Cu²⁺ ions in the active site of TYR.

Rigorous Free Energy Simulations in Virtual Screening

Virtual high throughput screening (vHTS) in drug discovery is a powerful approach to identify hits: when applied successfully, it can be much faster and cheaper than experimental high-throughput screening approaches. However, mainstream vHTS tools have significant limitations: ligand-based methods depend on knowledge of existing chemical matter, while structure-based tools such as docking involve significant approximations that limit their accuracy. Recent advances in scientific methods coupled with dramatic speedups in computational processing with GPUs make this an opportune time to consider the role of more rigorous methods that could improve the predictive power of vHTS workflows. In this Perspective, we assert that alchemical binding free energy methods using all-atom molecular dynamics simulations have matured to the point where they can be applied in virtual screening campaigns as a final scoring stage to prioritize the top molecules for experimental testing. Specifically, we propose that alchemical absolute binding free energy (ABFE) calculations offer the most direct and computationally efficient approach within a rigorous statistical thermodynamic framework for computing binding energies of diverse molecules, as is required for virtual screening. ABFE calculations are particularly attractive for drug discovery at this point in time, where the confluence of large-scale genomics data and insights from chemical biology have unveiled a large number of promising disease targets for which no small molecule binders are known, precluding ligand-based approaches, and where traditional docking approaches have foundered to find progressible chemical matter.

Multisecond ligand dissociation dynamics from atomistic simulations

Coarse-graining of fully atomistic molecular dynamics simulations is a long-standing goal in order to allow the description of processes occurring on biologically relevant timescales. For example, the prediction of pathways, rates and rate-limiting steps in protein-ligand unbinding is crucial for modern drug discovery. To achieve the enhanced sampling, we perform dissipation-corrected targeted molecular dynamics simulations, which yield free energy and friction profiles of molecular processes under consideration. Subsequently, we use these fields to perform temperature-boosted Langevin simulations which account for the desired kinetics occurring on multisecond timescales and beyond. Adopting the dissociation of solvated sodium chloride, trypsin-benzamidine and Hsp90-inhibitor protein-ligand complexes as test problems, we reproduce rates from molecular dynamics simulation and experiments within a factor of 2–20, and dissociation constants within a factor of 1–4. Analysis of friction profiles reveals that binding and unbinding dynamics are mediated by changes of the surrounding hydration shells in all investigated systems.

Protein-ligand unbinding processes are out of reach for atomistic simulations due to time-scale involved. Here the authors demonstrate an approach relying on dissipation-corrected targeted molecular dynamics that enables to provide binding and unbinding rates with a speed-up of several orders of magnitude.

Efficient Estimation of Absolute Binding Free Energy for a Homeodomain-DNA Complex from Nonequilibrium Pulling Simulations

Estimation of binding free energies is one of the central aims of simulations of biomolecular complexes. We explore the accuracy and efficiency of setups based on nonequilibrium pulling simulations applied to the estimation of binding affinities of DNA-binding proteins. Absolute binding free energies are calculated over a range of temperatures and compared to results obtained previously using an equilibrium method. We show that realistic binding affinities can be obtained with the presented nonequilibrium approach, which also entails lower computational requirements. Errors of the binding free energy estimates are investigated and are shown to be comparable to those observed previously. Bounds are provided on the convergence of the errors with respect to the number of pulling simulations performed and with respect to the applied pull rate.

Umbrella Visualization: A method of analysis dedicated to glycan flexibility with UnityMol

N-glycosylation is a post-translational modification heavily impacting protein functions. Some alterations of glycosylation, such as sialic acid hydrolysis, are related to protein dysfunction. Because of their high flexibility and the many reactive groups of the glycan chains, studying glycans with in vitro methods is a challenging task. Molecular dynamics is a useful tool and probably the only one in biology able to overcome this problem and gives access to conformational information through exhaustive sampling. To better decipher the impact of N-glycans, the analysis and visualization of their influence over time on protein structure is a prerequisite. We developed the Umbrella Visualization, a graphical method that assigns the glycan intrinsic flexibility during a molecular dynamics trajectory. The density plot generated by this method brought relevant informations regarding glycans dynamics and flexibility, but needs further development in order to integrate an accurate description of the protein topology and its interactions. We propose here to transform this analysis method into a visualization mode in UnityMol. UnityMol is a molecular editor, viewer and prototyping platform, coded in C#. The new representation of glycan chains presented in this study takes into account both the main positions adopted by each antenna of a glycan and their statistical relevance. By displaying the collected data on the protein surface, one is then able to investigate the protein/glycan interactions.