Deep Learning Model for Efficient Protein-Ligand Docking with Implicit Side-Chain Flexibility.

Protein-ligand docking is an essential tool in structure-based drug design with applications ranging from virtual high-throughput screening to pose prediction for lead optimization. Most docking programs for pose prediction are optimized for redocking to an existing cocrystallized protein structure, ignoring protein flexibility. In real-world drug design applications, however, protein flexibility is an essential feature of the ligand-binding process. Flexible protein-ligand docking still remains a significant challenge to computational drug design. To target this challenge, we present a deep learning (DL) model for flexible protein-ligand docking based on the prediction of an intermolecular Euclidean distance matrix (EDM), making the typical use of iterative search algorithms obsolete. The model was trained on a large-scale data set of protein-ligand complexes and evaluated on independent test sets. Our model generates high quality poses for a diverse set of protein and ligand structures and outperforms comparable docking methods.

Accurate Sampling of Macromolecular Conformations Using Adaptive Deep Learning and Coarse-Grained Representation.

Conformational sampling of protein structures is essential for understanding biochemical functions and for predicting thermodynamic properties such as free energies. Where previous approaches rely on sequential sampling procedures, recent developments in generative deep neural networks rendered possible the parallel, statistically independent sampling of molecular configurations. To be able to accurately generate samples of large molecular systems from a high-dimensional multimodal equilibrium distribution function, we developed a hierarchical approach based on expressive normalizing flows with rational quadratic neural splines and coarse-grained representation. Furthermore, system specific priors and adaptive and property-based controlled learning was designed to diminish the likelihood for the generation of high-energy structures during sampling. Finally, backmapping from a coarse-grained to fully atomistic representation is performed through an equivariant transformer model. We demonstrate the applicability of the method on the one-shot configurational sampling of a protein system with more than a hundred amino acids. The results show enhanced expressivity that diminish the invertibility constraints inherent in the normalizing flow framework. Moreover, the capacity of the hierarchical normalizing flow model was tested on a challenging case study of the folding/unfolding dynamics of the peptide chignolin.

Computational Assessment of Combination Therapy of Androgen Receptor-Targeting Compounds.

The ligand-binding domain of the androgen receptor (AR) is a target for drugs against prostate cancer and offers three distinct binding sites for small molecules. Drugs acting on the orthosteric hormone binding site suffer from resistance mechanisms that can, in the worst case, reverse their therapeutic effect. While many allosteric ligands targeting either the activation function-2 (AF-2) or the binding function-3 (BF-3) have been reported, their potential for simultaneous administration with currently prescribed antiandrogens was disregarded. Here, we report results of 60 µs molecular dynamics simulations to investigate combinations of orthosteric and allosteric AR antagonists. Our results suggest BF-3 inhibitors to be more suitable in combination with classical antiandrogens as opposed to AF-2 inhibitors based on binding free energies and binding modes. As a mechanistic explanation for these observations, we deduced a structural adaptation of helix-12 involved in the formation of the AF-2 site by classical AR antagonists. Additionally, the changes were accompanied by an expansion of the orthosteric binding site. Considering our predictions, the selective combination of AR-targeting compounds may improve the treatment of prostate cancer.

Conformational Changes of Thyroid Receptors in Response to Antagonists.

Thyroid hormone receptors (TRs) play a critical role in human development, growth, and metabolism. Antagonists of TRs offer an attractive strategy to treat hyperthyroidism without the disadvantage of a delayed onset of drug action. While it is challenging to examine the atomistic behavior of TRs in a laboratory setting, computational methods such as molecular dynamics (MD) simulations have proven their value to elucidate ligand-induced conformational changes in nuclear receptors. Here, we performed MD simulations of TRα and TRß complexed to their native ligand triiodothyronine (T3) as well as several antagonists. Based on the examination of 27 µs MD trajectories, we showed how binding of these compounds influences various structural features of the receptors including the helicity of helices 3 and 10 as well as the location of helix-12. Helices 3 and 12 are known to mediate coactivator association required for downstream signaling, suggesting these changes to be the molecular basis for TR antagonism. A mechanistic analysis of the trajectories revealed an allosteric pathway between H3 and H12 to be responsible for the conformational adaptations. Even though a mechanistic understanding of conformational adaptations triggered by TR antagonists is important for the development of novel therapeutics, they have not been previously examined in detail as it was done here.

