Calculating the binding free energies of charged species based on explicit-solvent simulations employing lattice-sum methods: an accurate correction scheme for electrostatic finite-size effects

BAT2: an Open-Source Tool for Flexible, Automated, and Low Cost Absolute Binding Free Energy Calculations

Absolute binding free energy (ABFE) calculations with all-atom molecular dynamics (MD) have the potential to greatly reduce costs in the first stages of drug discovery. Here, we introduce BAT2, the new version of the Binding Affinity Tool (BAT.py), designed to combine full automation of ABFE calculations with high-performance MD simulations, making it a potential tool for virtual screening. We describe and test several changes and new features that were incorporated into the code, such as relative restraints between the protein and the ligand instead of using fixed dummy atoms, support for the OpenMM simulation engine, a merged approach to the application/release of restraints, support for cobinders and proteins with multiple chains, and many others. We also reduced the simulation times for each ABFE calculation, assessing the effect on the expected robustness and accuracy of the calculations.

A Comprehensive Study of Factors Affecting the Prediction of the pKa Shift of Asp26 in Thioredoxin Protein

The stable protonation state of ionizable amino acids in a protein can be predicted by computing the pKa shift of that residue within the protein environment. Thermodynamic Integration (TI) is an ideal molecular dynamics-based approach for predicting the pKa shift of ionizable protein residues. Here, we probe TI-based simulation protocols for their ability to accurately predict the pKa shift of Asp26 in thioredoxin. While implicit solvent models can predict the pKa shift accurately, explicit solvent models result in substantial errors. To understand the underlying reason for this surprising discrepancy, we investigate the role of various factors such as solvent models, conformational sampling, background charges, and polarization.

Methodological and force field effects in the molecular dynamics-based prediction of binding free energies of host-guest systems

As a contribution to the understanding and rationalization of methodological and modeling effects in recent host-guest SAMPL challenges, using an alchemical molecular dynamics technique we have examined the impact of force field parameterization and ionic strength in connection with guest charge neutralization on computed dissociation free energies in two typical SAMPL heavily charged macrocyclic hosts encapsulating small protonated amines with disparate binding affinities. We have shown that the methodological treatment for host neutralization, with explicit ions or with the background neutralizing plasma in the context of alchemical calculations under periodic boundary conditions, has a moderate effect on the calculated affinities. On the other hand, we have shown that seemingly small differences in the force field parameterization in highly symmetric hosts can produce systematic effects on the structural features that can have a significant impact on the predicted binding affinities.

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.

Interactions of a Novel Anthracycline with Oligonucleotide DNA and Cyclodextrins in an Aqueous Environment

Berubicin, a chemotherapy medication belonging to the class of anthracyclines, is simulated in double-stranded DNA sequences and cyclodextrins in an aqueous environment via full-atom molecular dynamics simulations on the time scale of microseconds. The drug is studied in both the neutral and protonated states so as to better comprehend the role of its charge in the formed complexes. The noncovalent berubicin-DNA and berubicin-cyclodextrin complexes are investigated in detail, paying special attention to their thermodynamic description by employing the double decoupling method, the solvent balance method, the weighted solvent accessible surface model, and the linear interaction energy method. A novel approach for extracting the desolvation thermodynamics of the binding process is also presented. Both the binding and desolvation Gibbs energies are decomposed into entropic and enthalpic contributions so as to elucidate the nature of complexation and its driving forces. Selected structural and geometrical properties of all the complexes, which are all stable, are analyzed. Both cyclodextrins under consideration are widely utilized for drug delivery purposes, and a comparative investigation between their bound states with berubicin is carried out.

Multistate Method to Efficiently Account for Tautomerism and Protonation in Alchemical Free-Energy Calculations

The majority of drug-like molecules contain at least one ionizable group, and many common drug scaffolds are subject to tautomeric equilibria. Thus, these compounds are found in a mixture of protonation and/or tautomeric states at physiological pH. Intrinsically, standard classical molecular dynamics (MD) simulations cannot describe such equilibria between states, which negatively impacts the prediction of key molecular properties in silico. Following the formalism described by de Oliveira and co-workers (J. Chem. Theory Comput. 2019, 15, 424-435) to consider the influence of all states on the binding process based on alchemical free-energy calculations, we demonstrate in this work that the multistate method replica-exchange enveloping distribution sampling (RE-EDS) is well suited to describe molecules with multiple protonation and/or tautomeric states in a single simulation. We apply our methodology to a series of eight inhibitors of factor Xa with two protonation states and a series of eight inhibitors of glycogen synthase kinase 3β (GSK3β) with two tautomeric states. In particular, we show that given a sufficient phase-space overlap between the states, RE-EDS is computationally more efficient than standard pairwise free-energy methods.

Binding Sites of Bicarbonate in Phosphoenolpyruvate Carboxylase

Phosphoenolpyruvate carboxylase (PEPC) is used in plant metabolism for fruit maturation or seed development as well as in the C4 and crassulacean acid metabolism (CAM) mechanisms in photosynthesis, where it is used for the capture of hydrated CO2 (bicarbonate). To find the yet unknown binding site of bicarbonate in this enzyme, we have first identified putative binding sites with nonequilibrium molecular dynamics simulations and then ranked these sites with alchemical free energy calculations with corrections of computational artifacts. Fourteen pockets where bicarbonate could bind were identified, with three having realistic binding free energies with differences with the experimental value below 1 kcal/mol. One of these pockets is found far from the active site at 14 Å and predicted to be an allosteric binding site. In the two other binding sites, bicarbonate is in direct interaction with the magnesium ion; neither sequence alignment nor the study of mutant K606N allowed to discriminate between these two pockets, and both are good candidates as the binding site of bicarbonate in phosphoenolpyruvate carboxylase.

Quantitative Assessment of Energetic Contributions of Residues in a SARS-CoV-2 Viral Enzyme/Nanobody Interface

