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.

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.

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.

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

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.

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, 64Cu-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 ( NCF3Rpe) 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.

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.

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

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.

Microsecond simulations of mdm2 and its complex with p53 yield insight into force field accuracy and conformational dynamics.

The p53-MDM2 complex is both a major target for cancer drug development and a valuable model system for computational predictions of protein-ligand binding. To investigate the accuracy of molecular simulations of MDM2 and its complex with p53, we performed a number of long (200 ns to 1 µs) explicit-solvent simulations using a range of force fields. We systematically compared nine popular force fields (AMBER ff03, ff12sb, ff14sb, ff99sb, ff99sb-ildn, ff99sb-ildn-nmr, ff99sb-ildn-phi, CHARMM22*, and CHARMM36) against experimental chemical shift data, and found similarly accurate results, with microsecond simulations achieving better agreement compared to 200-ns trajectories. Although the experimentally determined apo structure has a closed binding cleft, simulations in all force fields suggest the apo state of MDM2 is highly flexible, and able to sample holo-like conformations, consistent with a conformational selection model. Initial structuring of the MDM2 lid region, known to competitively bind the binding cleft, is also observed in long simulations. Taken together, these results show molecular simulations can accurately sample conformations relevant for ligand binding. We expect this study to inform future computational work on folding and binding of MDM2 ligands.

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.

De novo structure prediction and experimental characterization of folded peptoid oligomers.

Peptoid molecules are biomimetic oligomers that can fold into unique three-dimensional structures. As part of an effort to advance computational design of folded oligomers, we present blind-structure predictions for three peptoid sequences using a combination of Replica Exchange Molecular Dynamics (REMD) simulation and Quantum Mechanical refinement. We correctly predicted the structure of a N-aryl peptoid trimer to within 0.2 Å rmsd-backbone and a cyclic peptoid nonamer to an accuracy of 1.0 Å rmsd-backbone. X-ray crystallographic structures are presented for a linear N-alkyl peptoid trimer and for the cyclic peptoid nonamer. The peptoid macrocycle structure features a combination of cis and trans backbone amides, significant nonplanarity of the amide bonds, and a unique "basket" arrangement of (S)-N(1-phenylethyl) side chains encompassing a bound ethanol molecule. REMD simulations of the peptoid trimers reveal that well folded peptoids can exhibit funnel-like conformational free energy landscapes similar to those for ordered polypeptides. These results indicate that physical modeling can successfully perform de novo structure prediction for small peptoid molecules.

Probing the nanosecond dynamics of a designed three-stranded beta-sheet with a massively parallel molecular dynamics simulation.

Recently a temperature-jump FTIR study of a designed three-stranded sheet showing a fast relaxation time of approximately 140 +/- 20 ns was published. We performed massively parallel molecular dynamics simulations in explicit solvent to probe the structural events involved in this relaxation. While our simulations produce similar relaxation rates, the structural ensemble is broad. We observe the formation of turn structure, but only very weak interaction in the strand regions, which is consistent with the lack of strong backbone-backbone NOEs in previous structural NMR studies. These results suggest that either (D)P(D)P-II folds at time scales longer than 240 ns, or that (D)P(D)P-II is not a well-defined three-stranded beta-sheet. This work also provides an opportunity to compare the performance of several popular forcefield models against one another.