Lambda-dynamics free energy simulation methods

How to Sample Dozens of Substitutions per Site with λ Dynamics

Alchemical free energy methods are useful in computer-aided drug design and computational protein design because they provide rigorous statistical mechanics-based estimates of free energy differences from molecular dynamics simulations. λ dynamics is a free energy method with the ability to characterize combinatorial chemical spaces spanning thousands of related systems within a single simulation, which gives it a distinct advantage over other alchemical free energy methods that are mostly limited to pairwise comparisons. Recently developed methods have improved the scalability of λ dynamics to perturbations at many sites; however, the size of chemical space that can be explored at each individual site has previously been limited to fewer than ten substituents. As the number of substituents increases, the volume of alchemical space corresponding to nonphysical alchemical intermediates grows exponentially relative to the size corresponding to the physical states of interest. Beyond nine substituents, λ dynamics simulations become lost in an alchemical morass of intermediate states. In this work, we introduce new biasing potentials that circumvent excessive sampling of intermediate states by favoring sampling of physical end points relative to alchemical intermediates. Additionally, we present a more scalable adaptive landscape flattening algorithm for these larger alchemical spaces. Finally, we show that this potential enables more efficient sampling in both protein and drug design test systems with up to 24 substituents per site, enabling, for the first time, simultaneous simulation of all 20 amino acids.

Replica Exchange of Expanded Ensembles: A Generalized Ensemble Approach with Enhanced Flexibility and Parallelizability

Generalized ensemble methods such as Hamiltonian replica exchange (HREX) and expanded ensemble (EE) have been shown effective in free energy calculations for various contexts, given their ability to circumvent free energy barriers via nonphysical pathways defined by states with different modified Hamiltonians. However, both HREX and EE methods come with drawbacks, such as limited flexibility in parameter specification or the lack of parallelizability for more complicated applications. To address this challenge, we present the method of replica exchange of expanded ensembles (REXEE), which integrates the principles of HREX and EE methods by periodically exchanging coordinates of EE replicas sampling different yet overlapping sets of alchemical states. With the solvation free energy calculation of anthracene and binding free energy calculation of the CB7-10 binding complex, we show that the REXEE method achieves the same level of accuracy in free energy calculations as the HREX and EE methods, while offering enhanced flexibility and parallelizability. Additionally, we examined REXEE simulations with various setups to understand how different exchange frequencies and replica configurations influence the sampling efficiency in the fixed-weight phase and the weight convergence in the weight-updating phase. The REXEE approach can be further extended to support asynchronous parallelization schemes, allowing looser communications between larger numbers of loosely coupled processors such as cloud computing and therefore promising much more scalable and adaptive executions of alchemical free energy calculations. All algorithms for the REXEE method are available in the Python package ensemble_md, which offers an interface for REXEE simulation management without modifying the source code in GROMACS.

A λ-Dynamics Investigation of Insulin Wakayama and Other A3 Variant Binding Affinities to the Insulin Receptor

Insulin Wakayama is a clinical insulin variant where a conserved valine at the third residue on insulin's A chain (ValA3) is replaced with a leucine (LeuA3), weakening insulin receptor (IR) binding by 140-500-fold. This severe impact on binding from a subtle modification has posed an intriguing problem for decades. Although experimental investigations of natural and unnatural A3 mutations have highlighted the sensitivity of insulin-IR binding at this site, atomistic explanations of these binding trends have remained elusive. We investigate this problem computationally using λ-dynamics free energy calculations to model structural changes in response to perturbations of the ValA3 side chain and to calculate associated relative changes in binding free energy (ΔΔGbind). The Wakayama LeuA3 mutation and seven other A3 substitutions were studied in this work. The calculated ΔΔGbind results showed high agreement compared to experimental binding potencies with a Pearson correlation of 0.88 and a mean unsigned error of 0.68 kcal/mol. Extensive structural analyses of λ-dynamics trajectories revealed that critical interactions were disrupted between insulin and the insulin receptor as a result of the A3 mutations. This investigation also quantifies the effect that adding an A3 Cδ atom or losing an A3 Cγ atom has on insulin's binding affinity to the IR. Thus, λ-dynamics was able to successfully model the effects of mutations to insulin's A3 side chain on its protein-protein interactions with the IR and shed new light on a decades-old mystery: the exquisite sensitivity of hormone-receptor binding to a subtle modification of an invariant insulin residue.

A λ-dynamics investigation of insulin Wakayama and other A3 variant binding affinities to the insulin receptor

Insulin Wakayama is a clinical insulin variant where a conserved valine at the third residue on insulin's A chain (ValA3) is replaced with a leucine (LeuA3), impairing insulin receptor (IR) binding by 140-500 fold. This severe impact on binding from such a subtle modification has posed an intriguing problem for decades. Although experimental investigations of natural and unnatural A3 mutations have highlighted the sensitivity of insulin-IR binding to minor changes at this site, an atomistic explanation of these binding trends has remained elusive. We investigate this problem computationally using λ-dynamics free energy calculations to model structural changes in response to perturbations of the ValA3 side chain and to calculate associated relative changes in binding free energy (ΔΔGbind). The Wakayama LeuA3 mutation and seven other A3 substitutions were studied in this work. The calculated ΔΔGbind results showed high agreement compared to experimental binding potencies with a Pearson correlation of 0.88 and a mean unsigned error of 0.68 kcal/mol. Extensive structural analyses of λ-dynamics trajectories revealed that critical interactions were disrupted between insulin and the insulin receptor as a result of the A3 mutations. This investigation also quantifies the effect that adding an A3 Cδ atom or losing an A3 Cγ atom has on insulin's binding affinity to the IR. Thus, λ-dynamics was able to successfully model the effects of subtle modifications to insulin's A3 side chain on its protein-protein interactions with the IR and shed new light on a decades-old mystery: the exquisite sensitivity of hormone-receptor binding to a subtle modification of an invariant insulin residue.

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.

