Accelerating Alchemical Free Energy Prediction Using a Multistate Method: Application to Multiple Kinases.

Alchemical free-energy methods based on molecular dynamics (MD) simulations have become important tools to identify modifications of small organic molecules that improve their protein binding affinity during lead optimization. The routine application of pairwise free-energy methods to rank potential binders from best to worst is impacted by the combinatorial increase in calculations to perform when the number of molecules to assess grows. To address this fundamental limitation, our group has developed replica-exchange enveloping distribution sampling (RE-EDS), a pathway-independent multistate method, enabling the calculation of alchemical free-energy differences between multiple ligands (N > 2) from a single MD simulation. In this work, we apply the method to a set of four kinases with diverse binding pockets and their corresponding inhibitors (42 in total), chosen to showcase the general applicability of RE-EDS in prospective drug design campaigns. We show that for the targets studied, RE-EDS is able to model up to 13 ligands simultaneously with high sampling efficiency, leading to a substantial decrease in computational cost when compared to pairwise methods.

Accounting for Solvation Correlation Effects on the Thermodynamics of Water Networks in Protein Cavities.

Macromolecular recognition and ligand binding are at the core of biological function and drug discovery efforts. Water molecules play a significant role in mediating the protein-ligand interaction, acting as more than just the surrounding medium by affecting the thermodynamics and thus the outcome of the binding process. As individual water contributions are impossible to measure experimentally, a range of computational methods have emerged to identify hydration sites in protein pockets and characterize their energetic contributions for drug discovery applications. Even though several methods model solvation effects explicitly, they focus on determining the stability of specific water sites independently and neglect solvation correlation effects upon replacement of clusters of water molecules, which typically happens in hit-to-lead optimization. In this work, we rigorously determine the conjoint effects of replacing all combinations of water molecules in protein binding pockets through the use of the RE-EDS multistate free-energy method, which combines Hamiltonian replica exchange (RE) and enveloping distribution sampling (EDS). Applications on the small bovine pancreatic trypsin inhibitor and four proteins of the bromodomain family illustrate the extent of solvation correlation effects on water thermodynamics, with the favorability of replacement of the water sites by pharmacophore probes highly dependent on the composition of the water network and the pocket environment. Given the ubiquity of water networks in biologically relevant protein targets, we believe our approach can be helpful for computer-aided drug discovery by providing a pocket-specific and a priori systematic consideration of solvation effects on ligand binding and selectivity.

Predicting Collagen Triple Helix Stability through Additive Effects of Terminal Residues and Caps.

Collagen model peptides (CMPs) consisting of proline-(2S,4R)-hydroxyproline-glycine (POG) repeats have provided a breadth of knowledge of the triple helical structure of collagen, the most abundant protein in mammals. Predictive tools for triple helix stability have, however, lagged behind since the effect of CMPs with different frames ([POG]n , [OGP]n , or [GPO]n ) and capped or uncapped termini have so far been underestimated. Here, we elucidated the impact of the frame, terminal functional group and its charge on the stability of collagen triple helices. Combined experimental and theoretical studies with frame-shifted, capped and uncapped CMPs revealed that electrostatic interactions, strand preorganization, interstrand H-bonding, and steric repulsion at the termini contribute to triple helix stability. We show that these individual contributions are additive and allow for the prediction of the melting temperatures of CMP trimers.

Frame Shifts Affect the Stability of Collagen Triple Helices.