The highly conserved protease enzyme from SARS-CoV-2 (MPro) is crucial for viral replication and is an attractive target for the design of novel inhibitory compounds. MPro is known to be conformationally flexible and has been stabilized in an extended conformation in a complex with a novel nanobody (NB2B4), which inhibits the dimerization of the enzyme via binding to an allosteric site. However, the energetic contributions of the nanobody residues stabilizing the MPro/nanobody interface remain unresolved. We probed these residues using all-atom MD simulations in combination with alchemical free energy calculations by studying the physical residue-residue interactions and discovered the role of hydrophobic and electrostatic interactions in stabilizing the complex. Specifically, we found via mutational analysis that three interfacial nanobody residues (Y59, R106, and L109) contributed significantly, two residues (L107 and P110) contributed moderately, and two residues (H112 and T113) contributed minimally to the overall binding affinity of the nanobody. We also discovered that the nanobody affinity could be enhanced via a charge-reversal mutation (D62R) that alters the local interfacial electrostatic environment of this residue in the complex. These findings are potentially useful in designing novel synthetic nanobodies as allosteric inhibitors of MPro.

Orientation and Membrane Partition Free Energy of PeT-Based Voltage-Sensitive Dyes from Molecular Simulations

Voltage measurement via small-molecule fluorescent indicators is a valuable approach in deciphering complex dynamics in electrically excitable cells. However, our understanding of various physicochemical properties governing the performance of fluorescent voltage sensors based on the photoinduced electron transfer (PeT) mechanism remains incomplete. Here, through extensive molecular dynamics and free energy calculations, we systematically examine the orientation and membrane partition of three PeT-based voltage-sensing VoltageFluor (VF) dyes in different lipid environment. We show that the symmetry of the molecular scaffold and the net charge of the hydrophilic headgroup of a given VF dye dominate its orientation and membrane partition, respectively. Our work provides a mechanistic understanding of the physical properties contributing to the voltage sensitivity, signal-to-noise ratio, as well as membrane distribution of VF dyes and sheds light onto rational design principles of PeT-based fluorescent probes in general.

An Efficient Approach to the Accurate Prediction of Mutational Effects in Antigen Binding to the MHC1

The major histocompatibility complex (MHC) can recognize and bind to external peptides to generate effective immune responses by presenting the peptides to T cells. Therefore, understanding the binding modes of peptide-MHC complexes (pMHC) and predicting the binding affinity of pMHCs play a crucial role in the rational design of peptide vaccines. In this study, we employed molecular dynamics (MD) simulations and free energy calculations with an Alanine Scanning with Generalized Born and Interaction Entropy (ASGBIE) method to investigate the protein-peptide interaction between HLA-A*02:01 and the G9209 peptide derived from the melanoma antigen gp100. The energy contribution of individual residue was calculated using alanine scanning, and hotspots on both the MHC and the peptides were identified. Our study shows that the pMHC binding is dominated by the van der Waals interactions. Furthermore, we optimized the ASGBIE method, achieving a Pearson correlation coefficient of 0.91 between predicted and experimental binding affinity for mutated antigens. This represents a significant improvement over the conventional MM/GBSA method, which yields a Pearson correlation coefficient of 0.22. The computational protocol developed in this study can be applied to the computational screening of antigens for the MHC1 as well as other protein-peptide binding systems.

What to Make of Zero: Resolving the Statistical Noise from Conformational Reorganization in Alchemical Binding Free Energy Estimates with Metadynamics Sampling

We introduce the self-relative binding free energy (self-RBFE) approach to evaluate the intrinsic statistical variance of dual-topology alchemical binding free energy estimators. The self-RBFE is the relative binding free energy between a ligand and a copy of the same ligand, and its true value is zero. Nevertheless, because the two copies of the ligand move independently, the self-RBFE value produced by a finite-length simulation fluctuates and can be used to measure the variance of the model. The results of this validation provide evidence that a significant fraction of the errors observed in benchmark studies reflect the statistical fluctuations of unconverged estimates rather than the models' accuracy. Furthermore, we find that ligand reorganization is a significant contributing factor to the statistical variance of binding free energy estimates and that metadynamics-accelerated conformational sampling of the torsional degrees of freedom of the ligand can drastically reduce the time to convergence.

Solvent effects on extractant conformational energetics in liquid-liquid extraction: a simulation study of molecular solvents and ionic liquids

Extractant design in liquid-liquid extraction (LLE) is a research frontier of metal ion separations that typically focuses on the direct extractant-metal interactions. However, a more detailed understanding of energetic drivers of separations beyond primary metal coordination is often lacking, including the role of solvent in the extractant phase. In this work, we propose a new mechanism for enhancing metal-complexant energetics with nanostructured solvents. Using molecular dynamics simulations with umbrella sampling, we find that the organic solvent can reshape the energetics of the extractant's intramolecular conformational landscape. We calculate free energy profiles of different conformations of a representative bidentate extractant, n-octyl(phenyl)-N,N-diisobutyl carbamoyl methyl phosphinoxide (CMPO), in four different solvents: dodecane, tributyl phosphate (TBP), and dry and wet ionic liquid (IL) 1-ethyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide ([EMIM][Tf2N]). By promoting reorganization of the extractant molecule into its binding conformation, our findings reveal how particular solvents can ameliorate this unfavorable step of the metal separation process. In particular, the charge alternating nanodomains formed in ILs substantially reduce the free energy penalty associated with extractant reorganization. Importantly, using alchemical free energy calculations, we find that this stabilization persists even when we explicitly include the extracted cation. These findings provide insight into the energetic drivers of metal ion separations and potentially suggest a new approach to designing effective separations using a molecular-level understanding of solvent effects.

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.

Performance and Analysis of the Alchemical Transfer Method for Binding-Free-Energy Predictions of Diverse Ligands

The Alchemical Transfer Method (ATM) is herein validated against the relative binding-free energies (RBFEs) of a diverse set of protein-ligand complexes. We employed a streamlined setup workflow, a bespoke force field, and AToM-OpenMM software to compute the RBFEs of the benchmark set prepared by Schindler and collaborators at Merck KGaA. This benchmark set includes examples of standard small R-group ligand modifications as well as more challenging scenarios, such as large R-group changes, scaffold hopping, formal charge changes, and charge-shifting transformations. The novel coordinate perturbation scheme and a dual-topology approach of ATM address some of the challenges of single-topology alchemical RBFE methods. Specifically, ATM eliminates the need for splitting electrostatic and Lennard-Jones interactions, atom mapping, defining ligand regions, and postcorrections for charge-changing perturbations. Thus, ATM is simpler and more broadly applicable than conventional alchemical methods, especially for scaffold-hopping and charge-changing transformations. Here, we performed well over 500 RBFE calculations for eight protein targets and found that ATM achieves accuracy comparable to that of existing state-of-the-art methods, albeit with larger statistical fluctuations. We discuss insights into the specific strengths and weaknesses of the ATM method that will inform future deployments. This study confirms that ATM can be applied as a production tool for RBFE predictions across a wide range of perturbation types within a unified, open-source framework.

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.

