Conventional Molecular Dynamics Research Articles

Unveiling dissociation mechanisms and binding patterns in the UHRF1-DPPA3 complex via multi-replica molecular dynamics simulations.

Ubiquitin-like with PHD and RING finger domain containing protein 1 (UHRF1) is responsible for preserving the stability of genomic methylation through the recruitment of DNA methyltransferase 1 (DNMT1). However, the interaction between Developmental pluripotency associated 3 (DPPA3) and the pre-PHD-PHD (PPHD) domain of UHRF1 hinders the nuclear localization of UHRF1. This disruption has implications for potential cancer treatment strategies. Drugs that mimic the binding pattern between DPPA3 and PPHD could offer a promising approach to cancer treatment. Our study reveals that DPPA3 undergoes dissociation from the C-terminal through three different modes of helix unfolding. Furthermore, we have identified key residue pairs involved in this dissociation process and potential drug-targeting residues. These findings offer valuable insights into the dissociation mechanism of DPPA3 from PPHD and have the potential to inform the design of novel drugs targeting UHRF1 for cancer therapy. To comprehend the dissociation process and binding patterns of PPHD-DPPA3, we employed enhanced sampling techniques, including steered molecular dynamics (SMD) and conventional molecular dynamics (cMD). Additionally, we utilized self-organizing maps (SOM) and time-resolved force distribution analysis (TRFDA) methodologies. The Gromacs software was used for performing molecular dynamics simulations, and the AMBER FF14SB force field was applied to the protein.

Two-Step Growth Mechanism of the Solid Electrolyte Interphase in Argyrodyte/Li-Metal Contacts.

The structure and growth of the solid electrolyte interphase (SEI) region between an electrolyte and an electrode is one of the most fundamental yet less well-understood phenomena in solid-state batteries. We present an atomistic simulation of the SEI growth for one of the currently promising solid electrolytes (Li6PS5Cl), based on ab initio-trained machine learning interatomic potentials, for over 30,000 atoms during 10 ns, well beyond the capabilities of conventional molecular dynamics. This unveils a two-step growth mechanism: a Li-argyrodite chemical reaction leading to the formation of an amorphous phase, followed by a kinetically slower crystallization of the reaction products into a 5Li2S·Li3P·LiCl solid solution. The simulation results support the recent, experimentally founded hypothesis of an indirect pathway of electrolyte reduction. These findings shed light on the intricate processes governing SEI evolution, providing a valuable foundation for the design and optimization of next-generation solid-state batteries.

Development of Membrane Permeability Coefficient by Means of Novel Molecular Dynamics Methods

The membrane permeability, and its evaluation, is crucial factor in the process of uptake of compounds from outside to inside the cell and in the inhibition of the activity of disease-causing target proteins. Although molecular dynamics (MD) simulations have been shown to be able to reproduce the conformational changes of compounds occurring during membrane permeation, it is still challenging to extract the membrane permeability at an affordable computational workload solely by conventional MD. Indeed, the time scale accessible by MD is far below the one characterizing the actual permeation process. Phenomena occurring in living organisms escaping the reach of standard MD are generally referred to as biological rare events, and the membrane permeation process is one of them. To overcome this time-scale problem, several enhanced sampling methods have been proposed over the years to improve conformational sampling. In this review, a hybrid sampling method that combines the parallel cascade selection MD (PaCS-MD) and the outlier flooding method (OFLOOD), introduced and developed by our group, is proposed as a tool to study the membrane permeation from structural sampling (rare-event sampling). The obtained trajectories are used to estimate the free energy profiles for the membrane permeation and to compute the membrane permeation coefficients. Moreover, we present an example of application of the free energy reaction network method as a versatile way for incorporating explicitly into reaction coordinates the degrees of freedom related to internal motion.

The Development of CDC25A-Derived Phosphoseryl Peptides That Bind 14-3-3ε with High Affinities.