Collagen model peptides (CMPs), composed of proline-(2S,4R)-hydroxyproline-glycine (POG) repeat units, have been extensively used to study the structure and stability of triple-helical collagen─the dominant structural protein in mammals─at the molecular level. Despite the more than 50-year history of CMPs and numerous studies on the relationship between the composition of single-stranded CMPs and the thermal stability of the assembled triple helices, little attention has been paid to the effects arising from their terminal residues. Here, we show that frame-shifted CMPs, which share POG repeat units but terminate with P, O, or G, form triple helices with vastly different thermal stabilities. A melting temperature difference as high as 16 °C was found for triple helices from 20-mers Ac-OG[POG]6-NH2 and Ac-[POG]6PO-NH2, and triple helices of the constitutional isomers Ac-[POG]7-NH2 and Ac-[GPO]7-NH2 melt 10 °C apart. A combination of thermal denaturation, circular dichroism and NMR spectroscopic studies, and molecular dynamics simulations revealed that the stability differences originate from the propensity of the peptide termini to preorganize into a polyproline-II helical structure. Our results advise that care must be taken when designing peptide mimics of structural proteins, as subtle changes in the terminal residues can significantly affect their properties. Our findings also provide a general and straightforward tool for tuning the stability of CMPs for applications as synthetic materials and biological probes.

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 as well as 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 in water and chloroform using both the GROMOS software package and a proof-of-concept implementation in OpenMM.

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

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

Relative free-energy calculations for scaffold hopping-type transformations with an automated RE-EDS sampling procedure.

The calculation of relative free-energy differences between different compounds plays an important role in drug design to identify potent binders for a given protein target. Most rigorous methods based on molecular dynamics simulations estimate the free-energy difference between pairs of ligands. Thus, the comparison of multiple ligands requires the construction of a "state graph", in which the compounds are connected by alchemical transformations. The computational cost can be optimized by reducing the state graph to a minimal set of transformations. However, this may require individual adaptation of the sampling strategy if a transformation process does not converge in a given simulation time. In contrast, path-free methods like replica-exchange enveloping distribution sampling (RE-EDS) allow the sampling of multiple states within a single simulation without the pre-definition of alchemical transition paths. To optimize sampling and convergence, a set of RE-EDS parameters needs to be estimated in a pre-processing step. Here, we present an automated procedure for this step that determines all required parameters, improving the robustness and ease of use of the methodology. To illustrate the performance, the relative binding free energies are calculated for a series of checkpoint kinase 1 inhibitors containing challenging transformations in ring size, opening/closing, and extension, which reflect changes observed in scaffold hopping. The simulation of such transformations with RE-EDS can be conducted with conventional force fields and, in particular, without soft bond-stretching terms.

Replica-Exchange Enveloping Distribution Sampling: Calculation of Relative Free Energies in GROMOS.

Molecular dynamics (MD) simulations have become an important tool to investigate biological systems. Free-energy calculations based on MD are playing an increasingly important role for computer-aided drug design and material discovery in recent years. Free-energy differences between pairs of end-states can be estimated using well-established methods such as thermodynamic integration (TI) or Bennett's acceptance ratio (BAR). An attractive alternative is the recently developed replica-exchange enveloping distribution sampling (RE-EDS) method, which enables estimating relative free-energy differences between multiple molecules from a single simulation. Here, we provide an introduction to the principles underlying RE-EDS and give an overview of the RE-EDS pipeline. In addition, we provide a description of the two complementary tools RestraintMaker and amber2gromos. We briefly discuss the findings of three recent applications of RE-EDS to calculate relative binding or hydration free energies. In all three studies, good agreement was found between the results obtained using RE-EDS and experimental values as well as values obtained using other free-energy methods.

Recent developments in multiscale free energy simulations.

Physics-based free energy simulations enable the rigorous calculation of properties, such as conformational equilibria, solvation or binding free energies. While historically most applications have occurred at the atomistic level of resolution, a range of advances in the past years make it possible now to reliably cross the temporal, spatial and theory scales for the modeling of complex systems or the efficient prediction of results at the accuracy level of expensive quantum-mechanical calculations. In this mini-review, we discuss recent methodological advances as well as opportunities opened up by the introduction of machine learning approaches, which tackle the diverse challenges across the different scales, improve the accuracy and feasibility, and push the boundaries of multiscale free energy simulations.

Incorporation of sensing modalities into de novo designed fluorescence-activating proteins.