Is the Functional Response of a Receptor Determined by the Thermodynamics of Ligand Binding?

For an effective drug, strong binding to the target protein is a prerequisite, but it is not enough. To produce a particular functional response, drugs need to either block the proteins' functions or modulate their activities by changing their conformational equilibrium. The binding free energy of a compound to its target is routinely calculated, but the timescales for the protein conformational changes are prohibitively long to be efficiently modeled via physics-based simulations. Thermodynamic principles suggest that the binding free energies of the ligands with different receptor conformations may infer their efficacy. However, this hypothesis has not been thoroughly validated. We present an actionable protocol and a comprehensive study to show that binding thermodynamics provides a strong predictor of the efficacy of a ligand. We apply the absolute binding free energy perturbation method to ligands bound to active and inactive states of eight G protein-coupled receptors and a nuclear receptor and then compare the resulting binding free energies. We find that carefully designed restraints are often necessary to efficiently model the corresponding conformational ensembles for each state. Our method achieves unprecedented performance in classifying ligands as agonists or antagonists across the various investigated receptors, all of which are important drug targets.

Alchemical Free Energy Calculations on Membrane-Associated Proteins

Membrane proteins have diverse functions within cells and are well-established drug targets. The advances in membrane protein structural biology have revealed drug and lipid binding sites on membrane proteins, while computational methods such as molecular simulations can resolve the thermodynamic basis of these interactions. Particularly, alchemical free energy calculations have shown promise in the calculation of reliable and reproducible binding free energies of protein-ligand and protein-lipid complexes in membrane-associated systems. In this review, we present an overview of representative alchemical free energy studies on G-protein-coupled receptors, ion channels, transporters as well as protein-lipid interactions, with emphasis on best practices and critical aspects of running these simulations. Additionally, we analyze challenges and successes when running alchemical free energy calculations on membrane-associated proteins. Finally, we highlight the value of alchemical free energy calculations calculations in drug discovery and their applicability in the pharmaceutical industry.

Alchemical Calculation of Relative Free Energies for Charge-Changing Mutations at Protein-Protein Interfaces Considering Fixed and Variable Protonation States

The calculation of relative free energies (ΔΔG) for charge-changing mutations at protein-protein interfaces through alchemical methods remains challenging due to variations in the system's net charge during charging steps, the possibility of mutated and contacting ionizable residues occurring in various protonation states, and undersampling issues. In this study, we present a set of strategies, collectively termed TIRST/TIRST-H+, to address some of these challenges. Our approaches combine thermodynamic integration (TI) with the prediction of pKa shifts to calculate ΔΔG values. Moreover, special sets of restraints are employed to keep the alchemically transformed molecules separated. The accuracy of the devised approaches was assessed on a large and diverse data set comprising 164 point mutations of charged residues (Asp, Glu, Lys, and Arg) to Ala at the protein-protein interfaces of complexes with known three-dimensional structures. Mean absolute and root-mean-square errors ranging from 1.38 to 1.66 and 1.89 to 2.44 kcal/mol, respectively, and Pearson correlation coefficients of ∼0.6 were obtained when testing the approaches on the selected data set using the GPU-TI module of Amber18 suite and the ff14SB force field. Furthermore, the inclusion of variable protonation states for the mutated acid residues improved the accuracy of the predicted ΔΔG values. Therefore, our results validate the use of TIRST/TIRST-H+ in prospective studies aimed at evaluating the impact of charge-changing mutations to Ala on the stability of protein-protein complexes.

Broadening the Scope of Binding Free Energy Calculations Using a Separated Topologies Approach

Binding free energy calculations predict the potency of compounds to protein binding sites in a physically rigorous manner and see broad application in prioritizing the synthesis of novel drug candidates. Relative binding free energy (RBFE) calculations have emerged as an industry-standard approach to achieve highly accurate rank-order predictions of the potency of related compounds; however, this approach requires that the ligands share a common scaffold and a common binding mode, restricting the methods' domain of applicability. This is a critical limitation since complex modifications to the ligands, especially core hopping, are very common in drug design. Absolute binding free energy (ABFE) calculations are an alternate method that can be used for ligands that are not congeneric. However, ABFE suffers from a known problem of long convergence times due to the need to sample additional degrees of freedom within each system, such as sampling rearrangements necessary to open and close the binding site. Here, we report on an alternative method for RBFE, called Separated Topologies (SepTop), which overcomes the issues in both of the aforementioned methods by enabling large scaffold changes between ligands with a convergence time comparable to traditional RBFE. Instead of only mutating atoms that vary between two ligands, this approach performs two absolute free energy calculations at the same time in opposite directions, one for each ligand. Defining the two ligands independently allows the comparison of the binding of diverse ligands without the artificial constraints of identical poses or a suitable atom-atom mapping. This approach also avoids the need to sample the unbound state of the protein, making it more efficient than absolute binding free energy calculations. Here, we introduce an implementation of SepTop. We developed a general and efficient protocol for running SepTop, and we demonstrated the method on four diverse, pharmaceutically relevant systems. We report the performance of the method, as well as our practical insights into the strengths, weaknesses, and challenges of applying this method in an industrial drug design setting. We find that the accuracy of the approach is sufficiently high to rank order ligands with an accuracy comparable to traditional RBFE calculations while maintaining the additional flexibility of SepTop.

Enhancing Hit Discovery in Virtual Screening through Absolute Protein-Ligand Binding Free-Energy Calculations

