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

Superimposing Ligands with a Ligand Overlay as an Alternate Topology Model for λ-Dynamics-Based Calculations

Alchemical free energy (AFE) calculations can predict binding affinity changes as a function of structural modifications and have become powerful tools for lead optimization and drug discovery. Central to the setup and performance of AFE calculations is the manner of mapping alchemical transformations, known as the topology model. Single, dual, and hybrid topology models have been used with various AFE methods in the field. In recent works, λ-dynamics (λD) free energy calculations, specifically, have preferred the use of a hybrid multiple topology (HMT) for sampling multiple ligand perturbations. In this work, we evaluate a new topology method called ligand overlay (LO) for use with λD-based calculations, including the recently introduced λ-dynamics with a bias-updated Gibbs sampling (LaDyBUGS) approach. LO is a full multiple topology model that allows entire ligands to be sampled and restrained within a λ-dynamics framework. Relative binding free energies were computed with HMT or LO topology models with LaDyBUGS for 45 ligands across five protein benchmark systems. An overall Pearson R correlation of 0.98 and mean unsigned error of 0.32 kcal/mol were observed, suggesting that LO is a viable alternative topology model for λD-based calculations. We discuss the merits of using an HMT or LO model for future ligand studies with λD or LaDyBUGS calculations.

CHARMM at 45: Enhancements in Accessibility, Functionality, and Speed

Since its inception nearly a half century ago, CHARMM has been playing a central role in computational biochemistry and biophysics. Commensurate with the developments in experimental research and advances in computer hardware, the range of methods and applicability of CHARMM have also grown. This review summarizes major developments that occurred after 2009 when the last review of CHARMM was published. They include the following: new faster simulation engines, accessible user interfaces for convenient workflows, and a vast array of simulation and analysis methods that encompass quantum mechanical, atomistic, and coarse-grained levels, as well as extensive coverage of force fields. In addition to providing the current snapshot of the CHARMM development, this review may serve as a starting point for exploring relevant theories and computational methods for tackling contemporary and emerging problems in biomolecular systems. CHARMM is freely available for academic and nonprofit research at https://academiccharmm.org/program.

Alchemical Transformations and Beyond: Recent Advances and Real-World Applications of Free Energy Calculations in Drug Discovery

Computational methods constitute efficient strategies for screening and optimizing potential drug molecules. A critical factor in this process is the binding affinity between candidate molecules and targets, quantified as binding free energy. Among various estimation methods, alchemical transformation methods stand out for their theoretical rigor. Despite challenges in force field accuracy and sampling efficiency, advancements in algorithms, software, and hardware have increased the application of free energy perturbation (FEP) calculations in the pharmaceutical industry. Here, we review the practical applications of FEP in drug discovery projects since 2018, covering both ligand-centric and residue-centric transformations. We show that relative binding free energy calculations have steadily achieved chemical accuracy in real-world applications. In addition, we discuss alternative physics-based simulation methods and the incorporation of deep learning into free energy calculations.

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.

A Membrane-Assisted Mechanism for the Release of Ceramide from the CERT START Domain

Ceramide transfer protein CERT is the mediator of nonvesicular transfer of ceramide from the ER to Golgi. In CERT, START is the domain responsible for the binding and transport of ceramide. A wealth of structural data has revealed a helix-grip fold surrounding a large hydrophobic cavity holding the ceramide. Yet, little is known about the mechanisms by which START releases the ceramide through the polar region and into the packed environment of cellular membranes. As such events do not lend themselves easily to experimental investigations, we used multiple unbiased microsecond-long molecular simulations. We propose a membrane-assisted mechanism in which the membrane acts as an allosteric effector initiating the release of ceramide and where the passage of the ceramide acyl chains is facilitated by the intercalation of a single phosphatidylcholine lipid in the cavity, practically greasing the ceramide way out. We verify using free energy calculation and experimental lipidomics data that CERT forms stable complexes with phosphatidylcholine lipids, in addition to ceramide, thus providing validation for the proposed mechanism.

Exploring the Limits of the Generalized CHARMM and AMBER Force Fields through Predictions of Hydration Free Energy of Small Molecules