Surface-Accelerated String Method for Locating Minimum Free Energy Paths

We present a surface-accelerated string method (SASM) to efficiently optimize low-dimensional reaction pathways from the sampling performed with expensive quantum mechanical/molecular mechanical (QM/MM) Hamiltonians. The SASM accelerates the convergence of the path using the aggregate sampling obtained from the current and previous string iterations, whereas approaches like the string method in collective variables (SMCV) or the modified string method in collective variables (MSMCV) update the path only from the sampling obtained from the current iteration. Furthermore, the SASM decouples the number of images used to perform sampling from the number of synthetic images used to represent the path. The path is optimized on the current best estimate of the free energy surface obtained from all available sampling, and the proposed set of new simulations is not restricted to being located along the optimized path. Instead, the umbrella potential placement is chosen to extend the range of the free energy surface and improve the quality of the free energy estimates near the path. In this manner, the SASM is shown to improve the exploration for a minimum free energy pathway in regions where the free energy surface is relatively flat. Furthermore, it improves the quality of the free energy profile when the string is discretized with too few images. We compare the SASM, SMCV, and MSMCV using 3 QM/MM applications: a ribozyme methyltransferase reaction using 2 reaction coordinates, the 2'-O-transphosphorylation reaction of Hammerhead ribozyme using 3 reaction coordinates, and a tautomeric reaction in B-DNA using 5 reaction coordinates. We show that SASM converges the paths using roughly 3 times less sampling than the SMCV and MSMCV methods. All three algorithms have been implemented in the FE-ToolKit package made freely available.

Comprehensive, Open-Source, and Automated Workflow for Multisite λ-Dynamics in Lead Optimization

Multisite λ-dynamics (MSLD) is a highly efficient binding free energy calculation method that samples multiple ligands in a single round by assigning different λ values to the alchemical part of each ligand. This method holds great promise for lead optimization (LO) in drug discovery. However, the complex data preparation and simulation process limits its widespread application in diverse protein-ligand systems. To address this challenge, we developed a comprehensive, open-source, and automated workflow for MSLD calculations based on the BLaDE dynamics engine. This workflow incorporates the Ligand Internal and Cartesian coordinate reconstruction-based alignment algorithm (LIC-align) and an optimized maximum common substructure (MCS) search algorithm to accurately generate MSLD multiple topologies with ideal perturbation patterns. Furthermore, our workflow is highly modularized, allowing straightforward integration and extension of various simulation techniques, and is highly accessible to nonexperts. This workflow was validated by calculating the relative binding free energies of large-scale congeneric ligands, many of which have large perturbing groups. The agreement between the calculations and experiments was excellent, with an average unsigned error of 1.08 ± 0.47 kcal/mol. More than 57.1% of the ligands had an error of less than 1.0 kcal/mol, and the perturbations of 6 targets were fully connected via the calculations, while those of 2 targets were connected via both calculations and experimental data. The Pearson correlation coefficient reached 0.88, indicating that the MSLD workflow provides accurate predictions that can guide lead optimization in drug discovery. We also examined the impact of single-site versus multisite perturbations, ligand grouping by perturbing group size, and the position of the anchor atom on the MSLD performance. By integrating our proposed LIC-align and optimized MCS search algorithm along with the coping strategies to handle challenging molecular substructures, our workflow can handle many realistic scenarios more reasonably than all previously published methods. Moreover, we observed that our MSLD workflow achieved similar accuracy to free energy perturbation (FEP) while improving computational efficiency by over 1 order of magnitude in speedup. These findings provide valuable insights and strategies for further MSLD development, making MSLD a competitive tool for lead optimization.

Investigating Polypharmacology through Targeting Known Human Neutrophil Elastase Inhibitors to Proteinase 3

Using a combination of multisite λ-dynamics (MSλD) together with in vitro IC50 assays, we evaluated the polypharmacological potential of a scaffold currently in clinical trials for inhibition of human neutrophil elastase (HNE), targeting cardiopulmonary disease, for efficacious inhibition of Proteinase 3 (PR3), a related neutrophil serine proteinase. The affinities we observe suggest that the dihydropyrimidinone scaffold can serve as a suitable starting point for the establishment of polypharmacologically targeting both enzymes and enhancing the potential for treatments addressing diseases like chronic obstructive pulmonary disease.

Biomolecular dynamics in the 21st century

The relevance of motions in biological macromolecules has been clear since the early structural analyses of proteins by X-ray crystallography. Computer simulations have been applied to provide a deeper understanding of the dynamics of biological macromolecules since 1976, and are now a standard tool in many labs working on the structure and function of biomolecules. In this mini-review we highlight some areas of current interest and active development for simulations, in particular all-atom molecular dynamics simulations.

In silico λ-dynamics predicts protein binding specificities to modified RNAs