Computational Selectivity Assessment of Protease Inhibitors against SARS-CoV-2.

The pandemic of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) poses a serious global health threat. Since no specific therapeutics are available, researchers around the world screened compounds to inhibit various molecular targets of SARS-CoV-2 including its main protease (Mpro) essential for viral replication. Due to the high urgency of these discovery efforts, off-target binding, which is one of the major reasons for drug-induced toxicity and safety-related drug attrition, was neglected. Here, we used molecular docking, toxicity profiling, and multiple molecular dynamics (MD) protocols to assess the selectivity of 33 reported non-covalent inhibitors of SARS-CoV-2 Mpro against eight proteases and 16 anti-targets. The panel of proteases included SARS-CoV Mpro, cathepsin G, caspase-3, ubiquitin carboxy-terminal hydrolase L1 (UCHL1), thrombin, factor Xa, chymase, and prostasin. Several of the assessed compounds presented considerable off-target binding towards the panel of proteases, as well as the selected anti-targets. Our results further suggest a high risk of off-target binding to chymase and cathepsin G. Thus, in future discovery projects, experimental selectivity assessment should be directed toward these proteases. A systematic selectivity assessment of SARS-CoV-2 Mpro inhibitors, as we report it, was not previously conducted.

Potential Inhibitors for Novel Coronavirus Protease Identified by Virtual Screening of 606 Million Compounds.

The rapid outbreak of the novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) in China followed by its spread around the world poses a serious global concern for public health. To this date, no specific drugs or vaccines are available to treat SARS-CoV-2 despite its close relation to the SARS-CoV virus that caused a similar epidemic in 2003. Thus, there remains an urgent need for the identification and development of specific antiviral therapeutics against SARS-CoV-2. To conquer viral infections, the inhibition of proteases essential for proteolytic processing of viral polyproteins is a conventional therapeutic strategy. In order to find novel inhibitors, we computationally screened a compound library of over 606 million compounds for binding at the recently solved crystal structure of the main protease (Mpro) of SARS-CoV-2. A screening of such a vast chemical space for SARS-CoV-2 Mpro inhibitors has not been reported before. After shape screening, two docking protocols were applied followed by the determination of molecular descriptors relevant for pharmacokinetics to narrow down the number of initial hits. Next, molecular dynamics simulations were conducted to validate the stability of docked binding modes and comprehensively quantify ligand binding energies. After evaluation of potential off-target binding, we report a list of 12 purchasable compounds, with binding affinity to the target protease that is predicted to be more favorable than that of the cocrystallized peptidomimetic compound. In order to quickly advise ongoing therapeutic intervention for patients, we evaluated approved antiviral drugs and other protease inhibitors to provide a list of nine compounds for drug repurposing. Furthermore, we identified the natural compounds (-)-taxifolin and rhamnetin as potential inhibitors of Mpro. Rhamnetin is already commercially available in pharmacies.

Modeling of Halogen-Protein Interactions in Co-Solvent Molecular Dynamics Simulations.

Co-solvent molecular dynamics (MD) simulations have recently become successful approaches in structure-based drug design but neglect important interactions such as halogen bonding. To be able to successfully model compound libraries containing halogenated ligands using co-solvent simulations, we investigate the use of halogenated benzene probes in co-solvent simulations on the two test systems human cathepsin L (hCatL) and the Y220C mutant of the tumor suppressor p53 (p53-Y220C). Our studies demonstrate that halogenated benzene probes indeed can unambiguously identify halogen-bonding interaction sites in the binding pocket and show superior correlation and ranking performance compared to standard co-solvent approaches.

