Avoiding False Positive Conclusions in Molecular Simulation: The Importance of Replicas

Identifying potential Alzheimer's disease therapeutics through GSK-3β inhibition: A molecular docking and dynamics approach

Emerging as a promising drug target for Alzheimer's disease (AD) therapy, glycogen synthase kinase 3β (GSK-3β) has garnered attention. This study sought to rigorously scrutinize a compendium of natural compounds retrieved from the ZINC database through pharmacodynamic experiments, employing a 1 H-indazole-3-carboxamide (INDZ) scaffold, to identify compounds capable of inhibiting the GSK-3β protein. Utilizing a multi-step approach, the study involved pharmacophore analysis, followed by molecular docking to select five promising ligands for further investigation. Subsequently, ESMACS simulations were employed to assess the stability of the ligand-protein interactions. Evaluation of the binding modes and free energy of the ligands revealed that five compounds (2a-6a) exhibited crucial interactions with the active site residues. Furthermore, various methodologies, including hydrogen bond and clustering analyses, were utilized to ascertain their inhibitory potential and elucidate the factors contributing to ligand binding in the protein's active site. The findings from MMPBSA/GBSA analysis indicated that these five selected small molecules closely approached the IC50 value of the reference ligand (OH8), yielding energy values of -34.85, -32.58, -31.71, and -30.39 kcal/mol, respectively. Additionally, an assessment of the interactions using hydrogen bond and dynamic analyses delineated the effective binding of the ligands with the binding pockets in the protein. Through computational analysis, we obtained valuable insights into the molecular mechanisms of GSK-3β, aiding in the development of more potent inhibitors.

Current State of Open Source Force Fields in Protein-Ligand Binding Affinity Predictions

In drug discovery, the in silico prediction of binding affinity is one of the major means to prioritize compounds for synthesis. Alchemical relative binding free energy (RBFE) calculations based on molecular dynamics (MD) simulations are nowadays a popular approach for the accurate affinity ranking of compounds. MD simulations rely on empirical force field parameters, which strongly influence the accuracy of the predicted affinities. Here, we evaluate the ability of six different small-molecule force fields to predict experimental protein-ligand binding affinities in RBFE calculations on a set of 598 ligands and 22 protein targets. The public force fields OpenFF Parsley and Sage, GAFF, and CGenFF show comparable accuracy, while OPLS3e is significantly more accurate. However, a consensus approach using Sage, GAFF, and CGenFF leads to accuracy comparable to OPLS3e. While Parsley and Sage are performing comparably based on aggregated statistics across the whole dataset, there are differences in terms of outliers. Analysis of the force field reveals that improved parameters lead to significant improvement in the accuracy of affinity predictions on subsets of the dataset involving those parameters. Lower accuracy can not only be attributed to the force field parameters but is also dependent on input preparation and sampling convergence of the calculations. Especially large perturbations and nonconverged simulations lead to less accurate predictions. The input structures, Gromacs force field files, as well as the analysis Python notebooks are available on GitHub.

6-Bromo quinazoline derivatives as cytotoxic agents: design, synthesis, molecular docking and MD simulation

Based on unselectively, several side effects and drug resistance of available anticancer agents, the development and research for novel anticancer agents is necessary. In this study, a new series of quinazoline-4(3H)-one derivatives having a thiol group at position 2 of the quinazoline ring (8a-8 h) were designed and synthesized as potential anticancer agents. The Chemical structures of all compounds were characterized by 1H-NMR, 13C-NMR, and Mass spectroscopy. The antiproliferative activity of all derivatives were determined against two cancer cell lines (MCF-7 and SW480) and one normal cell lines (MRC-5) by the MTT method. Cisplatin, Erlotinib and Doxorubicin were used as positive controls. The results of in vitro screening showed that 8a with an aliphatic linker to SH group was the most potent compound with IC50 values of 15.85 ± 3.32 and 17.85 ± 0.92 µM against MCF-7 and SW480 cell lines, respectively. 8a indicated significantly better potency compared to Erlotinib in the MCF-7 cell line. The cytotoxic results obtained from testing compound 8a on the normal cell line, revealing an IC50 value of 84.20 ± 1.72 µM, provide compelling evidence of its selectivity in distinguishing between tumorigenic and non-tumorigenic cell lines. Structure-activity relationship indicated that the variation in the anticancer activities of quinazoline-4(3H)-one derivatives was affected by different substitutions on the SH position. Molecular docking and MD simulation were carried out for consideration of the binding affinity of compounds against EGFR and EGFR-mutated. The binding energy of compounds 8a and 8c were calculated at -6.7 and - 5.3 kcal.mol- 1, respectively. Compounds 8a and 8c were found to establish hydrogen bonds and some other important interactions with key residue. The DFT analysis was also performed at the B3LYP/6-31 + G(d, p) level for compounds 8a, 8c and Erlotinib. Compound 8a was thermodynamically more stable than 8c. Also, the calculated theoretical and experimental data for the IR spectrum were in agreement. The obtained results delineated that the 8a can be considered an appropriate pharmacophore to develop as an anti-proliferative agent.

Biochemical, kinetic, and structural characterization of a Bacillus tequilensis nitroreductase

Nitroreductases (NRs) are NAD(P)H-dependent flavoenzymes that reduce nitro aromatic compounds to their corresponding arylamines via the nitroso and hydroxylamine intermediates. Because of their broad substrate scope and versatility, NRs have found application in multiple fields such as biocatalysis, bioremediation, cell-imaging and prodrug activation. However, only a limited number of members of the broad NR superfamily (> 24 000 sequences) have been experimentally characterized. Within this group of enzymes, only few are capable of amine synthesis, which is a fundamental chemical transformation for the pharmaceutical, agricultural, and textile industries. Herein, we provide a comprehensive description of a recently discovered NR from Bacillus tequilensis, named BtNR. This enzyme has previously been demonstrated to have the capability to fully convert nitro aromatic and heterocyclic compounds to their respective primary amines. In this study, we determined its biochemical, kinetic and structural properties, including its apparent melting temperature (Tm) of 59 °C, broad pH activity range (from pH 3 to 10) and a notably low redox potential (-236 ± 1 mV) in comparison to other well-known NRs. We also determined its steady-state and pre-steady-state kinetic parameters, which are consistent with other NRs. Additionally, we elucidated the crystal structure of BtNR, which resembles the well-characterized Escherichia coli oxygen-insensitive NAD(P)H nitroreductase (NfsB), and investigated the substrate binding in its active site through docking and molecular dynamics studies with four nitro aromatic substrates. Guided by these structural analyses, we probed the functional roles of active site residues by site-directed mutagenesis. Our findings provide valuable insights into the biochemical and structural properties of BtNR, as well as its potential applications in biotechnology.

A computational study to assess the pathogenicity of single or combinations of missense variants on respiratory complex I

Variants found in the respiratory complex I (CI) subunit genes encoded by mitochondrial DNA can cause severe genetic diseases. However, it is difficult to establish a priori whether a single or a combination of CI variants may impact oxidative phosphorylation. Here we propose a computational approach based on coarse-grained molecular dynamics simulations aimed at investigating new CI variants. One of the primary CI variants associated with the Leber hereditary optic neuropathy (m.14484T>C/MT-ND6) was used as a test case and was investigated alone or in combination with two additional rare CI variants whose role remains uncertain. We found that the primary variant positioned in the E-channel region, which is fundamental for CI function, stiffens the enzyme dynamics. Moreover, a new mechanism for the transition between π- and α-conformation in the helix carrying the primary variant is proposed. This may have implications for the E-channel opening/closing mechanism. Finally, our findings show that one of the rare variants, located next to the primary one, further worsens the stiffening, while the other rare variant does not affect CI function. This approach may be extended to other variants candidate to exert a pathogenic impact on CI dynamics, or to investigate the interaction of multiple variants.

Identifying Key Binding Interactions Between the Cardiac L-Type Calcium Channel and Calmodulin Using Molecular Dynamics Simulations