Overexpression of the 14-3-3ε protein is associated with suppression of apoptosis in cutaneous squamous cell carcinoma (cSCC). This antiapoptotic activity of 14-3-3ε is dependent on its binding to CDC25A; thus, inhibiting 14-3-3ε - CDC25A interaction is an attractive therapeutic approach to promote apoptosis in cSCC. In this regard, designing peptide inhibitors of 14-3-3ε - CDC25A interactions is of great interest. This work reports the rational design of peptide analogs of pS, a CDC25A-derived peptide that has been shown to inhibit 14-3-3ε-CDC25A interaction and promote apoptosis in cSCC with micromolar IC50. We designed new peptide analogs in silico by shortening the parent pS peptide from 14 to 9 amino acid residues; then, based on binding motifs of 14-3-3 proteins, we introduced modifications in the pS(174-182) peptide. We studied the binding of the peptides using conventional molecular dynamics (MD) and steered MD simulations, as well as biophysical methods. Our results showed that shortening the pS peptide from 14 to 9 amino acids reduced the affinity of the peptide. However, substituting Gln176 with either Phe or Tyr amino acids rescued the binding of the peptide. The optimized peptides obtained in this work can be candidates for inhibition of 14-3-3ε - CDC25A interactions in cSCC.

Uncertainty-biased molecular dynamics for learning uniformly accurate interatomic potentials

Efficiently creating a concise but comprehensive data set for training machine-learned interatomic potentials (MLIPs) is an under-explored problem. Active learning, which uses biased or unbiased molecular dynamics (MD) to generate candidate pools, aims to address this objective. Existing biased and unbiased MD-simulation methods, however, are prone to miss either rare events or extrapolative regions—areas of the configurational space where unreliable predictions are made. This work demonstrates that MD, when biased by the MLIP’s energy uncertainty, simultaneously captures extrapolative regions and rare events, which is crucial for developing uniformly accurate MLIPs. Furthermore, exploiting automatic differentiation, we enhance bias-forces-driven MD with the concept of bias stress. We employ calibrated gradient-based uncertainties to yield MLIPs with similar or, sometimes, better accuracy than ensemble-based methods at a lower computational cost. Finally, we apply uncertainty-biased MD to alanine dipeptide and MIL-53(Al), generating MLIPs that represent both configurational spaces more accurately than models trained with conventional MD.

Effect of quercetin on the protein-substrate interactions in SIRT6: Insight from MD simulations

SIRT6 is of interest for its promising effect in the treatment of aging-related diseases. Studies have shown quercetin (QUE) and its derivatives have varying degrees of effect on the catalytic effect of SIRT6. In the research, the effect of QUE on the protein-substrate interaction in the SIRT6-mediated mono-ADP ribosylation system was investigated by conventional molecular dynamics (MD) simulations combined with MM/PBSA binding free energy calculations. The results show that QUE can bind stably to SIRT6 with the binding energy of −22.8 kcal/mol and further affect the atomic interaction between SIRT6 and NAD+ (or H3K9), resulting in an increased affinity between SIRT6-NAD+ and decreased SIRT6-H3K9 binding capacity. At the same time, the binding of QUE can also alter some structural characteristics of the protein, with large shifts occurring in the residue regions involving the N-terminal (residues 1–27), Rossmann fold regions (residues 55–92), and ZBD (residues 164–179). Thus, QUE shows great potential as a scaffold for the design of novel potent SIRT6 modulators.

Understanding the Modulatory Role of E2012 on the γ-Secretase-Substrate Interaction.

Allosteric modulation plays a critical role in enzyme functionality and requires a deep understanding of the interactions between the active and allosteric sites. γ-Secretase (GS) is a key therapeutic target in the treatment of Alzheimer's disease (AD), through its role in the synthesis of amyloid β peptides that accumulate in AD patients. This study explores the structure and dynamic effects of GS modulation by E2012 binding, employing well-tempered metadynamics and conventional molecular dynamics simulations across three binding scenarios: (1) GS enzyme with and without L458 inhibitor, (2) the GS-substrate complex together with the modulator E2012 in two different binding modes, and (3) E2012 interacting with a C99 substrate fragment. Our findings reveal that the presence of L458 induces conformational changes that contribute to stabilization of the GS enzyme dynamics, previously reported as a key factor that allowed the resolution of the cryo-EM structure and the enhanced binding of E2012. Furthermore, we identified the most favorable binding site for E2012 within the GS-substrate complex, uncovering significant modulatory effects and a complex network of interactions that influence the position of the substrate for catalysis. In addition, we explore a potential substrate-modulator binding before the formation of the enzyme-substrate complex. The insights gained from our study emphasize the importance of these interactions in the development of potential therapeutic interventions that target the functionality of the GS enzyme in AD.