Through the efforts of many groups, a wide range of fluorescent protein reporters and sensors based on green fluorescent protein and its relatives have been engineered in recent years. Here we explore the incorporation of sensing modalities into de novo designed fluorescence-activating proteins, called mini-fluorescence-activating proteins (mFAPs), that bind and stabilize the fluorescent cis-planar state of the fluorogenic compound DFHBI. We show through further design that the fluorescence intensity and specificity of mFAPs for different chromophores can be tuned, and the fluorescence made sensitive to pH and Ca2+ for real-time fluorescence reporting. Bipartite split mFAPs enable real-time monitoring of protein-protein association and (unlike widely used split GFP reporter systems) are fully reversible, allowing direct readout of association and dissociation events. The relative ease with which sensing modalities can be incorporated and advantages in smaller size and photostability make de novo designed fluorescence-activating proteins attractive candidates for optical sensor engineering.

An integrated view of p53 dynamics, function, and reactivation.

The tumor suppressor p53 plays a vital role in responding to cell stressors such as DNA damage, hypoxia, and tumor formation by inducing cell-cycle arrest, senescence, or apoptosis. Expression level alterations and mutational frequency implicates p53 in most human cancers. In this review, we show how both computational and experimental methods have been used to provide an integrated view of p53 dynamics, function, and reactivation potential. We argue that p53 serves as an exceptional case study for developing methods in modeling intrinsically disordered proteins. We describe how these methods can be leveraged to improve p53 reactivation molecule design and other novel therapeutic modalities, such as PROteolysis TARgeting Chimeras (PROTACs).

The flexibility of ACE2 in the context of SARS-CoV-2 infection.

The coronavirus disease 2019 (COVID-19) pandemic has swept over the world in the past months, causing significant loss of life and consequences to human health. Although numerous drug and vaccine development efforts are underway, there are many outstanding questions on the mechanism of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) viral association to angiotensin-converting enzyme 2 (ACE2), its main host receptor, and host cell entry. Structural and biophysical studies indicate some degree of flexibility in the viral extracellular spike glycoprotein and at the receptor-binding domain (RBD)-receptor interface, suggesting a role in infection. Here, we perform explicitly solvated, all-atom, molecular dynamics simulations of the glycosylated, full-length, membrane-bound ACE2 receptor in both an apo and spike RBD-bound state to probe the intrinsic dynamics of the ACE2 receptor in the context of the cell surface. A large degree of fluctuation in the full-length structure is observed, indicating hinge bending motions at the linker region connecting the head to the transmembrane helix while still not disrupting the ACE2 homodimer or ACE2-RBD interfaces. This flexibility translates into an ensemble of ACE2 homodimer conformations that could sterically accommodate binding of the spike trimer to more than one ACE2 homodimer and suggests a mechanical contribution of the host receptor toward the large spike conformational changes required for cell fusion. This work presents further structural and functional insights into the role of ACE2 in viral infection that can potentially be exploited for the rational design of effective SARS-CoV-2 therapeutics.

AI-driven multiscale simulations illuminate mechanisms of SARS-CoV-2 spike dynamics.

We develop a generalizable AI-driven workflow that leverages heterogeneous HPC resources to explore the time-dependent dynamics of molecular systems. We use this workflow to investigate the mechanisms of infectivity of the SARS-CoV-2 spike protein, the main viral infection machinery. Our workflow enables more efficient investigation of spike dynamics in a variety of complex environments, including within a complete SARS-CoV-2 viral envelope simulation, which contains 305 million atoms and shows strong scaling on ORNL Summit using NAMD. We present several novel scientific discoveries, including the elucidation of the spike's full glycan shield, the role of spike glycans in modulating the infectivity of the virus, and the characterization of the flexible interactions between the spike and the human ACE2 receptor. We also demonstrate how AI can accelerate conformational sampling across different systems and pave the way for the future application of such methods to additional studies in SARS-CoV-2 and other molecular systems.

AI-Driven Multiscale Simulations Illuminate Mechanisms of SARS-CoV-2 Spike Dynamics.