RNA modifications shape gene expression through a smorgasbord of chemical changes to canonical RNA bases. Although numbering in the hundreds, only a few RNA modifications are well characterized, in part due to the absence of methods to identify modification sites. Antibodies remain a common tool to identify modified RNA and infer modification sites through straightforward applications. However, specificity issues can result in off-target binding and confound conclusions. This work utilizes in silico λ-dynamics to efficiently estimate binding free energy differences of modification-targeting antibodies between a variety of naturally occurring RNA modifications. Crystal structures of inosine and N6-methyladenosine (m6A) targeting antibodies bound to their modified ribonucleosides were determined and served as structural starting points. λ-Dynamics was utilized to predict RNA modifications that permit or inhibit binding to these antibodies. In vitro RNA-antibody binding assays supported the accuracy of these in silico results. High agreement between experimental and computed binding propensities demonstrated that λ-dynamics can serve as a predictive screen for antibody specificity against libraries of RNA modifications. More importantly, this strategy is an innovative way to elucidate how hundreds of known RNA modifications interact with biological molecules without the limitations imposed by in vitro or in vivo methodologies.

Direct Calculation of Electron Transfer Rates with the Binless Dynamic Histogram Analysis Method

Umbrella sampling molecular dynamics simulations are widely used to enhance sampling along the reaction coordinates of chemical reactions. The effect of the artificial bias can be removed using methods such as the dynamic weighted histogram analysis method (DHAM), which in addition to the global free energy profile also provides kinetic information about barrier-crossing rates directly from the Markov matrix. Here we present a binless formulation of DHAM that extends DHAM to high-dimensional and Hamiltonian-based biasing to allow the study of electron transfer (ET) processes, for which enhanced sampling is usually not possible based on simple geometric grounds. We show the capabilities of binless DHAM on examples such as aqueous ferrous-ferric ET and intramolecular ET in the radical anion of benzoquinone-tetrathiafulvalene-benzoquinone (Q-TTF-Q)-. From classical Hamiltonian-based umbrella sampling simulations and electronic coupling values from quantum chemistry calculations, binless DHAM provides ET rates for adiabatic and nonadiabatic ET reactions alike in excellent agreement with experimental results.

Computational investigation on lipid bilayer disruption induced by amphiphilic Janus nanoparticles: combined effect of Janus balance and charged lipid concentration

Janus nanoparticles (NPs) with charged/hydrophobic compartments have garnered attention for their potential antimicrobial activity. These NPs have been shown to disrupt lipid bilayers in experimental studies, yet the underlying mechanisms of this disruption at the particle-membrane interface remain unclear. To address this knowledge gap, the present study conducts a computational investigation to systematically examine the disruption of lipid bilayers induced by amphiphilic Janus NPs. The focus of this study is on the combined effects of the hydrophobicity of the Janus NP, referred to as the Janus balance, defined as the ratio of hydrophilic to hydrophobic surface coverage, and the concentration of charged phospholipids on the interactions between Janus NPs and lipid bilayers. Computational simulations were conducted using a coarse-grained molecular dynamics (MD) approach. The results of these MD simulations reveal that while the area change of the bilayer increases monotonically with the Janus balance, the effect of charged lipid concentration in the membrane is not easy to be predicted. Specifically, it was found that the concentration of negatively charged lipids is directly proportional to the intensity of membrane disruption. Conversely, positively charged lipids have a negligible effect on membrane defects. This study provides molecular insights into the significant role of Janus balance in the disruption of lipid bilayers by Janus NPs and supports the selectivity of Janus NPs for negatively charged lipid membranes. Furthermore, the anisotropic properties of Janus NPs were found to play a crucial role in their ability to disrupt the membrane via the combination of hydrophobic and electrostatic interactions. This finding is validated by testing the current Janus NP design on a bacterial membrane-mimicking model. This computational study may serve as a foundation for further studies aimed at optimizing the properties of Janus NPs for specific antimicrobial applications.

Distinct chemical factors in hydrolytic reactions catalyzed by metalloenzymes and metal complexes

The selective hydrolysis of the extremely stable phosphoester, peptide and ester bonds of molecules by bio-inspired metal-based catalysts (metallohydrolases) is required in a wide range of biological, biotechnological and industrial applications. Despite the impressive advances made in the field, the ultimate goal of designing efficient enzyme mimics for these reactions is still elusive. Its realization will require a deeper understanding of the diverse chemical factors that influence the activities of both natural and synthetic catalysts. They include catalyst-substrate complexation, non-covalent interactions and the electronic nature of the metal ion, ligand environment and nucleophile. Based on our computational studies, their roles are discussed for several mono- and binuclear metallohydrolases and their synthetic analogues. Hydrolysis by natural metallohydrolases is found to be promoted by a ligand environment with low basicity, a metal bound water and a heterobinuclear metal center (in binuclear enzymes). Additionally, peptide and phosphoester hydrolysis is dominated by two competing effects, i.e. nucleophilicity and Lewis acid activation, respectively. In synthetic analogues, hydrolysis is facilitated by the inclusion of a second metal center, hydrophobic effects, a biological metal (Zn, Cu and Co) and a terminal hydroxyl nucleophile. Due to the absence of the protein environment, hydrolysis by these small molecules is exclusively influenced by nucleophile activation. The results gleaned from these studies will enhance the understanding of fundamental principles of multiple hydrolytic reactions. They will also advance the development of computational methods as a predictive tool to design more efficient catalysts for hydrolysis, Diels-Alder reaction, Michael addition, epoxide opening and aldol condensation.

Alchemical Metadynamics: Adding Alchemical Variables to Metadynamics to Enhance Sampling in Free Energy Calculations