Defects in the binding of the calcium sensing protein calmodulin (CaM) to the L-type calcium channel (CaV1.2) or to the ryanodine receptor type 2 (RyR2) can lead to dangerous cardiac arrhythmias with distinct phenotypes, such as long-QT syndrome (LQTS) and catecholaminergic ventricular tachycardia (CPVT). Certain CaM mutations lead to LQTS while other mutations lead to CPVT, but the mechanisms by which a specific mutation can lead to each disease phenotype are not well-understood. In this study, we use long, 2 μs molecular dynamics simulations and a multitrajectory approach to identify the key binding interactions between the IQ domain of CaV1.2 and CaM. Five key interactions are found between CaV1.2 and CaM in the C-lobe, 1 in the central linker, and 2 in the N-lobe. In addition, while 5 key interactions appear between residues 120-149 in the C-lobe of CaM when it interacts with CaV1.2, only 1 key interaction is found within this region of CaM when it interacts with the RyR2. We show that this difference in the distribution of key interactions correlates with the known distribution of CaM mutations that lead to LQTS or CPVT. This correlation suggests that a disruption of key binding interactions is a plausible mechanism that can lead to these two different disease phenotypes.

Integrated Anchored Stapling and Hierarchical Dynamics: MSICDA-Driven CREBBP Bromodomain Inhibition

Despite recent success in the computational approaches of cyclic peptide design, current studies face challenges in modeling noncanonical amino acids and nonstandard cyclizations due to limited data. To address this challenge, we developed an integrated framework for the tailored design of stapled peptides (SPs) targeting the bromodomain of CREBBP (CREBBP-BrD). We introduce a powerful combination of anchored stapling and hierarchical molecular dynamics to design and optimize SPs by employing the MultiScale integrative conformational dynamics assessment (MSICDA) strategy, which involves an initial virtual screening of over 1.5 million SPs, followed by comprehensive simulations amounting to 154.54 μs across 5418 of instances. The MSICDA method provides a detailed and holistic stability view of peptide-protein interactions, systematically isolated optimized peptides and identified two leading candidates, DA#430 and DA#99409, characterized by their enhanced stability, optimized binding, and high affinity toward the CREBBP-BrD. In cell-free assays, DA#430 and DA#99409 exhibited 2- to 12-fold greater potency than inhibitor SGC-CBP30. Cell studies revealed higher peptide selectivity for cancerous versus normal cells over small molecules. DA#430 combined with (+)-JQ-1 showed promising synergistic effects. Our approach enables the identification of peptides with optimized binding, high affinity, and enhanced stability, leading to more precise and effective cyclic peptide design, thereby establishing MSICDA as a generalizable and transformative tool for uncovering novel targeted drug development in various therapeutic areas.

Modeling allosteric mechanisms of eukaryotic type II topoisomerases

Type II topoisomerases (TopoIIs) are ubiquitous enzymes that are involved in crucial nuclear processes such as genome organization, chromosome segregation, and other DNA metabolic processes. These enzymes function as large, homodimeric complexes that undergo a complex cycle of binding and hydrolysis of two ATP molecules in their ATPase domains, which regulates the capture and passage of one DNA double-helix through a second, cleaved DNA molecule. This process requires the transmission of information about the state of the bound nucleotide over vast ranges in the TopoII complex. How this information is transmitted at the molecular level to regulate TopoII functions and how protein substitutions disrupt these mechanisms remains largely unknown. Here, we employed extensive microsecond-scale molecular dynamics simulations of the yeast TopoII enzyme in multiple nucleotide-bound states and with amino acid substitutions near both the N and C termini of the complex. Simulation results indicate that the ATPase domains are remarkably flexible on the sub-microsecond timescale and that these dynamics are modulated by the identity of the bound nucleotides and both local and distant amino acid substitutions. Network analyses point toward specific allosteric networks that transmit information about the hydrolysis cycle throughout the complex, which include residues in both the protein and the bound DNA molecule. Amino acid substitutions weaken many of these pathways. Together, our results provide molecular level details on how the TopoII catalytic cycle is controlled through nucleotide binding and hydrolysis and how mutations may disrupt this process.

Directed Evolution of Protoglobin Optimizes the Enzyme Electric Field

To unravel why computational design fails in creating viable enzymes, while directed evolution (DE) succeeds, our research delves into the laboratory evolution of protoglobin. DE has adapted this protein to efficiently catalyze carbene transfer reactions. We show that the previously proposed enhanced substrate access and binding alone cannot account for increased yields during DE. The 3D electric field in the entire active site is tracked through protein dynamics, clustered using the affinity propagation algorithm, and subjected to principal component analysis. This analysis reveals notable changes in the electric field with DE, where distinct field topologies influence transition state energetics and mechanism. A chemically meaningful field component emerges and takes the lead during DE and facilitates crossing the barrier to carbene transfer. Our findings underscore intrinsic electric field dynamic's influence on enzyme function, the ability of the field to switch mechanisms within the same protein, and the crucial role of the field in enzyme design.

Illuminating ligand-induced dynamics in nuclear receptors through MD simulations

Nuclear receptors (NRs) regulate gene expression in critical physiological processes, with their functionality finely tuned by ligand-induced conformational changes. While NRs may sometimes undergo significant conformational motions in response to ligand-binding, these effects are more commonly subtle and challenging to study by traditional structural or biophysical methods. Molecular dynamics (MD) simulations are a powerful tool to bridge the gap between static protein-ligand structures and dynamical changes that govern NR function. Here, we summarize a handful of recent studies that apply MD simulations to study NRs. We present diverse methodologies for analyzing simulation data with a detailed examination of the information each method can yield. By delving into the strengths, limitations and unique contributions of these tools, this review provides guidance for extracting meaningful data from MD simulations to advance the goal of understanding the intricate mechanisms by which ligands orchestrate a range of functional outcomes in NRs.

Multistate Method to Efficiently Account for Tautomerism and Protonation in Alchemical Free-Energy Calculations

The majority of drug-like molecules contain at least one ionizable group, and many common drug scaffolds are subject to tautomeric equilibria. Thus, these compounds are found in a mixture of protonation and/or tautomeric states at physiological pH. Intrinsically, standard classical molecular dynamics (MD) simulations cannot describe such equilibria between states, which negatively impacts the prediction of key molecular properties in silico. Following the formalism described by de Oliveira and co-workers (J. Chem. Theory Comput. 2019, 15, 424-435) to consider the influence of all states on the binding process based on alchemical free-energy calculations, we demonstrate in this work that the multistate method replica-exchange enveloping distribution sampling (RE-EDS) is well suited to describe molecules with multiple protonation and/or tautomeric states in a single simulation. We apply our methodology to a series of eight inhibitors of factor Xa with two protonation states and a series of eight inhibitors of glycogen synthase kinase 3β (GSK3β) with two tautomeric states. In particular, we show that given a sufficient phase-space overlap between the states, RE-EDS is computationally more efficient than standard pairwise free-energy methods.

Optimal Molecular Dynamics System Size for Increased Precision and Efficiency for Epoxy Materials

Molecular dynamics (MD) simulation is an important tool for predicting thermo-mechanical properties of polymer resins at the nanometer length scale, which is particularly important for efficient computationally driven design of advanced composite materials and structures. Because of the statistical nature of modeling amorphous materials on the nanometer length scale, multiple MD models (replicates) are typically built and simulated for statistical sampling of predicted properties. Larger replicates generally provide higher precision in the predictions but result in higher simulation times. Unfortunately, there is insufficient information in the literature to establish guidelines between MD model size and the resulting precision in predicted thermo-mechanical properties. The objective of this study was to determine the optimal MD model size of epoxy resin to balance efficiency and precision. The results show that an MD model size of 15,000 atoms provides for the fastest simulations without sacrificing precision in the prediction of mass density, elastic properties, strength, and thermal properties of epoxy. The results of this study are important for efficient computational process modeling and integrated computational materials engineering (ICME) for the design of next-generation composite materials for demanding applications.

Machine Learning Integrating Protein Structure, Sequence, and Dynamics to Predict the Enzyme Activity of Bovine Enterokinase Variants

Despite recent advances in computational protein science, the dynamic behavior of proteins, which directly governs their biological activity, cannot be gleaned from sequence information alone. To overcome this challenge, we propose a framework that integrates the peptide sequence, protein structure, and protein dynamics descriptors into machine learning algorithms to enhance their predictive capabilities and achieve improved prediction of the protein variant function. The resulting machine learning pipeline integrates traditional sequence and structure information with molecular dynamics simulation data to predict the effects of multiple point mutations on the fold improvement of the activity of bovine enterokinase variants. This study highlights how the combination of structural and dynamic data can provide predictive insights into protein functionality and address protein engineering challenges in industrial contexts.