We develop a generalizable AI-driven workflow that leverages heterogeneous HPC resources to explore the time-dependent dynamics of molecular systems. We use this workflow to investigate the mechanisms of infectivity of the SARS-CoV-2 spike protein, the main viral infection machinery. Our workflow enables more efficient investigation of spike dynamics in a variety of complex environments, including within a complete SARS-CoV-2 viral envelope simulation, which contains 305 million atoms and shows strong scaling on ORNL Summit using NAMD. We present several novel scientific discoveries, including the elucidation of the spike's full glycan shield, the role of spike glycans in modulating the infectivity of the virus, and the characterization of the flexible interactions between the spike and the human ACE2 receptor. We also demonstrate how AI can accelerate conformational sampling across different systems and pave the way for the future application of such methods to additional studies in SARS-CoV-2 and other molecular systems.

Beyond Shielding: The Roles of Glycans in the SARS-CoV-2 Spike Protein.

The ongoing COVID-19 pandemic caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has resulted in more than 28,000,000 infections and 900,000 deaths worldwide to date. Antibody development efforts mainly revolve around the extensively glycosylated SARS-CoV-2 spike (S) protein, which mediates host cell entry by binding to the angiotensin-converting enzyme 2 (ACE2). Similar to many other viral fusion proteins, the SARS-CoV-2 spike utilizes a glycan shield to thwart the host immune response. Here, we built a full-length model of the glycosylated SARS-CoV-2 S protein, both in the open and closed states, augmenting the available structural and biological data. Multiple microsecond-long, all-atom molecular dynamics simulations were used to provide an atomistic perspective on the roles of glycans and on the protein structure and dynamics. We reveal an essential structural role of N-glycans at sites N165 and N234 in modulating the conformational dynamics of the spike's receptor binding domain (RBD), which is responsible for ACE2 recognition. This finding is corroborated by biolayer interferometry experiments, which show that deletion of these glycans through N165A and N234A mutations significantly reduces binding to ACE2 as a result of the RBD conformational shift toward the "down" state. Additionally, end-to-end accessibility analyses outline a complete overview of the vulnerabilities of the glycan shield of the SARS-CoV-2 S protein, which may be exploited in the therapeutic efforts targeting this molecular machine. Overall, this work presents hitherto unseen functional and structural insights into the SARS-CoV-2 S protein and its glycan coat, providing a strategy to control the conformational plasticity of the RBD that could be harnessed for vaccine development.

Beyond Shielding: The Roles of Glycans in SARS-CoV-2 Spike Protein.

The ongoing COVID-19 pandemic caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has resulted in more than 15,000,000 infections and 600,000 deaths worldwide to date. Antibody development efforts mainly revolve around the extensively glycosylated SARS-CoV-2 spike (S) protein, which mediates the host cell entry by binding to the angiotensin-converting enzyme 2 (ACE2). Similar to many other viruses, the SARS-CoV-2 spike utilizes a glycan shield to thwart the host immune response. Here, we built a full-length model of glycosylated SARS-CoV-2 S protein, both in the open and closed states, augmenting the available structural and biological data. Multiple microsecond-long, all-atom molecular dynamics simulations were used to provide an atomistic perspective on the roles of glycans, and the protein structure and dynamics. We reveal an essential structural role of N-glycans at sites N165 and N234 in modulating the conformational dynamics of the spike's receptor binding domain (RBD), which is responsible for ACE2 recognition. This finding is corroborated by biolayer interferometry experiments, which show that deletion of these glycans through N165A and N234A mutations significantly reduces binding to ACE2 as a result of the RBD conformational shift towards the "down" state. Additionally, end-to-end accessibility analyses outline a complete overview of the vulnerabilities of the glycan shield of SARS-CoV-2 S protein, which may be exploited by therapeutic efforts targeting this molecular machine. Overall, this work presents hitherto unseen functional and structural insights into the SARS-CoV-2 S protein and its glycan coat, providing a strategy to control the conformational plasticity of the RBD that could be harnessed for vaccine development.