Molecular Determinants of the Differential Modulation of Cav1.2 and Cav1.3 by Nifedipine and FPL 64176.

Nifedipine and FPL 64176 (FPL), which block and potentiate L-type voltage-gated Ca2+ channels, respectively, modulate Cav1.2 more potently than Cav1.3. To identify potential strategies for developing subtype-selective inhibitors, we investigated the role of divergent amino acid residues in transmembrane domains IIIS5 and the extracellular IIIS5-3P loop region in modulation of these channels by nifedipine and FPL. Insertion of the extracellular IIIS5-3P loop from Cav1.2 into Cav1.3 (Cav1.3+) reduced the IC50 of nifedipine from 289 to 101 nM, and substitution of S1100 with an A residue, as in Cav1.2, accounted for this difference. Substituting M1030 in IIIS5 to V in Cav1.3+ (Cav1.3+V) further reduced the IC50 of nifedipine to 42 nM. FPL increased current amplitude with an EC50 of 854 nM in Cav1.3, 103 nM in Cav1.2, and 99 nM in Cav1.3+V. In contrast to nifedipine block, substitution of M1030 to V in Cav1.3 had no effect on potency of FPL potentiation of current amplitude, but slowed deactivation in the presence and absence of 10 µM FPL. FPL had no effect on deactivation of Cav1.3/dihydropyridine-insensitive (DHPi), a channel with very low sensitivity to nifedipine block (IC50 ∼93 µM), but did shift the voltage-dependence of activation by ∼-10 mV. We conclude that the M/V variation in IIIS5 and the S/A variation in the IIIS5-3P loop of Cav1.2 and Cav1.3 largely determine the difference in nifedipine potency between these two channels, but the difference in FPL potency is determined by divergent amino acids in the IIIS5-3P loop.

Efficient and Accurate Hydration Site Profiling for Enclosed Binding Sites.

Molecular dynamics (MD) simulations allow for accurate prediction of the thermodynamic profile of binding-site water molecules critical for protein-ligand association. Whereas this hydration-site profiling converges rapidly for solvent-exposed sites independent of the initial water placement, an accurate and reliable placement is required for water molecules in occluded binding sites. Here, we present an accurate and efficient hydration-site prediction method for occluded binding sites combining water placement based on 3D-RISM and MD simulations using WATsite.

Molecular Modeling Evaluation of the Enantiomers of a Novel Adenylyl Cyclase 2 Inhibitor.

Adenylyl cyclase 2 (AC2) is one of nine membrane-bound isoforms of adenylyl cyclase that converts ATP into cyclic AMP (cAMP), an important second messenger molecule. Upregulation of AC2 is linked to cancers like pancreatic and small intestinal neuroendocrine tumors (NETs). The structures of the various isoforms of adenylyl cyclases are highly homologous, posing a significant challenge to drug discovery efforts for an effective, isoform-selective modulator of AC2. In a previous study, a screen identified a potential isoform-selective and noncompetitive inhibitor of AC2, SKF83566. In the present study, molecular modeling is used to explore the mode of inhibition of AC2 by SKF83566 and to investigate the active enantiomer of SKF83566. Homology models of hAC2 were built based on canine AC5-C1a and rat AC2-C2a templates. With these models, a combination of flexible docking, molecular dynamics simulations, and free energy calculations using the MM/GBSA methodology suggested an allosteric mechanism in which (S)-SKF83566 binds to an allosteric site near ATP and alters the protein conformation of the ATP binding site, potentially preventing the adenosine moiety of ATP from forming an archlike shape to form cAMP. The predicted binding preference for the (S)-SKF83566 enantiomer and the predicted free energy are consistent with the experimental data.

Ranking protein-protein docking results using steered molecular dynamics and potential of mean force calculations.

