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 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.

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 (k_cat/K_m) improves with increasing temperature and salinity, K_m is reduced with decreasing temperatures and increasing salinity, consistent with improved substrate binding at low temperatures. In contrast, k_cat 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.

Strategies capable of stapling unprotected peptides in a straightforward, chemoselective, and clean manner, as well as promoting cellular uptake are of great interest. Here the authors report a peptide macrocyclization and stapling strategy which satisfies those criteria, based on a fluorine thiol displacement reaction.

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].

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.

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.

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.

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.

Metal-Binding Q-Proline Macrocycles

We introduce the efficient Fmoc-SPPS and peptoid synthesis of Q-proline-based, metal-binding macrocycles (QPMs), which bind metal cations and display nine functional groups. Metal-free QPMs are disordered, evidenced by NMR and a crystal structure of QPM-3 obtained through racemic crystallization. Upon addition of metal cations, QPMs adopt ordered structures. Notably, the addition of a second functional group at the hydantoin amide position (R₂) converts the proline ring from Cγ-endo to Cγ-exo, due to steric interactions.

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.

SARS-CoV-2 Simulations Go Exascale to Capture Spike Opening and Reveal Cryptic Pockets Across the Proteome

SARS-CoV-2 has intricate mechanisms for initiating infection, immune evasion/suppression, and replication, which 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 an unprecedented 0.1 seconds of the viral proteome. Our simulations capture dramatic opening of the apo Spike complex, far beyond that seen experimentally, which explains and successfully predicts the existence of 'cryptic' epitopes. Different Spike homologues modulate the probabilities of open versus closed structures, balancing receptor binding and immune evasion. We also observe 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.

Microcanonical coarse-graining of the kinetic Ising model

We propose a scheme for coarse-graining the dynamics of the 2-D kinetic Ising model onto the microcanonical ensemble. At subcritical temperatures, 2-D and higher-dimensional Ising lattices possess two basins of attraction separated by a free energy barrier. Projecting onto the microcanonical ensemble has the advantage that the dependence of the crossing rate constant on environmental conditions can be obtained from a single Monte Carlo trajectory. Using various numerical methods, we computed the forward rate constants of coarse-grained representations of the Ising model and compared them with the true value obtained from brute force simulation. While coarse-graining preserves detailed balance, the computed rate constants for barrier heights between 5 kT and 9 kT were consistently 50% larger than the true value. Markovianity testing revealed loss of dynamical memory, which we propose accounts for coarse-graining error. Committor analysis did not support the alternative hypothesis that microcanonical projection is incompatible with an optimal reaction coordinate. The correct crossing rate constant was obtained by spectrally decomposing the diffusion coefficient near the free energy barrier and selecting the slowest (reactive) component. The spectral method also yielded the correct rate constant in the 3-D Ising lattice, where coarse-graining error was 6% and memory effects were diminished. We conclude that microcanonical coarse-graining supplemented by spectral analysis of short-term barrier fluctuations provides a comprehensive kinetic description of barrier crossing in a non-inertial continuous-time jump process.

Reconciling Simulated Ensembles of Apomyoglobin with Experimental Hydrogen/Deuterium Exchange Data Using Bayesian Inference and Multiensemble Markov State Models

Hydrogen/deuterium exchange (HDX) is a powerful technique to investigate protein conformational dynamics at amino acid resolution. Because HDX provides a measurement of solvent exposure of backbone hydrogens, ensemble-averaged over potentially slow kinetic processes, it has been challenging to use HDX protection factors to refine structural ensembles obtained from molecular dynamics simulations. This entails dual challenges: (1) identifying structural observables that best correlate with backbone amide protection from exchange and (2) restraining these observables in molecular simulations to model ensembles consistent with experimental measurements. Here, we make significant progress on both fronts. First, we describe an improved predictor of HDX protection factors from structural observables in simulated ensembles, parametrized from ultralong molecular dynamics simulation trajectory data, with a Bayesian inference approach used to retain the full posterior distribution of model parameters. We next present a new method for obtaining simulated ensembles in agreement with experimental HDX protection factors, in which molecular simulations are performed at various temperatures and restraint biases and used to construct multiensemble Markov State Models (MSMs). Finally, the BICePs (Bayesian Inference of Conformational Populations) algorithm is then used with our HDX protection factor predictor to infer which thermodynamic ensemble agrees best with the experiment and estimate populations of each conformational state in the MSM. To illustrate the approach, we use a combination of HDX protection factor restraints and chemical shift restraints to model the conformational ensemble of apomyoglobin at pH 6. The resulting ensemble agrees well with the experiment and gives insight into the all-atom structure of disordered helices F and H in the absence of heme.

