Efficient and Flexible Algorithm for Free Energy Calculations Using the λ-Dynamics Approach

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.

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.

Drug Discovery by Automated Adaptation of Chemical Structure and Identity

Computer-aided drug design offers the potential to dramatically reduce the cost and effort required for drug discovery. While screening-based methods are valuable in the early stages of hit identification, they are frequently succeeded by iterative, hypothesis-driven computations that require recurrent investment of human time and intuition. To increase automation, we introduce a computational method for lead refinement that combines concerted dynamics of the ligand/protein complex via molecular dynamics simulations with integrated Monte Carlo-based changes in the chemical formula of the ligand. This approach, which we refer to as ligand-exchange Monte Carlo molecular dynamics, accounts for solvent- and entropy-based contributions to competitive binding free energies by coupling the energetics of bound and unbound states during the ligand-exchange attempt. Quantitative comparison of relative binding free energies to reference values from free energy perturbation, conducted in vacuum, indicates that ligand-exchange Monte Carlo molecular dynamics simulations sample relevant conformational ensembles and are capable of identifying strongly binding compounds. Additional simulations demonstrate the use of an implicit solvent model. We speculate that the use of chemical graphs in which exchanges are only permitted between ligands with sufficient similarity may enable an automated search to capture some of the benefits provided by human intuition during hypothesis-guided lead refinement.

Addressing Intersite Coupling Unlocks Large Combinatorial Chemical Spaces for Alchemical Free Energy Methods

Alchemical free energy methods are playing a growing role in molecular design, both for computer-aided drug design of small molecules and for computational protein design. Multisite λ dynamics (MSλD) is a uniquely scalable alchemical free energy method that enables more efficient exploration of combinatorial alchemical spaces encountered in molecular design, but simulations have typically been limited to a few hundred ligands or sequences. Here, we focus on coupling between sites to enable scaling to larger alchemical spaces. We first discuss updates to the biasing potentials that facilitate MSλD sampling to include coupling terms and show that this can provide more thorough sampling of alchemical states. We then harness coupling between sites by developing a new free energy estimator based on the Potts models underlying direct coupling analysis, a method for predicting contacts from sequence coevolution, and find it yields more accurate free energies than previous estimators. The sampling requirements of the Potts model estimator scale with the square of the number of sites, a substantial improvement over the exponential scaling of the standard estimator. This opens up exploration of much larger alchemical spaces with MSλD for molecular design.

Exploring the trimerization process of a transmembrane helix with an ionizable residue by molecular dynamics simulations: a case study of transmembrane domain 5 of LMP-1

The oligomerization of membrane proteins is an important biological process that plays a critical role in the initialization of membrane protein receptor signaling. Unveiling how transmembrane segments oligomerize is critical for understanding the mechanism of membrane receptor signaling activation. Owing to the complicated membrane environment and the extraordinary dynamic properties of the ionizable residues in the transmembrane segment, it is extremely challenging to thoroughly understand the oligomerization process of the transmembrane domain. In this study, transmembrane domain 5 (TMD5) of latent membrane protein-1 from Epstein-Barr virus was used as a prototype model to investigate the trimerization process of the transmembrane segment with ionizable residues. The trimerization process of TMD5 was rebuilt and investigated via conventional molecular dynamics simulations and constant-pH molecular dynamics simulations. When TMD5s approached each other, the tilting angles of the TMD5 monomer decreased. TMD5s formed stable trimers until two interacting sites (D150s and Q139s) along each transmembrane helix were created to lock the TMD5s. The pK_a values of D150 shifted toward neutral states in the membrane environment. When TMD5s were monomers, the pK_a shift of D150 was mainly influenced by its microenvironment in the lipid bilayer. When TMD5s were moving close to each other, protein-protein interactions became the main contributing factor for the pK_a shift of D150s. Overall, this work elucidates the behavior of the TMD5 helix and the pK_a shift of ionizable residue D150 in the process of TMD5 oligomerization. This study may provide insight into the development of agents for targeting the oligomerization of membrane proteins.

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.

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.

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.

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.