Performing alchemical transformations, in which one molecular system is nonphysically changed to another system, is a popular approach adopted in performing free energy calculations associated with various biophysical processes, such as protein-ligand binding or the transfer of a molecule between environments. While the sampling of alchemical intermediate states in either parallel (e.g., Hamiltonian replica exchange) or serial manner (e.g., expanded ensemble) can bridge the high-probability regions in the configurational space between two end states of interest, alchemical methods can fail in scenarios where the most important slow degrees of freedom in the configurational space are, in large part, orthogonal to the alchemical variable, or if the system gets trapped in a deep basin extending in both the configurational and alchemical space. To alleviate these issues, we propose to use alchemical variables as an additional dimension in metadynamics, making it possible to both sample collective variables and to enhance sampling in free energy calculations simultaneously. In this study, we validate our implementation of "alchemical metadynamics" in PLUMED with test systems and alchemical processes with varying complexities and dimensionalities of collective variable space, including the interconversion between the torsional metastable states of a toy system and the methylation of a nucleoside both in the isolated form and in a duplex. We show that multidimensional alchemical metadynamics can address the challenges mentioned above and further accelerate sampling by introducing configurational collective variables. The method can trivially be combined with other metadynamics-based algorithms implemented in PLUMED. The necessary PLUMED code changes have already been released for general use in PLUMED 2.8.

Deconvoluting binding sites in amyloid nanofibrils using time-resolved spectroscopy

Steady-state fluorescence spectroscopy has a central role not only for sensing applications, but also in biophysics and imaging. Light switching probes, such as ruthenium dipyridophenazine complexes, have been used to study complex systems such as DNA, RNA, and amyloid fibrils. Nonetheless, steady-state spectroscopy is limited in the kind of information it can provide. In this paper, we use time-resolved spectroscopy for studying binding interactions between amyloid-β fibrillar structures and photoluminescent ligands. Using time-resolved spectroscopy, we demonstrate that ruthenium complexes with a pyrazino phenanthroline derivative can bind to two distinct binding sites on the surface of fibrillar amyloid-β, in contrast with previous studies using steady-state photoluminescence spectroscopy, which only identified one binding site for similar compounds. The second elusive binding site is revealed when deconvoluting the signals from the time-resolved decay traces, allowing the determination of dissociation constants of 3 and 2.2 μM. Molecular dynamic simulations agree with two binding sites on the surface of amyloid-β fibrils. Time-resolved spectroscopy was also used to monitor the aggregation of amyloid-β in real-time. In addition, we show that common polypyridine complexes can bind to amyloid-β also at two different binding sites. Information on how molecules bind to amyloid proteins is important to understand their toxicity and to design potential drugs that bind and quench their deleterious effects. The additional information contained in time-resolved spectroscopy provides a powerful tool not only for studying excited state dynamics but also for sensing and revealing important information about the system including hidden binding sites.

Reliable and Accurate Prediction of Single-Residue pKa Values through Free Energy Perturbation Calculations

Accurate prediction of the pK_a's of protein residues is crucial to many applications in biological simulation and drug discovery. Here, we present the use of free energy perturbation (FEP) calculations for the prediction of single-protein residue pK_a values. We begin with an initial set of 191 residues with experimentally determined pK_a values. To isolate sampling limitations from force field inaccuracies, we develop an algorithm to classify residues whose environments are significantly affected by crystal packing effects. We then report an approach to identify buried histidines that require significant sampling beyond what is achieved in typical FEP calculations. We therefore define a clean data set not requiring algorithms capable of predicting major conformational changes on which other pK_a prediction methods can be tested. On this data set, we report an RMSE of 0.76 pK_a units for 35 ASP residues, 0.51 pK_a units for 44 GLU residues, and 0.67 pK_a units for 76 HIS residues.

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.

Best Practices in Constant pH MD Simulations: Accuracy and Sampling

Various approaches have been proposed to include the effect of pH in molecular dynamics (MD) simulations. Among these, the λ-dynamics approach proposed by Brooks and co-workers [Kong, X.; Brooks III, C. L. J. Chem. Phys.1996, 105, 2414−2423] can be performed with little computational overhead and hfor each typeence be used to routinely perform MD simulations at microsecond time scales, as shown in the accompanying paper [Aho, N. et al. J. Chem. Theory Comput.2022, DOI: 10.1021/acs.jctc.2c00516]. At such time scales, however, the accuracy of the molecular mechanics force field and the parametrization becomes critical. Here, we address these issues and provide the community with guidelines on how to set up and perform long time scale constant pH MD simulations. We found that barriers associated with the torsions of side chains in the CHARMM36m force field are too high for reaching convergence in constant pH MD simulations on microsecond time scales. To avoid the high computational cost of extending the sampling, we propose small modifications to the force field to selectively reduce the torsional barriers. We demonstrate that with such modifications we obtain converged distributions of both protonation and torsional degrees of freedom and hence consistent pK_a estimates, while the sampling of the overall configurational space accessible to proteins is unaffected as compared to normal MD simulations. We also show that the results of constant pH MD depend on the accuracy of the correction potentials. While these potentials are typically obtained by fitting a low-order polynomial to calculated free energy profiles, we find that higher order fits are essential to provide accurate and consistent results. By resolving problems in accuracy and sampling, the work described in this and the accompanying paper paves the way to the widespread application of constant pH MD beyond pK_a prediction.

Alchemical Free Energy Methods Applied to Complexes of the First Bromodomain of BRD4