Accurate force field parameters, potential energy functions, and receptor-ligand models are essential for modeling the solvation and binding of drug-like molecules to a receptor. A large and ever-growing chemical space of medicinally relevant scaffolds has also required these factors, especially force field parameters, to be highly transferable. Generalized force fields such as the CHARMM General Force Field (CGenFF) and the generalized AMBER force field (GAFF) have accomplished this feat along with other contemporaneous ones like OPLS. Here, we analyze the limits in the parametrization of drug-like small molecules by CGenFF and GAFF in terms of the various functional groups represented within them. Specifically, we link the presence of specific functional groups to the error in the absolute hydration free energy of over 600 small molecules, predicted by alchemical free energy methods implemented in the CHARMM program. Our investigation reveals that molecules with (i) a nitro group in CGenFF and GAFF are, respectively, over- or undersolubilized in aqueous medium, (ii) amine groups are undersolubilized more so in CGenFF than in GAFF, and (iii) carboxyl groups are more oversolubilized in GAFF than in CGenFF. We present our analyses of the potential factors underlying these trends. We also showcase the use of a machine-learning-based approach combined with the SHapley Additive exPlanations framework to attribute these trends to specific functional groups, which can be easily adopted to explore the limits of other general force fields.

High Charge Density in Peptide Dendrimers is Required to Destabilize Membranes: Insights into Endosome Evasion

Peptide dendrimers are a type of branched, symmetric, and topologically well-defined molecule that have already been used as delivery systems for nucleic acid transfection. Several of the most promising sequences showed high efficiency in many key steps of transfection, namely, binding siRNA, entering cells, and evading the endosome. However, small changes to the peptide dendrimers, such as in the hydrophobic core, the amino acid chirality, or the total available charges, led to significantly different experimental results with unclear mechanistic insights. In this work, we built a computational model of several of those peptide dendrimers (MH18, MH13, and MH47) and some of their variants to study the molecular details of the structure and function of these molecules. We performed CpHMD simulations in the aqueous phase and in interaction with a lipid bilayer to assess how conformation and protonation are affected by pH in different environments. We found that while the different peptide dendrimer sequences lead to no substantial structural differences in the aqueous phase, the total charge and, more importantly, the total charge density are key for the capacity of the dendrimer to interact and destabilize the membrane. These dendrimers become highly charged when the pH changes from 7.5 to 4.5, and the presence of a high charge density, which is decreased for MH47 that has four fewer titratable lysines, is essential to trigger membrane destabilization. These findings are in excellent agreement with the experimental data and help us to understand the high efficiency of some dendrimers and why the dendrimer MH47 is unable to complete the transfection process. This evidence provides further understanding of the mode of action of these peptide dendrimers and will be pivotal for the future design of new sequences with improved transfection capabilities.

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.

Enhanced Sampling of Buried Charges in Free Energy Calculations Using Replica Exchange with Charge Tempering

Buried ionizable groups in proteins often play important structural and functional roles. However, it is generally challenging to study the detailed molecular mechanisms solely based on experimental measurements. Free energy calculations using atomistic simulations, on the other hand, complement experimental studies and can provide high temporal and spatial resolution information that can lead to mechanistic insights. Nevertheless, it is also well recognized that sufficient sampling of such atomistic simulations can be challenging, considering that structural changes related to the buried charges may be very slow. In the present study, we describe a simple but effective enhanced sampling technique called replica exchange with charge tempering (REChgT) with a novel free energy method, multisite λ dynamics (MSλD), to study two systems containing buried charges, pKa prediction of a small molecule, orotate, in complex with the dihydroorotate dehydrogenase, and relative stability of a Glu-Lys pair buried in the hydrophobic core of two variants of Staphylococcal nuclease. Compared to the original MSλD simulations, the usage of REChgT dramatically increases sampling in both conformational and alchemical spaces, which directly translates into a significant reduction of wall time to converge the free energy calculations. This study highlights the importance of sufficient sampling toward developing improved free energy methods.

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.

Selection pressures on evolution of ribonuclease H explored with rigorous free-energy-based design