Synthesis, NMR kinetics and dynamic structure of a 17-mer heptaloop RNA hairpin carrying a 3-N-methyluridine nucleotide residue in the loop region

A 17-mer RNA hairpin (5'GGGAGUXAGCGGCUCCC3') carrying 3-N-methyluridine (m3U) at position X (m3U7-RNA), designed to represent the anticodon stem-loop (ACSL) region of tRNAs to study an open loop state (O-state), was synthesized, purified by HPLC, and characterized by MALDI-ToF_MS and NMR methods. 1H-NMR data revealed primary (P-state in 56.1%), secondary (S-state in 43.9%) and tertiary (∼5-6%) ACSL conformations. Exchange rate constant (kex) for interconversion between P and S states is 112 sec-1 (<Δω ∼454 rad/sec), confirming a slow exchange regime between the two states. Forward (kPS) and backward (kSP) rate constants are 49.166 sec-1 and 62.792 sec-1, respectively, leading to a longer life-time (20.339 msec) for the P-state and a shorter life-time (15.926 msec) for the S-state. In accordance with conformational populations determined by 1H-NMR, dynamics of the P/S/tertiary states of m3U7-RNA and its wild-type counterpart (wt-RNA) were studied by three independent MD production simulations. Cluster analysis revealed that wt-RNA reflects the structural characteristics of the ACSL region of tRNAs. The P-state of m3U7-RNA was found to be structurally similar to wt-RNA but lacks an intraloop H-bond between m3U7 and C10 (U33 and nt36 in tRNAs). In the S-state of m3U7-RNA, m3U7 flips out of the loop region. O-state loop conformations of m3U7-RNA were also clustered (4.8%), wherein the loop nucleotides m3U7.A8.G9.C10.G11 stack one after another. We propose that the O-state of m3U7-RNA is the most suitable conformation that makes the loop accessible for complementary nucleotides and for non-enzymatic primordial replication of small circular RNAs.Communicated by Ramaswamy H. Sarma.

Comprehensive Comparison and Critical Assessment of RNA-Specific Force Fields

Molecular dynamics simulations play a pivotal role in elucidating the dynamic behaviors of RNA structures, offering a valuable complement to traditional methods such as nuclear magnetic resonance or X-ray. Despite this, the current precision of RNA force fields lags behind that of protein force fields. In this work, we systematically compared the performance of four RNA force fields (ff99bsc0χOL3, AMBERDES, ff99OL3_CMAP1, AMBERMaxEnt) across diverse RNA structures. Our findings highlight significant challenges in maintaining stability, particularly with regard to cross-strand and cross-loop hydrogen bonds. Furthermore, we observed the limitations in accurately describing the conformations of nonhelical structural motif, terminal nucleotides, and also base pairing and base stacking interactions by the tested RNA force fields. The identified deficiencies in existing RNA force fields provide valuable insights for subsequent force field development. Concurrently, these findings offer recommendations for selecting appropriate force fields in RNA simulations.

Unstructured linker regions play a role in the differential splicing activities of paralogous RNA binding proteins PTBP1 and PTBP2

RNA Binding Proteins regulate, in part, alternative pre-mRNA splicing and, in turn, gene expression patterns. Polypyrimidine tract binding proteins PTBP1 and PTBP2 are paralogous RNA binding proteins sharing 74% amino acid sequence identity. Both proteins contain four structured RNA-recognition motifs (RRMs) connected by linker regions and an N-terminal region. Despite their similarities, the paralogs have distinct tissue-specific expression patterns and can regulate discrete sets of target exons. How two highly structurally similar proteins can exert different splicing outcomes is not well understood. Previous studies revealed that PTBP2 is post-translationally phosphorylated in the unstructured N-terminal, Linker 1, and Linker 2 regions that share less sequence identity with PTBP1 signifying a role for these regions in dictating the paralog's distinct splicing activities. To this end, we conducted bioinformatics analysis to determine the evolutionary conservation of RRMs versus linker regions in PTBP1 and PTBP2 across species. To determine the role of PTBP2 unstructured regions in splicing activity, we created hybrid PTBP1-PTBP2 constructs that had counterpart PTBP1 regions swapped to an otherwise PTBP2 protein and assayed on differentially regulated exons. We also conducted molecular dynamics studies to investigate how negative charges introduced by phosphorylation in PTBP2 unstructured regions can alter their physical properties. Collectively, results from our studies reveal an important role for PTBP2 unstructured regions and suggest a role for phosphorylation in the differential splicing activities of the paralogs on certain regulated exons.

An Efficient Approach to the Accurate Prediction of Mutational Effects in Antigen Binding to the MHC1

The major histocompatibility complex (MHC) can recognize and bind to external peptides to generate effective immune responses by presenting the peptides to T cells. Therefore, understanding the binding modes of peptide-MHC complexes (pMHC) and predicting the binding affinity of pMHCs play a crucial role in the rational design of peptide vaccines. In this study, we employed molecular dynamics (MD) simulations and free energy calculations with an Alanine Scanning with Generalized Born and Interaction Entropy (ASGBIE) method to investigate the protein-peptide interaction between HLA-A*02:01 and the G9209 peptide derived from the melanoma antigen gp100. The energy contribution of individual residue was calculated using alanine scanning, and hotspots on both the MHC and the peptides were identified. Our study shows that the pMHC binding is dominated by the van der Waals interactions. Furthermore, we optimized the ASGBIE method, achieving a Pearson correlation coefficient of 0.91 between predicted and experimental binding affinity for mutated antigens. This represents a significant improvement over the conventional MM/GBSA method, which yields a Pearson correlation coefficient of 0.22. The computational protocol developed in this study can be applied to the computational screening of antigens for the MHC1 as well as other protein-peptide binding systems.

Inhibitory Activity of Flavonoid Scaffolds on SARS-CoV-2 3CLpro: Insights from the Computational and Experimental Investigations

The emergence of the COVID-19 situation has become a global issue due to the lack of effective antiviral drugs for treatment. Flavonoids are a class of plant secondary metabolites that have antiviral activity against SARS-CoV-2 through inhibition of the main protease (3CLpro). In this study, 22 flavonoids obtained from natural sources and semisynthetic approaches were investigated for their inhibitory activity against SARS-CoV-2 3CLpro, along with cytotoxicity on Vero cells. The protein-ligand interactions were examined using molecular dynamics simulation. Moreover, QSAR analysis was conducted to clarify the structural effects on bioactivity. Accordingly, the in vitro investigation demonstrated that four flavonoids, namely, tectochrysin (7), 6″,6″-dimethylchromeno-[2″,3″:7,8]-flavone (9), panduratin A (19), and genistein (20), showed higher protease inhibitory activity compared to the standard flavonoid baicalein. Finally, our finding suggests that genistein (20), an isoflavone discovered in Millettia brandisiana, has potential for further development as a SARS-CoV-2 3CLpro inhibitor.

Development of a Fluctuating Charge Model for Zinc-Containing Metalloproteins

Metalloproteins widely exist in biology and play important roles in various processes. To accurately simulate metalloprotein systems, modeling polarization and charge transfer effects is vital. The fluctuating charge (FQ) model can efficiently generate atomic charges and simulate the charge transfer effect; it has been developed for a wide range of applications, but few models have been specifically tailored for metalloproteins. In this study, we present a fluctuating charge model specifically for zinc-containing metalloproteins based on the extended charge equilibration (EQeq) scheme. Our model was parametrized to reproduce CM5 charges instead of RESP/CHELPG charges because the former is less dependent on the conformation or basis set, does not suffer from unphysical charges for buried atoms, and is still being able to well reproduce the molecular dipoles. During our study, we found that adding the Pauling-bond-order-like term (referred to as the "+C term" in a previous study) between the zinc ion and ligating atoms significantly improves the model's performance. Although our model was trained for four-coordinated zinc sites only, our results indicated it can well describe the atomic charges in diverse zinc sites. Morever, our model was used to generate partial charges for the metal sites in three different zinc-containing metalloproteins (with four-, five-, and six-coordinated metal sites, respectively). These charges exhibited performance comparable to that of the RESP charges in molecular dynamics (MD) simulations. Additional tests indicated our model could also well reproduce the CM5 charges when geometric changes were involved. Those results indicate that our model can efficiently calculate the atomic charges for metal sites and well simulate the charge transfer effect, which marks an important step toward developing versatile polarizable force fields for metalloproteins.

