Massively parallel free energy calculations for in silico affinity maturation of designed miniproteins.

Computational protein design efforts continue to make remarkable advances, yet the discovery of high-affinity binders typically requires large-scale experimental screening of site-saturated mutant (SSM) libraries. Here, we explore how massively parallel free energy methods can be used for in silico affinity maturation of de novo designed binding proteins. Using an expanded ensemble (EE) approach, we perform exhaustive relative binding free energy calculations for SSM variants of three miniproteins designed to bind influenza A H1 hemagglutinin by Chevalier et al. (2017). We compare our predictions to experimental ΔΔ G values inferred from a Bayesian analysis of the high-throughput sequencing data, and to state-of-the-art predictions made using the Flex ddG Rosetta protocol. A systematic comparison reveals prediction accuracies around 2 kcal/mol, and identifies net charge changes, large numbers of alchemical atoms, and slow side chain conformational dynamics as key contributors to the uncertainty of the EE predictions. Flex ddG predictions are more accurate on average, but highly conservative. In contrast, EE predictions can better classify stabilizing and destabilizing mutations. We also explored the ability of SSM scans to rationalize known affinity-matured variants containing multiple mutations, which are non-additive due to epistatic effects. Simple electrostatic models fail to explain non-additivity, but observed mutations are found at positions with higher Shannon entropies. Overall, this work suggests that simulation-based free energy methods can provide predictive information for in silico affinity maturation of designed miniproteins, with many feasible improvements to the efficiency and accuracy within reach.

Automatic Forward Model Parameterization with Bayesian Inference of Conformational Populations.

To quantify how well theoretical predictions of structural ensembles agree with experimental measurements, we depend on the accuracy of forward models. These models are computational frameworks that generate observable quantities from molecular configurations based on empirical relationships linking specific molecular properties to experimental measurements. Bayesian Inference of Conformational Populations (BICePs) is a reweighting algorithm that reconciles simulated ensembles with ensemble-averaged experimental observations, even when such observations are sparse and/or noisy. This is achieved by sampling the posterior distribution of conformational populations under experimental restraints as well as sampling the posterior distribution of uncertainties due to random and systematic error. In this study, we enhance the algorithm for the refinement of empirical forward model (FM) parameters. We introduce and evaluate two novel methods for optimizing FM parameters. The first method treats FM parameters as nuisance parameters, integrating over them in the full posterior distribution. The second method employs variational minimization of a quantity called the BICePs score that reports the free energy of "turning on" the experimental restraints. This technique, coupled with improved likelihood functions for handling experimental outliers, facilitates force field validation and optimization, as illustrated in recent studies (Raddi et al. 2023, 2024). Using this approach, we refine parameters that modulate the Karplus relation, crucial for accurate predictions of J -coupling constants based on dihedral angles ( ϕ ) between interacting nuclei. We validate this approach first with a toy model system, and then for human ubiquitin, predicting six sets of Karplus parameters for J H N H α 3 , J H α C ' 3 , J H N C ß 3 , J H N C ' 3 , J C ' C ß 3 , J C ' C ' 3 . This approach, which does not rely on any predetermined parameterization, enhances predictive accuracy and can be used for many applications.

Drug Discovery in Low Data Regimes: Leveraging a Computational Pipeline for the Discovery of Novel SARS-CoV-2 Nsp14-MTase Inhibitors.