Adaptive Markov state model estimation using short reseeding trajectories

In the last decade, advances in molecular dynamics (MD) and Markov State Model (MSM) methodologies have made possible accurate and efficient estimation of kinetic rates and reactive pathways for complex biomolecular dynamics occurring on slow time scales. A promising approach to enhanced sampling of MSMs is to use "adaptive" methods, in which new MD trajectories are "seeded" preferentially from previously identified states. Here, we investigate the performance of various MSM estimators applied to reseeding trajectory data, for both a simple 1D free energy landscape and mini-protein folding MSMs of WW domain and NTL9(1-39). Our results reveal the practical challenges of reseeding simulations and suggest a simple way to reweight seeding trajectory data to better estimate both thermodynamic and kinetic quantities.

Site-Specific Immuno-PET Tracer to Image PD-L1

The rapid ascension of immune checkpoint blockade treatments has placed an emphasis on the need for viable, robust, and noninvasive imaging methods for immune checkpoint proteins, which could be of diagnostic value. Immunoconjugate-based positron emission tomography (immuno-PET) allows for sensitive and quantitative imaging of target levels and has promising potential for the noninvasive evaluation of immune checkpoint proteins. However, the advancement of immuno-PET is currently limited by available imaging tools, which heavily rely on full-length IgGs with Fc-mediated effects and are heterogeneous mixtures upon random conjugation with chelators for imaging. Herein, we have developed a site-specific αPD-L1 Fab conjugate with the chelator 1,4,7-triazacyclononane- N, N', N″-triacetic acid (NOTA), enabling radiolabeling for PET imaging, using the amber suppression-mediated genetic incorporation of unnatural amino acid (UAA), p-azidophenylalanine. This Fab conjugate is homogeneous and demonstrated tight binding toward the PD-L1 antigen in vitro. The radiolabeled version, ⁶⁴Cu-NOTA-αPD-L1, has been employed in PET imaging to allow for effective visualization and mapping of the biodistribution of PD-L1 in two normal mouse models, including the capturing of different PD-L1 expression levels in the spleens of the different mouse types. Follow-up in vivo blocking studies and ex vivo fluorescent staining further validated specific tissue uptakes of the imaging agent. This approach illustrates the utility of UAA-based site-specific Fab conjugation as a general strategy for making sensitive PET imaging probes, which could facilitate the elucidation of the roles of a wide variety of immune checkpoint proteins in immunotherapy.

Fluorinated Aromatic Monomers as Building Blocks To Control α-Peptoid Conformation and Structure

Peptoids are peptidomimetics of interest in the fields of drug development and biomaterials. However, obtaining stable secondary structures is challenging, and designing these requires effective control of the peptoid tertiary amide cis/trans equilibrium. Herein, we report new fluorine-containing aromatic monomers that can control peptoid conformation. Specifically, we demonstrate that a fluoro-pyridine group can be used to circumvent the need for monomer chirality to control the cis/trans equilibrium. We also show that incorporation of a trifluoro-methyl group ( N^CF₃Rpe) rather than a methyl group ( NRpe) at the α-carbon of a monomer gives rise to a 5-fold increase in cis-isomer preference.

Exposing the Nucleation Site in α-Helix Folding: A Joint Experimental and Simulation Study