In the hit identification stage of drug discovery, a diverse chemical space needs to be explored to identify initial hits. Contrary to empirical scoring functions, absolute protein-ligand binding free-energy perturbation (ABFEP) provides a theoretically more rigorous and accurate description of protein-ligand binding thermodynamics and could, in principle, greatly improve the hit rates in virtual screening. In this work, we describe an implementation of an accurate and reliable ABFEP method in FEP+. We validated the ABFEP method on eight congeneric compound series binding to eight protein receptors including both neutral and charged ligands. For ligands with net charges, the alchemical ion approach is adopted to avoid artifacts in electrostatic potential energy calculations. The calculated binding free energies correlate with experimental results with a weighted average of R² = 0.55 for the entire dataset. We also observe an overall root-mean-square error (RMSE) of 1.1 kcal/mol after shifting the zero-point of the simulation data to match the average experimental values. Through ABFEP calculations using apo versus holo protein structures, we demonstrated that the protein conformational and protonation state changes between the apo and holo proteins are the main physical factors contributing to the protein reorganization free energy manifested by the overestimation of raw ABFEP calculated binding free energies using the holo structures of the proteins. Furthermore, we performed ABFEP calculations in three virtual screening applications for hit enrichment. ABFEP greatly improves the hit rates as compared to docking scores or other methods like metadynamics. The good performance of ABFEP in rank ordering compounds demonstrated in this work confirms it as a useful tool to improve the hit rates in virtual screening, thus facilitating hit discovery.

Accurate Host-Guest Binding Free Energies Using the AMOEBA Polarizable Force Field

A grand challenge of computational biophysics is accurate prediction of interactions between molecules. Molecular dynamics (MD) simulations have recently gained much interest as a tool to directly compute rigorous intermolecular binding affinities. The choice of a fixed point-charge or polarizable multipole force field used in MD is a topic of ongoing discussion. To compare alternative methods, we participated in the SAMPL7 and SAMPL8 Gibb octaacid host-guest challenges to assess the Atomic Multipole Optimized Energetics for Biomolecular Applications (AMOEBA) polarizable multipole force field. Advantages of AMOEBA over fixed charge models include improved representation of molecular electrostatic potentials and better description of water occupying the unligated host cavity. Prospective predictions for 26 host-guest systems exhibit a mean unsigned error vs experiment of 0.848 kcal/mol across all absolute binding free energies, demonstrating excellent agreement between computational and experimental results. In addition, we explore two topics related to the inclusion of ions in MD simulations: use of a neutral co-alchemical protocol and the effect of salt concentration on binding affinity. Use of the co-alchemical method minimally affects computed energies, but salt concentration significantly perturbs our binding results. Higher salt concentration strengthens binding through classical charge screening. In particular, added Na+ ions screen negatively charged carboxylate groups near the binding cavity, thereby diminishing repulsive coulomb interactions with negatively charged guests. Overall, the AMOEBA results demonstrate the accuracy available through a force field providing a detailed energetic description of the four octaacid hosts and 13 charged organic guests. Use of the AMOEBA polarizable atomic multipole force field in conjunction with an alchemical free energy protocol can achieve chemical accuracy in application to realistic molecular systems.

Validation of the Alchemical Transfer Method for the Estimation of Relative Binding Affinities of Molecular Series

The accurate prediction of protein-ligand binding affinities is crucial for drug discovery. Alchemical free energy calculations have become a popular tool for this purpose. However, the accuracy and reliability of these methods can vary depending on the methodology. In this study, we evaluate the performance of a relative binding free energy protocol based on the alchemical transfer method (ATM), a novel approach based on a coordinate transformation that swaps the positions of two ligands. The results show that ATM matches the performance of more complex free energy perturbation (FEP) methods in terms of Pearson correlation but with marginally higher mean absolute errors. This study shows that the ATM method is competitive compared to more traditional methods in speed and accuracy and offers the advantage of being applicable with any potential energy function.

Estimation of Nanoparticle's Surface Electrostatic Potential in Solution Using Acid-Base Molecular Probes I: In Silico Implementation for Surfactant Micelles

Surface electrostatic potential Ψ is a key characteristic of colloid particles. Since the surface of the particles adsorbs various compounds and facilitates chemical reactions between them, Ψ largely affects the properties of adsorbed reactants and governs the flow of chemical reactions occurring between them. One of the most popular methods for estimating Ψ in hydrophilic colloids, such as micellar surfactant solutions and related systems, is the application of molecular probes, predominantly acid-base indicator dyes. The Ψ value is calculated from the difference of the probe's indices of the apparent acidity constant between the examined colloid solution and, usually, some other colloid solution with noncharged particles. Here, we show how to implement this method in silico using alchemical free energy calculations within the framework of molecular dynamics simulations. The proposed implementation is tested on surfactant micelles and is shown to predict experimental Ψ values with quantitative accuracy depending on the kind of surfactant. The sources of errors in the method are discussed, and recommendations for its application are given.

The slow but steady rise of binding free energy calculations in drug discovery

Binding free energy calculations are increasingly used in drug discovery research to predict protein-ligand binding affinities and to prioritize candidate drug molecules accordingly. It has taken decades of collective effort to transform this academic concept into a technology adopted by the pharmaceutical and biotech industry. Having personally witnessed and taken part in this transformation, here I recount the (incomplete) list of problems that had to be solved to make this computational tool practical and suggest areas of future development.

Evolutionary divergence in the conformational landscapes of tyrosine vs serine/threonine kinases

Inactive conformations of protein kinase catalytic domains where the DFG motif has a "DFG-out" orientation and the activation loop is folded present a druggable binding pocket that is targeted by FDA-approved 'type-II inhibitors' in the treatment of cancers. Tyrosine kinases (TKs) typically show strong binding affinity with a wide spectrum of type-II inhibitors while serine/threonine kinases (STKs) usually bind more weakly which we suggest here is due to differences in the folded to extended conformational equilibrium of the activation loop between TKs vs. STKs. To investigate this, we use sequence covariation analysis with a Potts Hamiltonian statistical energy model to guide absolute binding free-energy molecular dynamics simulations of 74 protein-ligand complexes. Using the calculated binding free energies together with experimental values, we estimated free-energy costs for the large-scale (~17-20 Å) conformational change of the activation loop by an indirect approach, circumventing the very challenging problem of simulating the conformational change directly. We also used the Potts statistical potential to thread large sequence ensembles over active and inactive kinase states. The structure-based and sequence-based analyses are consistent; together they suggest TKs evolved to have free-energy penalties for the classical 'folded activation loop' DFG-out conformation relative to the active conformation, that is, on average, 4-6 kcal/mol smaller than the corresponding values for STKs. Potts statistical energy analysis suggests a molecular basis for this observation, wherein the activation loops of TKs are more weakly 'anchored' against the catalytic loop motif in the active conformation and form more stable substrate-mimicking interactions in the inactive conformation. These results provide insights into the molecular basis for the divergent functional properties of TKs and STKs, and have pharmacological implications for the target selectivity of type-II inhibitors.