The COVID-19 pandemic, caused by the SARS-CoV-2 virus, has led to significant global morbidity and mortality. A crucial viral protein, the non-structural protein 14 (nsp14), catalyzes the methylation of viral RNA and plays a critical role in viral genome replication and transcription. Due to the low mutation rate in the nsp region among various SARS-CoV-2 variants, nsp14 has emerged as a promising therapeutic target. However, discovering potential inhibitors remains a challenge. In this work, we introduce a computational pipeline for the rapid and efficient identification of potential nsp14 inhibitors by leveraging virtual screening and the NCI open compound collection, which contains 250,000 freely available molecules for researchers worldwide. The introduced pipeline provides a cost-effective and efficient approach for early-stage drug discovery by allowing researchers to evaluate promising molecules without incurring synthesis expenses. Our pipeline successfully identified seven promising candidates after experimentally validating only 40 compounds. Notably, we discovered NSC620333, a compound that exhibits a strong binding affinity to nsp14 with a dissociation constant of 427 ± 84 nM. In addition, we gained new insights into the structure and function of this protein through molecular dynamics simulations. We identified new conformational states of the protein and determined that residues Phe367, Tyr368, and Gln354 within the binding pocket serve as stabilizing residues for novel ligand interactions. We also found that metal coordination complexes are crucial for the overall function of the binding pocket. Lastly, we present the solved crystal structure of the nsp14-MTase complexed with SS148 (PDB:8BWU), a potent inhibitor of methyltransferase activity at the nanomolar level (IC50 value of 70 ± 6 nM). Our computational pipeline accurately predicted the binding pose of SS148, demonstrating its effectiveness and potential in accelerating drug discovery efforts against SARS-CoV-2 and other emerging viruses.

Markov State Model of Solvent Features Reveals Water Dynamics in Protein-Peptide Binding.

In this work, we investigate the role of solvent in the binding reaction of the p53 transactivation domain (TAD) peptide to its receptor MDM2. Previously, our group generated 831 µs of explicit-solvent aggregate molecular simulation trajectory data for the MDM2-p53 peptide binding reaction using large-scale distributed computing and subsequently built a Markov State Model (MSM) of the binding reaction (Zhou et al. 2017). Here, we perform a tICA analysis and construct an MSM with similar hyperparameters while using only solvent-based structural features. We find a remarkably similar landscape but accelerated implied timescales for the slowest motions. The solvent shells contributing most to the first tICA eigenvector are those centered on Lys24 and Thr18 of the p53 TAD peptide in the range of 3-6 Å. Important solvent shells were visualized to reveal solvation and desolvation transitions along the peptide-protein binding trajectories. Our results provide a solvent-centric view of the hydrophobic effect in action for a realistic peptide-protein binding scenario.

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

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

BICePs v2.0: Software for Ensemble Reweighting Using Bayesian Inference of Conformational Populations.

Bayesian Inference of Conformational Populations (BICePs) version 2.0 (v2.0) is a free, open-source Python package that reweights theoretical predictions of conformational state populations using sparse and/or noisy experimental measurements. In this article, we describe the implementation and usage of the latest version of BICePs (v2.0), a powerful, user-friendly and extensible package which makes several improvements upon the previous version. The algorithm now supports many experimental NMR observables (NOE distances, chemical shifts, J-coupling constants, and hydrogen-deuterium exchange protection factors), and enables convenient data preparation and processing. BICePs v2.0 can perform automatic analysis of the sampled posterior, including visualization, and evaluation of statistical significance and sampling convergence. We provide specific coding examples for these topics, and present a detailed example illustrating how to use BICePs v2.0 to reweight a theoretical ensemble using experimental measurements.

Folding@home: achievements from over twenty years of citizen science herald the exascale era.

Simulations of biomolecules have enormous potential to inform our understanding of biology but require extremely demanding calculations. For over twenty years, the Folding@home distributed computing project has pioneered a massively parallel approach to biomolecular simulation, harnessing the resources of citizen scientists across the globe. Here, we summarize the scientific and technical advances this perspective has enabled. As the project's name implies, the early years of Folding@home focused on driving advances in our understanding of protein folding by developing statistical methods for capturing long-timescale processes and facilitating insight into complex dynamical processes. Success laid a foundation for broadening the scope of Folding@home to address other functionally relevant conformational changes, such as receptor signaling, enzyme dynamics, and ligand binding. Continued algorithmic advances, hardware developments such as GPU-based computing, and the growing scale of Folding@home have enabled the project to focus on new areas where massively parallel sampling can be impactful. While previous work sought to expand toward larger proteins with slower conformational changes, new work focuses on large-scale comparative studies of different protein sequences and chemical compounds to better understand biology and inform the development of small molecule drugs. Progress on these fronts enabled the community to pivot quickly in response to the COVID-19 pandemic, expanding to become the world's first exascale computer and deploying this massive resource to provide insight into the inner workings of the SARS-CoV-2 virus and aid the development of new antivirals. This success provides a glimpse of what's to come as exascale supercomputers come online, and Folding@home continues its work.