One of the fundamental events in protein folding is α-helix formation, which involves sequential development of a series of helical hydrogen bonds between the backbone C═O group of residues i and the -NH group of residues i + 4. While we now know a great deal about α-helix folding dynamics, a key question that remains to be answered is where the productive helical nucleation event occurs. Statistically, a helical nucleus (or the first helical hydrogen-bond) can form anywhere within the peptide sequence in question; however, the one that leads to productive folding may only form at a preferred location. This consideration is based on the fact that the α-helical structure is inherently asymmetric, due to the specific alignment of the helical hydrogen bonds. While this hypothesis is plausible, validating it is challenging because there is not an experimental observable that can be used to directly pinpoint the location of the productive nucleation process. Therefore, in this study we combine several techniques, including peptide cross-linking, laser-induced temperature-jump infrared spectroscopy, and molecular dynamics simulations, to tackle this challenge. Taken together, our experimental and simulation results support an α-helix folding mechanism wherein the productive nucleus is formed at the N-terminus, which propagates toward the C-terminal end of the peptide to yield the folded structure. In addition, our results show that incorporation of a cross-linker can lead to formation of differently folded conformations, underscoring the need for all-atom simulations to quantitatively assess the proposed cross-linking design.

Simulations of the regulatory ACT domain of human phenylalanine hydroxylase (PAH) unveil its mechanism of phenylalanine binding

Phenylalanine hydroxylase (PAH) regulates phenylalanine (Phe) levels in mammals to prevent neurotoxicity resulting from high Phe concentrations as observed in genetic disorders leading to hyperphenylalaninemia and phenylketonuria. PAH senses elevated Phe concentrations by transient allosteric Phe binding to a protein-protein interface between ACT domains of different subunits in a PAH tetramer. This interface is present in an activated PAH (A-PAH) tetramer and absent in a resting-state PAH (RS-PAH) tetramer. To investigate this allosteric sensing mechanism, here we used the GROMACS molecular dynamics simulation suite on the Folding@home computing platform to perform extensive molecular simulations and Markov state model (MSM) analysis of Phe binding to ACT domain dimers. These simulations strongly implicated a conformational selection mechanism for Phe association with ACT domain dimers and revealed protein motions that act as a gating mechanism for Phe binding. The MSMs also illuminate a highly mobile hairpin loop, consistent with experimental findings also presented here that the PAH variant L72W does not shift the PAH structural equilibrium toward the activated state. Finally, simulations of ACT domain monomers are presented, in which spontaneous transitions between resting-state and activated conformations are observed, also consistent with a mechanism of conformational selection. These mechanistic details provide detailed insight into the regulation of PAH activation and provide testable hypotheses for the development of new allosteric effectors to correct structural and functional defects in PAH.

Elucidating the inhibition of peptidoglycan biosynthesis in Staphylococcus aureus by albocycline, a macrolactone isolated from Streptomyces maizeus

Antibiotic resistance is a serious threat to global public health, and methicillin-resistant Staphylococcus aureus (MRSA) is a poignant example. The macrolactone natural product albocycline, derived from various Streptomyces strains, was recently identified as a promising antibiotic candidate for the treatment of both MRSA and vancomycin-resistant S. aureus (VRSA), which is another clinically relevant and antibiotic resistant strain. Moreover, it was hypothesized that albocycline's antimicrobial activity was derived from the inhibition of peptidoglycan (i.e., bacterial cell wall) biosynthesis. Herein, preliminary mechanistic studies are performed to test the hypothesis that albocycline inhibits MurA, the enzyme that catalyzes the first step of peptidoglycan biosynthesis, using a combination of biological assays alongside molecular modeling and simulation studies. Computational modeling suggests albocycline exists as two conformations in solution, and computational docking of these conformations to an ensemble of simulated receptor structures correctly predicted preferential binding to S. aureus MurA-the enzyme that catalyzes the first step of peptidoglycan biosynthesis-over Escherichia coli (E. coli) MurA. Albocycline isolated from the producing organism (Streptomyces maizeus) weakly inhibited S. aureus MurA (IC50 of 480 μM) but did not inhibit E. coli MurA. The antimicrobial activity of albocycline against resistant S. aureus strains was superior to that of vancomycin, preferentially inhibiting Gram-positive organisms. Albocycline was not toxic to human HepG2 cells in MTT assays. While these studies demonstrate that albocycline is a promising lead candidate against resistant S. aureus, taken together they suggest that MurA is not the primary target, and further work is necessary to identify the major biological target.

Model Selection Using BICePs: A Bayesian Approach for Force Field Validation and Parameterization