Generation of conformational ensembles of small molecules via surrogate model-assisted molecular dynamics

The accurate prediction of thermodynamic properties is crucial in various fields such as drug discovery and materials design. This task relies on sampling from the underlying Boltzmann distribution, which is challenging using conventional approaches such as simulations. In this work, we introduce surrogate model-assisted molecular dynamics (SMA-MD), a new procedure to sample the equilibrium ensemble of molecules. First, SMA-MD leverages deep generative models to enhance the sampling of slow degrees of freedom. Subsequently, the generated ensemble undergoes statistical reweighting, followed by short simulations. Our empirical results show that SMA-MD generates more diverse and lower energy ensembles than conventional MD simulations. Furthermore, we showcase the application of SMA-MD for the computation of thermodynamical properties by estimating implicit solvation free energies.

Potential acridinedione derivatives for the development of a heterobifunctional PROTAC for targeted degradation of PCSK9 protein

Proprotein convertase substilisin-like/kexin type 9 (PCSK9) interacts with the low-density lipoprotein (LDL) receptor by protein-protein interaction (PPI). The inhibition of this PPI lowers the LDL cholesterol levels, thereby reducing the risks of cardiovascular diseases. We identified in-house synthesized acridinedione scaffolds as binders of PCSK9 through microsecond timescale conventional and biased molecular dynamics (MD) simulations. The free energy of binding was calculated by the Molecular Mechanics Poisson-Boltzmann Surface Area and umbrella sampling methods for acridinedione molecules and compared with reference molecules (co-crystallized inhibitors). Three acridinedione molecules (DSPD-2, DSPD-4, and DSPD-6) were identified as potential binders of the PCSK9 allosteric binding site from preliminary docking and steered MD simulations. The stability of molecules inside the allosteric pocket was established by the low values of root mean square deviations of backbone C-α atoms. The acridinedione derivatives also showed a shorter distance of hydrogen bonds as compared to the standard molecules. The molecule DSPD-6 showed the lowest binding free energy as compared to other selected and standard molecules. The binding of the selected molecules at the allosteric binding pocket had no effect on structural parameters of LDL receptor-binding site of PCSK9 protein. Finally, a proteolysis targeting chimera (PROTAC) was developed by attaching a proteasome recruiting tag to DSPD-6 (lowest binding affinity) to achieve ligand-induced protein degradation. These results qualify the development of the DSPD-6 molecule as a heterobifunctional PROTAC for targeted degradation of PCSK9 protein. Moreover, our computational strategy provides a framework for the identification of molecules against conventionally undruggable PPI targets.

Unraveling the Ligand-Binding Sites of CYP3A4 by Molecular Dynamics Simulations with Solvent Probes.

Cytochrome P450 3A4 (CYP3A4) is one of the most important drug-metabolizing enzymes in the human body and is well known for its complicated, atypical kinetic characteristics. The existence of multiple ligand-binding sites in CYP3A4 has been widely recognized as being capable of interfering with the active pocket through allosteric effects. The identification of ligand-binding sites other than the canonical active site above the heme is especially important for understanding the atypical kinetic characteristics of CYP3A4 and the intriguing association between the ligand and the receptor. In this study, we first employed mixed-solvent molecular dynamics (MixMD) simulations coupled with the online computational predictive tools to explore potential ligand-binding sites in CYP3A4. The MixMD approach demonstrates better performance in dealing with the receptor flexibility compared with other computational tools. From the sites identified by MixMD, we then picked out multiple sites for further exploration using ensemble docking and conventional molecular dynamics (cMD) simulations. Our results indicate that three extra sites are suitable for ligand binding in CYP3A4, including one experimentally confirmed site and two novel sites.

Understanding the Reversible Binding of a Multichain Protein-Protein Complex through Free-Energy Calculations.