Crystallization of protein-protein complexes can often be problematic and therefore computational structural models are often relied on. Such models are often generated using protein-protein docking algorithms, where one of the main challenges is selecting which of several thousand potential predictions represents the most near-native complex. We have developed a novel technique that involves the use of steered molecular dynamics (sMD) and umbrella sampling to identify near-native complexes among protein-protein docking predictions. Using this technique, we have found a strong correlation between our predictions and the interface RMSD (iRMSD) in ten diverse test systems. On two of the systems, we investigated if the prediction results could be further improved using potential of mean force calculations. We demonstrated that a near-native (<2.0 Å iRMSD) structure could be identified in the top-1 ranked position for both systems. © 2016 Wiley Periodicals, Inc.

Structural transitions and interactions in the early stages of human glucagon amyloid fibrillation.

A mechanistic understanding of the intermolecular interactions and structural changes during fibrillation is crucial for the design of safe and efficacious glucagon formulations. Amide hydrogen/deuterium exchange with mass spectrometric analysis was used to identify the interactions and amino acids involved in the initial stages of glucagon fibril formation at acidic pH. Kinetic measurements from intrinsic and thioflavin T fluorescence showed sigmoidal behavior. Secondary structural measurement of fibrillating glucagon using far-UV circular dichroism spectroscopy showed changes in structure from random coil → α-helix → ß-sheet, with increase in α-helix content during the lag phase followed by increase in ß-sheet content during the growth phase. Hydrogen/deuterium exchange with mass spectrometric analysis of fibrillating glucagon suggested that C-terminal residues 22-29 are involved in interactions during the lag phase, during which N-terminal residues 1-6 showed no changes. Molecular dynamics simulations of glucagon fragments showed C-terminal to C-terminal interactions with greater α-helix content for the 20-29 fragment, with hydrophobic and aromatic residues (Phe-22, Trp-25, Val-23, and Met-27) predominantly involved. Overall, the study shows that glucagon interactions during the early phase of fibrillation are mediated through C-terminal residues, which facilitate the formation of α-helix-rich oligomers, which further undergo structural rearrangement and elongation to form ß-sheet-rich mature fibrils.

Substrate tunnels in enzymes: structure-function relationships and computational methodology.

In enzymes, the active site is the location where incoming substrates are chemically converted to products. In some enzymes, this site is deeply buried within the core of the protein, and, in order to access the active site, substrates must pass through the body of the protein via a tunnel. In many systems, these tunnels act as filters and have been found to influence both substrate specificity and catalytic mechanism. Identifying and understanding how these tunnels exert such control has been of growing interest over the past several years because of implications in fields such as protein engineering and drug design. This growing interest has spurred the development of several computational methods to identify and analyze tunnels and how ligands migrate through these tunnels. The goal of this review is to outline how tunnels influence substrate specificity and catalytic efficiency in enzymes with buried active sites and to provide a brief summary of the computational tools used to identify and evaluate these tunnels.

Combining structure- and ligand-based approaches to improve site of metabolism prediction in CYP2C9 substrates.

PURPOSE: Predicting atoms in a potential drug compound that are susceptible to oxidation by cytochrome P450 (CYP) enzymes is of great interest to the pharmaceutical community. We aimed to develop a computational approach combining ligand- and structure-based design principles to accurately predict sites of metabolism (SoMs) in a series of CYP2C9 substrates. METHODS: We employed the reactivity model, SMARTCyp, ensemble docking, and pseudo-receptor modeling based on quantitative structure-activity relationships (QSAR) to account for influences of both the inherent reactivity of each atom and the physical structure of the CYP2C9 binding site. RESULTS: We tested ligand-based prediction alone (i.e. SMARTCyp), structure-based prediction alone (i.e. AutoDock Vina docking), the linear combination of the SMARTCYP and docking scores, and finally a pseudo-receptor QSAR model based on the docked compounds in combination with SMARTCyp. We found that by using the latter combined approach we were able to accurately predict 88% and 96% of the true SoMs, within the top-1 and top-2 predictions, respectively. CONCLUSIONS: We have outlined a novel combination approach for accurately predicting SoMs in CYP2C9 ligands. We believe that this method may be applied to other CYP2C9 ligands as well as to other CYP systems.

