Given the need for modern researchers to produce open, reproducible scientific output, the lack of standards and best practices for sharing data and workflows used to produce and analyze molecular dynamics (MD) simulations has become an important issue in the field. There are now multiple well-established packages to perform molecular dynamics simulations, often highly tuned for exploiting specific classes of hardware, each with strong communities surrounding them, but with very limited interoperability/transferability options. Thus, the choice of the software package often dictates the workflow for both simulation production and analysis. The level of detail in documenting the workflows and analysis code varies greatly in published work, hindering reproducibility of the reported results and the ability for other researchers to build on these studies. An increasing number of researchers are motivated to make their data available, but many challenges remain in order to effectively share and reuse simulation data. To discuss these and other issues related to best practices in the field in general, we organized a workshop in November 2018 (https://bioexcel.eu/events/workshop-on-sharing-data-from-molecular-simulations/). Here, we present a brief overview of this workshop and topics discussed. We hope this effort will spark further conversation in the MD community to pave the way toward more open, interoperable, and reproducible outputs coming from research studies using MD simulations.
Nanobody binding stabilizes G-protein-coupled receptors (GPCR) in a fully active state and modulates their affinity for bound ligands. However, the atomic-level basis for this allosteric regulation remains elusive. Here, we investigate the conformational changes induced by the binding of a nanobody (Nb80) on the active-like beta 2 adrenergic receptor (beta 2AR) via enhanced sampling molecular dynamics simulations. Dimensionality reduction analysis shows that Nb80 stabilizes structural features of the beta 2AR with an similar to 14 angstrom outward movement of transmembrane helix 6 and a close proximity of transmembrane (TM) helices 5 and 7, and favors the fully active-like conformation of the receptor, independent of ligand binding, in contrast to the conditions under which no intracellular binding partner is bound, in which case the receptor is only stabilized in an intermediateactive state. This activation is supported by the residues located at hotspots located on TMs 5, 6, and 7, as shown by supervised machine learning methods. Besides, ligand-specific subtle differences in the conformations assumed by intracellular loop 2 and extracellular loop 2 are captured from the trajectories of various ligand-bound receptors in the presence of Nb80. Dynamic network analysis further reveals that Nb80 binding triggers tighter and stronger local communication networks between the Nb80 and the ligand-binding sites, primarily involving residues around ICL2 and the intracellular end of TM3, TM5, TM6, as well as ECL2, ECL3, and the extracellular ends of TM6 and TM7. In particular, we identify unique allosteric signal transmission mechanisms between the Nb80-binding site and the extracellular domains in conformations modulated by a full agonist, BI167107, and a G-protein-biased partial agonist, salmeterol, involving mainly TM1 and TM2, and TM5, respectively. Altogether, our results provide insights into the effect of intracellular binding partners on the GPCR activation mechanism, which should be taken into account in structure-based drug discovery.
The farnesoid X receptor (FXR) is a bile acid-sensing transcription factor with indispensable roles in regulating metabolic processes. Nowadays, FXR has become a highly promising drug target for severe liver disorders, especially nonalcoholic steatohepatitis (NASH). A recent study showed that imatinib and its analogues were able to allosterically enhance agonist-induced FXR activation and its target gene expression. However, the allosteric modulation mechanism of FXR by these compounds remains unclear. In this work, the most effective imatinib analogue, P16, was used as a probe to explore this issue by computational approaches. Our results identified one potential allosteric site surrounded by residues Ile335, Phe336, Lys338, Glu339, Leu340, and Leu348, which could efficiently accommodate P16. In addition, the long-time molecular dynamics simulations indicated that the binding of P16 could significantly decrease the fluctuation of the co-activator and enhance the communications between the endogenous ligand chenodeoxycholic acid (CDCA) and FXR. By analyzing the residue interaction network, we observed two unique communication pathways connecting P16 and CDCA through three key residues, Arg331, Ser332, and Phe336. The communications of network organization in the P16-bound complex may allow the synergistic effect of the two compounds via robust signal transmission between the binding sites and global network bridges, which coordinate allosteric transitions and modulate the receptor activity. Our study offers insights into the allosteric modulation occurring in FXR and would be helpful for discovery of new allosteric modulators targeting FXR for further clinical research.
Constant pH molecular dynamics (MD) is a powerful technique that allows the protonation state of residues to change dynamically, thereby enabling the study of pH dependence in a manner that has not been possible before. Recently, a constant pH implementation was incorporated into the GROMACS MD package. Although this implementation provides good accuracy and performance, manual modification and the preparation of simulation input files are required, which can be complicated, tedious, and prone to errors. To simplify and automate the setup process, we present phbuilder, a tool that automatically prepares constant pH MD simulations for GROMACS by modifying the input structure and topology as well as generating the necessary parameter files. phbuilder can prepare constant pH simulations from both initial structures and existing simulation systems, and it also provides functionality for performing titrations and single-site parametrizations of new titratable group types. The tool is freely available at www.gitlab.com/gromacs-constantph. We anticipate that phbuilder will make constant pH simulations easier to set up, thereby making them more accessible to the GROMACS user community.
Thousands of anthropogenic chemicals are released into the environment each year, posing potential hazards to human and environmental health. Toxic chemicals may cause a variety of adverse health effects, triggering immediate symptoms or delayed effects over longer periods of time. It is thus crucial to develop methods that can rapidly screen and predict the toxicity of chemicals to limit the potential harmful impacts of chemical pollutants. Computational methods are being increasingly used in toxicity predictions. Here, the method of molecular docking is assessed for screening potential toxicity of a variety of xenobiotic compounds, including pesticides, pharmaceuticals, pollutants, and toxins derived from the chemical industry. The method predicts the binding energy of pollutants to a set of carefully selected receptors under the assumption that toxicity in many cases is related to interference with biochemical pathways. The strength of the applied method lies in its rapid generation of interaction maps between potential toxins and the targeted enzymes, which could quickly yield molecular-level information and insight into potential perturbation pathways, aiding in the prioritization of chemicals for further tests. Two scoring functions are compared: Autodock Vina and the machine-learning scoring function RF-Score-VS. The results are promising, although hampered by the accuracy of the scoring functions. The strengths and weaknesses of the docking protocol are discussed, as well as future directions for improving the accuracy for the purpose of toxicity predictions.
We define a molecular caging complex as a pair of molecules in which one molecule (the "host" or "cage") possesses a cavity that can encapsulate the other molecule (the "guest") and prevent it from escaping. Molecular caging complexes can be useful in applications such as molecular shape sorting, drug delivery, and molecular immobilization in materials science, to name just a few. However, the design and computational discovery of new caging complexes is a challenging task, as it is hard to predict whether one molecule can encapsulate another because their shapes can be quite complex. In this paper, we propose a computational screening method that predicts whether a given pair of molecules form a caging complex. Our method is based on a caging verification algorithm that was designed by our group for applications in robotic manipulation. We tested our algorithm on three pairs of molecules that were previously described in a pioneering work on molecular caging complexes and found that our results are fully consistent with the previously reported ones. Furthermore, we performed a screening experiment on a data set consisting of 46 hosts and four guests and used our algorithm to predict which pairs are likely to form caging complexes. Our method is computationally efficient and can be integrated into a screening pipeline to complement experimental techniques.
Human cytochrome P450 3A4 (CYP3A4) is responsible for the metabolism of similar to 50% clinically used drugs. Midazolam (MDZ) is a commonly used sedative drug and serves as a marker substrate for the CYP3A4 activity assessment. MDZ is metabolized by CYP3A4 to two hydroxylation products, 1'-OH-MDZ and 4-OH-MDZ. It has been reported that the ratio of 1'-OH-MDZ and 4-OH-MDZ is dependent on the MDZ concentration, which reflects the homotropic cooperative behavior in MDZ metabolism by CYP3A4. Here, we used quantum chemistry (QC), molecular docking, conventional molecular dynamics (cMD), and Gaussian accelerated molecular dynamics (GaMD) approaches to investigate the mechanism of the interactions between CYP3A4 and MDZ. QC calculations suggest that C1' is less reactive for hydroxylation than C4, which is a pro-chirality carbon. However, the 4-OH-MDZ product is likely to be racemic due to the chirality inversion in the rebound step. The MD simulation results indicate that MDZ at the peripheral allosteric site is not stable and the binding modes of the MDZ molecules at the productive site are in line with the experimental observations.
The plasticity of cytochromes P450 (P450s) is known to contribute significantly to their catalytic capacity of metabolizing various substrates. Although numerous studies have been performed, factors governing the plasticity and dynamics of P450s are still not fully understood. In this study, taking CYP2B4 as an example, we dissect the protein plasticity and dynamics in different environments. CYP2B4 is featured by a high degree of plasticity, which exhibits open, closed, and intermediate states. By analyzing the CYP2B4 crystal structures, we identified the structural features for the closed, open, and intermediate states. Interestingly, formation of the dimer structure was found in the open and intermediate states. The subsequent molecular dynamics (MD) simulations of the open structure in water confirmed the importance of the dimer form in stabilizing the open conformations. MD simulations of the closed and open structures in the membrane environment and the free energies for opening the F-G cassette obtained from the umbrella sampling calculations indicate that the membrane environment is important for stabilizing the F-G cassette. The dynamical network analysis indicates that Asp105 on the B-C loop plays an important role in transiting the structure from the open to the intermediate state. Our results thus unveil the mechanisms of dimer formation and open-to-intermediate transition for CYP2B4 in the water and membrane environments.
Sirtuins are a family of nicotinamide adenine dinucleotide (NAD(+))-dependent enzymes, which undergo robust deacetylase activity, resulting in the production of nicotinamide. It is well known that nicotinamide, which is one of the products, can also act as an inhibitor for further deacetylation process by forming NAD(+) again. Hence, the removal of nicotinamide from sirtuins is a demanding process, and the mechanistic understanding of the process remains elusive. In this investigation, we have made an attempt to unravel the unbinding pathways of nicotinamide from SIRT1, SIRT2, and SIRT3 (SIRT1-3) using Random Acceleration Molecular Dynamics (RAMD) Simulations, and we have successfully identified various unbinding channels. The selectivity of the egression channel is determined by using a thorough analysis of the frequency of egression trajectories. Similarly, various inhibitors have been docked with the active sites of SIRT1-3, and their egression pathways have been investigated to understand whether they follow the same egression pathway as that of nicotinamide. The residues that are responsible for the unbinding pathways have been determined from the analysis of RAMD trajectories. From these results, it is clear that phenylalanine and histidine residues play major roles in the egression of inhibitors. Additionally, the key residues Leu, Pro, Met, Phe, Tyr, and Ile are found to control the release by acting as gateway residues. The role of these residues from different egression channels has been studied by carrying out mutations with alanine residue. This is the first report on sirtuins, which demonstrates the novel unbinding pathways for nicotinamide/inhibitors. This work provides new insights for developing more promising SIRT1-3 inhibitors.
Monoamine oxidase B (MAO-B) is a potential biomarker for Parkinson's disease (PD), a neurodegenerative disease associated with the loss of motor activities in human subjects. The disease state is associated with dopamine deprival, and so the inhibitors of MAO-B can serve as therapeutic drugs for PD. Since the expression level of MAO-B directly correlates to the disease progress, the distribution and population of this enzyme can be employed to monitor disease development. One of the approaches available for estimating the population is two-photon imaging. The ligands used for two-photon imaging should have high binding affinity and binding specificity toward MAO-B along with significant two-photon absorption cross sections when they are bound to the target. In this article, we study using a multiscale modeling approach, the binding affinity and spectroscopic properties (one-and two-photon absorption) of three (Flu1, Flu2, Flu3) of the currently available probes for monitoring the MAO-B level. We report that the binding affinity of the probes can be explained using the molecular size and binding cavity volume. The experimentally determined one-photon absorption spectrum is well reproduced by the employed QM/MM approaches, and the most accurate spectral shifts, on passing from one probe to another, are obtained at the coupled-cluster (CC2) level of theory. An important conclusion from this study is also the demonstration that intrinsic molecular two-photon absorption strengths (delta(2PA)) increase in the order delta(2PA) (Flu1) > delta(2PA) (Flu2) > delta(2PA) (Flu3). This is in contrast with experimental data, which predict similar values of two-photon absorption cross sections for Flu1 and Flu3. We demontrate, based on the results of electronic-structure calculations for Flu1 that this discrepancy cannot be explained by an explicit account for neighboring residues (which could lead to charge transfer between a probe and neighboring aromatic amino acids thus boosting delta(2PA)). In summary, we show that the employed multiscale approach not only can optimize two-photon absorption properties and verify binding affinity, but it can also help in detailed analyses of experimental data.
Uncertainty was introduced to chemical descriptors of 16 publicly available data. sets to various degrees and in various-ways order to investigate the effect on the predictive performance Of the State of-the-art method decision tree ensembles. A number of strategies to handle uncertainty in:decision tree ensembles were evaluated. The main conclusion of the Study. is that uncertainty to a large extent may be introduced in chemical descriptors without. impairing the predictive performance of ensembles and without the predictive performance being significantly reduced from a practical point of view. The investigation. further showed that even When distributions of uncertain values were provided, the ensembles method could generate equally effective models from single-point samples from these distributions. Hence, there seems to be no advantage in using more., elaborate Methods for handling uncertainty in chemical descriptors when using decision tree ensembles as a modeling method for the considered types of introduced uncertainty.
Free fatty acid receptor 1 (FFAR1) is a potential therapeutic target for the treatment of type 2 diabetes (T2D). It has been validated that agonists targeting FFAR1 can achieve the initial therapeutic endpoints of T2D, and the epimer agonists (R,S) AM-8596 can activate FFAR1 differently, with one acting as a partial agonist and the other as a full agonist. Up to now, the origin of the stereoselectivity of FFAR1 agonists remains elusive. In this work, we used molecular simulation methods to elucidate the mechanism of the stereoselectivity of the FFAR1 agonists (R)-AM-8596 and (S)-AM-8596. We found that the full agonist (R)-AM-8596 disrupts the residue interaction network around the receptor binding pocket and promotes the opening of the binding site for the G-protein, thereby resulting in the full activation of FFAR1. In contrast, the partial agonist (S)-AM-8596 forms stable electrostatic interactions with FFAR1, which stabilizes the residue network and hinders the conformational transition of the receptor. Our work thus clarifies the selectivity and underlying molecular activation mechanism of FFAR1 agonists.
Our skin constitutes an effective permeability barrier that protects the body from exogenous substances but concomitantly severely limits the number of pharmaceutical drugs that can be delivered transdermally. In topical formulation design, chemical permeation enhancers (PEs) are used to increase drug skin permeability. In vitro skin permeability experiments can measure net effects of PEs on transdermal drug transport, but they cannot explain the molecular mechanisms of interactions between drugs, permeation enhancers, and skin structure, which limits the possibility to rationally design better new drug formulations. Here we investigate the effect of the PEs water, lauric acid, geraniol, stearic acid, thymol, ethanol, oleic acid, and eucalyptol on the transdermal transport of metronidazole, caffeine, and naproxen. We use atomistic molecular dynamics (MD) simulations in combination with developed molecular models to calculate the free energy difference between 11 PE-containing formulations and the skin’s barrier structure. We then utilize the results to calculate the final concentration of PEs in skin. We obtain an RMSE of 0.58 log units for calculated partition coefficients from water into the barrier structure. We then use the modified PE-containing barrier structure to calculate the PEs’ permeability enhancement ratios (ERs) on transdermal metronidazole, caffeine, and naproxen transport and compare with the results obtained from in vitro experiments. We show that MD simulations are able to reproduce rankings based on ERs. However, strict quantitative correlation with experimental data needs further refinement, which is complicated by significant deviations between different measurements. Finally, we propose a model for how to use calculations of the potential of mean force of drugs across the skin’s barrier structure in a topical formulation design.
Suicide inhibition of the CYP3A4 enzyme by a drug inactivates the enzyme in the drug biotransformation process and often shows safety concerns about the drug. Despite extensive experimental studies, the abnormal molecular mechanism of a suicide inhibitor that forms a covalent bond with the residue far away from the catalytically active center of CYP3A4 inactivating the enzyme remains elusive. Here, the authors used molecular simulation approaches to study in detail how diquinone methide (DQR), the metabolite product of raloxifene, unbinds from CYP3A4 and inactivates the enzyme at the atomistic level. The results dearly indicate that in one of the intermediate states formed in its unbinding process, DQR covalently binds to Cys239, a residue far away from the catalytically active center of CYP3A4, and hinders the substrate from entering or leaving the enzyme. This work therefore provides an unprecedented way of clarifying the abnormal mechanism of suicide inhibition of the CYP3A4 enzyme.
Understanding unbinding kinetics of protein-ligand systems is of great importance for the design of ligands with desired specificity and safety. In recent years, enhanced sampling techniques have emerged as effective tools for studying unbinding kinetics of protein-ligand systems at the atomistic level. However, in many protein-ligand systems, the ligand unbinding processes are strongly coupled to protein conformational changes and the disclosure of the hidden degrees of freedom closely related to the protein conformational changes so that sampling is enhanced over these degrees of freedom remains a great challenge. Here, we show how potential-scaled molecular dynamics (sMD) and infrequent metadynamics (InMetaD) simulation techniques can be combined to successfully reveal the unbinding mechanism of 3-(1,4-diazabicyclo[3.2.2]nonan-4-yl)-6-[F-18]fluorodibenzo[b,d]thiophen e 5,5-dioxide ([F-18]ASEM) from a chimera structure of the alpha 7-nicotinic acetylcholine receptor. By using sMD simulations, we disclosed that the "close to "open" conformational change of loop C plays a key role in the ASEM unbinding process. By carrying out InMetaD simulations with this conformational change taken into account as an additional collective variable, we further captured the key states in the unbinding process and clarified the unbinding mechanism of ASEM from the protein. Our work indicates that combining sMD and InMetaD simulation techniques can be an effective approach for revealing the unbinding mechanism of a protein-ligand system where protein conformational changes control the unbinding process.
Intrinsically disordered proteins (IDPs) exert their functions by binding to partner proteins via a complex process that includes coupled folding and binding. Because inhibiting the binding of the IDP p53 to its partner MDM2 has become a promising strategy for the design of anticancer drugs, we carried out metadynamics simulations to study the coupled folding and binding process linking the IDP p53 to MDM2 in atomic detail. Using bias-exchange metadynamics (BE-MetaD) and infrequent metadynamics (InMetaD), we estimated the binding free energy, the unbinding rate, and the binding rate. By analyzing the stable intermediates, we uncovered the role non-native interactions played in the p53-MDM2 binding/unbinding process. We used a three-state model to describe the whole binding/unbinding process and to obtain the corresponding rate constants. Our work shows that the binding of p53 favors an induced-fit mechanism which proceeds in a stepwise fashion. Our results can be helpful for gaining an in-depth understanding of the coupled folding and binding process needed for the design of MDM2 inhibitors.