Markov state models and NMR uncover an overlooked allosteric loop in p53.

The tumor suppressor p53 is the most frequently mutated gene in human cancer, and thus reactivation of mutated p53 is a promising avenue for cancer therapy. Analysis of wildtype p53 and the Y220C cancer mutant long-timescale molecular dynamics simulations with Markov state models and validation by NMR relaxation studies has uncovered the involvement of loop L6 in the slowest motions of the protein. Due to its distant location from the DNA-binding surface, the conformational dynamics of this loop has so far remained largely unexplored. We observe mutation-induced stabilization of alternate L6 conformations, distinct from all experimentally-determined structures, in which the loop is both extended and located further away from the DNA-interacting surface. Additionally, the effect of the L6-adjacent Y220C mutation on the conformational landscape of the functionally-important loop L1 suggests an allosteric role to this dynamic loop and the inactivation mechanism of the mutation. Finally, the simulations reveal a novel Y220C cryptic pocket that can be targeted for p53 rescue efforts. Our approach exemplifies the power of the MSM methodology for uncovering intrinsic dynamic and kinetic differences among distinct protein ensembles, such as for the investigation of mutation effects on protein function.

On the Interpretation of subtilisin Carlsberg Time-Resolved Fluorescence Anisotropy Decays: Modeling with Classical Simulations.

In this work, we discuss the challenging time-resolved fluorescence anisotropy of subtilisin Carlsberg (SC), which contains a single Trp residue and is a model fluorescence system. Experimental decay rates and quenching data suggest that the fluorophore should be exposed to water, but the Trp is partially buried in a hydrophobic pocket in the crystallographic structure. In order to study this inconsistency, molecular dynamics simulations were performed to predict the anisotropy decay rates and emission wavelengths of the Trp. We confirmed the inconsistency of the crystallographic structure with the experimentally observed fluorescence data and performed free energy calculations to show that the buried Trp conformation is 2 orders of magnitude (∼3 kcal/mol) more stable than the solvent-exposed one. However, molecular dynamics simulations in which the Trp side chain was restricted to solvent-exposed conformations displayed a maximum Trp emission wavelength shifted toward lower energies and decay rates compatible with the experimentally probed rates. Therefore, if the solvent-exposed conformations are the most important emitters, the experimental anisotropy can be compatibilized with the crystallographic structure. The most likely explanation is that the fluorescence of the most probable conformation in solution, observed in the crystal, is quenched, and this is consistent with the low quantum yield of Trp113 of SC. Additionally, some experiments might have probed denatured or lysed SC structures. SC anisotropy provides an interesting target for the study of fluorescence anisotropy using simulations, which can be used to test and exemplify how modeling can aid the interpretation of experimental data in a system where structure and solution experiments appear to be inconsistent.

Improving the Efficiency of Ligand-Binding Protein Design with Molecular Dynamics Simulations.

Custom-designed ligand-binding proteins represent a promising class of macromolecules with exciting applications toward the design of new enzymes or the engineering of antibodies and small-molecule recruited proteins for therapeutic interventions. However, several challenges remain in designing a protein sequence such that the binding site organization results in high affinity interaction with a bound ligand. Here, we study the dynamics of explicitly solvated designed proteins through all-atom molecular dynamics (MD) simulations to gain insight into the causes that lead to the low affinity or instability of most of these designs, despite the prediction of their success by the computational design methodology. Simulations ranging from 500 to 1000 ns per replicate were conducted on 37 designed protein variants encompassing two distinct folds and a range of ligand affinities, resulting in more than 180 µs of combined sampling. The simulations provide retrospective insights into the properties affecting ligand affinity that can prove useful in guiding further steps of design optimization. Features indicate that entropic components are particularly important for affinity, which are not easily incorporated in the empirical models often used in design protocols. Additionally, we demonstrate that the application of machine learning approaches built upon the output from the simulations can help discriminate between successful and failed binders, such that MD could act as a screening step in protein design, resulting in a more efficient process.