We demonstrate that the binding affinity of a multichain protein-protein complex, insulin dimer, can be accurately predicted using a streamlined route of standard binding free-energy calculations. We find that chains A and C, which do not interact directly during binding, stabilize the insulin monomer structures and reduce the binding affinity of the two monomers, therefore enabling their reversible association. Notably, we confirm that although classical methods can estimate the binding affinity of the insulin dimer, conventional molecular dynamics, enhanced sampling algorithms, and classical geometrical routes of binding free-energy calculations may not fully capture certain aspects of the role played by the noninteracting chains in the binding dynamics. Therefore, this study not only elucidates the role of noninteracting chains in the reversible binding of the insulin dimer but also offers a methodological guide for investigating the reversible binding of multichain protein-protein complexes utilizing streamlined free-energy calculations.

Base pair dynamics, electrostatics, and thermodynamics at the LTR-III quadruplex:duplex junction

G-quadruplexes (GQs) play key regulatory roles within the human genome and have also been identified to play similar roles in other eukaryotes, bacteria, archaea, and viruses. Human immunodeficiency virus 1, the etiological agent of acquired immunodeficiency syndrome, can form two GQs in its long terminal repeat (LTR) promoter region, each of which act to regulate viral gene expression in opposing manners. The major LTR GQ, called LTR-III, is a distinct hybrid GQ containing a 12-nucleotide duplex loop attached to the quadruplex motif. The resulting quadruplex:duplex junction (QDJ) has been hypothesized to serve as a selective drug targeting site. To better understand the dynamics of this QDJ, we performed conventional and enhanced-sampling molecular dynamics simulations using the Drude-2017 force field. We observed unbiased and reversible formation of additional base pairs in the QDJ, between Ade4:Thy14 and Gua3:Thy14. Both base pairs were electrostatically favored, but geometric constraints within the junction may drive the formation of, and preference for, the Ade4:Thy14 base pair. Finally, we demonstrated that the base pairs are separated only by small energy barriers that may enable transitions between both base-paired states. Together, these simulations provide new insights into the dynamics, electrostatics, and thermodynamics of the LTR-III QDJ.

Water-Hydrocarbon Interactions in Anionic Pyrene Monohydrate.

Interactions between water and polycyclic aromatic hydrocarbons are essential in many aspects of chemistry, from interstellar and atmospheric processes to interfacial hydrophobicity and wetting phenomena. Despite their growing importance, the intermolecular potentials of the water-hydrocarbon interactions are underdeveloped compared to the water-water potentials, and there are similarly few experimental probes that are sensitive to the details of the water-hydrocarbon potential. We present a combined experimental and computational study of anionic pyrene monohydrate, one of the simplest water/hydrocarbon clusters. The action spectrum in the OH region of the mass-selected cluster ion provides a rigorous benchmark for intermolecular potentials and computational methodologies. We identify missing intermolecular interactions and shortcomings in conventional dynamics calculations by comparing experimental data to density functional theory and classical molecular dynamics calculations. Kinetic trapping is prevalent, even for one water molecule and one pyrene molecule, leading to slow equilibration in conventional molecular dynamics calculations, even on nanosecond time scales and at low temperatures (50 K). At constant energy, temperature fluctuations for the pair of molecules are substantial. Immersing the system in a bath of soft spheres and employing parallel tempering alleviates kinetic trapping and dampens temperature fluctuations, bringing the system closer to the thermodynamic limit. With such augmented sampling, a simple, flexible water model reproduces the line width and the asymmetric broadening of the symmetric OH stretching mode, which we assign to spectral diffusion. In the OH stretching region, dynamics calculations predict a more intense antisymmetric peak than experiments observe but do not predict the bimodal split symmetric peak that the experiments show. Our work suggests that electronic polarization, missing in the empirical force field, is responsible for the first discrepancy and that quantum nuclear effects, captured neither in density functional theory nor in classical dynamics, may be responsible for the second.

MD investigation on the differences in the dynamic interactions between the specific ligand azamulin and two CYP3A isoforms, 3A4 and 3A5

The unmarked potential drug molecule azamulin has been found to be a specific inhibitor of CYP3A4 and CYP3A5 in recent years, but this molecule also shows different binding ability and affinity to the two CYP3A isoforms. In order to explore the microscopic mechanism, conventional molecular dynamics (MD) simulation methods were performed to study the dynamic interactions between two isoforms and azamulin. The simulation results show that the binding of the ligand leads to different structural properties of two CYP3A proteins. First of all, compared with apo-CYP3A4, the binding of the ligand azamulin can lead to changes in the structural flexibility of CYP3A4, i.e., holo-CYP3A4 is more flexible than apo-CYP3A4. The structural changes of CYP3A5 are just the opposite. The ligand binding increases the rigidity of CYP3A5. Furthermore, the representative structures of the production phase in the MD simulation were in details analyzed to obtain the microscopic interactions between the ligand azamulin and two CYP3A isoforms at the atomic level. It is speculated that the difference of composition and interaction of the active sites is the fundamental cause of the change of structural properties of the two proteins. Communicated by Ramaswamy H. Sarma

