Long-timescale Molecular Dynamics Simulations Research Articles

BaNDyT: Bayesian Network modeling of molecular Dynamics Trajectories.

Bayesian network modeling (BN modeling, or BNM) is an interpretable machine learning method for constructing probabilistic graphical models from the data. In recent years, it has been extensively applied to diverse types of biomedical datasets. Concurrently, our ability to perform long-timescale molecular dynamics (MD) simulations on proteins and other materials has increased exponentially. However, the analysis of MD simulation trajectories has not been data-driven but rather dependent on the user's prior knowledge of the systems, thus limiting the scope and utility of the MD simulations. Recently, we pioneered using BNM for analyzing the MD trajectories of protein complexes. The resulting BN models yield novel fully data-driven insights into the functional importance of the amino acid residues that modulate proteins' function. In this report, we describe the BaNDyT software package that implements the BNM specifically attuned to the MD simulation trajectories data. We believe that BaNDyT is the first software package to include specialized and advanced features for analyzing MD simulation trajectories using a probabilistic graphical network model. We describe here the software's uses, the methods associated with it, and a comprehensive Python interface to the underlying generalist BNM code. This provides a powerful and versatile mechanism for users to control the workflow. As an application example, we have utilized this methodology and associated software to study how membrane proteins, specifically the G protein-coupled receptors, selectively couple to G proteins. The software can be used for analyzing MD trajectories of any protein as well as polymeric materials.

Non-canonical amino acids uncover the significant impact of Tyr671 on Taq DNA polymerase catalytic activity.

Responsible for synthesizing the complementary strand of the DNA template, DNA polymerase is a crucial enzyme in DNA replication, recombination and repair. A highly conserved tyrosine (Tyr), located at the C-terminus of the O-helix in family A DNA polymerases, plays a critical role in enzyme activity and fidelity. Here, we combined the technology of genetic code extension to incorporate non-canonical amino acids and molecular dynamics (MD) simulations to uncover the mechanisms by which Tyr671 impacts substrate binding and conformation transitions in a DNA polymerase from Thermus aquaticus. Five non-canonical amino acids, namely l-3,4-dihydroxyphenylalanine (l-DOPA), p-aminophenylalanine (pAF), p-acetylphenylalanine (pAcF), p-cyanophenylalanine (pCNF) and p-nitrophenylalanine (pNTF), were individually incorporated at position 671. Strikingly, Y671pAF and Y671DOPA were active, but with lower activity compared to Y671F and wild-type. Y671pAF showed a higher fidelity than the Y671F, despite both possessing lower fidelity than the wild-type. Metadynamics and long-timescale MD simulations were carried out to probethe role of mutations in affecting protein structure, including open conformation, open-to-closed conformation transition, closed conformation, and closed-to-open conformation transition. The MD simulations clearly revealed that the size of the 671 amino acid residue and interactions with substrate or nearby residues were critical for Tyr671 to determine enzyme activity and fidelity.

Discovery of lirafugratinib (RLY-4008), a highly selective irreversible small-molecule inhibitor of FGFR2

Fibroblast growth factor receptor (FGFR) kinase inhibitors have been shown to be effective in the treatment of intrahepatic cholangiocarcinoma and other advanced solid tumors harboring FGFR2 alterations, but the toxicity of these drugs frequently leads to dose reduction or interruption of treatment such that maximum efficacy cannot be achieved. The most common adverse effects are hyperphosphatemia caused by FGFR1 inhibition and diarrhea due to FGFR4 inhibition, as current therapies are not selective among the FGFRs. Designing selective inhibitors has proved difficult with conventional approaches because the orthosteric sites of FGFR family members are observed to be highly similar in X-ray structures. In this study, aided by analysis of protein dynamics, we designed a selective, covalent FGFR2 inhibitor. In a key initial step, analysis of long-timescale molecular dynamics simulations of the FGFR1 and FGFR2 kinase domains allowed us to identify differential motion in their P-loops, which are located adjacent to the orthosteric site. Using this insight, we were able to design orthosteric binders that selectively and covalently engage the P-loop of FGFR2. Our drug discovery efforts culminated in the development of lirafugratinib (RLY-4008), a covalent inhibitor of FGFR2 that shows substantial selectivity over FGFR1 (~250-fold) and FGFR4 (~5,000-fold) in vitro, causes tumor regression in multiple FGFR2-altered human xenograft models, and was recently demonstrated to be efficacious in the clinic at doses that do not induce clinically significant hyperphosphatemia or diarrhea.