Understanding natural protein evolution and designing novel proteins are motivating interest in development of high-throughput methods to explore large sequence spaces. In this work, we demonstrate the application of multisite λ dynamics (MSλD), a rigorous free energy simulation method, and chemical denaturation experiments to quantify evolutionary selection pressure from sequence-stability relationships and to address questions of design. This study examines a mesophilic phylogenetic clade of ribonuclease H (RNase H), furthering its extensive characterization in earlier studies, focusing on E. coli RNase H (ecRNH) and a more stable consensus sequence (AncCcons) differing at 15 positions. The stabilities of 32,768 chimeras between these two sequences were computed using the MSλD framework. The most stable and least stable chimeras were predicted and tested along with several other sequences, revealing a designed chimera with approximately the same stability increase as AncCcons, but requiring only half the mutations. Comparing the computed stabilities with experiment for 12 sequences reveals a Pearson correlation of 0.86 and root mean squared error of 1.18 kcal/mol, an unprecedented level of accuracy well beyond less rigorous computational design methods. We then quantified selection pressure using a simple evolutionary model in which sequences are selected according to the Boltzmann factor of their stability. Selection temperatures from 110 to 168 K are estimated in three ways by comparing experimental and computational results to evolutionary models. These estimates indicate selection pressure is high, which has implications for evolutionary dynamics and for the accuracy required for design, and suggests accurate high-throughput computational methods like MSλD may enable more effective protein design.

Roles of Accelerated Molecular Dynamics Simulations in Predictions of Binding Kinetic Parameters

Rational predictions on binding kinetics parameters of drugs to targets play significant roles in future drug designs. Full conformational samplings of targets are requisite for accurate predictions of binding kinetic parameters. In this review, we mainly focus on the applications of enhanced sampling technologies in calculations of binding kinetics parameters and residence time of drugs. The methods involved in molecular dynamics simulations are applied to not only probe conformational changes of targets but also reveal calculations of residence time that is significant for drug efficiency. For this review, special attention are paid to accelerated molecular dynamics (aMD) and Gaussian aMD (GaMD) simulations that have been adopted to predict the association or disassociation rate constant. We also expect that this review can provide useful information for future drug design.

Fast free energy estimates from λ-dynamics with bias-updated Gibbs sampling

Relative binding free energy calculations have become an integral computational tool for lead optimization in structure-based drug design. Classical alchemical methods, including free energy perturbation or thermodynamic integration, compute relative free energy differences by transforming one molecule into another. However, these methods have high operational costs due to the need to perform many pairwise perturbations independently. To reduce costs and accelerate molecular design workflows, we present a method called λ-dynamics with bias-updated Gibbs sampling. This method uses dynamic biases to continuously sample between multiple ligand analogues collectively within a single simulation. We show that many relative binding free energies can be determined quickly with this approach without compromising accuracy. For five benchmark systems, agreement to experiment is high, with root mean square errors near or below 1.0 kcal mol-1. Free energy results are consistent with other computational approaches and within statistical noise of both methods (0.4 kcal mol-1 or less). Notably, large efficiency gains over thermodynamic integration of 18-66-fold for small perturbations and 100-200-fold for whole aromatic ring substitutions are observed. The rapid determination of relative binding free energies will enable larger chemical spaces to be more readily explored and structure-based drug design to be accelerated.

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

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

Enhanced Configurational Sampling Approaches to Alchemical Ligand Binding Free Energy Simulations: Current Status and Challenges

Ligand binding free energy simulations (LB-FES) have been routine tasks in modern drug discovery campaign. A long-standing challenge for LB-FES is the difficulty in adequately sampling nontrivial environmental reorganizations in response to ligand binding. Therefore, various enhanced configurational sampling (ECS) approaches were devised to speed up fluctuations of relevant slow degrees of freedom (SDOF) and ensure simulation convergence. However, in contrast to the achievements in parametrization, software performance, and workflow automation, efficient ECS methodology suitable for high throughput screening remains in an early stage of development. Here, a review of ECS developments with LB-FES is presented, revisiting current approaches and underlining the major technical pitfalls and challenges. This Perspective focuses on alchemical LB-FES on account of their predominant role in high throughput drug screening as well as the established partnership with ECS. The critical aspects of designing ECS approaches, from both theoretical and applied perspectives, are described. This work is intended to provide a contemporary review of the scientific, technical, and practical issues associated with the accelerating convergence of alchemical LB-FES.