An efficient method to predict protein thermostability in alanine mutation

The relationship between protein sequence and its thermodynamic stability is a critical aspect of computational protein design. In this work, we present a new theoretical method to calculate the free energy change (ΔΔG) resulting from a single-point amino acid mutation to alanine in a protein sequence. The method is derived based on physical interactions and is very efficient in estimating the free energy changes caused by a series of alanine mutations from just a single molecular dynamics (MD) trajectory. Numerical calculations are carried out on a total of 547 alanine mutations in 19 diverse proteins whose experimental results are available. The comparison between the experimental ΔΔG_exp and the calculated values shows a generally good correlation with a correlation coefficient of 0.67. Both the advantages and limitations of this method are discussed. This method provides an efficient and valuable tool for protein design and engineering.

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.

AMBER Drug Discovery Boost Tools: Automated Workflow for Production Free-Energy Simulation Setup and Analysis (ProFESSA)

We report an automated workflow for production free-energy simulation setup and analysis (ProFESSA) using the GPU-accelerated AMBER free-energy engine with enhanced sampling features and analysis tools, part of the AMBER Drug Discovery Boost package that has been integrated into the AMBER22 release. The workflow establishes a flexible, end-to-end pipeline for performing alchemical free-energy simulations that brings to bear technologies, including new enhanced sampling features and analysis tools, to practical drug discovery problems. ProFESSA provides the user with top-level control of large sets of free-energy calculations and offers access to the following key functionalities: (1) automated setup of file infrastructure; (2) enhanced conformational and alchemical sampling with the ACES method; and (3) network-wide free-energy analysis with the optional imposition of cycle closure and experimental constraints. The workflow is applied to perform absolute and relative solvation free-energy and relative ligand-protein binding free-energy calculations using different atom-mapping procedures. Results demonstrate that the workflow is internally consistent and highly robust. Further, the application of a new network-wide Lagrange multiplier constraint analysis that imposes key experimental constraints substantially improves binding free-energy predictions.

Electrostatic Contributions to the Binding Free Energy of Nicotine to the Acetylcholine Binding Protein

Biomolecular binding relies on specific attractive interactions between two partner molecules, including electrostatics, dispersion, hydrophobicity, and solvation. Assessing the contributions of electrostatic interactions to binding is key to the understanding of ligand binding mechanisms and the design of improved biomolecular binders. For example, nicotine is a well-known agonist of nicotinic acetylcholine receptors (nAChRs), but the molecular mechanisms for the differential action of nicotine on brain and muscle nAChRs remain elusive. In this work, we have chosen the acetylcholine binding protein (AChBP) in complex with nicotine as a model system to interrogate the electrostatic contributions to nicotine binding. Our absolute binding free energy simulations confirm that nicotine binds AChBP predominantly in its protonated (charged) form. By comparing energetic contributions from decomposed interactions for either neutral or charged nicotine, our calculations shed light on the nature of the binding of nicotine to the AChBP. The preferred binding of charged nicotine over neutral nicotine originates from its stronger electrostatic interactions with AChBP, a cation-π interaction to a tryptophan residue and a hydrogen bond between nicotine and the backbone carbonyl of the tryptophan, whereas the major force driving the binding process appears to be van der Waals interactions. The various nonelectrostatic terms can also indirectly modulate the electrostatic interactions through fine-tuning the binding pose of the ligand in the binding site, providing an explanation of why the binding specificity of nicotine to the brain versus muscle nAChRs is driven by electrostatic interaction, given that the immediate binding site residues, including the key tryptophan residue, are identical in the two receptors.

Meta-Analysis Reveals That Absolute Binding Free-Energy Calculations Approach Chemical Accuracy

Systematic and quantitative analysis of the reliability of formally exact methods that in silico calculate absolute protein-ligand binding free energies remains lacking. Here, we provide, for the first time, evidence-based information on the reliability of these methods by statistically studying 853 cases from 34 different research groups through meta-analysis. The results show that formally exact methods approach chemical accuracy (error = 1.58 kcal/mol), even if people are challenging difficult tasks such as blind drug screening in recent years. The geometrical-pathway-based methods prove to possess a better convergence ability than the alchemical ones, while the latter have a larger application range. We also reveal the importance of always using the latest force fields to guarantee reliability and discuss the pros and cons of turning to an implicit solvent model in absolute binding free-energy calculations. Moreover, based on the meta-analysis, an evidence-based guideline for in silico binding free-energy calculations is provided.

Best Practices in Constant pH MD Simulations: Accuracy and Sampling

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

Obtaining Consistent Free Energies for Ion Binding at Surfaces from Solution: Pathways versus Alchemy for Determining Kink Site Stability

Ion incorporation or removal from a solid at the interface with solution is a fundamental part of crystal growth. Despite this, there have been few quantitative determinations of the thermodynamics for such processes from atomistic molecular dynamics due to the associated technical challenges. In this study, we compute the free energies for ion removal from kink sites at the interface between NaCl and water as an illustrative example. To examine the influence of the free energy technique used, we compare methods that follow an explicit pathway for dissolution with those that focus on the thermodynamics of the initial and final states using metadynamics and free energy perturbation, respectively. While the initial results of the two approaches are found to be completely different, it is demonstrated that the thermodynamics can be reconciled with appropriate corrections for the standard states, thus illustrating the need for caution in interpreting raw free energy curves for ion binding as widely found in the literature. In addition, a new efficient approach is introduced to correct for the system size dependence of kink site energies both due to the periodic interaction of charges in an inhomogeneous dielectric system and due to the dipolar interactions between pairs of kinks along a row. Ultimately, it is shown that with suitable care, both classes of free energy techniques are capable of producing kink site stabilities that are consistent with the solubility of the underlying bulk solid. However, the precise values for individual kink sites exhibit a small systematic offset, which can be ascribed to the contribution of the interfacial potential to the pathway-based results. For the case of NaCl, the free energies of the kink sites relative to a 1 M aqueous solution for Na⁺ and Cl^- are found to be surprisingly different and of opposite sign, despite the ions having very similar hydration free energies.

Evaluating the use of absolute binding free energy in the fragment optimisation process