Molecular dynamics simulations of HIV-1 matrix-membrane interactions at different stages of viral maturation

Although the structural rearrangement of the membrane-bound matrix (MA) protein trimers upon HIV-1 maturation has been reported, the consequences of MA maturation on the MA-lipid interactions are not well understood. Long-timescale molecular dynamics simulations of the MA multimeric assemblies of immature and mature virus particles with our realistic asymmetric membrane model have explored MA-lipid interactions and lateral organization of lipids around MA complexes. The number of stable MA-phosphatidylserine and MA-phosphatidylinositol 4,5-bisphosphate (PIP2) interactions at the trimeric interface of the mature MA complex is observed to be greater compared to that of the immature MA complex. Our simulations identified an alternative PIP2-binding site in the immature MA complex where the multivalent headgroup of a PIP2 lipid with a greater negative charge binds to multiple basic amino acid residues such as ARG3 residues of both the MA monomers at the trimeric interface and highly basic region (HBR) residues (LYS29, LYS31) of one of the MA monomers. Our enhanced sampling simulations have explored the conformational space of phospholipids at different binding sites of the trimer-trimer interface of MA complexes that are not accessible by conventional unbiased molecular dynamics. Unlike the immature MA complex, the 2′ acyl tail of two PIP2 lipids at the trimeric interface of the mature MA complex is observed to sample stable binding pockets of MA consisting of helix-4 residues. Together, our results provide molecular-level insights into the interactions of MA trimeric complexes with membrane and different lipid conformations at the specific binding sites of MA protein before and after viral maturation.

In silico screening of natural products as uPAR inhibitors via multiple structure-based docking and molecular dynamics simulations

Cancer remains one of the most pressing challenges to global healthcare, exerting a significant impact on patient life expectancy. Cancer metastasis is a critical determinant of the lethality and treatment resistance of cancer. The urokinase-type plasminogen activator receptor (uPAR) shows great potential as a target for anticancer and antimetastatic therapies. In this work, we aimed to identify potential uPAR inhibitors by structural dynamics-based virtual screenings against a natural product library on four representative apo-uPAR structural models recently derived from long-timescale molecular dynamics (MD) simulations. Fifteen potential inhibitors (NP1-NP15) were initially identified through molecular docking, consensus scoring, and visual inspection. Subsequently, we employed MD-based molecular mechanics-generalized Born surface area (MM-GBSA) calculations to evaluate their binding affinities to uPAR. Structural dynamics analyses further indicated that all of the top 6 compounds exhibited stable binding to uPAR and interacted with the critical residues in the binding interface between uPAR and its endogenous ligand uPA, suggesting their potential as uPAR inhibitors by interrupting the uPAR-uPA interaction. We finally predicted the ADMET properties of these compounds. The natural products NP5, NP12, and NP14 with better binding affinities to uPAR than the uPAR inhibitors previously discovered by us were proven to be potentially orally active in humans. This work offers potential uPAR inhibitors that may contribute to the development of novel effective anticancer and antimetastatic therapeutics. Communicated by Ramaswamy H. Sarma

Tyr 118-Mediated Electron Transfer is Key to the Chlorite Decomposition in Heme-Dependent Chlorite Dismutase.

Chlorite dismutase (Cld) is a crucial enzyme that catalyzes the decomposition of chlorite ions into chloride ions (Cl-) and molecular oxygen (O2). Despite playing an important role in the detoxification of toxic chlorite ions, the mechanism of cleavage of the Cl-O bond by Cld remains highly debatable. The present study highlights the mechanism of such Cl-O bond cleavage in Cld using sophisticated computational tools such as hybrid quantum mechanical/molecular mechanical calculations and long-time scale molecular dynamics simulations. Here, we show that Cld forms a high spin ferric hexacoordinated substrate adduct in the presence of a chlorite ion, which subsequently reduces to a ferrous state. Our study shows a stepwise pathway with the homolytic cleavage of the Cl-O bond that produces a high spin Fe(III)-OH species and a diradicaloid species formed by the combination of a chlorine-based ClO• radical and a protein-based tyrosine118• radical. The findings provide significant insights into Cl-O bond cleavage and O2 formation which shows a crucial role of the tyrosine118 during the electron transfer process.