Accelerating equilibrium isotope effect calculations. II. Stochastic implementation of direct estimators

Path integral calculations of equilibrium isotope effects and isotopic fractionation are expensive due to the presence of path integral discretization errors, statistical errors, and thermodynamic integration errors. Whereas the discretization errors can be reduced by high-order factorization of the path integral and statistical errors by using centroid virial estimators, two recent papers proposed alternative ways to completely remove the thermodynamic integration errors: Cheng and Ceriotti [J. Chem. Phys. 141, 244112 (2015)] employed a variant of free-energy perturbation called "direct estimators," while Karandashev and Vaníček [J. Chem. Phys. 143, 194104 (2017)] combined the thermodynamic integration with a stochastic change of mass and piecewise-linear umbrella biasing potential. Here, we combine the former approach with the stochastic change in mass in order to decrease its statistical errors when applied to larger isotope effects and perform a thorough comparison of different methods by computing isotope effects first on a harmonic model and then on methane and methanium, where we evaluate all isotope effects of the form CH_4-xD_x/CH₄ and CH_5-xD_x ⁺/CH₅ ⁺, respectively. We discuss the reasons for a surprising behavior of the original method of direct estimators, which performed well for a much larger range of isotope effects than what had been expected previously, as well as some implications of our work for the more general problem of free energy difference calculations.

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.

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.

Investigating coordination flexibility of glycerophosphodiesterase (GpdQ) through interactions with mono-, di-, and triphosphoester (NPP, BNPP, GPE, and paraoxon) substrates

Interactions of the catalytically active binuclear form of glycerophosphodiesterase (GpdQ) with chemically diverse substrates, i.e. phosphomono-, phosphodi-, and phosphotriester have been investigated using molecular dynamics (MD) simulations.

A GPU-Accelerated Parameter Interpolation Thermodynamic Integration Free Energy Method

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

Gibbs Sampler-Based λ-Dynamics and Rao-Blackwell Estimator for Alchemical Free Energy Calculation

λ-dynamics is a generalized ensemble method for alchemical free energy calculations. In traditional λ-dynamics, the alchemical switch variable λ is treated as a continuous variable ranging from 0 to 1 and an empirical estimator is utilized to approximate the free energy. In the present article, we describe an alternative formulation of λ-dynamics that utilizes the Gibbs sampler framework, which we call Gibbs sampler-based λ-dynamics (GSLD). GSLD, like traditional λ-dynamics, can be readily extended to calculate free energy differences between multiple ligands in one simulation. We also introduce a new free energy estimator, the Rao-Blackwell estimator (RBE), for use in conjunction with GSLD. Compared with the current empirical estimator, the advantage of RBE is that RBE is an unbiased estimator and its variance is usually smaller than the current empirical estimator. We also show that the multistate Bennett acceptance ratio equation or the unbinned weighted histogram analysis method equation can be derived using the RBE. We illustrate the use and performance of this new free energy computational framework by application to a simple harmonic system as well as relevant calculations of small molecule relative free energies of solvation and binding to a protein receptor. Our findings demonstrate consistent and improved performance compared with conventional alchemical free energy methods.

Adaptive Landscape Flattening Accelerates Sampling of Alchemical Space in Multisite λ Dynamics

Multisite λ dynamics (MSλD) is a powerful emerging method in free energy calculation that allows prediction of relative free energies for a large set of compounds from very few simulations. Calculating free energy differences between substituents that constitute large volume or flexibility jumps in chemical space is difficult for free energy methods in general, and for MSλD in particular, due to large free energy barriers in alchemical space. This study demonstrates that a simple biasing potential can flatten these barriers and introduces an algorithm that determines system specific biasing potential coefficients. Two sources of error, deep traps at the end points and solvent disruption by hard-core potentials, are identified. Both scale with the size of the perturbed substituent and are removed by sharp biasing potentials and a new soft-core implementation, respectively. MSλD with landscape flattening is demonstrated on two sets of molecules: derivatives of the heat shock protein 90 inhibitor geldanamycin and derivatives of benzoquinone. In the benzoquinone system, landscape flattening leads to 2 orders of magnitude improvement in transition rates between substituents and robust solvation free energies. Landscape flattening opens up new applications for MSλD by enabling larger chemical perturbations to be sampled with improved precision and accuracy.