Application of information theory to a three-body coarse-grained representation of proteins in the PDB: insights into the structural and evolutionary roles of residues in protein structure.

Knowledge-based methods for analyzing protein structures, such as statistical potentials, primarily consider the distances between pairs of bodies (atoms or groups of atoms). Considerations of several bodies simultaneously are generally used to characterize bonded structural elements or those in close contact with each other, but historically do not consider atoms that are not in direct contact with each other. In this report, we introduce an information-theoretic method for detecting and quantifying distance-dependent through-space multibody relationships between the sidechains of three residues. The technique introduced is capable of producing convergent and consistent results when applied to a sufficiently large database of randomly chosen, experimentally solved protein structures. The results of our study can be shown to reproduce established physico-chemical properties of residues as well as more recently discovered properties and interactions. These results offer insight into the numerous roles that residues play in protein structure, as well as relationships between residue function, protein structure, and evolution. The techniques and insights presented in this work should be useful in the future development of novel knowledge-based tools for the evaluation of protein structure.

Including ligand-induced protein flexibility into protein tunnel prediction.

In proteins with buried active sites, understanding how ligands migrate through the tunnels that connect the exterior of the protein to the active site can shed light on substrate specificity and enzyme function. A growing body of evidence highlights the importance of protein flexibility in the binding site on ligand binding; however, the influence of protein flexibility throughout the body of the protein during ligand entry and egress is much less characterized. We have developed a novel tunnel prediction and evaluation method named IterTunnel, which includes the influence of ligand-induced protein flexibility, guarantees ligand egress, and provides detailed free energy information as the ligand proceeds along the egress route. IterTunnel combines geometric tunnel prediction with steered molecular dynamics in an iterative process to identify tunnels that open as a result of ligand migration and calculates the potential of mean force of ligand egress through a given tunnel. Applying this new method to cytochrome P450 2B6, we demonstrate the influence of protein flexibility on the shape and accessibility of tunnels. More importantly, we demonstrate that the ligand itself, while traversing through a tunnel, can reshape tunnels due to its interaction with the protein. This process results in the exposure of new tunnels and the closure of preexisting tunnels as the ligand migrates from the active site.

WATsite: hydration site prediction program with PyMOL interface.

Water molecules that mediate protein-ligand interactions or are released from the binding site on ligand binding can contribute both enthalpically and entropically to the free energy of ligand binding. To elucidate the thermodynamic profile of individual water molecules and their potential contribution to ligand binding, a hydration site analysis program WATsite was developed together with an easy-to-use graphical user interface based on PyMOL. WATsite identifies hydration sites from a molecular dynamics simulation trajectory with explicit water molecules. The free energy profile of each hydration site is estimated by computing the enthalpy and entropy of the water molecule occupying a hydration site throughout the simulation. The results of the hydration site analysis can be displayed in PyMOL. A key feature of WATsite is that it is able to estimate the protein desolvation free energy for any user specified ligand. The WATsite program and its PyMOL plugin are available free of charge from http://people.pnhs.purdue.edu/~mlill/software.

Analysis of factors influencing hydration site prediction based on molecular dynamics simulations.

Water contributes significantly to the binding of small molecules to proteins in biochemical systems. Molecular dynamics (MD) simulation based programs such as WaterMap and WATsite have been used to probe the locations and thermodynamic properties of hydration sites at the surface or in the binding site of proteins generating important information for structure-based drug design. However, questions associated with the influence of the simulation protocol on hydration site analysis remain. In this study, we use WATsite to investigate the influence of factors such as simulation length and variations in initial protein conformations on hydration site prediction. We find that 4 ns MD simulation is appropriate to obtain a reliable prediction of the locations and thermodynamic properties of hydration sites. In addition, hydration site prediction can be largely affected by the initial protein conformations used for MD simulations. Here, we provide a first quantification of this effect and further indicate that similar conformations of binding site residues (RMSD < 0.5 Å) are required to obtain consistent hydration site predictions.

