Deletions in the stalk of the influenza neuraminidase (NA) surface protein are associated with increased virulence, but the mechanisms responsible for this enhanced virulence are unclear. Here we use microsecond molecular dynamics simulations to explore the effect of stalk deletion on enzymatic activity, contrasting NA proteins from the A/swine/Shandong/N1/2009 strain both with and without a stalk deletion. By modeling and simulating neuraminidase apo glycoproteins embedded in complex-mixture lipid bilayers, we show that the geometry and dynamics of the neuraminidase enzymatic pocket may differ depending on stalk length, with possible repercussions on the binding of the endogenous sialylated-oligosaccharide receptors. We also use these simulations to predict previously unrecognized druggable "hotspots" on the neuraminidase surface that may prove useful for future efforts aimed at structure-based drug design.

The kinetic rate constants of binding were estimated for four biochemically relevant molecular systems by a method that uses milestoning theory to combine Brownian dynamics simulations with more detailed molecular dynamics simulations. The rate constants found using this method agreed well with experimentally and theoretically obtained values. We predicted the association rate of a small charged molecule toward both a charged and an uncharged spherical receptor and verified the estimated value with Smoluchowski theory. We also calculated the kon rate constant for superoxide dismutase with its natural substrate, O2-, in a validation of a previous experiment using similar methods but with a number of important improvements. We also calculated the kon for a new system: the N-terminal domain of Troponin C with its natural substrate Ca2+. The kon calculated for the latter two systems closely resemble experimentally obtained values. This novel multiscale approach is computationally cheaper and more parallelizable when compared to other methods of similar accuracy. We anticipate that this methodology will be useful for predicting kinetic rate constants and for understanding the process of binding between a small molecule and a protein receptor.

No vaccine exists against group A Streptococcus (GAS), a leading cause of worldwide morbidity and mortality. A severe hurdle is the hypervariability of its major antigen, the M protein, with >200 different M types known. Neutralizing antibodies typically recognize M protein hypervariable regions (HVRs) and confer narrow protection. In stark contrast, human C4b-binding protein (C4BP), which is recruited to the GAS surface to block phagocytic killing, interacts with a remarkably large number of M protein HVRs (apparently ∼90%). Such broad recognition is rare, and we discovered a unique mechanism for this through the structure determination of four sequence-diverse M proteins in complexes with C4BP. The structures revealed a uniform and tolerant 'reading head' in C4BP, which detected conserved sequence patterns hidden within hypervariability. Our results open up possibilities for rational therapies that target the M-C4BP interaction, and also inform a path towards vaccine design.

Ligand-induced protein allostery plays a central role in modulating cellular signalling pathways. Here using the conserved cyclic nucleotide-binding domain of protein kinase A's (PKA) regulatory subunit as a prototype signalling unit, we combine long-timescale, all-atom molecular dynamics simulations with Markov state models to elucidate the conformational ensembles of PKA's cyclic nucleotide-binding domain A for the cAMP-free (apo) and cAMP-bound states. We find that both systems exhibit shallow free-energy landscapes that link functional states through multiple transition pathways. This observation suggests conformational selection as the general mechanism of allostery in this canonical signalling domain. Further, we expose the propagation of the allosteric signal through key structural motifs in the cyclic nucleotide-binding domain and explore the role of kinetics in its function. Our approach integrates disparate lines of experimental data into one cohesive framework to understand structure, dynamics and function in complex biological systems.

The atomic-level mechanisms that coordinate ligand release from protein pockets are only known for a handful of proteins. Here, we report results from accelerated molecular dynamics simulations for benzene dissociation from the buried cavity of the T4 lysozyme Leu99Ala mutant (L99A). In these simulations, benzene is released through a previously characterized, sparsely populated room-temperature excited state of the mutant, explaining the coincidence for experimentally measured benzene off rate and apo protein slow-timescale NMR relaxation rates between ground and excited states. The path observed for benzene egress is a multistep ligand migration from the buried cavity to ultimate release through an opening between the F/G-, H-, and I-helices and requires a number of cooperative multiresidue and secondary-structure rearrangements within the C-terminal domain of L99A. These rearrangements are identical to those observed along the ground state to excited state transitions characterized by molecular dynamic simulations run on the Anton supercomputer. Analyses of the molecular properties of the residues lining the egress path suggest that protein surface electrostatic potential may play a role in the release mechanism. Simulations of wild-type T4 lysozyme also reveal that benzene-egress-associated dynamics in the L99A mutant are potentially exaggerations of the substrate-processivity-related dynamics of the wild type.