The Bayesian Inference of Conformational Populations (BICePs) algorithm reconciles theoretical predictions of conformational state populations with sparse and/or noisy experimental measurements. Among its key advantages is its ability to perform objective model selection through a quantity we call the BICePs score, which reflects the integrated posterior evidence in favor of a given model, computed through free energy estimation methods. Here, we explore how the BICePs score can be used for force field validation and parametrization. Using a 2D lattice protein as a toy model, we demonstrate that BICePs is able to select the correct value of an interaction energy parameter given ensemble-averaged experimental distance measurements. We show that if conformational states are sufficiently fine-grained, the results are robust to experimental noise and measurement sparsity. Using these insights, we apply BICePs to perform force field evaluations for all-atom simulations of designed β-hairpin peptides against experimental NMR chemical shift measurements. These tests suggest that BICePs scores can be used for model selection in the context of all-atom simulations. We expect this approach to be particularly useful for the computational foldamer design as a tool for improving general-purpose force fields given sparse experimental measurements.

Bridging Microscopic and Macroscopic Mechanisms of p53-MDM2 Binding with Kinetic Network Models

Under normal cellular conditions, the tumor suppressor protein p53 is kept at low levels in part due to ubiquitination by MDM2, a process initiated by binding of MDM2 to the intrinsically disordered transactivation domain (TAD) of p53. Many experimental and simulation studies suggest that disordered domains such as p53 TAD bind their targets nonspecifically before folding to a tightly associated conformation, but the microscopic details are unclear. Toward a detailed prediction of binding mechanisms, pathways, and rates, we have performed large-scale unbiased all-atom simulations of p53-MDM2 binding. Markov state models (MSMs) constructed from the trajectory data predict p53 TAD binding pathways and on-rates in good agreement with experiment. The MSM reveals that two key bound intermediates, each with a nonnative arrangement of hydrophobic residues in the MDM2 binding cleft, control the overall on-rate. Using microscopic rate information from the MSM, we parameterize a simple four-state kinetic model to 1) determine that induced-fit pathways dominate the binding flux over a large range of concentrations, and 2) predict how modulation of residual p53 helicity affects binding, in good agreement with experiment. These results suggest new ways in which microscopic models of peptide binding, coupled with simple few-state binding flux models, can be used to understand biological function in physiological contexts.

Precisely tuneable energy transfer system using peptoid helix-based molecular scaffold

The energy flow during natural photosynthesis is controlled by maintaining the spatial arrangement of pigments, employing helices as scaffolds. In this study, we have developed porphyrin-peptoid (pigment-helix) conjugates (PPCs) that can modulate the donor-acceptor energy transfer efficiency with exceptional precision by controlling the relative distance and orientation of the two pigments. Five donor-acceptor molecular dyads were constructed using zinc porphyrin and free base porphyrin (Zn(i + 2)–Zn(i + 6)), and highly efficient energy transfer was demonstrated with estimated efficiencies ranging from 92% to 96% measured by static fluorescence emission in CH₂Cl₂ and from 96.3% to 97.6% using femtosecond transient absorption measurements in toluene, depending on the relative spatial arrangement of the donor-acceptor pairs. Our results suggest that the remarkable precision and tunability exhibited by nature can be achieved by mimicking the design principles of natural photosynthetic proteins.

Computational and Experimental Evaluation of Designed β-Cap Hairpins Using Molecular Simulations and Kinetic Network Models

Molecular simulation has been used to model the detailed folding properties of peptides, yet prospective computational peptide design by such approaches remains challenging and nontrivial. To test the accuracy of simulation-based hairpin design, we characterized the folding properties of a series of so-called β-cap hairpin peptides designed to mimic a conserved hairpin of LapD, a bacterial intracellular signaling protein, both experimentally by NMR spectroscopy and computationally by implicit-solvent replica-exchange molecular dynamics using three different AMBER force fields (ff96, ff99sb-ildn, and ff99sb-ildn-NMR). A unique challenge presented by these designs is the presence of both a terminal Trp-Trp capping motif and a conserved GWxQ motif in the hairpin turn required for binding to LapG. Consistent with previous studies, we found AMBER ff96 to be the most accurate when used with the OBC GBSA implicit solvent model, despite its known bias toward β-sheet conformations when used in explicit-solvent simulations. To gain microscopic insight into the folding landscape of the hairpin designs, we additionally performed parallel simulations on the Folding@home distributed computing platform using AMBER ff99sb-ildn-NMR with TIP3P explicit solvent. Markov state models (MSMs) built from trajectory data reveal a number of non-native interactions between Trp and other amino acid side chains, creating potential problems in achieving well-folded hairpin structures in solution.