Biasing Potential Replica Exchange Multisite λ-Dynamics for Efficient Free Energy Calculations

Traditional free energy calculation methods are well-known for their drawbacks in scalability and speed in converging results particularly for calculations with large perturbations. In the present work, we report on the development of biasing potential replica exchange multisite λ-dynamics (BP-REX MSλD), which is a free energy method that is capable of performing simultaneous alchemical free energy transformations, including perturbations between flexible moieties. BP-REX MSλD and the original MSλD are applied to a series of symmetrical 2,5-benzoquinone derivatives covering a diverse chemical space and range of conformational flexibility. Improved λ-space sampling is observed for the BP-REX MSλD simulations, yielding a 2-5-fold increase in the number of transitions between substituents compared to traditional MSλD. We also demonstrate the efficacy of varying the value of c, the parameter that controls the ruggedness of the landscape mediating the sampling of λ-states, based on the flexibility of the fragment. Finally, we developed a protocol for maximizing the transition frequency between fragments. This protocol reduces the "kinetic barrier" for alchemically transforming fragments by grouping and ordering based on volume. These findings are applied to a challenging test set involving a series of geldanamycin-based inhibitors of heat shock protein 90 (Hsp90). Even though the perturbations span volume changes by as large as 60 Å(3), the values for the free energy change achieve an average unsigned error (AUE) of 1.5 kcal/mol relative to experimental Kd measurements with a reasonable correlation (R = 0.56). Our results suggest that the BP-REX MSλD algorithm is a highly efficient and scalable free energy method, which when utilized will enable routine calculations on the order of hundreds of compounds using only a few simulations.

Multistate λ-local-elevation umbrella-sampling (MS-λ-LEUS): method and application to the complexation of cations by crown ethers

An extension of the λ-local-elevation umbrella-sampling (λ-LEUS) scheme [ Bieler et al. J. Chem. Theory Comput. 2014 , 10 , 3006 ] is proposed to handle the multistate (MS) situation, i.e. the calculation of the relative free energies of multiple physical states based on a single simulation. The key element of the MS-λ-LEUS approach is to use a single coupling variable Λ controlling successive pairwise mutations between the states of interest in a cyclic fashion. The Λ variable is propagated dynamically as an extended-system variable, using a coordinate transformation with plateaus and a memory-based biasing potential as in λ-LEUS. Compared to other available MS schemes (one-step perturbation, enveloping distribution sampling and conventional λ-dynamics) the proposed method presents a number of important advantages, namely: (i) the physical states are visited explicitly and over finite time periods; (ii) the extent of unphysical space required to ensure transitions is kept minimal and, in particular, one-dimensional; (iii) the setup protocol solely requires the topologies of the physical states; and (iv) the method only requires limited modifications in a simulation code capable of handling two-state mutations. As an initial application, the absolute binding free energies of five alkali cations to three crown ethers in three different solvents are calculated. The results are found to reproduce qualitatively the main experimental trends and, in particular, the experimental selectivity of 18C6 for K(+) in water and methanol, which is interpreted in terms of opposing trends along the cation series between the solvation free energy of the cation and the direct electrostatic interactions within the complex.

Role of Deprotonation Free Energies in pKa Prediction and Molecule Ranking