Folding@home: Achievements from over 20 years of citizen science herald the exascale era.

Simulations of biomolecules have enormous potential to inform our understanding of biology but require extremely demanding calculations. For over 20 years, the Folding@home distributed computing project has pioneered a massively parallel approach to biomolecular simulation, harnessing the resources of citizen scientists across the globe. Here, we summarize the scientific and technical advances this perspective has enabled. As the project's name implies, the early years of Folding@home focused on driving advances in our understanding of protein folding by developing statistical methods for capturing long-timescale processes and facilitating insight into complex dynamical processes. Success laid a foundation for broadening the scope of Folding@home to address other functionally relevant conformational changes, such as receptor signaling, enzyme dynamics, and ligand binding. Continued algorithmic advances, hardware developments such as graphics processing unit (GPU)-based computing, and the growing scale of Folding@home have enabled the project to focus on new areas where massively parallel sampling can be impactful. While previous work sought to expand toward larger proteins with slower conformational changes, new work focuses on large-scale comparative studies of different protein sequences and chemical compounds to better understand biology and inform the development of small-molecule drugs. Progress on these fronts enabled the community to pivot quickly in response to the COVID-19 pandemic, expanding to become the world's first exascale computer and deploying this massive resource to provide insight into the inner workings of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) virus and aid the development of new antivirals. This success provides a glimpse of what is to come as exascale supercomputers come online and as Folding@home continues its work.

Effects of Salinity and Temperature on the Flexibility and Function of a Polyextremophilic Enzyme.

The polyextremophilic ß-galactosidase enzyme of the haloarchaeon Halorubrum lacusprofundi functions in extremely cold and hypersaline conditions. To better understand the basis of polyextremophilic activity, the enzyme was studied using steady-state kinetics and molecular dynamics at temperatures ranging from 10 °C to 50 °C and salt concentrations from 1 M to 4 M KCl. Kinetic analysis showed that while catalytic efficiency (kcat/Km) improves with increasing temperature and salinity, Km is reduced with decreasing temperatures and increasing salinity, consistent with improved substrate binding at low temperatures. In contrast, kcat was similar from 2-4 M KCl across the temperature range, with the calculated enthalpic and entropic components indicating a threshold of 2 M KCl to lower the activation barrier for catalysis. With molecular dynamics simulations, the increase in per-residue root-mean-square fluctuation (RMSF) was observed with higher temperature and salinity, with trends like those seen with the catalytic efficiency, consistent with the enzyme's function being related to its flexibility. Domain A had the smallest change in flexibility across the conditions tested, suggesting the adaptation to extreme conditions occurs via regions distant to the active site and surface accessible residues. Increased flexibility was most apparent in the distal active sites, indicating their importance in conferring salinity and temperature-dependent effects.

Oncogenic Mutations in the DNA-Binding Domain of FOXO1 that Disrupt Folding: Quantitative Insights from Experiments and Molecular Simulations.

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.

Estimation of binding rates and affinities from multiensemble Markov models and ligand decoupling.