Molecular simulations and free-energy calculations suggest conformation-dependent anion binding to a cytoplasmic site as a mechanism for Na+/K+-ATPase ion selectivity

Na⁺/K⁺-ATPase transports Na⁺ and K⁺ ions across the cell membrane via an ion-binding site becoming alternatively accessible to the intra- and extracellular milieu by conformational transitions that confer marked changes in ion-binding stoichiometry and selectivity. To probe the mechanism of these changes, we used molecular simulation and free-energy perturbation approaches to identify probable protonation states of Na⁺- and K⁺-coordinating residues in E1P and E2P conformations of Na⁺/K⁺-ATPase. Analysis of these simulations revealed a molecular mechanism responsible for the change in protonation state: the conformation-dependent binding of an anion (a chloride ion in our simulations) to a previously unrecognized cytoplasmic site in the loop between transmembrane helices 8 and 9, which influences the electrostatic potential of the crucial Na⁺-coordinating residue Asp⁹²⁶ This mechanistic model is consistent with experimental observations and provides a molecular-level picture of how E1P to E2P enzyme conformational transitions are coupled to changes in ion-binding stoichiometry and selectivity.

Diverted Total Synthesis of Carolacton-Inspired Analogs Yields Three Distinct Phenotypes in Streptococcus mutans Biofilms

The oral microbiome is a dynamic environment inhabited by both commensals and pathogens. Among these is Streptococcus mutans, the causative agent of dental caries, the most prevalent childhood disease. Carolacton has remarkably specific activity against S. mutans, causing acid-mediated cell death during biofilm formation; however, its complex structure limits its utility. Herein, we report the diverted total synthesis and biological evaluation of a rationally designed library of simplified analogs that unveiled three unique biofilm phenotypes further validating the role of natural product synthesis in the discovery of new biological phenomena.

Control of porphyrin interactions via structural changes of a peptoid scaffold

A template to control porphyrin interactions is constructed by displaying porphyrins at defined positions on a helical peptoid.

A Maximum-Caliber Approach to Predicting Perturbed Folding Kinetics Due to Mutations

We present a maximum-caliber method for inferring transition rates of a Markov state model (MSM) with perturbed equilibrium populations given estimates of state populations and rates for an unperturbed MSM. It is similar in spirit to previous approaches, but given the inclusion of prior information, it is more robust and simple to implement. We examine its performance in simple biased diffusion models of kinetics and then apply the method to predicting changes in folding rates for several highly nontrivial protein folding systems for which non-native interactions play a significant role, including (1) tryptophan variants of the GB1 hairpin, (2) salt-bridge mutations of the Fs peptide helix, and (3) MSMs built from ultralong folding trajectories of FiP35 and GTT variants of the WW domain. In all cases, the method correctly predicts changes in folding rates, suggesting the wide applicability of maximum-caliber approaches to efficiently predict how mutations perturb protein conformational dynamics.

Markov models of the apo-MDM2 lid region reveal diffuse yet two-state binding dynamics and receptor poses for computational docking

MDM2 is a negative regulator of p53 activity and an important target for cancer therapeutics. The N-terminal lid region of MDM2 modulates interactions with p53 via competition for its binding cleft, exchanging slowly between docked and undocked conformations in the absence of p53. To better understand these dynamics, we constructed Markov State Models (MSMs) from large collections of unbiased simulation trajectories of apo-MDM2, and find strong evidence for diffuse, yet two-state folding and binding of the N-terminal region to the p53 receptor site. The MSM also identifies holo-like receptor conformations highly suitable for computational docking, despite initiating trajectories from closed-cleft receptor structures unsuitable for docking. Fixed-anchor docking studies using a test set of high-affinity small molecules and peptides show simulated receptor ensembles achieve docking successes comparable to cross-docking studies using crystal structures of receptors bound by alternative ligands. For p53, the best-scoring receptor structures have the N-terminal region lid region bound in a helical conformation mimicking the bound structure of p53, suggesting lid region association induces receptor conformations suitable for binding. These results suggest that MD + MSM approaches can sample binding-competent receptor conformations suitable for computational peptidomimetic design, and that inclusion of disordered regions may be essential to capturing the correct receptor dynamics.