Molecular Simulations of Conformational Transitions within the Insulin Receptor Kinase Reveal Consensus Features in a Multistep Activation Pathway

Modulating the transitions between active and inactive conformations of protein kinases is the primary means of regulating their catalytic activity, achieved by phosphorylation of the activation loop (A-loop). To elucidate the mechanism of this conformational activation, we applied the string method to determine the conformational transition path of insulin receptor kinase between the active and inactive conformations and the corresponding free-energy profiles with and without A-loop phosphorylation. The conformational change was found to proceed in three sequential steps: first, the flipping of the DFG motif of the active site; second, rotation of the A-loop; finally, the inward movement of the αC helix. The main energetic bottleneck corresponds to the conformational change in the A-loop, while changes in the DFG motif and αC helix occur before and after A-loop conformational change, respectively. In accordance with this, two intermediate states are identified, the first state just after the DFG flipping and the second state after the A-loop rotation. These intermediates exhibit structural features characteristic of the corresponding inactive and active conformations of other protein kinases. To understand the impact of A-loop phosphorylation on kinase conformation, the free energies of A-loop phosphorylation were determined at several states along the conformational transition path using the free-energy perturbation simulations. The calculated free energies reveal that while the unphosphorylated kinase interconverts between the inactive and active conformations, A-loop phosphorylation restricts access to the inactive conformation, thereby increasing the active conformation population. Overall, this study suggests a consensus mechanism of conformational activation between different protein kinases.

pyCHARMM: Embedding CHARMM Functionality in a Python Framework

CHARMM is rich in methodology and functionality as one of the first programs addressing problems of molecular dynamics and modeling of biological macromolecules and their partners, e.g., small molecule ligands. When combined with the highly developed CHARMM parameters for proteins, nucleic acids, small molecules, lipids, sugars, and other biologically relevant building blocks, and the versatile CHARMM scripting language, CHARMM has been a trendsetting platform for modeling studies of biological macromolecules. To further enhance the utility of accessing and using CHARMM functionality in increasingly complex workflows associated with modeling biological systems, we introduce pyCHARMM, Python bindings, functions, and modules to complement and extend the extensive set of modeling tools and methods already available in CHARMM. These include access to CHARMM function-generated variables associated with the system (psf), coordinates, velocities and forces, atom selection variables, and force field related parameters. The ability to augment CHARMM forces and energies with energy terms or methods derived from machine learning or other sources, written in Python, CUDA, or OpenCL and expressed as Python callable routines is introduced together with analogous functions callable during dynamics calculations. Integration of Python-based graphical engines for visualization of simulation models and results is also accessible. Loosely coupled parallelism is available for workflows such as free energy calculations, using MBAR/TI approaches or high-throughput multisite λ-dynamics (MSλD) free energy methods, string path optimization calculations, replica exchange, and molecular docking with a new Python-based CDOCKER module. CHARMM accelerated platform kernels through the CHARMM/OpenMM API, CHARMM/DOMDEC, and CHARMM/BLaDE API are also readily integrated into this Python framework. We anticipate that pyCHARMM will be a robust platform for the development of comprehensive and complex workflows utilizing Python and its extensive functionality as well as an optimal platform for users to learn molecular modeling methods and practices within a Python-friendly environment such as Jupyter Notebooks.

Identification of non-histone substrates of the lysine methyltransferase PRDM9