A computationally efficient classical molecular simulation technique is derived for ranking the pKa values of a set of chemically similar congeneric molecules in an implicit solvent model of water. This uses the deprotonation free energy of the titratable group in the gas and aqueous phases obtained by thermodynamic integration (TI). For a series of alcohols and acids a strong linear correlation is demonstrated between the experimental pKa and the deprotonation free energy difference in the gas and liquid phases. These calculations also show that classical TI is more efficient than slow-growth TI in calculating deprotonation free energies for the series of molecules considered herein.

Membrane environment modulates the pKa values of transmembrane helices

In this work, we apply the recently developed constant pH molecular dynamics technique to study protonation equilibria of titratable side chains in the context of simple transmembrane (TM) helices and explore the effect of pH on their configurations in membrane bilayers. We observe that, despite a significant shift toward neutral states, considerable population of different side chains stay in the charged state that give rise to pK_a values around 9.6 for Asp and Glu and 4.5 to 6 for His and Lys side chains, respectively. These charged states are highly stabilized by favorable interactions between head groups, water molecules, and the charged side chains that are facilitated by substantial changes in the configuration of the peptides. The pH dependent configurations and the measured pK_a values are in good agreement with relatively recent solid state NMR measurements. Our results presented here demonstrate that all-atom constant pH molecular dynamics can be applied to membrane proteins and peptides to obtain reliable pK_a values and pH dependent behavior for these systems.

Distinguishing binders from false positives by free energy calculations: fragment screening against the flap site of HIV protease

Molecular docking is a powerful tool used in drug discovery and structural biology for predicting the structures of ligand–receptor complexes. However, the accuracy of docking calculations can be limited by factors such as the neglect of protein reorganization in the scoring function; as a result, ligand screening can produce a high rate of false positive hits. Although absolute binding free energy methods still have difficulty in accurately rank-ordering binders, we believe that they can be fruitfully employed to distinguish binders from nonbinders and reduce the false positive rate. Here we study a set of ligands that dock favorably to a newly discovered, potentially allosteric site on the flap of HIV-1 protease. Fragment binding to this site stabilizes a closed form of protease, which could be exploited for the design of allosteric inhibitors. Twenty-three top-ranked protein–ligand complexes from AutoDock were subject to the free energy screening using two methods, the recently developed binding energy analysis method (BEDAM) and the standard double decoupling method (DDM). Free energy calculations correctly identified most of the false positives (≥83%) and recovered all the confirmed binders. The results show a gap averaging ≥3.7 kcal/mol, separating the binders and the false positives. We present a formula that decomposes the binding free energy into contributions from the receptor conformational macrostates, which provides insights into the roles of different binding modes. Our binding free energy component analysis further suggests that improving the treatment for the desolvation penalty associated with the unfulfilled polar groups could reduce the rate of false positive hits in docking. The current study demonstrates that the combination of docking with free energy methods can be very useful for more accurate ligand screening against valuable drug targets.

Constant pH molecular dynamics of proteins in explicit solvent with proton tautomerism

pH is a ubiquitous regulator of biological activity, including protein-folding, protein-protein interactions, and enzymatic activity. Existing constant pH molecular dynamics (CPHMD) models that were developed to address questions related to the pH-dependent properties of proteins are largely based on implicit solvent models. However, implicit solvent models are known to underestimate the desolvation energy of buried charged residues, increasing the error associated with predictions that involve internal ionizable residue that are important in processes like hydrogen transport and electron transfer. Furthermore, discrete water and ions cannot be modeled in implicit solvent, which are important in systems like membrane proteins and ion channels. We report on an explicit solvent constant pH molecular dynamics framework based on multi-site λ-dynamics (CPHMD(MSλD)). In the CPHMD(MSλD) framework, we performed seamless alchemical transitions between protonation and tautomeric states using multi-site λ-dynamics, and designed novel biasing potentials to ensure that the physical end-states are predominantly sampled. We show that explicit solvent CPHMD(MSλD) simulations model realistic pH-dependent properties of proteins such as the Hen-Egg White Lysozyme (HEWL), binding domain of 2-oxoglutarate dehydrogenase (BBL) and N-terminal domain of ribosomal protein L9 (NTL9), and the pKa predictions are in excellent agreement with experimental values, with a RMSE ranging from 0.72 to 0.84 pKa units. With the recent development of the explicit solvent CPHMD(MSλD) framework for nucleic acids, accurate modeling of pH-dependent properties of both major class of biomolecules-proteins and nucleic acids is now possible.