Assessment of Kaistella jeonii esterase conformational dynamics in response to poly(ethylene terephthalate) binding

The pervasive presence of plastic in the environment has reached a concerning scale, being identified in many ecosystems. Bioremediation is the cheapest and most eco-friendly alternative to remove this polymer from affected areas. Recent work described that a novel cold-active esterase enzyme extracted from the bacteria Kaistella jeonii could promiscuously degrade PET. Compared to the well-known PETase from Ideonella sakaiensis, this novel esterase presents a low sequence identity yet has a remarkably similar folding. However, enzymatic assays demonstrated a lower catalytic efficiency. In this work, we employed a strict computational approach to investigate the binding mechanism between the esterase and PET. Understanding the underlying mechanism of binding can shed light on the evolutive mechanism of how enzymes have been evolving to degrade these artificial molecules and help develop rational engineering approaches to improve PETase-like enzymes. Our results indicate that this esterase misses a disulfide bridge, keeping the catalytic residues closer and possibly influencing its catalytic efficiency. Moreover, we describe the structural response to the interaction between enzyme and PET, indicating local and global effects. Our results aid in deepening the knowledge behind the mechanism of biological catalysis of PET degradation and as a base for the engineering of novel PETases.

Conformational dynamics of α-1 acid glycoprotein (AGP) in cancer: A comparative study of glycosylated and unglycosylated AGP

α-1 acid glycoprotein (AGP) is one of the most abundant plasma proteins. It fulfills two important functions: immunomodulation, and binding to various drugs and receptors. These different functions are closely associated and modulated via changes in glycosylation and cancer missense mutations. From a structural point of view, glycans alter the local biophysical properties of the protein leading to a diverse ligand-binding spectrum. However, glycans can typically not be observed in the resolved X-ray crystallography structure of AGP due to their high flexibility and microheterogeneity, so limiting our understanding of AGP's conformational dynamics 70 years after its discovery. We here investigate how mutations and glycosylation interfere with AGP's conformational dynamics changing its biophysical behavior, by using molecular dynamics (MD) simulations and sequence-based dynamics predictions. The MD trajectories show that glycosylation decreases the local backbone flexibility of AGP and increases the flexibility of distant regions through allosteric effects. We observe that mutations near the glycosylation site affect glycan's conformational preferences. Thus, we conclude that mutations control glycan dynamics which modulates the protein's backbone flexibility directly affecting its accessibility. These findings may assist in the drug design targeting AGP's glycosylation and mutations in cancer.

Aminopyrimidine Derivatives as Multiflavivirus Antiviral Compounds Identified from a Consensus Virtual Screening Approach

Around three billion people are at risk of infection by the dengue virus (DENV) and potentially other flaviviruses. Worldwide outbreaks of DENV, Zika virus (ZIKV), and yellow fever virus (YFV), the lack of antiviral drugs, and limitations on vaccine usage emphasize the need for novel antiviral research. Here, we propose a consensus virtual screening approach to discover potential protease inhibitors (NS3pro) against different flavivirus. We employed an in silico combination of a hologram quantitative structure-activity relationship (HQSAR) model and molecular docking on characterized binding sites followed by molecular dynamics (MD) simulations, which filtered a data set of 7.6 million compounds to 2,775 hits. Lastly, docking and MD simulations selected six final potential NS3pro inhibitors with stable interactions along the simulations. Five compounds had their antiviral activity confirmed against ZIKV, YFV, DENV-2, and DENV-3 (ranging from 4.21 ± 0.14 to 37.51 ± 0.8 μM), displaying aggregator characteristics for enzymatic inhibition against ZIKV NS3pro (ranging from 28 ± 7 to 70 ± 7 μM). Taken together, the compounds identified in this approach may contribute to the design of promising candidates to treat different flavivirus infections.

Insights into the selectivity of a brain-penetrant CDK4/6 vs CDK1/2 inhibitor for glioblastoma used in multiple replica molecular dynamics simulations

Cyclin dependent kinases (CDKs) play an important role in cell cycle regulation and their dysfunction is associated with many cancers. That is why CDKs have been attractive targets for the treatment of cancer. Glioblastoma is a cancer caused by the aberrant expression of CDK4/6, so exploring the mechanism of the selection of CDK4/6 toward inhibitors relative to the other family members CDK1/2 is essential. In this work, multiple replica molecular dynamics (MRMD) simulations, principal component analysis (PCA), free energy landscapes (FELs), molecular mechanics Poisson-Boltzmann/Generalized Born surface area (MM-PB/GBSA) and other methods were integrated to decipher the selectively binding mechanism of the inhibitor N1J to CDK4/6 and CDK1/2. Molecular electrostatic potential (MESP) analysis provides an explanation for the N1J selectivity. Residue-based free energy decomposition reveals that most of the hot residues are located at the same location of CDKs proteins, but the different types of residues in different proteins cause changes in binding energy, which is considered as a potential developmental direction to improve the selectivity of inhibitors to CDK4/6. These results provide insights into the source of inhibitor and CDK4/6 selectivity for the future development of more selective inhibitors.Communicated by Ramaswamy H. Sarma.

Identification of active compounds as novel dipeptidyl peptidase-4 inhibitors through machine learning and structure-based molecular docking simulations

The enzyme dipeptidyl peptidase 4 (DPP4) is a potential therapeutic target for type 2 diabetes (T2DM). Many synthetic anti-DPP4 medications are available to treat T2DM. The need for secure and efficient medicines has been unmet due to the adverse side effects of existing DPP4 medications. The present study implemented a combined approach to machine learning and structure-based virtual screening to identify DPP4 inhibitors. Two ML models were trained based on DPP4 IC50 datasets. The ML models random forest (RF) and multilayer perceptron (MLP) neural network showed good accuracy, with the area under the curve being 0.93 and 0.91, respectively. The natural compound library was screened through ML models, and 1% (217) of compounds were selected for further screening. Structure-based virtual screening was performed along with positive control sitagliptin to obtain more specific and selective leads for DPP4. Based on binding affinity, drug-likeness properties, and interaction with DPP4, Z-614 and Z-997 compounds showed high binding affinity and specificity in the catalytic pocket of DPP4. Finally, the stability conformation of the DPP4 enzyme complex was checked by a molecular dynamics (MD) simulation. The MD simulation showed that both compounds bind better in the catalytic pocket, but the Z-614 compound altered the DPP4 native conformation. Therefore, Z-614 showed a high deviation in the backbone. This combined approach (ML and structure-based) study reported that Z-997 binds most stably to DPP4 in their catalytic pocket with a binding free energy of -70.3 kJ/mol, suggesting its therapeutic potential as a treatment option for T2DM disease.Communicated by Ramaswamy H. Sarma.

When MELD Meets GaMD: Accelerating Biomolecular Landscape Exploration

We introduce Gaussian accelerated MELD (GaMELD) as a new method for exploring the energy landscape of biomolecules. GaMELD combines the strengths of Gaussian accelerated molecular dynamics (GaMD) and modeling employing limited data (MELD) to navigate complex energy landscapes. MELD uses replica-exchange molecular simulations to integrate limited and uncertain data into simulations via Bayesian inference. MELD has been successfully applied to problems of structure prediction like protein folding and complex structure prediction. However, the computational cost for MELD simulations has limited its broader applicability. The synergy of GaMD and MELD surmounts this limitation efficiently sampling the energy landscape at a lower computational cost (reducing the computational cost by a factor of 2 to six). Effectively, GaMD is used to shift energy distributions along replicas to increase the overlap in energy distributions across replicas, facilitating a random walk in replica space. We tested GaMELD on a benchmark set of 12 small proteins that have been previously studied through MELD and conventional MD. GaMELD consistently achieves accurate predictions with fewer replicas. By increasing the efficacy of replica exchange, GaMELD effectively accelerates convergence in the conformational space, enabling improved sampling across a diverse set of systems.