Lysine methylation is a dynamic, post-translational mark that regulates the function of histone and non-histone proteins. Many of the enzymes that mediate lysine methylation, known as lysine methyltransferases (KMTs), were originally identified to modify histone proteins but have also been discovered to methylate non-histone proteins. In this work, we investigate the substrate selectivity of the lysine methyltransferase PRDM9 to identify both potential histone and non-histone substrates. Though normally expressed in germ cells, PRDM9 is significantly upregulated across many cancer types. The methyltransferase activity of PRDM9 is essential for double-strand break formation during meiotic recombination. PRDM9 has been reported to methylate histone H3 at lysine residues 4 and 36; however, PRDM9 KMT activity had not previously been evaluated on non-histone proteins. Using lysine-oriented peptide (K-OPL) libraries to screen potential substrates of PRDM9, we determined that PRDM9 preferentially methylates peptide sequences not found in any histone protein. We confirmed PRDM9 selectivity through in vitro KMT reactions using peptides with substitutions at critical positions. A multisite λ-dynamics computational analysis provided a structural rationale for the observed PRDM9 selectivity. The substrate selectivity profile was then used to identify putative non-histone substrates, which were tested by peptide spot array, and a subset were further validated at the protein level by in vitro KMT assays on recombinant proteins. Finally, one of the non-histone substrates, CTNNBL1, was found to be methylated by PRDM9 in cells.

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.

GPU-Accelerated All-Atom Particle-Mesh Ewald Continuous Constant pH Molecular Dynamics in Amber

Constant pH molecular dynamics (MD) simulations sample protonation states on the fly according to the conformational environment and user specified pH conditions; however, the current accuracy is limited due to the use of implicit-solvent models or a hybrid solvent scheme. Here, we report the first GPU-accelerated implementation, parametrization, and validation of the all-atom continuous constant pH MD (CpHMD) method with particle-mesh Ewald (PME) electrostatics in the Amber22 pmemd.cuda engine. The titration parameters for Asp, Glu, His, Cys, and Lys were derived for the CHARMM c22 and Amber ff14sb and ff19sb force fields. We then evaluated the PME-CpHMD method using the asynchronous pH replica-exchange titration simulations with the c22 force field for six benchmark proteins, including BBL, hen egg white lysozyme (HEWL), staphylococcal nuclease (SNase), thioredoxin, ribonuclease A (RNaseA), and human muscle creatine kinase (HMCK). The root-mean-square deviation from the experimental pK_a's of Asp, Glu, His, and Cys is 0.76 pH units, and the Pearson's correlation coefficient for the pK_a shifts with respect to model values is 0.80. We demonstrated that a finite-size correction or much enlarged simulation box size can remove a systematic error of the calculated pK_a's and improve agreement with experiment. Importantly, the simulations captured the relevant biology in several challenging cases, e.g., the titration order of the catalytic dyad Glu35/Asp52 in HEWL and the coupled residues Asp19/Asp21 in SNase, the large pK_a upshift of the deeply buried catalytic Asp26 in thioredoxin, and the large pK_a downshift of the deeply buried catalytic Cys283 in HMCK. We anticipate that PME-CpHMD will offer proper pH control to improve the accuracies of MD simulations and enable mechanistic studies of proton-coupled dynamical processes that are ubiquitous in biology but remain poorly understood due to the lack of experimental tools and limitation of current MD simulations.

Constant pH molecular dynamics simulations: Current status and recent applications

Many important protein functions are carried out through proton-coupled conformational dynamics. Thus, the ability to accurately model protonation states dynamically has wide-ranging implications. Over the past two decades, two main types of constant pH methods (discrete and continuous) have been developed to enable proton-coupled molecular dynamics (MD) simulations. In this short review, we discuss the current status of the development and highlight recent applications that have advanced our understanding of protein structure-function relationships. We conclude the review by outlining the remaining challenges in the method development and projecting important areas for future applications.

Insights into the evolution of enzymatic specificity and catalysis: From Asgard archaea to human adenylate kinases