Role of Protonation States in the Stability of Molecular Dynamics Simulations of High-Resolution Membrane Protein Structures.

Classical molecular dynamics (MD) simulations provide unmatched spatial and time resolution of protein structure and function. However, the accuracy of MD simulations often depends on the quality of force field parameters and the time scale of sampling. Another limitation of conventional MD simulations is that the protonation states of titratable amino acid residues remain fixed during simulations, even though protonation state changes coupled to conformational dynamics are central to protein function. Due to the uncertainty in selecting protonation states, classical MD simulations are sometimes performed with all amino acids modeled in their standard charged states at pH 7. Here, we performed and analyzed classical MD simulations on high-resolution cryo-EM structures of two large membrane proteins that transfer protons by catalyzing protonation/deprotonation reactions. In simulations performed with titratable amino acids modeled in their standard protonation (charged) states, the structure diverges far from its starting conformation. In comparison, MD simulations performed with predetermined protonation states of amino acid residues reproduce the structural conformation, protein hydration, and protein-water and protein-protein interactions of the structure much better. The results support the notion that it is crucial to perform basic protonation state calculations, especially on structures where protonation changes play an important functional role, prior to the launch of any conventional MD simulations. Furthermore, the combined approach of fast protonation state prediction and MD simulations can provide valuable information about the charge states of amino acids in the cryo-EM sample. Even though accurate prediction of protonation states in proteinaceous environments currently remains a challenge, we introduce an approach of combining pKa prediction with cryo-EM density map analysis that helps in improving not only the protonation state predictions but also the atomic modeling of density data.

The Role of Force Fields and Water Models in Protein Folding and Unfolding Dynamics.

Protein folding is a fascinating, not fully understood phenomenon in biology. Molecular dynamics (MD) simulations are an invaluable tool to study conformational changes in atomistic detail, including folding and unfolding processes of proteins. However, the accuracy of the conformational ensembles derived from MD simulations inevitably relies on the quality of the underlying force field in combination with the respective water model. Here, we investigate protein folding, unfolding, and misfolding of fast-folding proteins by examining different force fields with their recommended water models, i.e., ff14SB with the TIP3P model and ff19SB with the OPC model. To this end, we generated long conventional MD simulations highlighting the perks and pitfalls of these setups. Using Markov state models, we defined kinetically independent conformational substates and emphasized their distinct characteristics, as well as their corresponding state probabilities. Surprisingly, we found substantial differences in thermodynamics and kinetics of protein folding, depending on the combination of the protein force field and water model, originating primarily from the different water models. These results emphasize the importance of carefully choosing the force field and the respective water model as they determine the accuracy of the observed dynamics of folding events. Thus, the findings support the hypothesis that the water model is at least equally important as the force field and hence needs to be considered in future studies investigating protein dynamics and folding in all areas of biophysics.

Size dependent lithium-ion conductivity of solid electrolytes in machine learning molecular dynamics simulations

Solid-state electrolytes are key ingredients in next-generation devices for energy storage and release. Machine learning molecular dynamics (MLMD) has shown great promise in studying the diffusivity of mobile ions in solid-state electrolytes, with much higher efficiency than conventional ab initio molecular dynamics (AIMD). In this work, we combine an efficient embedded atom neural network (EANN) approach and an uncertainty-driven active learning algorithm that optimally selects data points from high-temperature AIMD trajectories to construct ML potentials for solid-state electrolytes and validate this strategy in a benchmark system, Li3YCl6, for which several controversy theoretical results exist. Through systematic MLMD simulations, we find that a typically used small supercell in AIMD simulations fails to predict the supersonic transition at a critical temperature, leading to a significant overestimation of the Li+ conductivity in Li3YCl6 at room temperature. Fortunately, thanks to the scalability of the EANN potential, extended MLMD simulations in a sufficiently large cell does yield a notable change of temperature-dependence in conductivity at ∼420 K and a much lower room-temperature conductivity in excellent with experiment. Interestingly, our results are all based on a semi-local PBE density functional, which was argued unable to predict the superionic transition. We analyze possible reasons of the seemingly inconsistent MLMD results reported in literature with different ML potentials. This work paves the way of simply using high-temperature AIMD data to generate more reliable MLMD results of low-temperature ionic conductivities in solid-state electrolytes.