Accurate and efficient simulation of the thermodynamics and kinetics of protein-ligand interactions is crucial for computational drug discovery. Multiensemble Markov Model (MEMM) estimators can provide estimates of both binding rates and affinities from collections of short trajectories but have not been systematically explored for situations when a ligand is decoupled through scaling of non-bonded interactions. In this work, we compare the performance of two MEMM approaches for estimating ligand binding affinities and rates: (1) the transition-based reweighting analysis method (TRAM) and (2) a Maximum Caliber (MaxCal) based method. As a test system, we construct a small host-guest system where the ligand is a single uncharged Lennard-Jones (LJ) particle, and the receptor is an 11-particle icosahedral pocket made from the same atom type. To realistically mimic a protein-ligand binding system, the LJ ϵ parameter was tuned, and the system was placed in a periodic box with 860 TIP3P water molecules. A benchmark was performed using over 80 µs of unbiased simulation, and an 18-state Markov state model was used to estimate reference binding affinities and rates. We then tested the performance of TRAM and MaxCal when challenged with limited data. Both TRAM and MaxCal approaches perform better than conventional Markov state models, with TRAM showing better convergence and accuracy. We find that subsampling of trajectories to remove time correlation improves the accuracy of both TRAM and MaxCal and that in most cases, only a single biased ensemble to enhance sampled transitions is required to make accurate estimates.

Unprotected peptide macrocyclization and stapling via a fluorine-thiol displacement reaction.

We report the discovery of a facile peptide macrocyclization and stapling strategy based on a fluorine thiol displacement reaction (FTDR), which renders a class of peptide analogues with enhanced stability, affinity, cellular uptake, and inhibition of cancer cells. This approach enabled selective modification of the orthogonal fluoroacetamide side chains in unprotected peptides in the presence of intrinsic cysteines. The identified benzenedimethanethiol linker greatly promoted the alpha helicity of a variety of peptide substrates, as corroborated by molecular dynamics simulations. The cellular uptake of benzenedimethanethiol stapled peptides appeared to be universally enhanced compared to the classic ring-closing metathesis (RCM) stapled peptides. Pilot mechanism studies suggested that the uptake of FTDR-stapled peptides may involve multiple endocytosis pathways in a distinct pattern in comparison to peptides stapled by RCM. Consistent with the improved cell permeability, the FTDR-stapled lead Axin and p53 peptide analogues demonstrated enhanced inhibition of cancer cells over the RCM-stapled analogues and the unstapled peptides.

Expanded Ensemble Methods Can be Used to Accurately Predict Protein-Ligand Relative Binding Free Energies.

Alchemical free energy methods have become indispensable in computational drug discovery for their ability to calculate highly accurate estimates of protein-ligand affinities. Expanded ensemble (EE) methods, which involve single simulations visiting all of the alchemical intermediates, have some key advantages for alchemical free energy calculation. However, there have been relatively few examples published in the literature of using expanded ensemble simulations for free energies of protein-ligand binding. In this paper, as a test of expanded ensemble methods, we compute relative binding free energies using the Open Force Field Initiative force field (codename "Parsley") for 24 pairs of Tyk2 inhibitors derived from a congeneric series of 16 compounds. The EE predictions agree well with the experimental values (root-mean-square error (RMSE) of 0.94 ± 0.13 kcal mol-1 and mean unsigned error (MUE) of 0.75 ± 0.12 kcal mol-1). We find that while increasing the number of alchemical intermediates can improve the phase space overlap, faster convergence can be obtained with fewer intermediates, as long as acceptance rates are sufficient. We also find that convergence can be improved using more aggressive updating of biases, and that estimates can be improved by performing multiple independent EE calculations. This work demonstrates that EE is a viable option for alchemical free energy calculation. We discuss the implications of these findings for rational drug design, as well as future directions for improvement.

Stacking Gaussian processes to improve [Formula: see text] predictions in the SAMPL7 challenge.