Enzymatic catalysis is critically dependent on selectivity, active site architecture, and dynamics. To contribute insights into the interplay of these properties, we established an approach with NMR, crystallography, and MD simulations focused on the ubiquitous phosphotransferase adenylate kinase (AK) isolated from Odinarchaeota (OdinAK). Odinarchaeota belongs to the Asgard archaeal phylum that is believed to be the closest known ancestor to eukaryotes. We show that OdinAK is a hyperthermophilic trimer that, contrary to other AK family members, can use all NTPs for its phosphorylation reaction. Crystallographic structures of OdinAK-NTP complexes revealed a universal NTP-binding motif, while 19F NMR experiments uncovered a conserved and rate-limiting dynamic signature. As a consequence of trimerization, the active site of OdinAK was found to be lacking a critical catalytic residue and is therefore considered to be "atypical." On the basis of discovered relationships with human monomeric homologs, our findings are discussed in terms of evolution of enzymatic substrate specificity and cold adaptation.

A Guide to the Continuous Constant pH Molecular Dynamics Methods in Amber and CHARMM [Article v1.0]

Like temperature and pressure, solution pH is an important environmental variable in biomolecular simulations. Virtually all proteins depend on pH to maintain their structure and function. In conventional molecular dynamics (MD) simulations of proteins, pH is implicitly accounted for by assigning and fixing protonation states of titratable sidechains. This is a significant limitation, as the assigned protonation states may be wrong and they may change during dynamics. In this tutorial, we guide the reader in learning and using the various continuous constant pH MD methods in Amber and CHARMM packages, which have been applied to predict pK a values and elucidate proton-coupled conformational dynamics of a variety of proteins including enzymes and membrane transporters.

Leveraging the Sampling Efficiency of RE-EDS in OpenMM Using a Shifted Reaction-Field With an Atom-Based Cutoff

Replica-exchange enveloping distribution sampling (RE-EDS) is a pathway-independent multistate free-energy method, currently implemented in the GROMOS software package for molecular dynamics (MD) simulations. It has a high intrinsic sampling efficiency as the interactions between the unperturbed particles have to be calculated only once for multiple end-states. As a result, RE-EDS is an attractive method for the calculation of relative solvation and binding free energies. An essential requirement for reaching this high efficiency is the separability of the nonbonded interactions into solute-solute, solute-environment, and environment-environment contributions. Such a partitioning is trivial when using a Coulomb term with a reaction-field (RF) correction to model the electrostatic interactions, but not when using lattice- sum schemes. To avoid cutoff artifacts, the RF correction is typically used in combination with a charge-group based cutoff, which is not supported by most small-molecule force fields and other MD engines. To address this issue, we investigate the combination of RE-EDS simulations with a recently introduced RF scheme including a shifting function that enables the rigorous calculation of RF electrostatics with atom-based cutoffs. The resulting approach is validated by calculating solvation free energies with the generalized AMBER force field (GAFF) in water and chloroform using both the GROMOS software package and a proof-of-concept implementation in OpenMM.

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

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

Molecular Dynamics Simulations Establish the Molecular Basis for the Broad Allostery Hotspot Distributions in the Tetracycline Repressor

It is imperative to identify the network of residues essential to the allosteric coupling for the purpose of rationally engineering allostery in proteins. Deep mutational scanning analysis has emerged as a function-centric approach for identifying such allostery hotspots in a comprehensive and unbiased fashion, leading to observations that challenge our understanding of allostery at the molecular level. Specifically, a recent deep mutational scanning study of the tetracycline repressor (TetR) revealed an unexpectedly broad distribution of allostery hotspots throughout the protein structure. Using extensive molecular dynamics simulations (up to 50 μs) and free energy computations, we establish the molecular and energetic basis for the strong anticooperativity between the ligand and DNA binding sites. The computed free energy landscapes in different ligation states illustrate that allostery in TetR is well described by a conformational selection model, in which the apo state samples a broad set of conformations, and specific ones are selectively stabilized by either ligand or DNA binding. By examining a range of structural and dynamic properties of residues at both local and global scales, we observe that various analyses capture different subsets of experimentally identified hotspots, suggesting that these residues modulate allostery in distinct ways. These results motivate the development of a thermodynamic model that qualitatively explains the broad distribution of hotspot residues and their distinct features in molecular dynamics simulations. The multifaceted strategy that we establish here for hotspot evaluations and our insights into their mechanistic contributions are useful for modulating protein allostery in mechanistic and engineering studies.

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.