As ever larger and more complex biological systems are modeled in silico, approximating physiological lipid bilayers with simple planar models becomes increasingly unrealistic. In order to build accurate large-scale models of subcellular environments, models of lipid membranes with carefully considered, biologically relevant curvature will be essential. In the current work, we present a multi-scale utility called LipidWrapper capable of creating curved membrane models with geometries derived from various sources, both experimental and theoretical. To demonstrate its utility, we use LipidWrapper to examine an important mechanism of influenza virulence. A copy of the program can be downloaded free of charge under the terms of the open-source FreeBSD License from http://nbcr.ucsd.edu/lipidwrapper. LipidWrapper has been tested on all major computer operating systems.

Allostery can occur by way of subtle cooperation among protein residues (e.g., amino acids) even in the absence of large conformational shifts. Dynamical network analysis has been used to model this cooperation, helping to computationally explain how binding to an allosteric site can impact the behavior of a primary site many ångstroms away. Traditionally, computational efforts have focused on the most optimal path of correlated motions leading from the allosteric to the primary active site. We present a program called Weighted Implementation of Suboptimal Paths (WISP) capable of rapidly identifying additional suboptimal pathways that may also play important roles in the transmission of allosteric signals. Aside from providing signal redundancy, suboptimal paths traverse residues that, if disrupted through pharmacological or mutational means, could modulate the allosteric regulation of important drug targets. To demonstrate the utility of our program, we present a case study describing the allostery of HisH-HisF, an amidotransferase from T. maritima thermotiga. WISP and its VMD-based graphical user interface (GUI) can be downloaded from http://nbcr.ucsd.edu/wisp.

The recently discovered 150-cavity in the active site of group-1 influenza A neuraminidase (NA) proteins provides a target for rational structure-based drug development to counter the increasing frequency of antiviral resistance in influenza. Surprisingly, the 2009 H1N1 pandemic virus (09N1) neuraminidase was crystalized without the 150-cavity characteristic of group-1 NAs. Here we demonstrate, through a total sum of 1.6 μs of biophysical simulations, that 09N1 NA exists in solution preferentially with an open 150-cavity. Comparison with simulations using avian N1, human N2 and 09N1 with a I149V mutation and an extensive bioinformatics analysis suggests that the conservation of a key salt bridge is crucial in the stabilization of the 150-cavity across both subtypes. This result provides an atomic-level structural understanding of the recent finding that antiviral compounds designed to take advantage of contacts in the 150-cavity can inactivate both 2009 H1N1 pandemic and avian H5N1 viruses.

Chlamydia trachomatis is the most prevalent cause of bacterial sexually transmitted diseases and the leading cause of preventable blindness worldwide. Global control of Chlamydia will best be achieved with a vaccine, a primary target for which is the major outer membrane protein, MOMP, which comprises ~60% of the outer membrane protein mass of this bacterium. In the absence of experimental structural information on MOMP, three previously published topology models presumed a16-stranded barrel architecture. Here, we use the latest β-barrel prediction algorithms, previous 2D topology modeling results, and comparative modeling methodology to build a 3D model based on the 16-stranded, trimeric assumption. We find that while a 3D MOMP model captures many structural hallmarks of a trimeric 16-stranded β-barrel porin, and is consistent with most of the experimental evidence for MOMP, MOMP residues 320-334 cannot be modeled as β-strands that span the entire membrane, as is consistently observed in published 16-stranded β-barrel crystal structures. Given the ambiguous results for β-strand delineation found in this study, recent publications of membrane β-barrel structures breaking with the canonical rule for an even number of β-strands, findings of β-barrels with strand-exchanged oligomeric conformations, and alternate folds dependent upon the lifecycle of the bacterium, we suggest that although the MOMP porin structure incorporates canonical 16-stranded conformations, it may have novel oligomeric or dynamic structural changes accounting for the discrepancies observed.

The tumor suppressor protein p53 can lose its function upon single-point missense mutations in the core DNA-binding domain ("cancer mutants"). Activity can be restored by second-site suppressor mutations ("rescue mutants"). This paper relates the functional activity of p53 cancer and rescue mutants to their overall molecular dynamics (MD), without focusing on local structural details. A novel global measure of protein flexibility for the p53 core DNA-binding domain, the number of clusters at a certain RMSD cutoff, was computed by clustering over 0.7 µs of explicitly solvated all-atom MD simulations. For wild-type p53 and a sample of p53 cancer or rescue mutants, the number of clusters was a good predictor of in vivo p53 functional activity in cell-based assays. This number-of-clusters (NOC) metric was strongly correlated (r(2) = 0.77) with reported values of experimentally measured ΔΔG protein thermodynamic stability. Interpreting the number of clusters as a measure of protein flexibility: (i) p53 cancer mutants were more flexible than wild-type protein, (ii) second-site rescue mutations decreased the flexibility of cancer mutants, and (iii) negative controls of non-rescue second-site mutants did not. This new method reflects the overall stability of the p53 core domain and can discriminate which second-site mutations restore activity to p53 cancer mutants.