Effect of AQP4 and its palmitoylation on the permeability of exogenous reactive oxygen species: Insights from computational study

Aquaporin 4 (AQP4) facilitates the transport of reactive oxygen species (ROS). Both cancer cells and the ionizing radiation microenvironment can induce posttranslational modifications (PTMs) in AQP4, which may affect its permeability to ROS. Because this ROS diffusion process is rapid, microscopic, and instantaneous within and outside cells, conventional experimental methods are inadequate for elucidating the molecular mechanisms involved. In this study, computational methods were employed to investigate the permeability of exogenous ROS mediated by radiation in AQP4 at a molecular scale. We constructed a simulation system incorporating AQP4 and AQP4-Cysp13 in a complex lipid environment with ROS. Long-timescale molecular dynamics simulations were conducted to assess the structural stability of both AQP4 and AQP4-Cysp13. Free energy calculations were utilized to determine the ROS transport capability of the two AQP4 proteins. Computational electrophysiology and channel structural analysis quantitatively evaluated changes in ROS transport capacity under various radiation-induced transmembrane voltage microenvironments. Our findings demonstrate the distinct transport capabilities of AQP4 channels for water molecules and various types of ROS and reveal a decrease in transport efficiency when AQP4 undergoes palmitoylation modification. In addition, we have simulated the radiation-induced alteration of cell membrane voltage, which significantly affected the ROS transport capacity. We propose that this research will enhance the understanding of the molecular mechanisms governing the transport of exogenous ROS by AQP4 and elucidate the influence of palmitoylation on ROS transport. This study will also help clarify how different structural features of AQP4 affect the transport of exogenous ROS mediated by radiotherapy, thereby providing a theoretical molecular basis for the development of new treatment strategies that combine with radiotherapy.

Elucidating protein-ligand binding kinetics based on returning probability theory.

The returning probability (RP) theory, a rigorous diffusion-influenced reaction theory, enables us to analyze the binding process systematically in terms of thermodynamics and kinetics using molecular dynamics (MD) simulations. Recently, the theory was extended to atomistically describe binding processes by adopting the host-guest interaction energy as the reaction coordinate. The binding rate constants can be estimated by computing the thermodynamic and kinetic properties of the reactive state existing in the binding processes. Here, we propose a methodology based on the RP theory in conjunction with the energy representation theory of solution, applicable to complex binding phenomena, such as protein-ligand binding. The derived scheme of calculating the equilibrium constant between the reactive and dissociate states, required in the RP theory, can be used for arbitrary types of reactive states. We apply the present method to the bindings of small fragment molecules [4-hydroxy-2-butanone (BUT) and methyl methylthiomethyl sulphoxide (DSS)] to FK506 binding protein (FKBP) in an aqueous solution. Estimated binding rate constants are consistent with those obtained from long-timescale MD simulations. Furthermore, by decomposing the rate constants to the thermodynamic and kinetic contributions, we clarify that the higher thermodynamic stability of the reactive state for DSS causes the faster binding kinetics compared with BUT.

Functional dynamics and allosteric modulation of TRPA1

The cation channel TRPA1 is a potentially important drug target, and characterization of TRPA1 functional dynamics might help guide structure-based drug design. Here, we present results from long-timescale molecular dynamics simulations of TRPA1 with an allosteric activator, allyl isothiocyanate (AITC), in which we observed spontaneous transitions from a closed, non-conducting channel conformation into an open, conducting conformation. Based on these transitions, we propose a gating mechanism in which movement of a regulatory TRP-like domain allosterically translates into pore opening in a manner reminiscent of pore opening in voltage-gated ion channels. In subsequent experiments, we found that mutations that disrupt packing of the S4-S5 linker-TRP-like domain and the S5 and S6 helices also affected channel activity. In simulations, we also observed A-967079, a known allosteric inhibitor, binding between helices S5 and S6, suggesting that A-967079 may suppress activity by stabilizing a non-conducting pore conformation-a finding consistent with our proposed gating mechanism.

Advanced virtual screening enables the discovery of a host-targeting and broad-spectrum antiviral agent

