Alchemical Free Energy Research Articles

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.

Predicting binding affinity changes from long-distance mutations using molecular dynamics simulations and Rosetta.

Computationally modeling how mutations affect protein-protein binding not only helps uncover the biophysics of protein interfaces, but also enables the redesign and optimization of protein interactions. Traditional high-throughput methods for estimating binding free energy changes are currently limited to mutations directly at the interface due to difficulties in accurately modeling how long-distance mutations propagate their effects through the protein structure. However, the modeling and design of such mutations is of substantial interest as it allows for greater control and flexibility in protein design applications. We have developed a method that combines high-throughput Rosetta-based side-chain optimization with conformational sampling using classical molecular dynamics simulations, finding significant improvements in our ability to accurately predict long-distance mutational perturbations to protein binding. Our approach uses an analytical framework grounded in alchemical free energy calculations while enabling exploration of a vastly larger sequence space. When comparing to experimental data, we find that our method can predict internal long-distance mutational perturbations with a level of accuracy similar to that of traditional methods in predicting the effects of mutations at the protein-protein interface. This work represents a new and generalizable approach to optimize protein free energy landscapes for desired biological functions.

Routine Molecular Dynamics Simulations Including Nuclear Quantum Effects: From Force Fields to Machine Learning Potentials.

We report the implementation of a multi-CPU and multi-GPU massively parallel platform dedicated to the explicit inclusion of nuclear quantum effects (NQEs) in the Tinker-HP molecular dynamics (MD) package. The platform, denoted Quantum-HP, exploits two simulation strategies: the Ring-Polymer Molecular Dynamics (RPMD) that provides exact structural properties at the cost of a MD simulation in an extended space of multiple replicas and the adaptive Quantum Thermal Bath (adQTB) that imposes the quantum distribution of energy on a classical system via a generalized Langevin thermostat and provides computationally affordable and accurate (though approximate) NQEs. We discuss some implementation details, efficient numerical schemes, and parallelization strategies and quickly review the GPU acceleration of our code. Our implementation allows an efficient inclusion of NQEs in MD simulations for very large systems, as demonstrated by scaling tests on water boxes with more than 200,000 atoms (simulated using the AMOEBA polarizable force field). We test the compatibility of the approach with Tinker-HP's recently introduced Deep-HP machine learning potentials module by computing water properties using the DeePMD potential with adQTB thermostatting. Finally, we show that the platform is also compatible with the alchemical free energy estimation capabilities of Tinker-HP and fast enough to perform simulations. Therefore, we study how NQEs affect the hydration free energy of small molecules solvated with the recently developed Q-AMOEBA water force field. Overall, the Quantum-HP platform allows users to perform routine quantum MD simulations of large condensed-phase systems and will help to shed new light on the quantum nature of important interactions in biological matter.

Recent Advances in Alchemical Binding Free Energy Calculations for Drug Discovery.

Rigorous physics-based methods to calculate binding free energies of protein-ligand complexes have become a valued component of structure-based drug design. Relative and absolute binding free energy calculations have been deployed prospectively in support of solving diverse drug discovery challenges. Here we review recent applications of binding free energy calculations to fragment growing and linking, scaffold hopping, binding pose validation, virtual screening, covalent enzyme inhibition, and positional analogue scanning. Furthermore, we discuss the merits of using protein models and highlight recent efforts to replace costly binding free energy calculations with predictions from machine learning models trained on a limited number of free energy perturbation or thermodynamic integration calculations thereby allowing for extended chemical space exploration.

Alchemical free energy to probe RNA/ligand interactions.