Controlling Uncertainty of Empirical First-Passage Times in the Small-Sample Regime

We derive general bounds on the probability that the empirical first-passage time τ[over ¯]_{n}≡∑_{i=1}^{n}τ_{i}/n of a reversible ergodic Markov process inferred from a sample of n independent realizations deviates from the true mean first-passage time by more than any given amount in either direction. We construct nonasymptotic confidence intervals that hold in the elusive small-sample regime and thus fill the gap between asymptotic methods and the Bayesian approach that is known to be sensitive to prior belief and tends to underestimate uncertainty in the small-sample setting. We prove sharp bounds on extreme first-passage times that control uncertainty even in cases where the mean alone does not sufficiently characterize the statistics. Our concentration-of-measure-based results allow for model-free error control and reliable error estimation in kinetic inference, and are thus important for the analysis of experimental and simulation data in the presence of limited sampling.

A Degron Blocking Strategy Towards Improved CRL4CRBN Recruiting PROTAC Selectivity

Small molecules inducing protein degradation are important pharmacological tools to interrogate complex biology and are rapidly translating into clinical agents. However, to fully realise the potential of these molecules, selectivity remains a limiting challenge. Herein, we addressed the issue of selectivity in the design of CRL4CRBN recruiting PROteolysis TArgeting Chimeras (PROTACs). Thalidomide derivatives used to generate CRL4CRBN recruiting PROTACs have well described intrinsic monovalent degradation profiles by inducing the recruitment of neo-substrates, such as GSPT1, Ikaros and Aiolos. We leveraged structural insights from known CRL4CRBN neo-substrates to attenuate and indeed remove this monovalent degradation function in well-known CRL4CRBN molecular glues degraders, namely CC-885 and Pomalidomide. We then applied these design principles on a previously published BRD9 PROTAC (dBRD9-A) and generated an analogue with improved selectivity profile. Finally, we implemented a computational modelling pipeline to show that our degron blocking design does not impact PROTAC-induced ternary complex formation. We believe that the tools and principles presented in this work will be valuable to support the development of targeted protein degradation.

The structure of a Lactobacillus helveticus chlorogenic acid esterase and the dynamics of its insertion domain provide insights into substrate binding

Chlorogenic acid esterases (ChlEs) are a useful class of enzymes that hydrolyze chlorogenic acid (CGA) into caffeic and quinic acids. ChlEs can break down CGA in foods to improve their sensory properties and release caffeic acid in the digestive system to improve the absorption of bioactive compounds. This work presents the structure, molecular dynamics, and biochemical characterization of a ChlE from Lactobacillus helveticus (Lh). Molecular dynamics simulations suggest that substrate access to the active site of LhChlE is modulated by two hairpin loops above the active site. Docking simulations and mutational analysis suggest that two residues within the loops, Gln145 and Lys164 , are important for CGA binding. Lys164 provides a slight substrate preference for CGA, whereas Gln145 is required for efficient turnover. This work is the first to examine the dynamics of a bacterial ChlE and provides insights on substrate binding preference and turnover in this type of enzyme.

Ensemble-Based Approaches Ensure Reliability and Reproducibility

It is increasingly widely recognized that ensemble-based approaches are required to achieve reliability, accuracy, and precision in molecular dynamics calculations. The purpose of the present article is to address a frequently raised question: what is the optimal way to perform ensemble simulation to calculate quantities of interest?

Accelerating Alchemical Free Energy Prediction Using a Multistate Method: Application to Multiple Kinases

Alchemical free-energy methods based on molecular dynamics (MD) simulations have become important tools to identify modifications of small organic molecules that improve their protein binding affinity during lead optimization. The routine application of pairwise free-energy methods to rank potential binders from best to worst is impacted by the combinatorial increase in calculations to perform when the number of molecules to assess grows. To address this fundamental limitation, our group has developed replica-exchange enveloping distribution sampling (RE-EDS), a pathway-independent multistate method, enabling the calculation of alchemical free-energy differences between multiple ligands (N > 2) from a single MD simulation. In this work, we apply the method to a set of four kinases with diverse binding pockets and their corresponding inhibitors (42 in total), chosen to showcase the general applicability of RE-EDS in prospective drug design campaigns. We show that for the targets studied, RE-EDS is able to model up to 13 ligands simultaneously with high sampling efficiency, leading to a substantial decrease in computational cost when compared to pairwise methods.

Comparison of Equilibrium and Nonequilibrium Approaches for Relative Binding Free Energy Predictions

Alchemical relative binding free energy calculations have recently found important applications in drug optimization. A series of congeneric compounds are generated from a preidentified lead compound, and their relative binding affinities to a protein are assessed in order to optimize candidate drugs. While methods based on equilibrium thermodynamics have been extensively studied, an approach based on nonequilibrium methods has recently been reported together with claims of its superiority. However, these claims pay insufficient attention to the basis and reliability of both methods. Here we report a comparative study of the two approaches across a large data set, comprising more than 500 ligand transformations spanning in excess of 300 ligands binding to a set of 14 diverse protein targets. Ensemble methods are essential to quantify the uncertainty in these calculations, not only for the reasons already established in the equilibrium approach but also to ensure that the nonequilibrium calculations reside within their domain of validity. If and only if ensemble methods are applied, we find that the nonequilibrium method can achieve accuracy and precision comparable to those of the equilibrium approach. Compared to the equilibrium method, the nonequilibrium approach can reduce computational costs but introduces higher computational complexity and longer wall clock times. There are, however, cases where the standard length of a nonequilibrium transition is not sufficient, necessitating a complete rerun of the entire set of transitions. This significantly increases the computational cost and proves to be highly inconvenient during large-scale applications. Our findings provide a key set of recommendations that should be adopted for the reliable implementation of nonequilibrium approaches to relative binding free energy calculations in ligand-protein systems.

Molecular docking and dynamics guided approach to identify potential anti-inflammatory molecules as NRF2 activator to protect against drug-induced liver injury (DILI): a computational study

Inflammation and oxidative stress can contribute to the etiology of metabolic and chronic illnesses. The ability to prevent oxidative stress induced diseases such as cancer, cardiovascular disease, Alzheimer's disease, and others has been the subject of global research. Drug-induced liver injury (DILI) pathogenesis can be either due to oxidative stress or inflammatory response elicited by the drug, its metabolite, or herbal supplements. Our present research uses computational studies to identify a molecule with anti-inflammatory properties that can operate as an NRF2 activator. Acquiring and preparing the KEAP1-NRF2 Protein (PDB: 4L7D) with Schrodinger Suite was followed by developing a ligand library (Anti-inflammatory library downloaded from ChemDiv database). Molecular docking studies were performed in HTVS, SP, and XP modes, respectively. Based on the docking score, interaction, ADMET and binding free energy, the top ten compounds were selected and subjected to induced-fit docking (IFD) analysis for further study. The top three molecules were chosen for a molecular dynamics (MD) simulation study. Using the Desmond module of the Schrodinger Suite, the stability of the protein-ligand complex and protein-ligand contact throughout 100ns were evaluated during the MD simulation study. In our study, it was observed that three compounds exhibit exceptional stability and retain the essential interaction throughout the studies, and it is anticipated that these compounds may act as effective NRF2 activators. Further in vitro and in vivo assessments can be conducted to determine its potential to prevent DILI via acting as an NRF2 activator for future drug development.

Deciphering the Catalytic Mechanism of Virginiamycin B Lyase with Multiscale Methods and Molecular Dynamics Simulations

Due to the emergence of antibiotic resistance, the need to explore novel antibiotics and/or novel strategies to counter antibiotic resistance is of utmost importance. In this work, we explored the molecular and mechanistic details of the degradation of a streptogramin B antibiotic by virginiamycin B (Vgb) lyase of Staphylococcus aureus using classical molecular dynamics simulations and multiscale quantum mechanics/molecular mechanics methods. Our results were in line with available experimental kinetic information. Although we were able to identify a stepwise mechanism, in the wild-type enzyme, the intermediate is short-lived, showing a small barrier to decay to the product state. The impact of point mutations on the reaction was also assessed, showing not only the importance of active site residues to the reaction catalyzed by Vgb lyase but also of near positive and negative residues surrounding the active site. Using molecular dynamics simulations, we also predicted the most likely protonation state of the 3-hydroxypicolinic moiety of the antibiotic and the impact of mutants on antibiotic binding. All this information will expand our understanding of linearization reactions of cyclic antibiotics, which are crucial for the development of novel strategies that aim to tackle antibiotic resistance.