Implementation and Validation of an OpenMM Plugin for the Deep Potential Representation of Potential Energy.

Machine learning potentials, particularly the deep potential (DP) model, have revolutionized molecular dynamics (MD) simulations, striking a balance between accuracy and computational efficiency. To facilitate the DP model's integration with the popular MD engine OpenMM, we have developed a versatile OpenMM plugin. This plugin supports a range of applications, from conventional MD simulations to alchemical free energy calculations and hybrid DP/MM simulations. Our extensive validation tests encompassed energy conservation in microcanonical ensemble simulations, fidelity in canonical ensemble generation, and the evaluation of the structural, transport, and thermodynamic properties of bulk water. The introduction of this plugin is expected to significantly expand the application scope of DP models within the MD simulation community, representing a major advancement in the field.

Structure Comparison of Beta Amyloid Peptide Aβ 1-42 Isoforms. Molecular Dynamics Modeling.

Beta amyloid peptide Aβ 1-42 (Aβ42) has a unique dual role in the human organism, as both the peptide with an important physiological function and one of the most toxic biological compounds provoking Alzheimer's disease (AD). There are several known Aβ42 isoforms that we discuss here that are highly neurotoxic and lead to the early onset of AD. Aβ42 is an intrinsically disordered protein with no experimentally solved structure under physiological conditions. The objective of this research was to establish the appropriate molecular dynamics (MD) methodology and model a uniform set of structures for the Aβ42 isoforms that form the core of this study. For that purpose, force field selection and verification including convergence testing for MD simulations was made. Replica exchange MD and conventional MD modeling of several Aβ42 and Aβ16 isoforms that have neurotoxic and amyloidogenic effects impacting the severity of Alzheimer's disease were carried out with the optimal force field and solvent parameters. A standardized ensemble of structures for the Aβ42 and Aβ16 isoforms covering 30-50% of the conformational ensembles extracted from the free energy minima was calculated from MD trajectories. The resulting data set of modeled structures includes Aβ42 wild type, isoD7, pS8, D7H, and H6R-Aβ42 and Aβ16 wild type, isoD7, pS8, D7H, and H6R-Aβ16. The representative structures are given in the Supporting Information; they are open for public access. In the study, we also evaluated the differences between the structures of Aβ42 isoforms and speculate on their possible relevance to the known functions. Utilizing several representative structures for a single disordered protein for docking, with their subsequent averaging by conformations, would markedly increase the reliability of docking results.

Unraveling the allosteric mechanisms of prolyl endopeptidases for celiac disease therapy: Insights from molecular dynamics simulations

Prolyl endopeptidases (PEP) from Sphingomonas capsulata (sc) and Myxococcus xanthus (mx) selectively degrade gluten peptides in vitro, offering a potential therapeutic strategy for celiac disease. However, the mechanisms governing the interaction of these enzymes with their substrates remain unclear. In this study, conventional molecular dynamics simulations with a microsecond timescale and targeted molecular dynamics simulations were performed to investigate the native states of mxPEP and scPEP enzymes, as well as their allosteric binding with a representative substrate, namely, Z-Ala-Pro-p-nitroanilide (pNA). The simulations reveal that the native scPEP is in an open state, while the native mxPEP is in a closed state. When pNA approaches a closed mxPEP, it binds to an allosteric pocket located at the first and second β-sheet of the β-propeller domain, inducing the opening of this enzyme. Neither enzyme is active in the open or partly-open states. Enzymatic activity is enabled only when the catalytic pocket in the closed state fully accommodates the substrates. The internal capacity of the catalytic pocket of PEP in the closed state determines the maximum size of the gluten peptides that the enzymes can catalyze. The present work provides essential molecular dynamics information for the redesign or engineering of PEP enzymes.