We employed an advanced virtual screening (AVS) approach to identify potential inhibitors of human dihydroorotate dehydrogenase (DHODH), a validated target for development of broad-spectrum antivirals. We screened a library of 495118 compounds and identified 495 compounds that exhibited better binding scores than the reference ligands involved in the screening. From the top 100 compounds, we selected 28 based on their consensus docking scores and structural novelty. Then, we conducted in vitro experiments to investigate the antiviral activity of selected compounds on HSV-1 infection, which is susceptible to DHODH inhibitors. Among the tested compounds, seven displayed statistically significant antiviral effects, with Comp 19 being the most potent inhibitor. We found that Comp 19 exerted its antiviral effect in a dose-dependent manner (IC50 = 1.1 μM) and exhibited the most significant antiviral effect when added before viral infection. In the biochemical assay, Comp 19 inhibited human DHODH in a dose-dependent manner with the IC50 value of 7.3 μM. Long-timescale molecular dynamics simulations (1000 ns) revealed that Comp 19 formed a very stable complex with human DHODH. Comp 19 also displayed broad-spectrum antiviral activity and suppressed cytokine production in THP-1 cells. Overall, our study provides evidence that AVS could be successfully implemented to discover novel DHODH inhibitors with broad-spectrum antiviral activity.

A Conserved Local Structural Motif Controls the Kinetics of PTP1B Catalysis.

Protein tyrosine phosphatase 1B (PTP1B) is a negative regulator of the insulin and leptin signaling pathways, making it a highly attractive target for the treatment of type II diabetes. For PTP1B to perform its enzymatic function, a loop referred to as the "WPD loop" must transition between open (catalytically incompetent) and closed (catalytically competent) conformations, which have both been resolved by X-ray crystallography. Although prior studies have established this transition as the rate-limiting step for catalysis, the transition mechanism for PTP1B and other PTPs has been unclear. Here we present an atomically detailed model of WPD loop transitions in PTP1B based on unbiased, long-timescale molecular dynamics simulations and weighted ensemble simulations. We found that a specific WPD loop region─the PDFG motif─acted as the key conformational switch, with structural changes to the motif being necessary and sufficient for transitions between long-lived open and closed states of the loop. Simulations starting from the closed state repeatedly visited open states of the loop that quickly closed again unless the infrequent conformational switching of the motif stabilized the open state. The functional importance of the PDFG motif is supported by the fact that it is well conserved across PTPs. Bioinformatic analysis shows that the PDFG motif is also conserved, and adopts two distinct conformations, in deiminases, and the related DFG motif is known to function as a conformational switch in many kinases, suggesting that PDFG-like motifs may control transitions between structurally distinct, long-lived conformational states in multiple protein families.

Discovery and Validation of the Binding Poses of Allosteric Fragment Hits to Protein Tyrosine Phosphatase 1b: From Molecular Dynamics Simulations to X-ray Crystallography.

Fragment-based drug discovery has led to six approved drugs, but the small sizes of the chemical fragments used in such methods typically result in only weak interactions between the fragment and its target molecule, which makes it challenging to experimentally determine the three-dimensional poses fragments assume in the bound state. One computational approach that could help address this difficulty is long-timescale molecular dynamics (MD) simulations, which have been used in retrospective studies to recover experimentally known binding poses of fragments. Here, we present the results of long-timescale MD simulations that we used to prospectively discover binding poses for two series of fragments in allosteric pockets on a difficult and important pharmaceutical target, protein tyrosine phosphatase 1b (PTP1b). Our simulations reversibly sampled the fragment association and dissociation process. One of the binding pockets found in the simulations has not to our knowledge been previously observed with a bound fragment, and the other pocket adopted a very rare conformation. We subsequently obtained high-resolution crystal structures of members of each fragment series bound to PTP1b, and the experimentally observed poses confirmed the simulation results. To the best of our knowledge, our findings provide the first demonstration that MD simulations can be used prospectively to determine fragment binding poses to previously unidentified pockets.

Computational Insights into Prostaglandin E2 Ligand Binding and Activation of G-Protein-Coupled Receptors.