pH-sensitive residues in the p19 RNA silencing suppressor protein from carnation Italian ringspot virus affect siRNA binding stability

Tombusviruses, such as Carnation Italian ringspot virus (CIRV), encode a protein homodimer called p19 that is capable of suppressing RNA silencing in their infected hosts by binding to and sequestering short-interfering RNA (siRNA) away from the RNA silencing pathway. P19 binding stability has been shown to be sensitive to changes in pH but the specific amino acid residues involved have remained unclear. Using constant pH molecular dynamics simulations, we have identified key pH-dependent residues that affect CIRV p19-siRNA binding stability at various pH ranges based on calculated changes in the free energy contribution from each titratable residue. At high pH, the deprotonation of Lys60, Lys67, Lys71, and Cys134 has the largest effect on the binding stability. Similarly, deprotonation of several acidic residues (Asp9, Glu12, Asp20, Glu35, and/or Glu41) at low pH results in a decrease in binding stability. At neutral pH, residues Glu17 and His132 provide a small increase in the binding stability and we find that the optimal pH range for siRNA binding is between 7.0 and 10.0. Overall, our findings further inform recent experiments and are in excellent agreement with data on the pH-dependent binding profile.

pH-dependent dynamics of complex RNA macromolecules

The role of pH-dependent protonation equilibrium in modulating RNA dynamics and function is one of the key unanswered questions in RNA biology. Molecular dynamics (MD) simulations can provide insight into the mechanistic roles of protonated nucleotides, but it is only capable of modeling fixed protonation states and requires prior knowledge of the key residue's protonation state. Recently, we developed a framework for constant pH molecular dynamics simulations (CPHMDMSλD) of nucleic acids, where the nucleotides' protonation states are modeled as dynamic variables that are coupled to the structural dynamics of the RNA. In the present study, we demonstrate the application of CPHMDMSλD to the lead-dependent ribozyme; establishing the validity of this approach for modeling complex RNA structures. We show that CPHMDMSλD accurately predicts the direction of the pKa shifts and reproduces experimentally-measured microscopic pKa values with an average unsigned error of 1.3 pKa units. The effects of coupled titration states in RNA structures are modeled, and the importance of conformation sampling is highlighted. The general accuracy of CPHMDMSλD simulations in reproducing pH-dependent observables reported in this work demonstrates that constant pH simulations provides a powerful tool to investigate pH-dependent processes in nucleic acids.

An introduction to best practices in free energy calculations

Free energy calculations are extremely useful for investigating small-molecule biophysical properties such as protein-ligand binding affinities and partition coefficients. However, these calculations are also notoriously difficult to implement correctly. In this chapter, we review standard methods for computing free energy via simulation, discussing current best practices and examining potential pitfalls for computational researchers performing them for the first time. We include a variety of examples and tips for how to set up and conduct these calculations, including applications to relative binding affinities and small-molecule solvation free energies.

Insights from free-energy calculations: protein conformational equilibrium, driving forces, and ligand-binding modes

Accurate free-energy calculations provide mechanistic insights into molecular recognition and conformational equilibrium. In this work, we performed free-energy calculations to study the thermodynamic properties of different states of molecular systems in their equilibrium basin, and obtained accurate absolute binding free-energy calculations for protein-ligand binding using a newly developed M2 algorithm. We used a range of Asp-Phe-Gly (DFG)-in/out p38α mitogen-activated protein kinase inhibitors as our test cases. We also focused on the flexible DFG motif, which is closely connected to kinase activation and inhibitor binding. Our calculations explain the coexistence of DFG-in and DFG-out states of the loop and reveal different components (e.g., configurational entropy and enthalpy) that stabilize the apo p38α conformations. To study novel ligand-binding modes and the key driving forces behind them, we computed the absolute binding free energies of 30 p38α inhibitors, including analogs with unavailable experimental structures. The calculations revealed multiple stable, complex conformations and changes in p38α and inhibitor conformations, as well as balance in several energetic terms and configurational entropy loss. The results provide relevant physics that can aid in designing inhibitors and understanding protein conformational equilibrium. Our approach is fast for use with proteins that contain flexible regions for structure-based drug design.