Alchemical free energy calculations are increasingly being applied in academic research as well as in the pharmaceutical industry. These methods originate in theories of Kirkwood and Zwanzig, where Kirkwood introduced a coupling parameter approach, and Zwanzig showed that the free energy difference between two states can be computed via an appropriate exponential average of energy differences over an ensemble of configurations. In recent years, alchemical free energy calculations have benefited from improved force fields, new sampling algorithms, and the emergence of parallel computing, which have been shown to yield accurate results. These calculations are routinely applied to large-scale free energy calculations for ligand modification, and for testing mutations in amino acids and nucleotides. Here, we report alchemical calculations on RNA/peptide and RNA/small-molecule complexes, which are conducted using free energy perturbation (FEP) and thermodynamic integration (TI) methods. We discovered that the peptides employ base-specific hydrogen bonding, sugar-phosphate backbone interactions, and cation-pi interactions for the recognition of RNA molecules. In contrast, the small molecules recognize the RNA via base-specific and pi-pi stacking interactions.

Absolute binding free energy calculations for molecules binding to the GPCR GCGR at the protein/bilayer/water interface.

Numerous structures have been solved of small molecules binding at the protein-lipid interface of membrane proteins, such as GPCRs and ion channels. Allosteric modulators that bind at these lipid-exposed sites provide a novel way to target these therapeutically important proteins with greater specificity. Computational tools play a major role in modern drug discovery pipelines, but these tools must be validated, extended, and developed before being applied prospectively to lipid-exposed sites. The creation of a small molecule with a high affinity for its target is a key step in the drug discovery process and the accurate prediction of affinity remains a central challenge in computational chemistry. Alchemical absolute binding free energy (ABFE) calculations have been shown to be highly accurate in predicting the binding affinity of small molecules to soluble proteins. However, ABFE calculations remain challenging, with little work assessing their application to membrane proteins. Our work has focused on extending ABFE calculations to lipid exposed sites. The binding free energy calculations for these sites are complicated by the fact the binding site is ‘solvated’ by lipids. We present a study on the GPCR GCGR in complex with the small molecule MK-0893. MK-0893 sits at the protein/bilayer/water interface. This, coupled with the size and charge of the molecule, makes ABFE calculations for this system particularly challenging. We discuss the adaptations we made to our ABFE workflow to overcome these challenges. Particularly, we make suggestions as to how to monitor the generation of a ‘correct’ apo state for the non-interacting ligand, with regards to the occupancy of the binding site and side-chain rotamers. These adaptations significantly improved our results and show a good agreement with experimental values from a publicly available SAR series.

Application of molecular dynamics to elucidate key binding interactions for the development of therapeutics against opioid overdose.

Opioid overdose has been a long-standing public health issue in the United States. This crisis has only been exacerbated by the Covid-19 pandemic, with overdose deaths involving opioids increasing from an estimated 70,029 in 2020 to 80,816 in 2021. While there are treatments for opioid overdose, there are currently no therapeutic tools to prevent it. Moreover, these treatments are limited to emergency scenarios, most notably due to induction of withdrawal and loss of analgesia. Fatal opioid overdoses are primarily attributed to opioid-induced respiratory depression (OIRD). As part of an ongoing collaboration, a class of cysteine esters have been identified that reverse OIRD without blocking the analgesic effects of the opioid or inducing withdrawal. Our working hypothesis is that these esters function by binding β-arrestin and altering either its interaction with the opioid receptor, downstream proteins, or both. Accordingly, we have used “flooding” molecular dynamics simulations to identify candidate binding sites for several related ligands on the surface of β-arrestin 1 and 2, which we then validate using alchemical ligand binding free energy calculations. The results will be used to suggest mutations that would disrupt binding and to better understand the mechanism by which these compounds combat opioid overdose.

Absolute binding affinities of phospholipids to Erwinia ligand-gated ion channel (ELIC) by streamlined alchemical free energy perturbations (SAFEP).