Key to the fragment optimisation process within drug design is the need to accurately capture the changes in affinity that are associated with a given set of chemical modifications. Due to the weakly binding nature of fragments, this has proven to be a challenging task, despite recent advancements in leveraging experimental and computational methods. In this work, we evaluate the use of Absolute Binding Free Energy (ABFE) calculations in guiding fragment optimisation decisions, retrospectively calculating binding free energies for 59 ligands across 4 fragment elaboration campaigns. We first demonstrate that ABFEs can be used to accurately rank fragment-sized binders with an overall Spearman's r of 0.89 and a Kendall τ of 0.67, although often deviating from experiment in absolute free energy values with an RMSE of 2.75 kcal/mol. We then also show that in several cases, retrospective fragment optimisation decisions can be supported by the ABFE calculations. Comparing against cheaper endpoint methods, namely Nwat-MM/GBSA, we find that ABFEs offer better ranking power and correlation metrics. Our results indicate that ABFE calculations can usefully guide fragment elaborations to maximise affinity.

Absolute binding free energies of mucroporin and its analog mucroporin-M1 with the heptad repeat 1 domain and RNA-dependent RNA polymerase of SARS-CoV-2

The peptide Mucroporin and its analog Mucroporin-M1 were studied using the molecular docking and molecular dynamics simulation of their complexation with two protein targets, the Heptad Repeat 1 (HR1) domain and RNA-dependent RNA polymerase (RdRp) of SARS-CoV-2. The molecular docking of the peptide-protein complexes was performed using the glowworm swarm optimization algorithm. The lowest energy poses were submitted to molecular dynamics simulation. Then, the binding free energies of Mucroporin and its analog Mucroporin-M1 with these two protein targets were calculated using the Multistate Bennett Acceptance Ratio (MBAR) method. It was verified that the peptides/HR1 domain complex showed stability in the interaction site determined by molecular docking. It was also found that Mucroporin-M1 has a much higher affinity than Mucroporin to the HR1 protein target. The peptides showed similar stability and affinity at the NTP binding site in the RdRp protein. Additional experimental studies are needed to confirm the antiviral activity of Mucroporin-M1 and a possible mechanism of action against SARS-CoV-2. However, here we indicate that Mucroporin-M1 may have potential antiviral activity against the HR1 domain with the possibility for further peptide optimization.Communicated by Ramaswamy H. Sarma.

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.

Impacts of targeting different hydration free energy references on the development of ion potentials

Hydration free energy (HFE) as the most important solvation parameter is often targeted in ion model development, even though the reported values differ by dozens of kcal mol^-1 mainly due to the experimentally undetermined HFE of the proton ΔG°(H⁺). The choice of ΔG°(H⁺) obviously affects the hydration of single ions and the relative HFE between the ions with different (magnitude or sign) charges, and the impacts of targeted HFEs on the ion solvation and ion-ion interactions are largely unrevealed. Here we designed point charge models of K⁺, Mg²⁺, Al³⁺, and Cl^- ions targeting a variety of HFE references and then investigated the HFE influences on the simulations of dilute and concentrated ion solutions and of the salt ion pairs in gas, liquid, and solid phases. Targeting one more property of ion-water oxygen distances (IOD) leaves the ion-water binding distance invariant, while the binding strength increases with the decreasing (more negative) HFE of ions as a result of a decrease in ΔG°(H⁺) for the cation and an increase in ΔG°(H⁺) for the anion. The increase in ΔG°(H⁺) leads to strengthened cation-anion interactions and thus to close ion-ion contacts, low osmotic pressures, and small activity derivatives in concentrated ion solutions as well as too stable ion pairs of the salts in different phases. The ion diffusivity and water exchange rates around the ions are simply not HFE dependent but rather more complex. Targeting both the aqueous IOD and salt crystal properties of KCl was also attempted and the comparison between different models indicates the complexity and challenge in obtaining a balanced performance between different phases using classical force fields. Our results also support that a real ΔG°(H⁺) value of -259.8 kcal mol^-1 recommended by Hünenberger and Reif guides ion models to reproduce ion-water and ion-ion interactions reasonably at relatively low salt concentrations. Simulations of a metalloprotein show that a relatively more positive ΔG°(H⁺) for Mg²⁺ model is better for a reasonable description of the metal binding network.

Improving the Potential of Mean Force and Nonequilibrium Pulling Simulations by Simultaneous Alchemical Modifications

We present an approach combining alchemical modifications and physical-pathway methods to calculate absolute binding free energies. The employed physical-pathway method is either a stratified umbrella sampling to calculate a potential of mean force or nonequilibrium pulling. We devised two basic approaches: the simultaneous approach (S-approach), where, along the physical unbinding pathway, an alchemical transformation of ligand-protein interactions is installed and deinstalled, and the prior-plus-simultaneous approach (PPS-approach), where, prior to the physical-pathway simulation, an alchemical transformation of ligand-protein interactions is installed in the binding site and deinstalled during the physical-pathway simulation. Using a mutant of T4 lysozyme with a benzene ligand as an example, we show that installation and deinstallation of soft-core interactions concurrent with physical ligand unbinding (S-approach) allow successful potential of mean force calculations and nonequilibrium pulling simulations despite the problems posed by the occluded nature of the lysozyme binding pocket. Good agreement between the potential of the mean-force-based S-approach and double decoupling simulations as well as a remarkable efficiency and accuracy of the nonequilibrium-pulling-based S-approach is found. The latter turned out to be more compute-efficient than the potential of mean force calculation by approximately 70%. Furthermore, we illustrate the merits of reducing ligand-protein interactions prior to potential of mean force calculations using the murine double minute homologue protein MDM2 with a p53-derived peptide ligand (PPS-approach). Here, the problem of breaking strong interactions in the binding pocket is transferred to a prior alchemical transformation that reduces the free-energy barrier between the bound and unbound state in the potential of mean force. Besides, disentangling physical ligand displacement from the deinstallation of ligand-protein interactions was seen to allow a more uniform sampling of distance histograms in the umbrella sampling. In the future, physical ligand unbinding combined with simultaneous alchemical modifications may prove useful in the calculation of protein-protein binding free energies, where sampling problems posed by multiple, possibly sticky interactions and potential steric clashes can thus be reduced.

How Peptides Bind to PSD-95/Discs-Large/ZO-1 Domains