Cardiolipin occupancy profiles of YidC paralogs reveal the significance of respective TM2 helix residues in determining paralog-specific phenotypes

YidC belongs to an evolutionarily conserved family of insertases, YidC/Oxa1/Alb3, in bacteria, mitochondria, and chloroplasts, respectively. Unlike Gram-negative bacteria, Gram-positives including Streptococcus mutans harbor two paralogs of YidC. The mechanism for paralog-specific phenotypes of bacterial YidC1 versus YidC2 has been partially attributed to the differences in their cytoplasmic domains. However, we previously identified a W138R gain-of-function mutation in the YidC1 transmembrane helix 2. YidC1W138R mostly phenocopied YidC2, yet the mechanism remained unknown. Primary sequence comparison of streptococcal YidCs led us to identify and mutate the YidC1W138 analog, YidC2S152 to W/A, which resulted in a loss of YidC2- and acquisition of YidC1-like phenotype. The predicted lipid-facing side chains of YidC1W138/YidC2S152 led us to propose a role for membrane phospholipids in specific-residue dependent phenotypes of S. mutans YidC paralogs. Cardiolipin (CL), a prevalent phospholipid in the S. mutans cytoplasmic membrane during acid stress, is encoded by a single gene, cls. We show a concerted mechanism for cardiolipin and YidC2 under acid stress based on similarly increased promoter activities and similar elimination phenotypes. Using coarse grain molecular dynamics simulations with the Martini2.2 Forcefield, YidC1 and YidC2 wild-type and mutant interactions with CL were assessed in silico. We observed substantially increased CL interaction in dimeric versus monomeric proteins, and variable CL occupancy in YidC1 and YidC2 mutant constructs that mimicked characteristics of the other wild-type paralog. Hence, paralog-specific amino acid- CL interactions contribute to YidC1 and YidC2-associated phenotypes that can be exchanged by point mutation at positions 138 or 152, respectively.

The nuclear receptor LRH-1 discriminates between ligands using distinct allosteric signaling circuits

Nuclear receptors (NRs) are transcription factors that regulate essential biological processes in response to cognate ligands. An important part of NR function involves ligand-induced conformational changes that recruit coregulator proteins to the activation function surface (AFS), ~15 Å away from the ligand-binding pocket. Ligands must communicate with the AFS to recruit appropriate coregulators and elicit different transcriptional outcomes, but this communication is poorly understood. These studies illuminate allosteric communication networks underlying activation of liver receptor homolog-1 (LRH-1), a NR that regulates development, metabolism, cancer progression, and intestinal inflammation. Using >100 μs of all-atom molecular dynamics simulations involving 74 LRH-1 complexes, we identify distinct signaling circuits used by active and inactive ligands for AFS communication. Inactive ligands communicate via strong, coordinated motions along paths through the receptor to the AFS. Activating ligands disrupt the "inactive" circuit and induce connectivity with a second allosteric site. Ligand-contacting residues in helix 7 help mediate the switch between circuits, suggesting new avenues for developing LRH-1-targeted therapeutics. We also elucidate aspects of coregulator signaling, showing that localized, destabilizing fluctuations are induced by inappropriate ligand-coregulator pairings. These studies have uncovered novel features of LRH-1 allostery, and the quantitative approach used to analyze many simulations provides a framework to study allosteric signaling in other receptors.

CHAPERONg: A tool for automated GROMACS-based molecular dynamics simulations and trajectory analyses

Molecular dynamics (MD) simulation is a powerful computational tool used in biomolecular studies to investigate the dynamics, energetics, and interactions of a wide range of biological systems at the atomic level. GROMACS is a widely used free and open-source biomolecular MD simulation software recognized for its efficiency, accuracy, and extensive range of simulation options. However, the complexity of setting up, running, and analyzing MD simulations for diverse systems often poses a significant challenge, requiring considerable time, effort, and expertise. Here, we introduce CHAPERONg, a tool that automates the GROMACS MD simulation pipelines for protein and protein-ligand systems. CHAPERONg also integrates seamlessly with GROMACS modules and third-party tools to provide comprehensive analyses of MD simulation trajectories, offering up to 20 post-simulation processing and trajectory analyses. It also streamlines and automates established pipelines for conducting and analyzing biased MD simulations via the steered MD-umbrella sampling workflow. Thus, CHAPERONg makes MD simulations more accessible to beginner GROMACS users whilst empowering experts to focus on data interpretation and other less programmable aspects of MD simulation workflows. CHAPERONg is written in Bash and Python, and the source code is freely available at https://github.com/abeebyekeen/CHAPERONg. Detailed documentation and tutorials are available online at dedicated web pages accessible via https://abeebyekeen.com/chaperong-online.

In Silico Prediction and Biophysical Validation of Novel 14-3-3σ Homodimer Stabilizers

14-3-3σ plays an important role in controlling tumor metabolic reprogramming and cancer cell growth. However, its function is often compromised in many cancers due to its downregulation. Previous studies found that homodimerization of 14-3-3σ is critical for its activity. However, to date, it is not known if stabilization of 14-3-3σ homodimers can improve its activity or prevent its degradation. In our previous work, we have showed that GCP-Lys-OMe is a potential 14-3-3σ homodimer stabilizer. However, its stabilizing effect was not experimentally validated. Therefore, in this study, we have attempted to predict few potential peptides that can stabilize the dimeric form of 14-3-3σ using similar in silico techniques as described previously for GCP-Lys-OMe. Subsequent [1H]-CPMG NMR experiments confirmed the binding of the peptides (peptides 3, 5, 9, and 16) on 14-3-3σ, with peptide 3 showing the strongest binding. Competitive [1H]-CPMG assays further revealed that while peptide 3 does not compete with a 14-3-3σ binding peptide (ExoS) for the protein's amphipathic groove, it was found to improve ExoS binding on 14-3-3σ. When 14-3-3σ was subjected to dynamic light scattering experiments, the 14-3-3σ homodimer was found to undergo dissociation into monomers prior to aggregation. Intriguingly, the presence of peptide 3 increased 14-3-3σ stability against aggregation. Overall, our findings suggest that (1) docking accompanied by MD simulations can be used to identify potential homodimer stabilizing compounds of 14-3-3σ and (2) peptide 3 can slow down 14-3-3σ aggregation (presumably by preventing its dissociation into monomers), as well as improving the binding of 14-3-3σ to ExoS protein.

Putative Pore Structures of Amyloid β 25-35 in Lipid Bilayers

The amyloid β peptide aggregates to form extracellular plaques in the brains of Alzheimer's disease patients. Certain of its fragments have been found to have similar properties to those of the full-length peptide. The best-studied of these is 25-35, which aggregates into fibrils, is toxic to neurons, and forms ion channels in synthetic lipid bilayers. Here, we investigate possible pore-forming structures of oligomers of this peptide in a POPC/POPG membrane. We consider octameric and decameric β-barrels of different topology, strand orientation, and shear, evaluate their stability in an implicit membrane model, and subject the best models to multimicrosecond all-atom molecular dynamics simulations. We find two decameric structures that are kinetically stable in membranes on this time scale: an imperfectly closed antiparallel β-barrel with K28 in the pore lumen and a short parallel β-barrel with K28 toward the membrane interface. Both structures exhibit dehydrated gaps in the pore lumen, which are larger for the antiparallel barrel. Based on these results, the experimental cation selectivity, the dependence of ion channel activity on voltage direction, and certain mutation data, the parallel model seems more compatible with experimental data.

Design, synthesis, computational study and cytotoxic evaluation of some new quinazoline derivatives containing pyrimidine moiety