Accurate and rapid predictions of the binding affinity of a compound to a target are one of the ultimate goals of computer aided drug design. Alchemical approaches to free energy estimations follow the path from an initial state of the system to the final state through alchemical changes of the energy function during a molecular dynamics simulation. Herein, we explore the accuracy and efficiency of two such techniques: relative free energy perturbation (FEP) and multisite lambda dynamics (MSλD). These are applied to a series of inhibitors for the bromodomain-containing protein 4 (BRD4). We demonstrate a procedure for obtaining accurate relative binding free energies using MSλD when dealing with a change in the net charge of the ligand. This resulted in an impressive comparison with experiment, with an average difference of 0.4 ± 0.4 kcal mol-1. In a benchmarking study for the relative FEP calculations, we found that using 20 lambda windows with 0.5 ns of equilibration and 1 ns of data collection for each window gave the optimal compromise between accuracy and speed. Overall, relative FEP and MSλD predicted binding free energies with comparable accuracy, an average of 0.6 kcal mol-1 for each method. However, MSλD makes predictions for a larger molecular space over a much shorter time scale than relative FEP, with MSλD requiring a factor of 18 times less simulation time for the entire molecule space.

Optimizing Multisite λ-Dynamics Throughput with Charge Renormalization

With the ability to sample combinations of alchemical perturbations at multiple sites off a small molecule core, multisite λ-dynamics (MSλD) has become an attractive alternative to conventional alchemical free energy methods for exploring large combinatorial chemical spaces. However, current software implementations dictate that combinatorial sampling with MSλD must be performed with a multiple topology model (MTM), which is nontrivial to create by hand, especially for a series of ligand analogues which may have diverse functional groups attached. This work introduces an automated workflow, referred to as msld_py_prep, to assist in the creation of a MTM for use with MSλD. One approach for partitioning partial atomic charges between ligands to create a MTM, called charge renormalization, is also presented and rigorously evaluated. We find that msld_py_prep greatly accelerates the preparation of MSλD ready-to-use files and that charge renormalization can provide a successful approach for MTM generation, as long as bookending calculations are applied to correct small differences introduced by charge renormalization. Charge renormalization also facilitates the use of many different force field parameters with MSλD, broadening the applicability of MSλD for computer-aided drug design.

Organic Host Encapsulation Effects on Nitrosobenzene Monomer-Dimer Distribution and C-NO Bond Rotation in an Aqueous Solution

The intermolecular (monomer-dimer equilibrium) and intramolecular (C-NO and C-NMe₂ rotations) dynamics of 4-nitrosocumene (1a) and 4-(N,N-dimethylamino)nitrosobenzene (1b), respectively, were found to be controlled by the medium (water) and the host environment (organic capsules and cavitands). The ability of water to shift the equilibrium toward the dimer appears to result from dipolar stabilization of the polar dimer structure and has a resemblance to water's known ability to favor organic cycloaddition reactions. In an aqueous medium, a range of organic hosts selectively include only the nitrosocumene monomer 1a. Encapsulation in the octa acid duplex (OA₂) selects two 1a monomers rather than a dimer structure. Octa acid encapsulation also results in more restricted intramolecular C-N rotations of the guest 1b.

Interactions of Urea-Based Inhibitors with Prostate-Specific Membrane Antigen for Boron Neutron Capture Therapy

In this study, molecular interactions of prostate-specific membrane antigen (PSMA) with five chemically distinct urea-based boron-containing inhibitors have been investigated at the atomic level using molecular docking and molecular dynamics simulations. The PSMA-inhibitor complexations have been analyzed by comparing their binding modes, secondary structures, root-mean-square deviations, noncovalent interactions, principal components, and binding free energies. PSMA is a cell surface glycoprotein upregulated in cancerous cells and can be targeted by boron-labeled inhibitors for boron neutron capture therapy (BNCT). The effective BNCT requires the selective boron delivery to the tumor area and highly specific PSMA-mediated cellular uptake by tumor. Thus, a potent inhibitor must exhibit both high binding affinity and high boron density. The computational results suggest that the chemical nature of inhibitors affects the binding mode and their association with PSMA is primarily dominated by hydrogen bonding, salt bridge, electrostatic, and π-π interactions. The binding free energies (-28.0, -15.2, -43.9, -23.2, and -38.2 kcal/mol) calculated using λ-dynamics for all inhibitors (In1-5) predict preferential binding that is in accordance with experimental data. Among all inhibitors, In5 was found to be the best candidate for BNCT. The binding of this inhibitor to PSMA preserved its overall secondary structure. These results provide computational insights into the coordination flexibility of PSMA and its interaction with various inhibitors. They can be used for the design and synthesis of efficient BNCT agents with improved drug selectivity and high boron percentage.

Pandemic drugs at pandemic speed: infrastructure for accelerating COVID-19 drug discovery with hybrid machine learning- and physics-based simulations on high-performance computers

The race to meet the challenges of the global pandemic has served as a reminder that the existing drug discovery process is expensive, inefficient and slow. There is a major bottleneck screening the vast number of potential small molecules to shortlist lead compounds for antiviral drug development. New opportunities to accelerate drug discovery lie at the interface between machine learning methods, in this case, developed for linear accelerators, and physics-based methods. The two in silico methods, each have their own advantages and limitations which, interestingly, complement each other. Here, we present an innovative infrastructural development that combines both approaches to accelerate drug discovery. The scale of the potential resulting workflow is such that it is dependent on supercomputing to achieve extremely high throughput. We have demonstrated the viability of this workflow for the study of inhibitors for four COVID-19 target proteins and our ability to perform the required large-scale calculations to identify lead antiviral compounds through repurposing on a variety of supercomputers.

BLaDE: A Basic Lambda Dynamics Engine for GPU-Accelerated Molecular Dynamics Free Energy Calculations