Accurate predictions of acid dissociation constants are essential to rational molecular design in the pharmaceutical industry and elsewhere. There has been much interest in developing new machine learning methods that can produce fast and accurate pKa predictions for arbitrary species, as well as estimates of prediction uncertainty. Previously, as part of the SAMPL6 community-wide blind challenge, Bannan et al. approached the problem of predicting [Formula: see text]s by using a Gaussian process regression to predict microscopic [Formula: see text]s, from which macroscopic [Formula: see text] values can be analytically computed (Bannan et al. in J Comput-Aided Mol Des 32:1165-1177). While this method can make reasonably quick and accurate predictions using a small training set, accuracy was limited by the lack of a sufficiently broad range of chemical space in the training set (e.g., the inclusion of polyprotic acids). Here, to address this issue, we construct a deep Gaussian Process (GP) model that can include more features without invoking the curse of dimensionality. We trained both a standard GP and a deep GP model using a database of approximately 3500 small molecules curated from public sources, filtered by similarity to targets. We tested the model on both the SAMPL6 and more recent SAMPL7 challenge, which introduced a similar lack of ionizable sites and/or environments found between the test set and the previous training set. The results show that while the deep GP model made only minor improvements over the standard GP model for SAMPL6 predictions, it made significant improvements over the standard GP model in SAMPL7 macroscopic predictions, achieving a MAE of 1.5 [Formula: see text].

Metal Cation-Binding Mechanisms of Q-Proline Peptoid Macrocycles in Solution.

The rational design of foldable and functionalizable peptidomimetic scaffolds requires the concerted application of both computational and experimental methods. Recently, a new class of designed peptoid macrocycle incorporating spiroligomer proline mimics (Q-prolines) has been found to preorganize when bound by monovalent metal cations. To determine the solution-state structure of these cation-bound macrocycles, we employ a Bayesian inference method (BICePs) to reconcile enhanced-sampling molecular simulations with sparse ROESY correlations from experimental NMR studies to predict and design conformational and binding properties of macrocycles as functional scaffolds for peptidomimetics. Conformations predicted to be most populated in solution were then simulated in the presence of explicit cations to yield trajectories with observed binding events, revealing a highly preorganized all-trans amide conformation, whose formation is likely limited by the slow rate of cis/trans isomerization. Interestingly, this conformation differs from a racemic crystal structure solved in the absence of cation. Free energies of cation binding computed from distance-dependent potentials of mean force suggest Na+ has a higher affinity to the macrocycle than K+, with both cations binding much more strongly in acetonitrile than water. The simulated affinities are able to correctly rank the extent to which different macrocycle sequences exhibit preorganization in the presence of different metal cations and solvents, suggesting our approach is suitable for solution-state computational design.

Assigning confidence to molecular property prediction.

Introduction: Computational modeling has rapidly advanced over the last decades. Recently, machine learning has emerged as a powerful and cost-effective strategy to learn from existing datasets and perform predictions on unseen molecules. Accordingly, the explosive rise of data-driven techniques raises an important question: What confidence can be assigned to molecular property predictions and what techniques can be used?Areas covered: The authors discuss popular strategies for predicting molecular properties, their corresponding uncertainty sources and methods to quantify uncertainty. First, the authors' considerations for assessing confidence begin with dataset bias and size, data-driven property prediction and feature design. Next, the authors discuss property simulation via computations of binding affinity in detail. Lastly, they investigate how these uncertainties propagate to generative models, as they are usually coupled with property predictors.Expert opinion: Computational techniques are paramount to reduce the prohibitive cost of brute-force experimentation during exploration. The authors believe that assessing uncertainty in property prediction models is essential whenever closed-loop drug design campaigns relying on high-throughput virtual screening are deployed. Accordingly, considering sources of uncertainty leads to better-informed validations, more reliable predictions and more realistic expectations of the entire workflow. Overall, this increases confidence in the predictions and, ultimately, accelerates drug design.

SARS-CoV-2 simulations go exascale to predict dramatic spike opening and cryptic pockets across the proteome.