Best practices in free energy calculations for drug design

Free energy calculations are increasingly of interest for computing biophysical properties of novel small molecules of interest in drug design, such as protein-ligand binding affinities and small molecule partition coefficients. However, these calculations are also notoriously difficult to implement correctly. In this article, we review standard methods for computing free energy differences via simulation, discuss current best practices, and examine potential pitfalls for computational researchers without extensive experience in such calculations. We include a variety of examples and tips for how to set up and conduct these calculations, including applications to relative binding affinities and absolute binding free energies.

Constant pH Molecular Dynamics Simulations of Nucleic Acids in Explicit Solvent

The nucleosides of adenine and cytosine have pKa values of 3.50 and 4.08, respectively, and are assumed to be unprotonated under physiological conditions. However, evidence from recent NMR and X-Ray crystallography studies has revealed the prevalence of protonated adenine and cytosine in RNA macromolecules. Such nucleotides with elevated pKa values may play a role in stabilizing RNA structure and participate in the mechanism of ribozyme catalysis. With the work presented here, we establish the framework and demonstrate the first constant pH MD simulations (CPHMD) for nucleic acids in explicit solvent in which the protonation state is coupled to the dynamical evolution of the RNA system via λ-dynamics. We adopt the new functional form λ(Nexp) for λ that was recently developed for Multi-Site λ-Dynamics (MSλD) and demonstrate good sampling characteristics in which rapid and frequent transitions between the protonated and unprotonated states at pH = pKa are achieved. Our calculated pKa values of simple nucleotides are in a good agreement with experimentally measured values, with a mean absolute error of 0.24 pKa units. This work demonstrates that CPHMD can be used as a powerful tool to investigate pH-dependent biological properties of RNA macromolecules.

Applying efficient implicit nongeometric constraints in alchemical free energy simulations

Several strategies have been developed for satisfying bond lengths, angle, and other geometric constraints in molecular dynamics simulations. Advanced variations of alchemical free energy perturbation simulations, however, also require nongeometric constraints. In our recently developed multisite λ-dynamics simulation method, the conventional λ parameters that are associated with the progress variables in alchemical transformations are treated as dynamic variables and are constrained such that: 0 ≤ λ(i) ≤ 1 and Σ(i = 1)(N) λ(i) = 1. Here, we present four functional forms of λ that implicitly satisfy these nongeometric constraints, whose values and forces are facile to compute and that yield stable simulations using a 2 fs integration timestep. Using model systems, we present the sampling characteristics of these functional forms and demonstrate the enhanced sampling profiles and improved convergence rates that are achieved by the functional form: λ(i) = e(c sinθ(i))/Σ(j = 1)(N) e(c sinθ(j)) that oscillates between λ(i) = 0 and λ(i) = 1 and has relatively steep transitions between these endpoints.

λ-Meta Dynamics Approach To Compute Absolute Solvation Free Energy

We present a new approach to combine λ dynamics with meta-dynamics (named λ-meta dynamics) to compute free energy surface with respect to λ. Particularly, the λ-meta dynamics method extends meta-dynamics to a single virtual variable λ, i.e., the coupling parameter between solute and solvent, to compute absolute solvation free energy as an exemplary application. We demonstrate that λ-meta dynamics simulations can recover the accurate potential of mean force surface with respect to λ compared to the benchmark results from traditional λ-dynamics with umbrella sampling. The solvation free energy results for five small organic molecules from λ-meta dynamics simulations using the same filling scheme show that the statistical errors are within ±0.5 kcal/mol. The new λ-meta dynamics method is general and other variables such as order parameters to describe conformational changes can be easily combined with λ-meta dynamics. This should allow for efficient samplings on high-dimension free energy landscapes.