There is an accelerating interest in practical applications of alchemical free energy methods to problems in protein design, constant pH simulations, and especially computer-aided drug design. In the present paper, we describe a basic lambda dynamics engine (BLaDE) that enables alchemical free energy simulations, including multisite λ dynamics (MSλD) simulations, on graphical processor units (GPUs). We find that BLaDE is 5 to 8 times faster than the current GPU implementation of MSλD-based free energy calculations in CHARMM. We also demonstrate that BLaDE running standard molecular dynamics attains a performance competitive with and sometimes exceeding that of the highly optimized OpenMM GPU code. BLaDE is available as a standalone program and through an API in CHARMM.

Generalizing the Discrete Gibbs Sampler-Based λ-Dynamics Approach for Multisite Sampling of Many Ligands

In this work, the discrete λ variant of the Gibbs sampler-based λ-dynamics (d-GSλD) method is developed to enable multiple functional group perturbations to be investigated at one or more sites of substitution off a common ligand core. The theoretical framework and special considerations for constructing discrete λ states for multisite d-GSλD are presented. The precision and accuracy of the d-GSλD method is evaluated with three test cases of increasing complexity. Specifically, methyl → methyl symmetric perturbations in water, 1,4-benzene hydration free energies and protein-ligand binding affinities for an example HIV-1 reverse transcriptase inhibitor series are computed with d-GSλD. Complementary MSλD calculations were also performed to compare with d-GSλD's performance. Excellent agreement between d-GSλD and MSλD is observed, with mean unsigned errors of 0.12 and 0.22 kcal/mol for computed hydration and binding free energy test cases, respectively. Good agreement with experiment is also observed, with errors of 0.5-0.7 kcal/mol. These findings support the applicability of the d-GSλD free energy method for a variety of molecular design problems, including structure-based drug design. Finally, a discussion of d-GSλD versus MSλD approaches is presented to compare and contrast features of both methods.

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.

Hybrid MC/MD for protein design

Computational protein design relies on simulations of a protein structure, where selected amino acids can mutate randomly, and mutations are selected to enhance a target property, such as stability. Often, the protein backbone is held fixed and its degrees of freedom are modeled implicitly to reduce the complexity of the conformational space. We present a hybrid method where short molecular dynamics (MD) segments are used to explore conformations and alternate with Monte Carlo (MC) moves that apply mutations to side chains. The backbone is fully flexible during MD. As a test, we computed side chain acid/base constants or pK_a's in five proteins. This problem can be considered a special case of protein design, with protonation/deprotonation playing the role of mutations. The solvent was modeled as a dielectric continuum. Due to cost, in each protein we allowed just one side chain position to change its protonation state and the other position to change its type or mutate. The pK_a's were computed with a standard method that scans a range of pH values and with a new method that uses adaptive landscape flattening (ALF) to sample all protonation states in a single simulation. The hybrid method gave notably better accuracy than standard, fixed-backbone MC. ALF decreased the computational cost a factor of 13.

Accounting for the Central Role of Interfacial Water in Protein-Ligand Binding Free Energy Calculations

Rigorous binding free energy methods in drug discovery are growing in popularity because of a combination of methodological advances, improvements in computer hardware, and workflow automation. These calculations typically use molecular dynamics (MD) to sample from the Boltzmann distribution of conformational states. However, when part or all of the binding sites is inaccessible to the bulk solvent, the time needed for water molecules to equilibrate between bulk solvent and the binding site can be well beyond what is practical with standard MD. This sampling limitation is problematic in relative binding free energy calculations, which compute the reversible work of converting ligand 1 to ligand 2 within the binding site. Thus, if ligand 1 is smaller and/or more polar than ligand 2, the perturbation may allow additional water molecules to occupy a region of the binding site. However, this change in hydration may not be captured by standard MD simulations and may therefore lead to errors in the computed free energy. We recently developed a hybrid Monte Carlo/MD (MC/MD) method, which speeds up the equilibration of water between bulk solvent and buried cavities, while sampling from the intended distribution of states. Here, we report on the use of this approach in the context of alchemical binding free energy calculations. We find that using MC/MD markedly improves the accuracy of the calculations and also reduces hysteresis between the forward and reverse perturbations, relative to matched calculations using only MD with or without the crystallographic water molecules. The present method is available for use in AMBER simulation software.

Accounting for the Central Role of Interfacial Water in Protein-Ligand Binding Free Energy Calculations

Rigorous binding free energy methods in drug discovery are growing in popularity because of a combination of methodological advances, improvements in computer hardware, and workflow automation. These calculations typically use molecular dynamics (MD) to sample from the Boltzmann distribution of conformational states. However, when part or all of the binding sites is inaccessible to the bulk solvent, the time needed for water molecules to equilibrate between bulk solvent and the binding site can be well beyond what is practical with standard MD. This sampling limitation is problematic in relative binding free energy calculations, which compute the reversible work of converting ligand 1 to ligand 2 within the binding site. Thus, if ligand 1 is smaller and/or more polar than ligand 2, the perturbation may allow additional water molecules to occupy a region of the binding site. However, this change in hydration may not be captured by standard MD simulations and may therefore lead to errors in the computed free energy. We recently developed a hybrid Monte Carlo/MD (MC/MD) method, which speeds up the equilibration of water between bulk solvent and buried cavities, while sampling from the intended distribution of states. Here, we report on the use of this approach in the context of alchemical binding free energy calculations. We find that using MC/MD markedly improves the accuracy of the calculations and also reduces hysteresis between the forward and reverse perturbations, relative to matched calculations using only MD with or without the crystallographic water molecules. The present method is available for use in AMBER simulation software.

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.

An Approach to Computing Solvent Reorganization Energy