Biosynthesis of the bacterial antibiotic 3,7-dihydroxytropolone through enzymatic salvaging of catabolic shunt products.

The non-benzenoid aromatic tropone ring is a structural motif of numerous microbial and plant natural products with potent bioactivities. In bacteria, tropone biosynthesis involves early steps of the widespread CoA-dependent phenylacetic acid (paa) catabolon, from which a shunt product is sequestered and surprisingly further utilized as a universal precursor for structurally and functionally diverse tropone derivatives such as tropodithietic acid or (hydroxy)tropolones. Here, we elucidate the biosynthesis of the antibiotic 3,7-dihydroxytropolone in Actinobacteria by in vitro pathway reconstitution using paa catabolic enzymes as well as dedicated downstream tailoring enzymes, including a thioesterase (TrlF) and two flavoprotein monooxygenases (TrlCD and TrlE). We furthermore mechanistically and structurally characterize the multifunctional key enzyme TrlE, which mediates an unanticipated ipso-substitution involving a hydroxylation and subsequent decarboxylation of the CoA-freed side chain, followed by ring oxidation to afford tropolone. This study showcases a remarkably efficient strategy for 3,7-dihydroxytropolone biosynthesis and illuminates the functions of the involved biosynthetic enzymes.

Stabilization of an unusual salt bridge in ubiquitin by the extra C-terminal domain of the proteasome-associated deubiquitinase UCH37 as a mechanism of its exo specificity.

Ubiquitination is countered by a group of enzymes collectively called deubiquitinases (DUBs); ∼100 of them can be found in the human genome. One of the most interesting aspects of these enzymes is the ability of some members to selectively recognize specific linkage types between ubiquitin in polyubiquitin chains and their endo and exo specificity. The structural basis of exo-specific deubiquitination catalyzed by a DUB is poorly understood. UCH37, a cysteine DUB conserved from fungi to humans, is a proteasome-associated factor that regulates the proteasome by sequentially cleaving polyubiquitin chains from their distal ends, i.e., by exo-specific deubiquitination. In addition to the catalytic domain, the DUB features a functionally uncharacterized UCH37-like domain (ULD), presumed to keep the enzyme in an inhibited state in its proteasome-free form. Herein we report the crystal structure of two constructs of UCH37 from Trichinella spiralis in complex with a ubiquitin-based suicide inhibitor, ubiquitin vinyl methyl ester (UbVME). These structures show that the ULD makes direct contact with ubiquitin stabilizing a highly unusual intramolecular salt bridge between Lys48 and Glu51 of ubiquitin, an interaction that would be favored only with the distal ubiquitin but not with the internal ones in a Lys48-linked polyubiquitin chain. An inspection of 39 DUB-ubiquitin structures in the Protein Data Bank reveals the uniqueness of the salt bridge in ubiquitin bound to UCH37, an interaction that disappears when the ULD is deleted, as revealed in the structure of the catalytic domain alone bound to UbVME. The structural data are consistent with previously reported mutational data on the mammalian enzyme, which, together with the fact that the ULD residues that bind to ubiquitin are conserved, points to a similar mechanism behind the exo specificity of the human enzyme. To the best of our knowledge, these data provide the only structural example so far of how the exo specificity of a DUB can be determined by its noncatalytic domain. Thus, our data show that, contrary to its proposed inhibitory role, the ULD actually contributes to substrate recognition and could be a major determinant of the proteasome-associated function of UCH37. Moreover, our structures show that the unproductively oriented catalytic cysteine in the free enzyme is aligned correctly when ubiquitin binds, suggesting a mechanism for ubiquitin selectivity.