PSD-95/discs-large/ZO-1 (PDZ) domains form a large family of adaptor proteins that bind to the C-terminal tails of their binding partner proteins. Via extensive molecular dynamics simulations and alchemical free energy calculations, we characterized the binding modi of phosphorylated and unphosphorylated EQVSAV peptides and of a EQVEAV phosphate mimic to the hPTP1E PDZ2 and MAGI1 PDZ1 domains. The simulations reproduced the well-known binding characteristics such as tight coordination of the peptidic carboxyl tail and pronounced hydrogen bonding between the peptide backbone and the backbone atoms of a β-sheet in PDZ. Overall, coordination by hPTP1E PDZ2 appeared tighter than by MAGI1 PDZ1. Simulations of wild-type PDZ and arginine mutants suggest that contacts with Arg79/85 in hPTP1E/MAGI1 are more important for the EQVEAV peptide than for EQVSAV. Alchemical free energy calculations and PaCS-MD simulations could well reproduce the difference in binding free energy between unphosphorylated EQVSAV and EQVEAV peptides and the absolute binding free energy of EQVSAV. However, likely due to small force field inaccuracies, the simulations erroneously favored binding of the phosphorylated peptide instead of its unphosphorylated counterpart, which is in contrast to the experiment.

A Computational Model for the PLP-Dependent Enzyme Methionine γ-Lyase

Pyridoxal-5′-phosphate (PLP) is a cofactor in the reactions of over 160 enzymes, several of which are implicated in diseases. Methionine γ-lyase (MGL) is of interest as a therapeutic protein for cancer treatment. It binds PLP covalently through a Schiff base linkage and digests methionine, whose depletion is damaging for cancer cells but not normal cells. To improve MGL activity, it is important to understand and engineer its PLP binding. We develop a simulation model for MGL, starting with force field parameters for PLP in four main states: two phosphate protonation states and two tautomeric states, keto or enol for the Schiff base moiety. We used the force field to simulate MGL complexes with each form, and showed that those with a fully-deprotonated PLP phosphate, especially keto, led to the best agreement with MGL structures in the PDB. We then confirmed this result through alchemical free energy simulations that compared the keto and enol forms, confirming a moderate keto preference, and the fully-deprotonated and singly-protonated phosphate forms. Extensive simulations were needed to adequately sample conformational space, and care was needed to extrapolate the protonation free energy to the thermodynamic limit of a macroscopic, dilute protein solution. The computed phosphate pK_a was 5.7, confirming that the deprotonated, −2 form is predominant. The PLP force field and the simulation methods can be applied to all PLP enzymes and used, as here, to reveal fine details of structure and dynamics in the active site.

A deep potential model with long-range electrostatic interactions

Machine learning models for the potential energy of multi-atomic systems, such as the deep potential (DP) model, make molecular simulations with the accuracy of quantum mechanical density functional theory possible at a cost only moderately higher than that of empirical force fields. However, the majority of these models lack explicit long-range interactions and fail to describe properties that derive from the Coulombic tail of the forces. To overcome this limitation, we extend the DP model by approximating the long-range electrostatic interaction between ions (nuclei + core electrons) and valence electrons with that of distributions of spherical Gaussian charges located at ionic and electronic sites. The latter are rigorously defined in terms of the centers of the maximally localized Wannier distributions, whose dependence on the local atomic environment is modeled accurately by a deep neural network. In the DP long-range (DPLR) model, the electrostatic energy of the Gaussian charge system is added to short-range interactions that are represented as in the standard DP model. The resulting potential energy surface is smooth and possesses analytical forces and virial. Missing effects in the standard DP scheme are recovered, improving on accuracy and predictive power. By including long-range electrostatics, DPLR correctly extrapolates to large systems the potential energy surface learned from quantum mechanical calculations on smaller systems. We illustrate the approach with three examples: the potential energy profile of the water dimer, the free energy of interaction of a water molecule with a liquid water slab, and the phonon dispersion curves of the NaCl crystal.

Correction Schemes for Absolute Binding Free Energies Involving Lipid Bilayers

Absolute binding free-energy (ABFE) calculations are playing an increasing role in drug design, especially as they can be performed on a range of disparate compounds and direct comparisons between them can be made. It is, however, especially important to ensure that they are as accurate as possible, as unlike relative binding free-energy (RBFE) calculations, one does not benefit as much from a cancellation of errors during the calculations. In most modern implementations of ABFE calculations, a particle mesh Ewald scheme is typically used to treat the electrostatic contribution to the free energy. A central requirement of such schemes is that the box preserves neutrality throughout the calculation. There are many ways to deal with this problem that have been discussed over the years ranging from a neutralizing plasma with a post hoc correction term through to a simple co-alchemical ion within the same box. The post hoc correction approach is the most widespread. However, the vast majority of these studies have been applied to a soluble protein in a homogeneous solvent (water or salt solution). In this work, we explore which of the more common approaches would be the most suitable for a simulation box with a lipid bilayer within it. We further develop the idea of the so-called Rocklin correction for lipid-bilayer systems and show how such a correction could work. However, we also show that it will be difficult to make this generalizable in a practical way and thus we conclude that the use of a “co-alchemical ion” is the most useful approach for simulations involving lipid membrane systems.

Exploring the Effect of Enhanced Sampling on Protein Stability Prediction

Changes in protein stability due to side-chain mutations are evaluated by alchemical free-energy calculations based on classical molecular dynamics (MD) simulations in explicit solvent using the GROMOS force field. Three proteins of different complexity with a total number of 93 single-point mutations are analyzed, and the relative free-energy differences are discussed with respect to configurational sampling and (dis)agreement with experimental data. For the smallest protein studied, a 34-residue WW domain, the starting structure dependence of the alchemical free-energy changes, is discussed in detail. Deviations from previous simulations for the two other proteins are shown to result from insufficient sampling in the earlier studies. Hamiltonian replica exchange in combination with multiple starting structures and sufficient sampling time of more than 100 ns per intermediate alchemical state is required in some cases to reach convergence.

Computational Tools for Accurate Binding Free-Energy Prediction

A quantitative thermodynamic understanding of the noncovalent association of (bio)molecules is of central importance in molecular life sciences. An important quantity characterizing (bio)molecular association is the binding affinity or absolute binding free energy. In recent years, the computational prediction of absolute binding free energies has evolved considerably in terms of accuracy, computational speed, and user-friendliness. In this chapter, we first give an overview of how absolute free energies are defined and how they can be determined with computational means. We proceed with an outline of the theoretical basis of the two most reliable methods, potential of mean force, and double decoupling calculations. In particular, we describe how the sampling problem can be alleviated by application of restraints. Finally, we provide step-by-step instructions of how to set up corresponding molecular simulations with a commonly employed molecular dynamics simulation engine.