Conventional drug design embraces the "one gene, one drug, one disease" philosophy. Polypharmacology, which focuses on multi-target drugs, has emerged as a new paradigm in drug discovery. The rational design of drugs that act via polypharmacological mechanisms can produce compounds that exhibit increased therapeutic potency and against which resistance is less likely to develop. Additionally, identifying multiple protein targets is also critical for side-effect prediction. One third of potential therapeutic compounds fail in clinical trials or are later removed from the market due to unacceptable side effects often caused by off-target binding. In the current work, we introduce a multidimensional strategy for the identification of secondary targets of known small-molecule inhibitors in the absence of global structural and sequence homology with the primary target protein. To demonstrate the utility of the strategy, we identify several targets of 4,5-dihydroxy-3-(1-naphthyldiazenyl)-2,7-naphthalenedisulfonic acid, a known micromolar inhibitor of Trypanosoma brucei RNA editing ligase 1. As it is capable of identifying potential secondary targets, the strategy described here may play a useful role in future efforts to reduce drug side effects and/or to increase polypharmacology.

Hemagglutinins (HA's) from duck, swine, and human influenza viruses have previously been shown to prefer avian and human glycan receptor analogues with distinct topological profiles, pentasaccharides LSTa (alpha-2,3 linkage) and LSTc (alpha-2,6 linkage), in comparative molecular dynamics studies. On the basis of detailed analyses of the dynamic motions of the receptor binding domains (RBDs) and interaction energy profiles with individual glycan residues, we have identified approximately 30 residue positions in the RBD that present distinct profiles with the receptor analogues. Glycan binding constrained the conformational space sampling by the HA. Electrostatic steering appeared to play a key role in glycan binding specificity. The complex dynamic behaviors of the major SSE and trimeric interfaces with or without bound glycans suggested that networks of interactions might account for species specificity in these low affinity and high avidity (multivalent) interactions between different HA and glycans. Contact frequency, energetic decomposition, and H-bond analyses revealed species-specific differences in HA-glycan interaction profiles, not readily discernible from crystal structures alone. Interaction energy profiles indicated that mutation events at the set of residues such as 145, 156, 158, and 222 would favor human or avian receptor analogues, often through interactions with distal asialo-residues. These results correlate well with existing experimental evidence, and suggest new opportunities for simulation-based vaccine and drug development.

The interactions among associating (macro)molecules are dynamic, which adds to the complexity of molecular recognition. While ligand flexibility is well accounted for in computational drug design, the effective inclusion of receptor flexibility remains an important challenge. The relaxed complex scheme (RCS) is a promising computational methodology that combines the advantages of docking algorithms with dynamic structural information provided by molecular dynamics (MD) simulations, therefore explicitly accounting for the flexibility of both the receptor and the docked ligands. Here, we briefly review the RCS and discuss new extensions and improvements of this methodology in the context of ligand binding to two example targets: kinetoplastid RNA editing ligase 1 and the W191G cavity mutant of cytochrome c peroxidase. The RCS improvements include its extension to virtual screening, more rigorous characterization of local and global binding effects, and methods to improve its computational efficiency by reducing the receptor ensemble to a representative set of configurations. The choice of receptor ensemble, its influence on the predictive power of RCS, and the current limitations for an accurate treatment of the solvent contributions are also briefly discussed. Finally, we outline potential methodological improvements that we anticipate will assist future development.

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

The discovery of new drugs is a time consuming and expensive process. Methods such as virtual screening, which can filter out ineffective compounds from drug libraries prior to expensive experimental study, have become popular research topics. As the computational drug discovery community has grown, in order to benchmark the various advances in methodology, organizations such as the Drug Design Data Resource have begun hosting blinded grand challenges seeking to identify the best methods for ligand pose-prediction, ligand affinity ranking, and free energy calculations. Such open challenges offer a unique opportunity for researchers to partner with junior students (e.g., high school and undergraduate) to validate basic yet fundamental hypotheses considered to be uninteresting to domain experts. Here, we, a group of high school-aged students and their mentors, present the results of our participation in Grand Challenge 4 where we predicted ligand affinity rankings for the Cathepsin S protease, an important protein target for autoimmune diseases. To investigate the effect of incorporating receptor dynamics on ligand affinity rankings, we employed the Relaxed Complex Scheme, a molecular docking method paired with molecular dynamics-generated receptor conformations. We found that Cathepsin S is a difficult target for molecular docking and we explore some advanced methods such as distance-restrained docking to try to improve the correlation with experiments. This project has exemplified the capabilities of high school students when supported with a rigorous curriculum, and demonstrates the value of community-driven competitions for beginners in computational drug discovery.

Mucin-domain glycoproteins are densely O-glycosylated and play critical roles in a host of biological functions. In particular, the T cell immunoglobulin and mucin-domain containing family of proteins (TIM-1, -3, -4) decorate immune cells and act as key checkpoint inhibitors in cancer. However, their dense O-glycosylation remains enigmatic both in terms of glycoproteomic landscape and structural dynamics, primarily due to the challenges associated with studying mucin domains. Here, we present a mucinase (SmE) and demonstrate its ability to selectively cleave along the mucin glycoprotein backbone, similar to others of its kind. Unlike other mucinases, though, SmE harbors the unique ability to cleave at residues bearing extremely complex glycans which enabled improved mass spectrometric analysis of several mucins, including the entire TIM family. With this information in-hand, we performed molecular dynamics (MD) simulations of TIM-3 and -4 to demonstrate how glycosylation affects structural features of these proteins. Overall, we present a powerful workflow to better understand the detailed molecular structures of the mucinome.

Influenza A has caused several deadly pandemics throughout human history. The virus is often resistant to developed treatments because of its genetic drift or shift property. Broad-spectrum antibodies show a promising potential to overcome the resistance of influenza viruses. In silico studies on broad-reactive antibodies and their interactions with hemagglutinins might shed light on the rational design of a universal vaccine. In this study, 11 broad-spectrum antibodies (or antigen-binding fragments) and 14 hemagglutinins of H3N2 and H5N1 strains were docked and analyzed to provide information about the construction of the scaffold for using universal antibodies against the influenza A virus. Antigen-binding fragments that have high number of appearances in the top 3 within each H3 and H5 subtypes were chosen for protein-protein interaction analysis. The results show that while the hydrogen bond is important for Ab/Fab binding to H3, the H5-Ab/Fab system may need cation-pi interaction for a strong interaction.

High-fidelity replication of the large RNA genome of coronaviruses (CoVs) is mediated by a 3'-to-5' exoribonuclease (ExoN) in non-structural protein 14 (nsp14), which excises nucleotides including antiviral drugs mis-incorporated by the low-fidelity viral RNA-dependent RNA polymerase (RdRp) and has also been implicated in viral RNA recombination and resistance to innate immunity. Here we determined a 1.6-Å resolution crystal structure of SARS-CoV-2 ExoN in complex with its essential co-factor, nsp10. The structure shows a highly basic and concave surface flanking the active site, comprising several Lys residues of nsp14 and the N-terminal amino group of nsp10. Modeling suggests that this basic patch binds to the template strand of double-stranded RNA substrates to position the 3' end of the nascent strand in the ExoN active site, which is corroborated by mutational and computational analyses. Molecular dynamics simulations further show remarkable flexibility of multi-domain nsp14 and suggest that nsp10 stabilizes ExoN for substrate RNA-binding to support its exoribonuclease activity. Our high-resolution structure of the SARS-CoV-2 ExoN-nsp10 complex serves as a platform for future development of anti-coronaviral drugs or strategies to attenuate the viral virulence.

Recent outbreaks of highly pathogenic and occasional drug-resistant influenza strains have highlighted the need to develop novel anti-influenza therapeutics. Here, we report computational and experimental efforts to identify influenza neuraminidase inhibitors from among the 3000 natural compounds in the Malaysian-Plants Natural-Product (NADI) database. These 3000 compounds were first docked into the neuraminidase active site. The five plants with the largest number of top predicted ligands were selected for experimental evaluation. Twelve specific compounds isolated from these five plants were shown to inhibit neuraminidase, including two compounds with IC50 values less than 92 μM. Furthermore, four of the 12 isolated compounds had also been identified in the top 100 compounds from the virtual screen. Together, these results suggest an effective new approach for identifying bioactive plant species that will further the identification of new pharmacologically active compounds from diverse natural-product resources.

APOBEC3 family DNA cytosine deaminases provide overlapping defenses against pathogen infections. However, most viruses have elaborate evasion mechanisms such as the HIV-1 Vif protein, which subverts cellular CBF-β and a polyubiquitin ligase complex to neutralize these enzymes. Despite advances in APOBEC3 and Vif biology, a full understanding of this direct host-pathogen conflict has been elusive. We combine virus adaptation and computational studies to interrogate the APOBEC3F-Vif interface and build a robust structural model. A recurring compensatory amino acid substitution from adaptation experiments provided an initial docking constraint, and microsecond molecular dynamic simulations optimized interface contacts. Virus infectivity experiments validated a long-lasting electrostatic interaction between APOBEC3F E289 and HIV-1 Vif R15. Taken together with mutagenesis results, we propose a wobble model to explain how HIV-1 Vif has evolved to bind different APOBEC3 enzymes and, more generally, how pathogens may evolve to escape innate host defenses.