Using Kinetic Network Models To Probe Non-Native Salt-Bridge Effects on α-Helix Folding

Salt-bridge interactions play an important role in stabilizing many protein structures, and have been shown to be designable features for protein design. In this work, we study the effects of non-native salt bridges on the folding of a soluble alanine-based peptide (Fs peptide) using extensive all-atom molecular dynamics simulations performed on the Folding@home distributed computing platform. Using Markov State Models, we show how non-native salt-bridges affect the folding kinetics of Fs peptide by perturbing specific conformational states. Furthermore, we present methods for the automatic detection and analysis of such states. These results provide insight into helix folding mechanisms and useful information to guide simulation-based computational protein design.

Insights into Peptoid Helix Folding Cooperativity from an Improved Backbone Potential

Peptoids (N-substituted oligoglycines) are biomimetic polymers that can fold into a variety of unique structural scaffolds. Peptoid helices, which result from the incorporation of bulky chiral side chains, are a key peptoid structural motif whose formation has not yet been accurately modeled in molecular simulations. Here, we report that a simple modification of the backbone φ-angle potential in GAFF is able to produce well-folded cis-amide helices of (S)-N-(1-phenylethyl)glycine (Nspe), consistent with experiment. We validate our results against both QM calculations and NMR experiments. For this latter task, we make quantitative comparisons to sparse NOE data using the Bayesian Inference of Conformational Populations (BICePs) algorithm, a method we have recently developed for this purpose. We then performed extensive REMD simulations of Nspe oligomers as a function of chain length and temperature to probe the molecular forces driving cooperative helix formation. Analysis of simulation data by Lifson-Roig helix-coil theory show that the modified potential predicts much more cooperative folding for Nspe helices. Unlike peptides, per-residue entropy changes for helix nucleation and extension are mostly positive, suggesting that steric bulk provides the main driving force for folding. We expect these results to inform future work aimed at predicting and designing peptoid peptidomimetics and tertiary assemblies of peptoid helices.

Molecular Simulation of Conformational Pre-Organization in Cyclic RGD Peptides

To test the ability of molecular simulations to accurately predict the solution-state conformational properties of peptidomimetics, we examined a test set of 18 cyclic RGD peptides selected from the literature, including the anticancer drug candidate cilengitide, whose favorable binding affinity to integrin has been ascribed to its pre-organization in solution. For each design, we performed all-atom replica-exchange molecular dynamics simulations over several microseconds and compared the results to extensive published NMR data. We find excellent agreement with experimental NOE distance restraints, suggesting that molecular simulation can be a useful tool for the computational design of pre-organized solution-state structure. Moreover, our analysis of conformational populations estimates that, despite the potential for increased flexibility due to backbone amide isomerizaton, N-methylation provides about 0.5 kcal/mol of reduced conformational entropy to cyclic RGD peptides. The combination of pre-organization and binding-site compatibility explains the strong binding affinity of cilengitide to integrin.

A molecular interpretation of 2D IR protein folding experiments with Markov state models

The folding mechanism of the N-terminal domain of ribosomal protein L9 (NTL91-39) is studied using temperature-jump (T-jump) amide I' two-dimensional infrared (2D IR) spectroscopy in combination with spectral simulations based on a Markov state model (MSM) built from millisecond-long molecular dynamics trajectories. The results provide evidence for a compact well-structured folded state and a heterogeneous fast-exchanging denatured state ensemble exhibiting residual secondary structure. The folding rate of 26.4 μs(-1) (at 80°C), extracted from the T-jump response of NTL91-39, compares favorably with the 18 μs(-1) obtained from the MSM. Structural decomposition of the MSM and analysis along the folding coordinate indicates that helix-formation nucleates the global folding. Simulated difference spectra, corresponding to the global folding transition of the MSM, are in qualitative agreement with measured T-jump 2D IR spectra. The experiments demonstrate the use of T-jump 2D IR spectroscopy as a valuable tool for studying protein folding, with direct connections to simulations. The results suggest that in addition to predicting the correct native structure and folding time constant, molecular dynamics simulations carried out with modern force fields provide an accurate description of folding mechanisms in small proteins.