G-protein coupled receptors (GPCRs) are eukaryotic integral membrane proteins that regulate signal transduction cascade pathways implicated in a variety of human diseases and are consequently of interest as drug targets. For this reason, it is of interest to investigate the way in which specific ligands bind and trigger conformational changes in the receptor during activation and how this in turn modulates intracellular signaling. In the present study, we investigate the way in which the ligand Prostaglandin E2 interacts with three GPCRs in the E-prostanoid family: EP1, EP2, and EP3. We examine information transfer pathways based on long-time scale molecular dynamics simulations using transfer entropy and betweenness centrality to measure the physical transfer of information among residues in the system. We monitor specific residues involved in binding to the ligand and investigate how the information transfer behavior of these residues changes upon ligand binding. Our results provide key insights that enable a deeper understanding of EP activation and signal transduction functioning pathways at the molecular level, as well as enabling us to make some predictions about the activation pathway for the EP1 receptor, for which little structural information is currently available. Our results should advance ongoing efforts in the development of potential therapeutics targeting these receptors.

Desensitization dynamics of the AMPA receptor

To perform their physiological functions, amino methyl propionic acid receptors (AMPARs) cycle through active, resting, and desensitized states, and dysfunction in AMPAR activity is associated with various neurological disorders. Transitions among AMPAR functional states, however, are largely uncharacterized at atomic resolution and are difficult to examine experimentally. Here, we report long-timescale molecular dynamics simulations of dimerized AMPAR ligand-binding domains (LBDs), whose conformational changes are tightly coupled to changes in AMPAR functional states, in which we observed LBD dimer activation and deactivation upon ligand binding and unbinding at atomic resolution. Importantly, we observed the ligand-bound LBD dimer transition from the active conformation to several other conformations, which may correspond with distinct desensitized conformations. We also identified a linker region whose structural rearrangements heavily affected the transitions to and among these putative desensitized conformations, and confirmed, using electrophysiology experiments, the importance of the linker region in these functional transitions.

Structural Dynamics-Driven Discovery of Anticancer and Antimetastatic Effects of Diltiazem and Glibenclamide Targeting Urokinase Receptor.

Diltiazem and glibenclamide are commonly used hypotensive and antidiabetic drugs. This study reports the discovery of the potential antitumor and antimetastatic effects of these two drugs using a structural dynamics-driven virtual screening targeting urokinase receptor (uPAR). Owing to uPAR's high flexibility, currently resolved crystal structures of uPAR, all in ligand-bound states, provide limited representations of its physiological conformation. To improve the accuracy of screening, we performed a long-timescale molecular dynamics simulation and obtained the representative conformations of apo-uPAR as the targets for our screening. Experimentally, we demonstrated that diltiazem and glibenclamide bound uPAR with KD values in the micromolar range. In addition, both compounds effectively suppressed tumor growth and metastasis in a uPAR-dependent manner in vitro and in vivo. This work not only provides two potent uPAR inhibitors but also reports a proof-of-concept study on the potential off-label antitumor and antimetastatic uses of diltiazem and glibenclamide.

Elucidating the origins of acyltransferase selectivity in primary and secondary metabolism.

Carrier-protein-S-acyltransferases (ATs) are enzymes that load acyl groups onto carrier proteins (CP). These catalysts can exist as components of larger megasynth(et)ases that act in-cis (i.e., intramolecularly) or as so-called trans-ATs that act “in-trans”, forming intermolecular complexes with a corresponding CP or CP domain. A growing number of crystal structures of trans-ATs, including have been resolved in recent years. In many cases, these structures feature a CP covalently tethered to its cognate AT, revealing the specific protein-protein interactions and conformational reorganization that mediate CP-AT binding and catalysis. To complement insights offered by these experimental structures, we have employed long-timescale molecular dynamics simulations to provide a dynamical understanding of biomolecular recognition and acyltransferase function. These simulations reveal key differences that may explain the distinct specificities of trans-AT that operate in primary and secondary metabolism, while also demonstrating that certain internal rigid body motions of the AT that modulate active site access are common among all trans-ATs examined herein.

Gating and modulation of an inward-rectifier potassium channel.