An increasing number of cryo-EM structures of membrane proteins contain putative lipid binding sites with partially resolved lipid densities. In principle, these densities might be identified by computing the binding free energies of candidate ligands by alchemical free energy perturbation (FEP). In practice, however, these systems have several challenging features including lipid flexibility, system heterogeneity, and low binding affinity. Each of these confounds the distinction between the bound and unbound state: flexible molecules have large accessible volumes in conformation space; heterogeneous environments complicate the unbound state; and open, low-affinity sites mean that a “bound” ligand may still have considerably higher entropy than more in a more traditional, high-affinity binding pocket. Our SAFEP method seeks to address these problems by defining a single collective variable to define the bound state. This distance-from-bound configuration (DBC) metric reliably classifies ligand-protein conformations as “bound” or “unbound.” SAFEP has been successfully applied to the computation of absolute binding affinities of cholesterol to three GPCRs and relative binding affinities of phospholipids to two mutants of ELIC based on structural data. We have extended SAFEP to the calculation of absolute binding affinities of phospholipids. This has required additional sampling of each state and refinement to the collective variables used, especially in the unbound state. While this is somewhat more computationally expensive than computing relative free energies, absolute affinities make functional questions more accessible and mitigate the error propagation inherent to relative binding affinities. It is our hope that these methods will aid in the interpretation of similar systems.

Comparison of cyclic AMP binding activity within HCN channel subfamily using atomistic molecular dynamics simulations of isolated cyclic nucleotide binding domains.

Hyperpolarization-activated cyclic nucleotide-gated (HCN) channels are a subfamily of homotetrameric, voltage gated potassium/sodium channels whose activity is modulated by cyclic AMP (cAMP). The binding of cAMP to an intracellular cyclic nucleotide binding domain (CNBD) promotes an increase in channel activity through a depolarizing shift in the voltage threshold necessary for channel activation, as well as kinetically favoring the active channel conformation, leading to improved conductance. While the precise mechanisms of this change are still being investigated, the shift in the CNDB region which is induced by the binding of cAMP is able to be propagated to the remainder of the channel through an adjacent C-linker, causing an alteration in the channel that affects conductance. Although various static structures of the four HCN isoforms with and without bound cyclic nucleotide have been solved and deposited on the Protein Data Bank, there is still a need to understand the precise binding interactions which lead to different sensitivities and selectivities for cyclic nucleotides within the subfamily. Using atomistic molecular dynamics, we have simulated the isolated CNBD of HCN1-4 bound to cAMP and have used both equilibrium simulations and alchemical free energy perturbation to characterize and compare the binding activities of the different isoforms. These simulations highlight the differential cyclic nucleotide binding behavior between HCN channels.

Oncogenic mutations in the DNA-binding domain of FOXO1 that disrupt folding: Quantitative insights from experiments and molecular simulation.

FOXO1, a member of the family of winged-helix motif Forkhead box (FOX) transcription factors, is the most abundantly expressed FOXO member in mature B cells. Sequencing of diffuse large B-cell lymphoma (DLBCL) tumors and cell lines identified specific mutations in the forkhead domain linked to loss of function. Differential scanning calorimetry and thermal shift assays were used to characterize how eight of these mutations affect the stability of the FOX domain. Mutations L183P and L183R were found to be particularly destabilizing. Electrophoresis mobility shift assays show these same mutations also disrupt FOXO1 binding to their canonical DNA sequences, suggesting that the loss of function is due to destabilization of the folded structure. Computational modeling of the effect of mutations on FOXO1 folding was performed using alchemical free energy perturbation (FEP), and a Markov model of the entire folding reaction was constructed from massively parallel molecular simulations, which predicts folding pathways involving the late folding of helix α3. Although FEP can qualitatively predict the destabilization from L183 mutations, we find that a simple hydrophobic transfer model, combined with estimates of unfolded-state solvent-accessible surface areas from molecular simulations, is able to more accurately predict changes in folding free energies due to mutations. These results suggest that the atomic detail provided by simulations is important for the accurate prediction of mutational effects on folding stability. Corresponding disease-associated mutations in other FOX family members support further experimental and computational studies of the folding mechanism of FOX domains.

Regioselective metabolism of cannabinoids by cytochrome P450 2J2.