Relative Binding Free Energy Calculations for Ligands with Diverse Scaffolds with the Alchemical Transfer Method

We present an extension of the alchemical transfer method (ATM) for the estimation of relative binding free energies of molecular complexes applicable to conventional, as well as scaffold-hopping, alchemical transformations. Named ATM-RBFE, the method is implemented in the free and open-source OpenMM molecular simulation package and aims to provide a simpler and more generally applicable route to the calculation of relative binding free energies than what is currently available. ATM-RBFE is based on sound statistical mechanics theory and a novel coordinate perturbation scheme designed to swap the positions of a pair of ligands such that one is transferred from the bulk solvent to the receptor binding site while the other moves simultaneously in the opposite direction. The calculation is conducted directly in a single solvent box with a system prepared with conventional setup tools, without splitting of electrostatic and nonelectrostatic transformations, and without pairwise soft-core potentials. ATM-RBFE is validated here against the absolute binding free energies of the SAMPL8 GDCC host-guest benchmark set and against protein-ligand benchmark sets that include complexes of the estrogen receptor ERα and those of the methyltransferase EZH2. In each case the method yields self-consistent and converged relative binding free energy estimates in agreement with absolute binding free energies and reference literature values, as well as experimental measurements.

Rational Design of Nonbonded Point Charge Models for Monovalent Ions with Lennard-Jones 12-6 Potential

Ions are of central importance in nature, and a variety of potential models was proposed to model ions in different phases for an in-depth exploration of ion-related systems. Here, we developed point charge models of 14 monovalent ions with the traditional 12-6 Lennard-Jones (LJ) potential for use in conjunction with 11 water models of TIP3P, OPC3, SPC/E, SPC/E_b, TIP3P-FB, a99SB-disp, TIP4P-Ew, OPC, TIP4P/2005, TIP4P-D, and TIP4P-FB. The designed models reproduced the real hydration free energy (HFE) of ions and the ion-oxygen distance (IOD) in the first hydration shell accurately and simultaneously, a performance similar to the previously reported 12-6-4 LJ-type ion models (12-6 LJ plus an attractive C₄ term for cations or a repulsive one for anions). This work, along with our previous work on di-, tri-, and tetravalent metal cations (J. Chem. Inf. Model. 2021, 61, 4031-4044; J. Chem. Inf. Model. 2021, 61, 4613-4629), demonstrates the feasibility of the simple 12-6 LJ potential in ion modeling. In order for the 12-6 LJ potential to reproduce both the HFE and IOD, the LJ R parameters need to be close to Shannon's ionic radii for the highly charged cations and to the Stokes's van der Waals (vdW) radii for the monovalent ions. With an additional C₄ term, the R parameters of 12-6-4 LJ ion models agree well with the Stokes's vdW radii and have a more physical meaning. It appears that the C₄ term can be merged into the 12-6 LJ potential by a rational tuning of R and the LJ well depth. Simulations of the osmotic coefficients of alkali chloride solutions and the properties of gaseous and solid alkali halides indicate the necessity of further optimizing ion-ion interactions via, for instance, targeting more properties or using a more physical (polarizable) model.

Studying the Effect of Amino Acid Substitutions in the M2 Ion Channel of the Influenza Virus on the Antiviral Activity of the Aminoadamantane Derivative In Vitro and In Silico

Purpose: The aminoadamantane derivative of L-histidyl-1-adamantayl ethylamine hydrochloride (HCl*H-His-Rim) has showed a high inhibition level against influenza A virus strains in vitro. The aim of this work is to search and establish evidence of the direct effect of the drug on influenza A virus proton channel M2.

Methods: The compound HCl*H-His-Rim was obtained by classical peptide synthesis methods. Influenza A virus mutants of A/PuertoRico/8/34(H1N1) strain were obtained by reverse genetics methods. The mutant samples of the virus were cultured on chicken embryos with a virus titer in the hemagglutination test. ELISA was carried out on Madin-Darby canine kidney (MDCK) monolayer cells when multiplying the virus 10^-4-10^-6. The binding stability of HCl*H-His-Rim was compared to those of M2 (S31N) and M2 (S31N_A30T) channels by molecular dynamic (MD) modeling. The calculation was performed taking into account the interaction with the model lipid bilayer (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) in the presence of water molecules in accordance with the three-center model.

Results: It was found that HCl*H-His-Rim is a direct action drug against influenza A. The most likely conformation of drug binding to target protein has been shown. It has been found that the A30T mutation reduces the binding energy of the drug, and the results obtained in vitro have confirmed the data calculated in silico.

Conclusion: The mechanism of action of HCl*H-His-Rim is directly related to the suppression of the function of the proton channel M2 of influenza A virus.

Quantum simulations of SARS-CoV-2 main protease Mpro enable high-quality scoring of diverse ligands

The COVID-19 pandemic has led to unprecedented efforts to identify drugs that can reduce its associated morbidity/mortality rate. Computational chemistry approaches hold the potential for triaging potential candidates far more quickly than their experimental counterparts. These methods have been widely used to search for small molecules that can inhibit critical proteins involved in the SARS-CoV-2 replication cycle. An important target is the SARS-CoV-2 main protease Mpro, an enzyme that cleaves the viral polyproteins into individual proteins required for viral replication and transcription. Unfortunately, standard computational screening methods face difficulties in ranking diverse ligands to a receptor due to disparate ligand scaffolds and varying charge states. Here, we describe full density functional quantum mechanical (DFT) simulations of Mpro in complex with various ligands to obtain absolute ligand binding energies. Our calculations are enabled by a new cloud-native parallel DFT implementation running on computational resources from Amazon Web Services (AWS). The results we obtain are promising: the approach is quite capable of scoring a very diverse set of existing drug compounds for their affinities to M pro and suggest the DFT approach is potentially more broadly applicable to repurpose screening against this target. In addition, each DFT simulation required only ~ 1 h (wall clock time) per ligand. The fast turnaround time raises the practical possibility of a broad application of large-scale quantum mechanics in the drug discovery pipeline at stages where ligand diversity is essential.