Quinazoline derivatives, as an important category of heterocyclic compounds, have received much attention for the design and development of new drugs due to their various pharmacological properties. Besides, there is a great deal of evidence showing pyrimidine analogs as anticancer agents. Thus, in the present study, for the design of new target compounds with cytotoxic activity, we focused on various quinazolinone and pyrimidine hybrids. A new series of quinazoline-pyrimidine hybrid derivatives (6a-6n) have been designed and synthesized as novel antiproliferative agents. All the synthesized compounds characterized based on their IR, NMR and Mass spectroscopic data. Antiproliferative activities of the new compounds were evaluated against three human cancer cell lines (MCF-7, A549, SW-480). The compounds were found to have appropriate potential with IC50 values ranging from 2.3 ± 5.91 to 176.5 ± 0.7 μM against the tested cell lines. Compound 6n exerted the highest antiproliferative activity with IC50 values of 5.9 ± 1.69 μM, 2.3 ± 5.91 μM and 5.65 ± 2.33 μM against A549, SW-480 and MCF-7 respectively. The results indicated that 6n could induce apoptosis in A549 cell line in a dose dependent manner and arrest in the S phase of cell cycle. Docking studies were also done to investigate the detailed binding pattern of the synthesized compounds against EGFR. Furthermore, molecular dynamic simulation and binding free energy calculation have been done to rescore initial docking pose of the synthesized compounds using ensemble-based MMGB/PBSA free energy method. According to the results, free energy calculation confirmed biological activity of compounds and also, Arg 817 and Lys 721 residues had the pivotal role in the high potency of 6n. Finally, the drug likeness and in silico ADME study were also predicted.

3DDPDs: describing protein dynamics for proteochemometric bioactivity prediction. A case for (mutant) G protein-coupled receptors

Proteochemometric (PCM) modelling is a powerful computational drug discovery tool used in bioactivity prediction of potential drug candidates relying on both chemical and protein information. In PCM features are computed to describe small molecules and proteins, which directly impact the quality of the predictive models. State-of-the-art protein descriptors, however, are calculated from the protein sequence and neglect the dynamic nature of proteins. This dynamic nature can be computationally simulated with molecular dynamics (MD). Here, novel 3D dynamic protein descriptors (3DDPDs) were designed to be applied in bioactivity prediction tasks with PCM models. As a test case, publicly available G protein-coupled receptor (GPCR) MD data from GPCRmd was used. GPCRs are membrane-bound proteins, which are activated by hormones and neurotransmitters, and constitute an important target family for drug discovery. GPCRs exist in different conformational states that allow the transmission of diverse signals and that can be modified by ligand interactions, among other factors. To translate the MD-encoded protein dynamics two types of 3DDPDs were considered: one-hot encoded residue-specific (rs) and embedding-like protein-specific (ps) 3DDPDs. The descriptors were developed by calculating distributions of trajectory coordinates and partial charges, applying dimensionality reduction, and subsequently condensing them into vectors per residue or protein, respectively. 3DDPDs were benchmarked on several PCM tasks against state-of-the-art non-dynamic protein descriptors. Our rs- and ps3DDPDs outperformed non-dynamic descriptors in regression tasks using a temporal split and showed comparable performance with a random split and in all classification tasks. Combinations of non-dynamic descriptors with 3DDPDs did not result in increased performance. Finally, the power of 3DDPDs to capture dynamic fluctuations in mutant GPCRs was explored. The results presented here show the potential of including protein dynamic information on machine learning tasks, specifically bioactivity prediction, and open opportunities for applications in drug discovery, including oncology.

Pentagalloyl Glucose-Targeted Inhibition of P-Glycoprotein and Re-Sensitization of Multidrug-Resistant Leukemic Cells (K562/ADR) to Doxorubicin: In Silico and Functional Studies

Combining phytochemicals with chemotherapeutic drugs has demonstrated the potential to surmount drug resistance. In this paper, we explore the efficacy of pentagalloyl glucose (PGG) in modulating P-gp and reversing multidrug resistance (MDR) in drug-resistant leukemic cells (K562/ADR). The cytotoxicity of PGG was evaluated using a CCK-8 assay, and cell apoptosis was assessed using flow cytometry. Western blotting was used to analyze protein expression levels. P-glycoprotein (P-gp) activity was evaluated by monitoring the kinetics of P-gp-mediated efflux of pirarubicin (THP). Finally, molecular docking, molecular dynamics simulation, and molecular mechanics with generalized Born and surface area solvation (MM-GBSA) calculation were conducted to investigate drug-protein interactions. We found that PGG selectively induced cytotoxicity in K562/ADR cells and enhanced sensitivity to doxorubicin (DOX), indicating its potential as a reversal agent. PGG reduced the expression of P-gp and its gene transcript levels. Additionally, PGG inhibited P-gp-mediated efflux and increased intracellular drug accumulation in drug-resistant cells. Molecular dynamics simulations and MM-GBSA calculation provided insights into the binding affinity of PGG to P-gp, suggesting that PGG binds tightly to both the substrate and the ATP binding sites of P-gp. These findings support the potential of PGG to target P-gp, reverse drug resistance, and enhance the efficacy of anticancer therapies.

Can Free Energy Perturbation Simulations Coupled with Replica-Exchange Molecular Dynamics Study Ligands with Distributed Binding Sites?

Free energy perturbation coupled with replica exchange with solute tempering (FEP/REST) offers a rigorous approach to compute relative free energy changes for ligands. To determine the applicability of FEP/REST for the ligands with distributed binding poses, we considered two alchemical transformations involving three putative inhibitors I0, I1, and I2 of the Venezuelan equine encephalitis virus nuclear localization signal sequence binding to the importin-α (impα) transporter protein. I0 → I1 and I0 → I2 transformations, respectively, increase or decrease the polarity of the parent molecule. Our objective was three-fold─(i) to verify FEP/REST technical performance and convergence, (ii) to estimate changes in binding free energy ΔΔG, and (iii) to determine the utility of FEP/REST simulations for conformational binding analysis. Our results are as follows. First, our FEP/REST implementation properly follows FEP/REST formalism and produces converged ΔΔG estimates. Due to ligand inherent unbinding, the better FEP/REST strategy lies in performing multiple independent trajectories rather than extending their length. Second, I0 → I1 and I0 → I2 transformations result in overall minor changes in inhibitor binding free energy, slightly strengthening the affinity of I1 and weakening that of I2. Electrostatic interactions dominate binding interactions, determining the enthalpic changes. The two transformations cause opposite entropic changes, which ultimately govern binding affinities. Importantly, we confirm the validity of FEP/REST free energy estimates by comparing them with our previous REST simulations, directly probing binding of three ligands to impα. Third, we established that FEP/REST simulations can sample binding ensembles of ligands. Thus, FEP/REST can be applied (i) to study the energetics of the ligand binding without defined poses and showing minor differences in affinities |ΔΔG| ≲ 0.5 kcal/mol and (ii) to collect ligand binding conformational ensembles.

A dynamic biomimetic model of the membrane-bound CD4-CD3-TCR complex during pMHC disengagement

The coordinated (dis)engagement of the membrane-bound T cell receptor (TCR)-CD3-CD4 complex from the peptide-major histocompatibility complex (pMHC) is fundamental to TCR signal transduction and T cell effector function. As such, an atomic-scale understanding would not only enhance our basic understanding of the adaptive immune response but would also accelerate the rational design of TCRs for immunotherapy. In this study, we explore the impact of the CD4 coreceptor on the TCR-pMHC (dis)engagement by constructing a molecular-level biomimetic model of the CD3-TCR-pMHC and CD4-CD3-TCR-pMHC complexes within a lipid bilayer. After allowing the system complexes to equilibrate (engage), we use steered molecular dynamics to dissociate (disengage) the pMHC. We find that 1) the CD4 confines the pMHC closer to the T cell by 1.8 nm at equilibrium; 2) CD4 confinement shifts the TCR along the MHC binding groove engaging a different set of amino acids and enhancing the TCR-pMHC bond lifetime; 3) the CD4 translocates under load increasing the interaction strength between the CD4-pMHC, CD4-TCR, and CD4-CD3; and 4) upon dissociation, the CD3-TCR complex undergoes structural oscillation and increased energetic fluctuation between the CD3-TCR and CD3-lipids. These atomic-level simulations provide mechanistic insight on how the CD4 coreceptor impacts TCR-pMHC (dis)engagement. More specifically, our results provide further support (enhanced bond lifetime) for a force-dependent kinetic proofreading model and identify an alternate set of amino acids in the TCR that dominate the TCR-pMHC interaction and could thus impact the design of TCRs for immunotherapy.