Cannabichromene (CBC) is a major non-psychoactive phytocannabinoid extracted from the cannabis plant. CBC has been shown to have anti-inflammatory effects; hence it can serve as a non-narcotic alternative to tetrahydrocannabinol (THC). Nevertheless, the metabolic pathways and metabolites of CBC have not been fully characterized. Metabolism of phytocannabinoids typically occurs via a superfamily of hemoproteins called cytochrome P450s (CYP450s). Our liquid chromatography-mass spectrometry (LC/MS) experiments show that CBC is metabolized to three different products by various human CYP450s. We identified that CYP450 2J2 (CYP2J2) predominantly produces 8'-hydroxy-CBC, an allylic hydroxide, out of the three metabolites. To investigate the CBC metabolic regioselectivity in CYP2J2, we leveraged atomistic molecular dynamics (MD) and modeling of the system. Using 200 unique conformations of CYP2J2 obtained from MD, we performed ensemble molecular docking of CBC to reveal three putative binding modes. We then assessed the stability of the obtained CBC binding modes by simulating them individually for 1 μs. Our simulations captured an average distance of 4.8 Å between the CYP2J2 heme and the 8’ carbon atom on CBC for one of the predominant binding modes, at which the allylic hydroxide metabolite occurs experimentally. At the regioselective binding mode, we observe conserved hydrogen bonding networks, mainly between the CBC hydroxyl group and the backbone oxygen of I487, as well as N379, stabilized by a water molecule. We further corroborated the stability of the regioselective binding modes by calculating relative binding free energies using alchemical free energy perturbation calculations. In conjunction, our computational and experimental methodologies elucidate a likely mode for regioselective binding and metabolism of CBC by CYP2J2, furthering the way towards accessible plant-based non-psychoactive anti-inflammatory therapeutics.

Complementary in vitro functional studies and in silico alchemical free energy simulations of the sarco-endoplasmic reticulum calcium pump (SERCA) and phospholamban (PLN) complex.

SERCA is ubiquitously expressed in all eukaryotic cells and is responsible for actively transporting calcium ions from the cytoplasm into the SR/ER lumen. PLN is a small transmembrane peptide that regulates SERCA's calcium transport activity in cardiac muscle. Human genetic variants of PLN have been implicated in cardiomyopathies. PLN is a 52 amino acid peptide with an N-terminal cytoplasmic helix (residues 1-17), a short linker region (residues 18-27) and a C-terminal transmembrane helix (residues 28-52). Previous structural studies have identified the binding groove for PLN formed by transmembrane segments M2, M6, and M9 of SERCA (PDB: 4KYT). The structure of the SERCA-PLN complex consists of SERCA in a calcium-free, E1-like state and includes the linker region and transmembrane helix of PLN (residues 20-52). Here, we carry out in vitro functional studies and in silico free energy calculations to examine the effects of alanine scanning mutations in PLN (residues 20-52) on SERCA regulation. From the functional studies, changes in kinetics parameters such as calcium affinity (KCa) and maximal activity (VMax) of SERCA are determined for each of the alanine substitutions. From the in silico alchemical free energy simulations (using Compute Canada supercomputing resources), each of the residues in PLN are alchemically transformed to alanine to determine the change in Gibbs free energy (Δ;Δ;G) relative to the wild-type SERCA-PLN complex. Together, the complementary in vitro and in silico mutational studies provide novel insights into a membrane protein-peptide complex. The PLN variants are quantified in terms of their physiological significance – gain-of-function, loss-of-function, or neutral – and the molecular determinants that control binding specificity.

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.

Exploring PROTAC Cooperativity with Coarse-Grained Alchemical Methods.