A method for accurate calculation of reorganization free energy of an explicit solvent solvating a solute molecule is presented. The method relies on the knowledge of solvation free energy known either from experimental measurement or from accurate calculation. An important part of this approach lies in the calculation of entropy in solute-solvent interaction free energy using the interaction entropy method, combined with the calculation of enthalpy of solvent reorganization based on a finite number of solvent molecules near the solute molecule. This interaction entropy-solvent reorganization or IESR method enables one to compute both enthalpy and entropy components of the solvent reorganization energy from MD simulation of the solvated system. The calculated results are determined by the molecular force field only without any empirical parameter. The current method is applied to computing the reorganization energies of water solvent solvating a wide range of solute molecules including both hydrophilic and hydrophobic ones. The accuracy of the approach is indirectly verified by the excellent agreement of individual enthalpies and entropies of the solvation energies between theory and experiment.

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

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

Exploring pH Dependent Host/Guest Binding Affinities

When the electrostatic environment surrounding binding partners changes between unbound and bound states, the net uptake or release of a proton is possible by either binding partner. This process is pH-dependent in that the free energy required to uptake or release the proton varies with pH. This pH-dependence is typically not considered in conventional free energy methods where the use of fixed protonation states is the norm. In the present paper, we apply a simple two-step approach to calculate the pH-dependent binding free energy of a model cucubit[7]uril host/guest system. By use of λ-dynamics with an enhanced sampling protocol, adaptive landscape flattening, pK_a shifts and reference binding free energies upon complexation were determined. This information enables the construction of pH-dependent binding profiles that accurately capture the pK_a shifts and reproduce binding free energies at the different pH conditions that were observed experimentally. Our calculations illustrate a general framework for computing pH-dependent binding free energies but also point to some issues in modeling the molecular charge distributions within this series of molecules with CGenFF. However, by introducing some minor charge modifications to the CGenFF force field, we saw significant improvement in accuracy of the calculated pK_a shifts.

Computationally predicting binding affinity in protein-ligand complexes: free energy-based simulations and machine learning-based scoring functions

Accurately predicting protein-ligand binding affinities can substantially facilitate the drug discovery process, but it remains as a difficult problem. To tackle the challenge, many computational methods have been proposed. Among these methods, free energy-based simulations and machine learning-based scoring functions can potentially provide accurate predictions. In this paper, we review these two classes of methods, following a number of thermodynamic cycles for the free energy-based simulations and a feature-representation taxonomy for the machine learning-based scoring functions. More recent deep learning-based predictions, where a hierarchy of feature representations are generally extracted, are also reviewed. Strengths and weaknesses of the two classes of methods, coupled with future directions for improvements, are comparatively discussed.

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.

Fast, Accurate, and Reliable Protocols for Routine Calculations of Protein-Ligand Binding Affinities in Drug Design Projects Using AMBER GPU-TI with ff14SB/GAFF

Accurate prediction of the absolute or relative protein-ligand binding affinity is one of the major tasks in computer-aided drug design projects, especially in the stage of lead optimization. In principle, the alchemical free energy (AFE) methods such as thermodynamic integration (TI) or free-energy perturbation (FEP) can fulfill this task, but in practice, a lot of hurdles prevent them from being routinely applied in daily drug design projects, such as the demanding computing resources, slow computing processes, unavailable or inaccurate force field parameters, and difficult and unfriendly setting up and post-analysis procedures. In this study, we have exploited practical protocols of applying the CPU (central processing unit)-TI and newly developed GPU (graphic processing unit)-TI modules and other tools in the AMBER software package, combined with ff14SB/GAFF1.8 force fields, to conduct efficient and accurate AFE calculations on protein-ligand binding free energies. We have tested 134 protein-ligand complexes in total for four target proteins (BACE, CDK2, MCL1, and PTP1B) and obtained overall comparable performance with the commercial Schrodinger FEP+ program (WangJ. Am. Chem. Soc.2015, 137, 2695-2703). The achieved accuracy fits within the requirements for computations to generate effective guidance for experimental work in drug lead optimization, and the needed wall time is short enough for practical application. Our verified protocol provides a practical solution for routine AFE calculations in real drug design projects.

M2 amphipathic helices facilitate pH-dependent conformational transition in influenza A virus

The matrix-2 (M2) protein from influenza A virus is a tetrameric, integral transmembrane (TM) protein that plays a vital role in viral replication by proton flux into the virus. The His37 tetrad is a pH sensor in the center of the M2 TM helix that activates the channel in response to the low endosomal pH. M2 consists of different regions that are believed to be involved in membrane targeting, packaging, nucleocapsid binding, and proton transport. Although M2 has been the target of many experimental and theoretical studies that have led to significant insights into its structure and function under differing conditions, the main mechanism of proton transport, its conformational dynamics, and the role of the amphipathic helices (AHs) on proton conductance remain elusive. To this end, we have applied explicit solvent constant pH molecular dynamics using the multisite λ-dynamics approach (CpHMD^MSλD) to investigate the buried ionizable residues comprehensively and to elucidate their effect on the conformational transition. Our model recapitulates the pH-dependent conformational transition of M2 from closed to open state when the AH domain is included in the M2 construct, revealing the role of the amphipathic helices on this transition and shedding light on the proton-transport mechanism. This work demonstrates the importance of including the amphipathic helices in future experimental and theoretical studies of ion channels. Finally, our work shows that explicit solvent CpHMD^MSλD provides a realistic pH-dependent model for membrane proteins.

In-silico profiling and structural insights into the impact of nSNPs in the P. falciparum acetyl-CoA transporter gene to understand the mechanism of drug resistance in malaria