Multi-Site λ-dynamics for simulated Structure-Activity Relationship studies

Multi-Site λ-dynamics (MSλD) is a new free energy simulation method that is based on λ-dynamics. It has been developed to enable multiple substituents at multiple sites on a common ligand core to be modeled simultaneously and their free energies assessed. The efficacy of MSλD for estimating relative hydration free energies and relative binding affinties is demonstrated using three test systems. Model compounds representing multiple identical benzene, dihydroxybenzene and dimethoxybenzene molecules show total combined MSλD trajectory lengths of ~1.5 ns are sufficient to reliably achieve relative hydration free energy estimates within 0.2 kcal/mol and are less sensitive to the number of trajectories that are used to generate these estimates for hybrid ligands that contain up to ten substituents modeled at a single site or five substituents modeled at each of two sites. Relative hydration free energies among six benzene derivatives calculated from MSλD simulations are in very good agreement with those from alchemical free energy simulations (with average unsigned differences of 0.23 kcal/mol and R(2)=0.991) and experiment (with average unsigned errors of 1.8 kcal/mol and R(2)=0.959). Estimates of the relative binding affinities among 14 inhibitors of HIV-1 reverse transcriptase obtained from MSλD simulations are in reasonable agreement with those from traditional free energy simulations and experiment (average unsigned errors of 0.9 kcal/mol and R(2)=0.402). For the same level of accuracy and precision MSλD simulations are achieved ~20-50 times faster than traditional free energy simulations and thus with reliable force field parameters can be used effectively to screen tens to hundreds of compounds in structure-based drug design applications.

Path Integral Computation of Quantum Free Energy Differences Due to Alchemical Transformations Involving Mass and Potential

Thermodynamic integration, perturbation theory, and λ-dynamics methods were applied to path integral molecular dynamics calculations to investigate free energy differences due to "alchemical" transformations. Several estimators were formulated to compute free energy differences in solvable model systems undergoing changes in mass and/or potential. Linear and nonlinear alchemical interpolations were used for the thermodynamic integration. We find improved convergence for the virial estimators, as well as for the thermodynamic integration over nonlinear interpolation paths. Numerical results for the perturbative treatment of changes in mass and electric field strength in model systems are presented. We used thermodynamic integration in ab initio path integral molecular dynamics to compute the quantum free energy difference of the isotope transformation in the Zundel cation. The performance of different free energy methods is discussed.

Extending fragment-based free energy calculations with library Monte Carlo simulation: annealing in interaction space

Pre-calculated libraries of molecular fragment configurations have previously been used as a basis for both equilibrium sampling (via library-based Monte Carlo) and for obtaining absolute free energies using a polymer-growth formalism. Here, we combine the two approaches to extend the size of systems for which free energies can be calculated. We study a series of all-atom poly-alanine systems in a simple dielectric solvent and find that precise free energies can be obtained rapidly. For instance, for 12 residues, less than an hour of single-processor time is required. The combined approach is formally equivalent to the annealed importance sampling algorithm; instead of annealing by decreasing temperature, however, interactions among fragments are gradually added as the molecule is grown. We discuss implications for future binding affinity calculations in which a ligand is grown into a binding site.

Lambda-dynamics free energy simulation methods

Free energy calculations are fundamental to obtaining accurate theoretical estimates of many important biological phenomena including hydration energies, protein-ligand binding affinities and energetics of conformational changes. Unlike traditional free energy perturbation and thermodynamic integration methods, lambda-dynamics treats the conventional "lambda" as a dynamic variable in free energy simulations and simultaneously evaluates thermodynamic properties for multiple states in a single simulation. In the present article, we provide an overview of the theory of lambda-dynamics, including the use of biasing and restraining potentials to facilitate conformational sampling. We review how lambda-dynamics has been used to rapidly and reliably compute relative hydration free energies and binding affinities for series of ligands, to accurately identify crystallographically observed binding modes starting from incorrect orientations, and to model the effects of mutations upon protein stability. Finally, we suggest how lambda-dynamics may be extended to facilitate modeling efforts in structure-based drug design.