SARS-CoV-2 has intricate mechanisms for initiating infection, immune evasion/suppression and replication that depend on the structure and dynamics of its constituent proteins. Many protein structures have been solved, but far less is known about their relevant conformational changes. To address this challenge, over a million citizen scientists banded together through the Folding@home distributed computing project to create the first exascale computer and simulate 0.1 seconds of the viral proteome. Our adaptive sampling simulations predict dramatic opening of the apo spike complex, far beyond that seen experimentally, explaining and predicting the existence of 'cryptic' epitopes. Different spike variants modulate the probabilities of open versus closed structures, balancing receptor binding and immune evasion. We also discover dramatic conformational changes across the proteome, which reveal over 50 'cryptic' pockets that expand targeting options for the design of antivirals. All data and models are freely available online, providing a quantitative structural atlas.

Reconciling Simulations and Experiments With BICePs: A Review.

Bayesian Inference of Conformational Populations (BICePs) is an algorithm developed to reconcile simulated ensembles with sparse experimental measurements. The Bayesian framework of BICePs enables population reweighting as a post-simulation processing step, with several advantages over existing methods, including the proper use of reference potentials, and the estimation of a Bayes factor-like quantity called the BICePs score for model selection. Here, we summarize the theory underlying this method in context with related algorithms, review the history of BICePs applications to date, and discuss current shortcomings along with future plans for improvement.

Solution-State Preorganization of Cyclic ß-Hairpin Ligands Determines Binding Mechanism and Affinities for MDM2.

Understanding mechanisms of protein folding and binding is crucial to designing their molecular function. Molecular dynamics (MD) simulations and Markov state model (MSM) approaches provide a powerful way to understand complex conformational change that occurs over long time scales. Such dynamics are important for the design of therapeutic peptidomimetic ligands, whose affinity and binding mechanism are dictated by a combination of folding and binding. To examine the role of preorganization in peptide binding to protein targets, we performed massively parallel explicit-solvent MD simulations of cyclic ß-hairpin ligands designed to mimic the p53 transactivation domain and competitively bind mouse double minute 2 homologue (MDM2). Disrupting the MDM2-p53 interaction is a therapeutic strategy to prevent degradation of the p53 tumor suppressor in cancer cells. MSM analysis of over 3 ms of aggregate trajectory data enabled us to build a detailed mechanistic model of coupled folding and binding of four cyclic peptides which we compare to experimental binding affinities and rates. The results show a striking relationship between the relative preorganization of each ligand in solution and its affinity for MDM2. Specifically, changes in peptide conformational populations predicted by the MSMs suggest that entropy loss upon binding is the main factor influencing affinity. The MSMs also enable detailed examination of non-native interactions which lead to misfolded states and comparison of structural ensembles with experimental NMR measurements. In contrast to an MSM study of p53 transactivation domain (TAD) binding to MDM2, MSMs of cyclic ß-hairpin binding show a conformational selection mechanism. Finally, we make progress toward predicting accurate off rates of cyclic peptides using multiensemble Markov models (MEMMs) constructed from unbiased and biased simulated trajectories.

Markov State Models to Elucidate Ligand Binding Mechanism.

Molecular dynamics simulations can now routinely access the microsecond timescale, making feasible direct sampling of ligand association events. While Markov State Model (MSM) approaches offer a useful framework for analyzing such trajectory data to gain insight into binding mechanisms, accurate modeling of ligand association pathways and kinetics must be done carefully. We describe methods and good practices for constructing MSMs of ligand binding from unbiased trajectory data and discuss how to use time-lagged independent component analysis (tICA) to build informative models, using as an example recent simulation work to model the binding of phenylalanine to the regulatory ACT domain dimer of phenylalanine hydroxylase. We describe a variety of methods for estimating association rates from MSMs and discuss how to distinguish between conformational selection and induced-fit mechanisms using MSMs. In addition, we review some examples of MSMs constructed to elucidate the mechanisms by which p53 transactivation domain (TAD) and related peptides bind the oncoprotein MDM2.