Analysis of the Distribution and Influencing Factors of Diffusion Coefficient Model Parameters Based on Molecular Dynamics Simulations

The establishment of mathematical models to predict the diffusion coefficients of gas and liquid systems have important theoretical significance and practical value. In this work, based on the previously proposed diffusion coefficient model D_LV, the distribution and influencing factors of the model parameters characteristic length (L) and diffusion velocity (V) were further investigated using molecular dynamics simulations. The statistical analysis of L and V for 10 gas systems and 10 liquid systems was presented in the paper. New distribution functions were established to describe the probability distributions of molecular motion L and V. The mean values of the correlation coefficients were 0.98 and 0.99, respectively. Meanwhile, the effects of molecular molar mass and system temperature on the molecular diffusion coefficients were discussed. The results show that the effect of molecular molar mass on the diffusion coefficient mainly affects the molecular motion L, and the effect of system temperature on the diffusion coefficient mainly affects V. For the gas system, the average relative deviation of D_LV and D_MSD is 10.73% and that of D_LV and experimental value is 12.63%; for the solution system, the average relative deviation of D_LV and D_MSD is 12.93% and that of D_LV and experimental value is 18.86%, which indicates the accuracy of the new model results. The new model reveals the potential mechanism of molecular motion and provides a theoretical basis for further study of the diffusion process.

Binding Mechanism of Inhibitors to Heat Shock Protein 90 Investigated by Multiple Independent Molecular Dynamics Simulations and Prediction of Binding Free Energy

The heat shock protein (HSP90) has been an import target of drug design in the treatment of human disease. An exploration of the conformational changes in HSP90 can provide useful information for the development of efficient inhibitors targeting HSP90. In this work, multiple independent all-atom molecular dynamics (AAMD) simulations followed by calculations of the molecular mechanics generalized Born surface area (MM-GBSA) were performed to explore the binding mechanism of three inhibitors (W8Y, W8V, and W8S) to HSP90. The dynamics analyses verified that the presence of inhibitors impacts the structural flexibility, correlated movements, and dynamics behavior of HSP90. The results of the MM-GBSA calculations suggest that the selection of GB models and empirical parameters has important influences on the predicted results and verify that van der Waals interactions are the main forces that determine inhibitor-HSP90 binding. The contributions of separate residues to the inhibitor-HSP90 binding process indicate that hydrogen-bonding interactions (HBIs) and hydrophobic interactions play important roles in HSP90-inhibitor identifications. Moreover, residues L34, N37, D40, A41, D79, I82, G83, M84, F124, and T171 are recognized as hot spots of inhibitor-HSP90 binding and provide significant target sites of for the design of drugs related to HSP90. This study aims to contribute to the development of efficient inhibitors that target HSP90 by providing an energy-based and theoretical foundation.

Identification Mechanism of BACE1 on Inhibitors Probed by Using Multiple Separate Molecular Dynamics Simulations and Comparative Calculations of Binding Free Energies

β-amyloid cleaving enzyme 1 (BACE1) is regarded as an important target of drug design toward the treatment of Alzheimer's disease (AD). In this study, three separate molecular dynamics (MD) simulations and calculations of binding free energies were carried out to comparatively determine the identification mechanism of BACE1 for three inhibitors, 60W, 954 and 60X. The analyses of MD trajectories indicated that the presence of three inhibitors influences the structural stability, flexibility and internal dynamics of BACE1. Binding free energies calculated by using solvated interaction energy (SIE) and molecular mechanics generalized Born surface area (MM-GBSA) methods reveal that the hydrophobic interactions provide decisive forces for inhibitor-BACE1 binding. The calculations of residue-based free energy decomposition suggest that the sidechains of residues L91, D93, S96, V130, Q134, W137, F169 and I179 play key roles in inhibitor-BACE1 binding, which provides a direction for future drug design toward the treatment of AD.

Multiscale models of integrins and cellular adhesions

Computational models of integrin-based adhesion complexes have revealed important insights into the mechanisms by which cells establish connections with their external environment. However, how changes in conformation and function of individual adhesion proteins regulate the dynamics of whole adhesion complexes remains largely elusive. This is because of the large separation in time and length scales between the dynamics of individual adhesion proteins (nanoseconds and nanometers) and the emergent dynamics of the whole adhesion complex (seconds and micrometers), and the limitations of molecular simulation approaches in extracting accurate free energies, conformational transitions, reaction mechanisms, and kinetic rates, that can inform mechanisms at the larger scales. In this review, we discuss models of integrin-based adhesion complexes and highlight their main findings regarding: (i) the conformational transitions of integrins at the molecular and macromolecular scales and (ii) the molecular clutch mechanism at the mesoscale. Lastly, we present unanswered questions in the field of modeling adhesions and propose new ideas for future exciting modeling opportunities.

All-Atom Simulations Reveal the Intricacies of Signal Transduction upon Binding of the HLA-E Ligand to the Transmembrane Inhibitory CD94/NKG2A Receptor

Natural killer (NK) cells play an important role in the innate immune response against tumors and various pathogens such as viruses and bacteria. Their function is controlled by a wide array of activating and inhibitory receptors, which are expressed on their cell surface. Among them is a dimeric NKG2A/CD94 inhibitory transmembrane (TM) receptor which specifically binds to the non-classical MHC I molecule HLA-E, which is often overexpressed on the surface of senescent and tumor cells. Using the Alphafold 2 artificial intelligence system, we constructed the missing segments of the NKG2A/CD94 receptor and generated its complete 3D structure comprising extracellular (EC), TM, and intracellular regions, which served as a starting point for the multi-microsecond all-atom molecular dynamics simulations of the receptor with and without the bound HLA-E ligand and its nonameric peptide. The simulated models revealed that an intricate interplay of events is taking place between the EC and TM regions ultimately affecting the intracellular immunoreceptor tyrosine-based inhibition motif (ITIM) regions that host the point at which the signal is transmitted further down the inhibitory signaling cascade. Signal transduction through the lipid bilayer was also coupled with the changes in the relative orientation of the NKG2A/CD94 TM helices in response to linker reorganization, mediated by fine-tuned interactions in the EC region of the receptor, taking place after HLA-E binding. This research provides atomistic details of the cells' protection mechanism against NK cells and broadens the knowledge regarding the TM signaling of ITIM-bearing receptors.

Deciphering the intrinsic dynamics of unphosphorylated IRAK4 kinase bound to type I and type II inhibitors

Interleukin-1 receptor-associated kinase 4 (IRAK4) is a vital protein involved in Toll-like and interleukin-1 receptor signal transduction. Several studies have reported regarding the crystal structure, dynamic properties, and interactions with inhibitors of the phosphorylated form of IRAK4. However, no dynamic properties of inhibitor-bound unphosphorylated IRAK4 have been previously studied. Herein, we report the intrinsic dynamics of unphosphorylated IRAK4 (uIRAK4) bound to type I and type II inhibitors. The corresponding apo and inhibitor-bound forms of uIRAK4 were subjected to three independent simulations of 500 ns (total 1.5 μs) each, and their trajectories were analyzed. The results indicated that all three systems were relatively stable, except for the type II inhibitor-bound form of uIRAK4, which exhibited less compact folding and higher solvent surface area. The intra-hydrogen bonds corroborated the structural deformation of the type-II inhibitor-bound complex, which could be attributed to the long molecular structure of the type-II inhibitor. Moreover, the type II inhibitor bound to uIRAK4 showed higher binding free energy with uIRAK4 than the type I inhibitor. The free energy landscape analysis showed a reorientation of Phe330 side chain from the DFG motif at different metastable states for all the systems. The intra-residual distance between residues Lys213, Glu233, Tyr262, and Phe330 suggests a functional interplay when the inhibitors are bound to uIRAK4, thereby hinting at their crucial role in the inhibition mechanism. Ultimately, the intrinsic dynamics study observed between type I/II inhibitor-bound forms of uIRAK4 may assist in better understanding the enzyme and designing therapeutic compounds.