Replica exchange with dynamical scaling

A replica exchange method is presented which requires fewer replicas and is designed to be used for large systems. In this method, dynamically scaled replicas are placed between conventional replicas at broadly spaced temperatures. The potential of the scaled replicas is linearly scaled by a dynamical variable which varies between 0 and 1. When the variable is near either end point the replica can undergo exchanges with one of its neighboring replicas. Two different versions of the method are presented for a model system of a small peptide in water. The scaled replica can replace many replicas and the method can be up to ten times more efficient than conventional replica exchange.

Molecular grand-canonical ensemble density functional theory and exploration of chemical space

We present a rigorous description of chemical space within a molecular grand-canonical ensemble multi-component density functional theory framework. A total energy density functional for chemical compounds in contact with an electron and a proton bath is introduced using Lagrange multipliers which correspond to the energetic response to changes of the elementary particle densities. From a generalized Gibbs-Duhem equation analog, reactivity indices such as the nuclear hardness and a molecular Fukui function, which couples the grand-canonical electronic and nuclear degrees of freedom, are obtained. Maxwell relations between composition particles, ionic displacements, and the external potential are discussed. Numerical results for the molecular Fukui function are presented as well as finite temperature estimates for the oxidation of ammonia.

Efficient and precise solvation free energies via alchemical adiabatic molecular dynamics

A new molecular dynamics method for calculating free energies associated with transformations of the thermodynamic state or chemical composition of a system (also known as alchemical transformations) is presented. The new method extends the adiabatic dynamics approach recently introduced by Rosso et al. [J. Chem. Phys. 116, 4389 (2002)] and is based on the use of an additional degree of freedom, lambda, that is used as a switching parameter between the potential energy functions that characterize the two states. In the new method, the coupling parameter lambda is introduced as a fictitious dynamical variable in the Hamiltonian, and a system of switching functions is employed that leads to a barrier in the lambda free energy profile between the relevant thermodynamic end points. The presence of such a barrier, therefore, enhances sampling in the end point (lambda = 0 and lambda = 1) regions which are most important for computing relevant free energy differences. In order to ensure efficient barrier crossing, a high temperature T(lambda) is assigned to lambda and a fictitious mass m(lambda) is introduced as a means of creating an adiabatic separation between lambda and the rest of the system. Under these conditions, it is shown that the lambda free energy profile can be directly computed from the adiabatic probability distribution function of lambda without any postprocessing or unbiasing of the output data. The new method is illustrated on two model problems and in the calculation of the solvation free energy of amino acid side-chain analogs in TIP3P water. Comparisons to previous work using thermodynamic integration and free energy perturbation show that the new lambda adiabatic free energy dynamics method results in very precise free energy calculations using significantly shorter trajectories.

Conformational sampling via a self-regulating effective energy surface

The difficulty of efficiently sampling the phase space of complex systems with rough energy surfaces is well known. Typical solutions to the problem involve accelerating the crossing of barriers, but such methods often have the secondary problem that the low-energy states of interest are inadequately sampled, unless the parameters of the search algorithm are modified as the system evolves. A method is presented to improve the sampling with particular emphasis on the low-energy conformations, which make the most important contributions to the thermodynamics of the system. The algorithm proposed here samples the details of the minima, while easily surmounting barriers. This is achieved by introducing a self-regulating sampling variable which depends on the current state of the system. Two replicas of the system are introduced and the sampling variable is treated as a particle coupled to the physical system. The method is illustrated with a simple model system and is applied to the realistic example of barrier crossing in a protein-ligand complex.