Proteolysis targeting chimera (PROTAC) is a novel drug modality that facilitates the degradation of a target protein by inducing proximity with an E3 ligase. In this work, we present a new computational framework to model the cooperativity between PROTAC-E3 binding and PROTAC-target binding principally through protein-protein interactions (PPIs) induced by the PROTAC. Due to the scarcity and low resolution of experimental measurements, the physical and chemical drivers of these non-native PPIs remain to be elucidated. We develop a coarse-grained (CG) approach to model interactions in the target-PROTAC-E3 complexes, which enables converged thermodynamic estimations using alchemical free energy calculation methods despite an unconventional scale of perturbations. With minimal parametrization, we successfully capture fundamental principles of cooperativity, including the optimality of intermediate PROTAC linker lengths that originates from configurational entropy. We qualitatively characterize the dependency of cooperativity on PROTAC linker lengths and protein charges and shapes. Minimal inclusion of sequence- and conformation-specific features in our current force field, however, limits quantitative modeling to reproduce experimental measurements, but further development of the CG model may allow for efficient computational screening to optimize PROTAC cooperativity.

Taming multiple binding poses in alchemical binding free energy prediction: the β-cyclodextrin host-guest SAMPL9 blinded challenge.

We apply the Alchemical Transfer Method (ATM) and a bespoke fixed partial charge force field to the SAMPL9 bCD host-guest binding free energy prediction challenge that comprises a combination of complexes formed between five phenothiazine guests and two cyclodextrin hosts. Multiple chemical forms, competing binding poses, and computational modeling challenges pose significant obstacles to obtaining reliable computational predictions for these systems. The phenothiazine guests exist in solution as racemic mixtures of enantiomers related by nitrogen inversions that bind the hosts in various binding poses, each requiring an individual free energy analysis. Due to the large size of the guests and the conformational reorganization of the hosts, which prevent a direct absolute binding free energy route, binding free energies are obtained by a series of absolute and relative binding alchemical steps for each chemical species in each binding pose. Metadynamics-accelerated conformational sampling was found to be necessary to address the poor convergence of some numerical estimates affected by conformational trapping. Despite these challenges, our blinded predictions quantitatively reproduced the experimental affinities for the β-cyclodextrin host and, to a lesser extent, those with a methylated derivative. The work illustrates the challenges of obtaining reliable free energy data in in silico drug design for even seemingly simple systems and introduces some of the technologies available to tackle them.

High-Throughput Prediction of the Impact of Genetic Variability on Drug Sensitivity and Resistance Patterns for Clinically Relevant Epidermal Growth Factor Receptor Mutations from Atomistic Simulations

Mutations in the kinase domain of the epidermal growth factor receptor (EGFR) can be drivers of cancer and also trigger drug resistance in patients receiving chemotherapy treatment based on kinase inhibitors. A priori knowledge of the impact of EGFR variants on drug sensitivity would help to optimize chemotherapy and design new drugs that are effective against resistant variants before they emerge in clinical trials. To this end, we explored a variety of in silico methods, from sequence-based to "state-of-the-art" atomistic simulations. We did not find any sequence signal that can provide clues on when a drug-related mutation appears or the impact of such mutations on drug activity. Low-level simulation methods provide limited qualitative information on regions where mutations are likely to cause alterations in drug activity, and they can predict around 70% of the impact of mutations on drug efficiency. High-level simulations based on nonequilibrium alchemical free energy calculations show predictive power. The integration of these "state-of-the-art" methods into a workflow implementing an interface for parallel distribution of the calculations allows its automatic and high-throughput use, even for researchers with moderate experience in molecular simulations.

Open-channel structure of a pentameric ligand-gated ion channel reveals a mechanism of leaflet-specific phospholipid modulation

Pentameric ligand-gated ion channels (pLGICs) mediate synaptic transmission and are sensitive to their lipid environment. The mechanism of phospholipid modulation of any pLGIC is not well understood. We demonstrate that the model pLGIC, ELIC (Erwinia ligand-gated ion channel), is positively modulated by the anionic phospholipid, phosphatidylglycerol, from the outer leaflet of the membrane. To explore the mechanism of phosphatidylglycerol modulation, we determine a structure of ELIC in an open-channel conformation. The structure shows a bound phospholipid in an outer leaflet site, and structural changes in the phospholipid binding site unique to the open-channel. In combination with streamlined alchemical free energy perturbation calculations and functional measurements in asymmetric liposomes, the data support a mechanism by which an anionic phospholipid stabilizes the activated, open-channel state of a pLGIC by specific, state-dependent binding to this site.