The continuous emergence of resistance to the available drugs poses major constraints in the development of effective therapeutics against malaria. Malaria drug resistance has been attributed to be the manifestation of numerous factors. For example, mutations in the parasite transporter protein acetyl-CoA transporter (Pfact) can remarkably affect its uptake affinity for a drug molecule against malaria, and hence enhance its susceptibility to resistance. To identify major contributors to its loss of functionality, we have thoroughly scrutinized eight such recently reported resistant mutants, via in-silico tools in terms of alterations in different properties. We performed molecular dynamics simulations of the selected Pfact mutants to gain deeper insights into the structural perturbation and dynamicity. Comparison of residue interaction network map of mutants with that of Wild type (WT) protein suggests structural changes as a result of the mutation(s) that translate into the weakening of intra-protein interactions, especially around the drug binding pocket. This, in turn, diminishes the affinity of drug molecules towards the binding site, which was validated by docking analysis. Finally, collating all the observations, we have delineated R108K mutant to deviate the most from WT protein, which, intriguingly suggests that replacing an amino acid with another of similar nature can even translate into greater functional effects as those with dissimilar substitutions.Communicated by Ramaswamy H. Sarma.

Approaching protein design with multisite λ dynamics: Accurate and scalable mutational folding free energies in T4 lysozyme

The estimation of changes in free energy upon mutation is central to the problem of protein design. Modern protein design methods have had remarkable success over a wide range of design targets, but are reaching their limits in ligand binding and enzyme design due to insufficient accuracy in mutational free energies. Alchemical free energy calculations have the potential to supplement modern design methods through more accurate molecular dynamics based prediction of free energy changes, but suffer from high computational cost. Multisite λ dynamics (MSλD) is a particularly efficient and scalable free energy method with potential to explore combinatorially large sequence spaces inaccessible with other free energy methods. This work aims to quantify the accuracy of MSλD and demonstrate its scalability. We apply MSλD to the classic problem of calculating folding free energies in T4 lysozyme, a system with a wealth of experimental measurements. Single site mutants considering 32 mutations show remarkable agreement with experiment with a Pearson correlation of 0.914 and mean unsigned error of 1.19 kcal/mol. Multisite mutants in systems with up to five concurrent mutations spanning 240 different sequences show comparable agreement with experiment. These results demonstrate the promise of MSλD in exploring large sequence spaces for protein design.

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

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

Overcoming Challenging Substituent Perturbations with Multisite λ-Dynamics: A Case Study Targeting β-Secretase 1

Alchemical free energy calculations have made a dramatic impact upon the field of structure-based drug design by allowing functional group modifications to be explored computationally prior to experimental synthesis and assay evaluation, thereby informing and directing synthetic strategies. In furthering the advancement of this area, a series of 21 β-secretase 1 (BACE1) inhibitors developed by Janssen Pharmaceuticals were examined to evaluate the ability to explore large substituent perturbations, some of which contain scaffold modifications, with multisite λ-dynamics (MSλD), an innovative alchemical free energy framework. Our findings indicate that MSλD is able to efficiently explore all structurally diverse ligand end-states simultaneously within a single MD simulation with a high degree of precision and with reduced computational costs compared to the widely used approach TI/MBAR. Furthermore, computational predictions were shown to be accurate to within 0.5-0.8 kcal/mol when CM1A partial atomic charges were combined with CHARMM or OPLS-AA-based force fields, demonstrating that MSλD is force field independent and a viable alternative to FEP or TI approaches for drug design.

Space Constrained Stereoselective Geometric Isomerization of 1,2-Diphenylcyclopropane and Stilbenes in an Aqueous Medium

Confinement provided by the reaction space alters the photostationary state isomer distribution during the geometric isomerization of excited 1,2-diphenylcyclopropane and stilbenes. The selectivity in 1,2-diphenylcyclopropane is suggested to result from the supramolecular steric hindrance exerted by the medium for the rotational motion. The alteration in the selectivity between a dimethyl and n-propyl substituted stilbenes is attributed to the medium influence on the location of the transition state on the ground state surface.

Template-Based Method for Conformation Generation and Scoring for Congeneric Series of Ligands

Physics-based prediction of protein-ligand binding affinities for a congeneric series of ligands in lead optimization requires their geometries as a first step. In this paper, we report a method that uses the 3D conformation of a lead compound in complex with a protein as a template to generate conformations of a series of related analog compounds. The method uses the Maximal Common Substructure (MCS) computed between lead and analog ligands to assign coordinates for the atoms shared between the ligands. For the differing atoms, a conformation generation procedure is implemented that results in a diversity of conformations. The generated conformations are sorted using a score based on the Molecular Mechanics and Generalized Born with Solvent Accessible Surface Area contribution (MM-GBSA) method. The accuracy of the generated conformations is tested retrospectively using a cross-validation approach applied to four data sets obtained from the Drug Design Data Resource (D3R) by measuring the RMSD of the top scored conformation with respect to the crystallographic pose. The scoring ability of the method is independently assessed using data for the same protein targets to test the rank ordering ability and separating active and inactive ligands. We tested the effect of protein flexibility during structural optimization and scoring approaches with and without strain energies. Retrospective validation on data sets comprising 4 targets shows that the method outperforms random selection for all targets and outperforms a molecular weight-based null model in 3 out of 4 targets in separating active and inactive compounds. Therefore, the presented method is expected to be of utility in lead optimization for rapidly screening analog ligands and generating initial conformations for use in more detailed physics-based binding affinity prediction methods.

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.