Inward-rectifier potassium channels (Kirs) are lipid-gated ion channels that differ from other K+ channels in that they allow K+ ions to flow more easily into, rather than out of, the cell. Inward rectification is known to result from endogenous magnesium ions or polyamines (e.g., spermine) binding to Kirs, resulting in a block of outward potassium currents, but questions remain regarding the structural and dynamic basis of the rectification process and lipid-dependent channel activation. Here, we present the results of long-timescale molecular dynamics simulations starting from a crystal structure of phosphatidylinositol 4,5-bisphosphate (PIP2)-bound chicken Kir2.2 with a non-conducting pore. After introducing a mutation (G178R) that is known to increase the open probability of a homologous channel, we were able to observe transitions to a stably open, ion-conducting pore, during which key conformational changes occurred in the main activation gate and the cytoplasmic domain. PIP2 binding appeared to increase stability of the pore in its open and conducting state, as PIP2 removal resulted in pore closure, with a median closure time about half of that with PIP2 present. To investigate structural details of inward rectification, we simulated spermine binding to and unbinding from the open pore conformation at positive and negative voltages, respectively, and identified a spermine-binding site located near a previously hypothesized site between the pore cavity and the selectivity filter. We also studied the effects of long-range electrostatics on conduction and spermine binding by mutating charged residues in the cytoplasmic domain and found that a finely tuned charge density, arising from basic and acidic residues within the cytoplasmic domain, modulated conduction and rectification.

Toward long, all-atom simulation of the E. coli 70S ribosome

The precise machinery of large ribosomal complexes underlies the fundamental process of translation from RNA to protein in biology. In order for an amino acid to be incorporated into a growing polypeptide chain in both an accurate and efficient manner, the various components of the ribosome must coordinate in a dynamic interplay that moves the ribosome through the states of the elongation cycle. A deeper understanding of these processes at atomic resolution can potentially be achieved through the use of long-timescale molecular dynamics simulations.

A comprehensive evaluation of the potential binding poses of fentanyl and its analogs at the µ-opioid receptor

Fentanyl and its analogs are selective agonists of the µ-opioid receptor (MOR). Among novel synthetic opioids (NSOs), they dominate the recreational drug market and are the main culprits for the opioid crisis, which has been exacerbated by the COVID-19 pandemic. By taking advantage of the crystal structures of the MOR, several groups have investigated the binding mechanism of fentanyl, but have not reached a consensus, in terms of both the binding orientation and the fentanyl conformation. Thus, the binding mechanism of fentanyl at the MOR remains an unsolved and challenging question. Here, we carried out a systematic computational study to investigate the preferred fentanyl conformations, and how these conformations are being accommodated in the MOR binding pocket. We characterized the free energy landscape of fentanyl conformations with metadynamics simulations, and compared and evaluated several possible fentanyl binding conditions in the MOR with long-timescale molecular dynamics simulations. Our results indicate that the most preferred binding pose in the MOR binding pocket corresponds well with the global minimum on the energy landscape of fentanyl in the absence of the receptor, while the energy landscape can be reconfigured by modifying the fentanyl scaffold. The interactions with the receptor may stabilize a slightly unfavored fentanyl conformation in an alternative binding pose. By extending similar investigations to fentanyl analogs, our findings establish a structure–activity relationship of fentanyl binding at the MOR. In addition to providing a structural basis to understand the potential toxicity of the emerging NSOs, such insights will contribute to developing new, safer analgesics.

Probing the formation, structure and free energy relationships of M protein dimers of SARS-CoV-2

The M protein of the novel coronavirus 2019 (SARS-CoV-2) is the major structural component of the viral envelope and is also the minimum requirement for virus particle budding. M proteins generally exist as dimers. In virus assembly, they are the main driving force for envelope formation through lateral interactions and interactions with other viral structural proteins that play a central role. We built 100 candidate models and finally analyzed the six most convincing structural features of the SARS-CoV-2 M protein dimer based on long-timescale molecular dynamics (MD) simulations, multiple free energy analyses (potential mean force (PMF) and molecular mechanics Poisson-Boltzmann surface area (MMPBSA)) and principal component analysis (PCA) to obtain the most reasonable structure. The dimer stability was found to depend on the Leu-Ile zipper motif and aromatic amino acids in the transmembrane domain (TMD). Furthermore, the C-terminal domain (CTD) effects were relatively small. These results highlight a model in which there is sufficient binding affinity between the TMDs of M proteins to form dimers through the residues at the interface of the three transmembrane helices (TMHs). This study aims to help find more effective inhibitors of SARS-CoV-2 M dimers and to develop vaccines based on structural information.