High-Throughput Molecular Dynamics-Based Alchemical Free Energy Calculations for Predicting the Binding Free Energy Change Associated with the Selected Omicron Mutations in the Spike Receptor-Binding Domain of SARS-CoV-2.

The ongoing pandemic caused by SARS-CoV-2 has gone through various phases. Since the initial outbreak, the virus has mutated several times, with some lineages showing even stronger infectivity and faster spread than the original virus. Among all the variants, omicron is currently classified as a variant of concern (VOC) by the World Health Organization, as the previously circulating variants have been replaced by it. In this work, we have focused on the mutations observed in omicron sub lineages BA.1, BA.2, BA.4 and BA.5, particularly at the receptor-binding domain (RBD) of the spike protein that is responsible for the interactions with the host ACE2 receptor and binding of antibodies. Studying such mutations is particularly important for understanding the viral infectivity, spread of the disease and for tracking the escape routes of this virus from antibodies. Molecular dynamics (MD) based alchemical free energy calculations have been shown to be very accurate in predicting the free energy change, due to a mutation that could have a deleterious or a stabilizing effect on either the protein itself or its binding affinity to another protein. Here, we investigated the significance of five spike RBD mutations on the stability of the spike protein binding to ACE2 by free energy calculations using high throughput MD simulations. For comparison, we also used conventional MD simulations combined with a Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) based approach, and compared our results with the available experimental data. Overall, the alchemical free energy calculations performed far better than the MM-GBSA approach in predicting the individual impact of the mutations. When considering the experimental variation, the alchemical free energy method was able to produce a relatively accurate prediction for N501Y, the mutant that has previously been reported to increase the binding affinity to hACE2. On the other hand, the other individual mutations seem not to have a significant effect on the spike RBD binding affinity towards hACE2.

Free energy perturbation calculations of tetrahydroquinolines complexed to the first bromodomain of BRD4

Alchemical free energy perturbation (FEP) theory is widely used nowadays to calculate protein–ligand binding energies, often in support of drug discovery endeavours. We assess the accuracy and sensitivity of absolute FEP binding energies with respect to the CHARMM/CGenFF and the AMBER/GAFF force field parameterisations for a set of tetrahydroquinoline inhibitors of the first bromodomain of BRD4, a target of keen interest for the development of anti-cancer drugs. We find that AMBER/GAFF is better able than CHARMM/CGenFF to cover the range of and to distinguish between the relative binding energies of the 16 ligands.

A New Formulation for the Concerted Alchemical Calculation of van der Waals and Coulomb Components of Solvation Free Energies.

Alchemical free energy calculations via molecular dynamics have been widely used to obtain thermodynamic properties related to protein-ligand binding and solute-solvent interactions. Although soft-core modeling is the most common approach, the linear basis function (LBF) methodology [Naden, L. N.; et al. J. Chem. Theory Comput.2014, 10 (3), 1128; 2015, 11 (6), 2536] has emerged as a suitable alternative. It overcomes the end-point singularity of the scaling method while maintaining essential advantages such as ease of implementation and high flexibility for postprocessing analysis. In the present work, we propose a simple LBF variant and formulate an efficient protocol for evaluating van der Waals and Coulomb components of an alchemical transformation in tandem, in contrast to the prevalent sequential evaluation mode. To validate our proposal, which results from a careful optimization study, we performed solvation free energy calculations and obtained octanol-water partition coefficients of small organic molecules. Comparisons with results obtained via the sequential mode using either another LBF approach or the soft-core model attest to the effectiveness and correctness of our method. In addition, we show that a reaction field model with an infinite dielectric constant can provide very accurate hydration free energies when used instead of a lattice-sum method to model solute-solvent electrostatics.