Analyzing and Biasing Simulations with PLUMED

Dynamic interchange between two protonation states is characteristic of active sites of cholinesterases

Cholinesterases are well-known and widely studied enzymes crucial to human health and involved in neurology, Alzheimer's, and lipid metabolism. The protonation pattern of active sites of cholinesterases influences all the chemical processes within, including reaction, covalent inhibition by nerve agents, and reactivation. Despite its significance, our comprehension of the fine structure of cholinesterases remains limited. In this study, we employed enhanced-sampling quantum-mechanical/molecular-mechanical calculations to show that cholinesterases predominantly operate as dynamic mixtures of two protonation states. The proton transfer between two non-catalytic glutamate residues follows the Grotthuss mechanism facilitated by a mediator water molecule. We show that this uncovered complexity of active sites presents a challenge for classical molecular dynamics simulations and calls for special treatment. The calculated proton transfer barrier of 1.65 kcal/mol initiates a discussion on the potential existence of two coupled low-barrier hydrogen bonds in the inhibited form of butyrylcholinesterase. These findings expand our understanding of structural features expressed by highly evolved enzymes and guide future advances in cholinesterase-related protein and drug design studies.

Mechanism of Abnormal Activation of MEK1 Induced by Dehydroalanine Modification

Mitogen-activated protein kinase kinase 1 (MAPK kinase 1, MEK1) is a key kinase in the mitogen-activated protein kinase (MAPK) signaling pathway. MEK1 mutations have been reported to lead to abnormal activation that is closely related to the malignant growth and spread of various tumors, making it an important target for cancer treatment. Targeting MEK1, four small-molecular drugs have been approved by the FDA, including Trametinib, Cobimetinib, Binimetinib, and Selumetinib. Recently, a study showed that modification with dehydroalanine (Dha) can also lead to abnormal activation of MEK1, which has the potential to promote tumor development. In this study, we used molecular dynamics simulations and metadynamics to explore the mechanism of abnormal activation of MEK1 caused by the Dha modification and predicted the inhibitory effects of four FDA-approved MEK1 inhibitors on the Dha-modified MEK1. The results showed that the mechanism of abnormal activation of MEK1 caused by the Dha modification is due to the movement of the active segment, which opens the active pocket and exposes the catalytic site, leading to sustained abnormal activation of MEK1. Among four FDA-approved inhibitors, only Selumetinib clearly blocks the active site by changing the secondary structure of the active segment from α-helix to disordered loop. Our study will help to explain the mechanism of abnormal activation of MEK1 caused by the Dha modification and provide clues for the development of corresponding inhibitors.

Peptide Drug Design Using Alchemical Free Energy Calculation: An Application and Validation on Agonists of Ghrelin Receptor

With recent large-scale applications and validations, the relative binding free energy (RBFE) calculated using alchemical free energy methods has been proven to be an accurate measure to probe the binding of small-molecule drug candidates. On the other hand, given the flexibility of peptides, it is of great interest to find out whether sufficient sampling could be achieved within the typical time scale of such calculation, and a similar level of accuracy could be reached for peptide drugs. However, the systematic evaluation of such calculations on protein-peptide systems has been less reported. Most reported studies of peptides were restricted to a limited number of data points or lacking experimental support. To demonstrate the applicability of the alchemical free energy method for protein-peptide systems in a typical real-world drug discovery project, we report an application of the thermodynamic integration (TI) method to the RBFE calculation of ghrelin receptor and its peptide agonists. Along with the calculation, the synthesis and in vitro EC50 activity of relamorelin and 17 new peptide derivatives were also reported. A cost-effective criterion to determine the data collection time was proposed for peptides in the TI simulation. The average of three TI repeats yielded a mean absolute error of 0.98 kcal/mol and Pearson's correlation coefficient (R) of 0.77 against the experimental free energy derived from the in vitro EC50 activity, showing good repeatability of the proposed method and a slightly better agreement than the results obtained from the arbitrary time frames up to 20 ns. Although it is limited by having one target and a deduced binding pose, we hope that this study can add some insights into alchemical free energy calculation of protein-peptide systems, providing theoretical assistance to the development of peptide drugs.

Accurate Characterization of Binding Kinetics and Allosteric Mechanisms for the HSP90 Chaperone Inhibitors Using AI-Augmented Integrative Biophysical Studies

The binding kinetics of drugs to their targets are gradually being recognized as a crucial indicator of the efficacy of drugs in vivo, leading to the development of various computational methods for predicting the binding kinetics in recent years. However, compared with the prediction of binding affinity, the underlying structure and dynamic determinants of binding kinetics are more complicated. Efficient and accurate methods for predicting binding kinetics are still lacking. In this study, quantitative structure-kinetics relationship (QSKR) models were developed using 132 inhibitors targeting the ATP binding domain of heat shock protein 90α (HSP90α) to predict the dissociation rate constant (koff), enabling a direct assessment of the drug-target residence time. These models demonstrated good predictive performance, where hydrophobic and hydrogen bond interactions significantly influence the koff prediction. In subsequent applications, our models were used to assist in the discovery of new inhibitors for the N-terminal domain of HSP90α (N-HSP90α), demonstrating predictive capabilities on an experimental validation set with a new scaffold. In X-ray crystallography experiments, the loop-middle conformation of apo N-HSP90α was observed for the first time (previously, the loop-middle conformation had only been observed in holo-N-HSP90α structures). Interestingly, we observed different conformations of apo N-HSP90α simultaneously in an asymmetric unit, which was also observed in a holo-N-HSP90α structure, suggesting an equilibrium of conformations between different states in solution, which could be one of the determinants affecting the binding kinetics of the ligand. Different ligands can undergo conformational selection or alter the equilibrium of conformations, inducing conformational rearrangements and resulting in different effects on binding kinetics. We then used molecular dynamics simulations to describe conformational changes of apo N-HSP90α in different conformational states. In summary, the study of the binding kinetics and molecular mechanisms of N-HSP90α provides valuable information for the development of more targeted therapeutic approaches.

Mechanistic insight into the competition between interfacial and bulk reactions in microdroplets through N2O5 ammonolysis and hydrolysis

Reactive uptake of dinitrogen pentaoxide (N2O5) into aqueous aerosols is a major loss channel for NOx in the troposphere; however, a quantitative understanding of the uptake mechanism is lacking. Herein, a computational chemistry strategy is developed employing high-level quantum chemical methods; the method offers detailed molecular insight into the hydrolysis and ammonolysis mechanisms of N2O5 in microdroplets. Specifically, our calculations estimate the bulk and interfacial hydrolysis rates to be (2.3 ± 1.6) × 10-3 and (6.3 ± 4.2) × 10-7 ns-1, respectively, and ammonolysis competes with hydrolysis at NH3 concentrations above 1.9 × 10-4 mol L-1. The slow interfacial hydrolysis rate suggests that interfacial processes have negligible effect on the hydrolysis of N2O5 in liquid water. In contrast, N2O5 ammonolysis in liquid water is dominated by interfacial processes due to the high interfacial ammonolysis rate. Our findings and strategy are applicable to high-chemical complexity microdroplets.

Multiscale In Silico Study of the Mechanism of Activation of the RtcB Ligase by the PTP1B Phosphatase

Inositol-requiring enzyme 1 (IRE1) is a transmembrane sensor that is part of a trio of sensors responsible for controlling the unfolded protein response within the endoplasmic reticulum (ER). Upon the accumulation of unfolded or misfolded proteins in the ER, IRE1 becomes activated and initiates the cleavage of a 26-nucleotide intron from human X-box-containing protein 1 (XBP1). The cleavage is mediated by the RtcB ligase enzyme, which splices together two exons, resulting in the formation of the spliced isoform XBP1s. The XBP1s isoform activates the transcription of genes involved in ER-associated degradation to maintain cellular homeostasis. The catalytic activity of RtcB is controlled by the phosphorylation and dephosphorylation of three tyrosine residues (Y306, Y316, and Y475), which are regulated by the ABL1 tyrosine kinase and PTP1B phosphatase, respectively. This study focuses on investigating the mechanism by which the PTP1B phosphatase activates the RtcB ligase using a range of advanced in silico methods. Protein-protein docking identified key interacting residues between RtcB and PTP1B. Notably, the phosphorylated Tyr306 formed hydrogen bonds and salt bridge interactions with the "gatekeeper" residues Arg47 and Lys120 of the inactive PTP1B. Classical molecular dynamics simulation emphasized the crucial role of Asp181 in the activation of PTP1B, driving the conformational change from an open to a closed state of the WPD-loop. Furthermore, QM/MM-MD simulations provided insights into the free energy landscape of the dephosphorylation reaction mechanism of RtcB, which is mediated by the PTP1B phosphatase.

Understanding the free-energy landscape of phase separation in lipid bilayers using molecular dynamics

Liquid-liquid phase separation inside the cell often results in biological condensates that can critically affect cell homeostasis. Such phase separation events occur in multiple parts of cells, including the cell membranes, where the "lipid raft" hypothesis posits the formation of ordered domains floating in a sea of disordered lipids. The resulting lipid domains often have functional roles. However, the thermodynamics of lipid phase separation and their resulting mechanistic effects on cell function and dysfunction are poorly understood. Understanding such complex phenomena in cell membranes, with their diverse lipid compositions, is exceptionally difficult. For these reasons, simple model systems that can recapitulate similar behavior are widely used to study this phenomenon. Despite these simplifications, the timescale and length scales of domain formation pose a challenge for molecular dynamics (MD) simulations. Thus, most MD studies focus on spontaneous lipid phase separation-essentially measuring the sign (but not the amplitude) of the free-energy change upon separation-rather than directly interrogating the thermodynamics. Here, we propose a proof-of-concept pipeline that can directly measure this free energy by combining coarse-grained MD with enhanced sampling protocols using a novel collective variable. This approach will be a useful tool to help connect the thermodynamics of phase separation with the mechanistic insights already available from MD simulations.

Accelerated Photolysis of H2O2 at the Air-Water Interface of a Microdroplet

Photochemical homolysis of hydrogen peroxide (H2O2) occurs widely in nature and is a key source of hydroxyl radicals (·OH). The kinetics of H2O2 photolysis play a pivotal role in determining the efficiency of ·OH production, which is currently mainly investigated in bulk systems. Here, we report considerably accelerated H2O2 photolysis at the air-water interface of microdroplets, with a rate 1.9 × 103 times faster than that in bulk water. Our simulations show that due to the trans quasiplanar conformational preference of H2O2 at the air-water interface compared to the bulk or gas phase, the absorption peak in the spectrum of H2O2 is significantly redshifted by 45 nm, corresponding to greater absorbance of photons in the sunlight spectrum and faster photolysis of H2O2. This discovery has great potential to solve current problems associated with ·OH-centered heterogeneous photochemical processes in aerosols. For instance, we show that accelerated H2O2 photolysis in microdroplets could lead to markedly enhanced oxidation of SO2 and volatile organic compounds.

Temperature sensitive contact modes allosterically gate TRPV3

TRPV Ion channels are sophisticated molecular sensors designed to respond to distinct temperature thresholds. The recent surge in cryo-EM structures has provided numerous insights into the structural rearrangements accompanying their opening and closing; however, the molecular mechanisms by which TRPV channels establish precise and robust temperature sensing remain elusive. In this work we employ molecular simulations, multi-ensemble contact analysis, graph theory, and machine learning techniques to reveal the temperature-sensitive residue-residue interactions driving allostery in TRPV3. We find that groups of residues exhibiting similar temperature-dependent contact frequency profiles cluster at specific regions of the channel. The dominant mode clusters on the ankyrin repeat domain and displays a linear melting trend while others display non-linear trends. These modes describe the residue-level temperature response patterns that underlie the channel's functional dynamics. With network analysis, we find that the community structure of the channel changes with temperature. And that a network of high centrality contacts connects distant regions of the protomer to the gate, serving as a means for the temperature-sensitive contact modes to allosterically regulate channel gating. Using a random forest model, we show that the contact states of specific temperature-sensitive modes are indeed predictive of the channel gate's state. Supporting the physical validity of these modes and networks are several residues identified with our analyses that are reported in literature to be functionally critical. Our results offer high resolution insight into thermo-TRP channel function and demonstrate the utility of temperature-sensitive contact analysis.

Structure and Dynamics of DNA and RNA Double Helices Formed by d(CTG), d(GTC), r(CUG), and r(GUC) Trinucleotide Repeats and Associated DNA-RNA Hybrids

Myotonic dystrophy type 1 is the most frequent form of muscular dystrophy in adults caused by an abnormal expansion of the CTG trinucleotide. Both the expanded DNA and the expanded CUG RNA transcript can fold into hairpins. Co-transcriptional formation of stable RNA·DNA hybrids can also enhance the instability of repeat tracts. We performed molecular dynamics simulations of homoduplexes associated with the disease, d(CTG)n and r(CUG)n, and their corresponding r(CAG)n:d(CTG)n and r(CUG)n:d(CAG)n hybrids that can form under bidirectional transcription and of non-pathological d(GTC)n and d(GUC)n homoduplexes. We characterized their conformations, stability, and dynamics and found that the U·U and T·T mismatches are dynamic, favoring anti-anti conformations inside the helical core, followed by anti-syn and syn-syn conformations. For DNA, the secondary minima in the non-expanding d(GTC)n helices are deeper, wider, and longer-lived than those in d(CTG)n, which constitutes another biophysical factor further differentiating the expanding and non-expanding sequences. The hybrid helices are closer to A-RNA, with the A-T and A-U pairs forming two stable Watson-Crick hydrogen bonds. The neutralizing ion distribution around the non-canonical pairs is also described.

Kinetics from Metadynamics: Principles, Applications, and Outlook

Metadynamics is a popular enhanced sampling algorithm for computing the free energy landscape of rare events by using molecular dynamics simulation. Ten years ago, Tiwary and Parrinello introduced the infrequent metadynamics approach for calculating the kinetics of transitions across free energy barriers. Since then, metadynamics-based methods for obtaining rate constants have attracted significant attention in computational molecular science. Such methods have been applied to study a wide range of problems, including protein-ligand binding, protein folding, conformational transitions, chemical reactions, catalysis, and nucleation. Here, we review the principles of elucidating kinetics from metadynamics-like approaches, subsequent methodological developments in this area, and successful applications on chemical, biological, and material systems. We also highlight the challenges of reconstructing accurate kinetics from enhanced sampling simulations and the scope of future developments.

Temperature-Sensitive Contact Modes Allosterically Gate TRPV3

TRPV Ion channels are sophisticated molecular sensors designed to respond to distinct temperature thresholds. The recent surge in cryo-EM structures has provided numerous insights into the structural rearrangements accompanying their opening and closing; however, the molecular mechanisms by which TRPV channels establish precise and robust temperature sensing remain elusive. In this work we employ molecular simulations, multi-ensemble contact analysis, graph theory, and machine learning techniques to reveal the temperature-sensitive residue-residue interactions driving allostery in TRPV3. We find that groups of residues exhibiting similar temperature-dependent contact frequency profiles cluster at specific regions of the channel. The dominant mode clusters on the ankyrin repeat domain and displays a linear melting trend while others display non-linear trends. These modes describe the residue-level temperature response patterns that underlie the channel's functional dynamics. With network analysis, we find that the community structure of the channel changes with temperature. And that a network of high centrality contacts connects distant regions of the protomer to the gate, serving as a means for the temperature-sensitive contact modes to allosterically regulate channel gating. Using a random forest model, we show that the contact states of specific temperature-sensitive modes are indeed predictive of the channel gate's state. Supporting the physical validity of these modes and networks are several residues identified with our analyses that are reported in literature to be functionally critical. Our results offer high resolution insight into thermo-TRP channel function and demonstrate the utility of temperature-sensitive contact analysis.

Combining Experiments and Simulations to Examine the Temperature-Dependent Behavior of a Disordered Protein

Intrinsically disordered proteins are a class of proteins that lack stable folded conformations and instead adopt a range of conformations that determine their biochemical functions. The temperature-dependent behavior of such disordered proteins is complex and can vary depending on the specific protein and environment. Here, we have used molecular dynamics simulations and previously published experimental data to investigate the temperature-dependent behavior of histatin 5, a 24-residue-long polypeptide. We examined the hypothesis that histatin 5 undergoes a loss of polyproline II (PPII) structure with increasing temperature, leading to more compact conformations. We found that the conformational ensembles generated by the simulations generally agree with small-angle X-ray scattering data for histatin 5, but show some discrepancies with the hydrodynamic radius as probed by pulsed-field gradient NMR spectroscopy, and with the secondary structure information derived from circular dichroism. We attempted to reconcile these differences by reweighting the conformational ensembles against the scattering and NMR data. By doing so, we were in part able to capture the temperature-dependent behavior of histatin 5 and to link the observed decrease in hydrodynamic radius with increasing temperature to a loss of PPII structure. We were, however, unable to achieve agreement with both the scattering and NMR data within experimental errors. We discuss different possible reasons for this including inaccuracies in the force field, differences in conditions of the NMR and scattering experiments, and issues related to the calculation of the hydrodynamic radius from conformational ensembles. Our study highlights the importance of integrating multiple types of experimental data when modeling conformational ensembles of disordered proteins and how environmental factors such as the temperature influence them.

Binding pocket dynamics along the recovery stroke of human β-cardiac myosin

The druggability of small-molecule binding sites can be significantly affected by protein motions and conformational changes. Ligand binding, protein dynamics and protein function have been shown to be closely interconnected in myosins. The breakthrough discovery of omecamtiv mecarbil (OM) has led to an increased interest in small molecules that can target myosin and modulate its function for therapeutic purposes (myosin modulators). In this work, we use a combination of computational methods, including steered molecular dynamics, umbrella sampling and binding pocket tracking tools, to follow the evolution of the OM binding site during the recovery stroke transition of human β-cardiac myosin. We found that steering two internal coordinates of the motor domain can recapture the main features of the transition and in particular the rearrangements of the binding site, which shows significant changes in size, shape and composition. Possible intermediate conformations were also identified, in remarkable agreement with experimental findings. The differences in the binding site properties observed along the transition can be exploited for the future development of conformation-selective myosin modulators.

Sub-millisecond conformational dynamics of the A2A adenosine receptor revealed by single-molecule FRET

The complex pharmacology of G-protein-coupled receptors (GPCRs) is defined by their multi-state conformational dynamics. Single-molecule Förster Resonance Energy Transfer (smFRET) is well suited to quantify dynamics for individual protein molecules; however, its application to GPCRs is challenging. Therefore, smFRET has been limited to studies of inter-receptor interactions in cellular membranes and receptors in detergent environments. Here, we performed smFRET experiments on functionally active human A_2A adenosine receptor (A_2AAR) molecules embedded in freely diffusing lipid nanodiscs to study their intramolecular conformational dynamics. We propose a dynamic model of A_2AAR activation that involves a slow (>2 ms) exchange between the active-like and inactive-like conformations in both apo and antagonist-bound A_2AAR, explaining the receptor's constitutive activity. For the agonist-bound A_2AAR, we detected faster (390 ± 80 µs) ligand efficacy-dependent dynamics. Our work establishes a general smFRET platform for GPCR investigations that can potentially be used for drug screening and/or mechanism-of-action studies.

Free Energy Differences from Molecular Simulations: Exact Confidence Intervals from Transition Counts

Here, we demonstrate a method to estimate the uncertainty (confidence intervals and standard errors) of free energy differences calculated by molecular simulations. The widths of confidence intervals and standard errors can be calculated solely from temperature and the number of transitions between states. Uncertainty (95% confidence interval) lower than ±1 kcal/mol can be achieved by a simulation with four forward and four reverse transitions. For a two-state Markovian system, the confidence interval is exact, regardless the number of transitions.

Modern non-polarizable force fields diverge in modeling the enzyme-substrate complex of a canonical serine protease

Classical molecular dynamics simulation is a powerful and established method of modern computational chemistry. Being able to obtain accurate information on molecular behavior is crucial to get valuable insights into structure-function relationships that translate into fundamental findings and practical applications. Active sites of enzymes are known to be particularly intricate, therefore, simpler non-polarizable force fields may provide an inaccurate description. In this work, we addressed this hypothesis in a case of a canonical serine triad protease trypsin in its complex with a substrate-mimicking inhibitor. We tested six modern and popular force fields to find that significantly diverging results may be obtained. Amber FB-15 and OPLS-AA/M turned out to model the active site incorrectly. Amber ff19sb and ff15ipq demonstrated mixed performance. The best performing force fields were CHARMM36m and Amber ff99sb-ildn, therefore, they are recommended for use with this and related systems. We speculate that a similar lack of cross-force field convergence may be characteristic of other enzymatic systems. Therefore, we advocate for careful consideration of different force fields in any study within the field of computational enzymology.

Funnel metadynamics and behavioral studies reveal complex effect of D2AAK1 ligand on anxiety-like processes

Anxiety is a troublesome symptom for many patients, especially those suffering from schizophrenia. Its regulation involves serotonin receptors, targeted e.g. by antipsychotics or psychedelics such as LSD. 5-HT_2A receptors are known for an extremely long LSD residence time, enabling minute doses to exert a long-lasting effect. In this work, we explore the changes in anxiety-like processes induced by the previously reported antipsychotic, D2AAK1. In vivo studies revealed that the effect of D2AAK1 on the anxiety is mediated through serotonin 5-HT_1A and 5-HT_2A receptors, and that it is time-dependent (anxiogenic after 30 min, anxiolytic after 60 min) and dose-dependent. The funnel metadynamics simulations suggest complicated ligand-5HT_2AR interactions, involving an allosteric site located under the third extracellular loop, which is a possible explanation of the time-dependency. The binding of D2AAK1 at the allosteric site results in a broader opening of the extracellular receptor entry, possibly altering the binding kinetics of orthosteric ligands.

Temperature-sensitive contacts in disordered loops tune enzyme I activity

Homologous enzymes with identical folds often exhibit different thermal and kinetic behaviors. Understanding how an enzyme sequence encodes catalytic activity at functionally optimal temperatures is a fundamental problem in biophysics. Recently it was shown that the residues that tune catalytic activities of thermophilic/mesophilic variants of the C-terminal domain of bacterial enzyme I (EIC) are largely localized within disordered loops, offering a model system with which to investigate this phenomenon. In this work, we use molecular dynamics simulations and mutagenesis experiments to reveal a mechanism of sequence-dependent activity tuning of EIC homologs. We find that a network of contacts in the catalytic loops is particularly sensitive to changes in temperature, with some contacts exhibiting distinct linear or nonlinear temperature-dependent trends. Moreover, these trends define structurally clustered dynamical modes and can distinguish regions that tend toward order or disorder at higher temperatures. Assaying several thermophilic EIC mutants, we show that complementary mesophilic mutations to the most temperature-sensitive positions exhibit the most enhanced activity, while mutations to relatively temperature insensitive positions exhibit the least enhanced activities. These results provide a mechanistic explanation of sequence-dependent temperature tuning and offer a computational method for rational enzyme modification.

Computer simulation of the early stages of self-assembly and thermal decomposition of ZIF-8

We employ all-atom well-tempered metadynamics simulations to study the mechanistic details of both the early stages of nucleation and crystal decomposition for the benchmark metal–organic framework (MOF) ZIF-8. To do so, we developed and validated a force field that reliably models the modes of coordination bonds via a Morse potential functional form and employs cationic and anionic dummy atoms to capture coordination symmetry. We also explored a set of physically relevant collective variables and carefully selected an appropriate subset for our problem at hand. After a rapid increase of the Zn–N connectivity, we observe the evaporation of small clusters in favor of a few large clusters, which leads to the formation of an amorphous highly connected aggregate. [Formula: see text] and [Formula: see text] complexes are observed with lifetimes in the order of a few picoseconds, while larger structures, such as four-, five-, and six-membered rings, have substantially longer lifetimes of a few nanoseconds. The free ligands act as “templating agents” for the formation of sodalite cages. ZIF-8 crystal decomposition results in the formation of a vitreous phase. Our findings contribute to a fundamental understanding of MOF’s synthesis that paves the way to controlling synthesis products. Furthermore, our developed force field and methodology can be applied to model solution processes that require coordination bond reactivity for other ZIFs besides ZIF-8.

Between Protein Fold and Nucleophile Identity: Multiscale Modeling of the TEV Protease Enzyme-Substrate Complex

The cysteine protease from the tobacco etch virus (TEVp) is a well-known and widely utilized enzyme. TEVp's chymotrypsin-like fold is generally associated with serine catalytic triads that differ in terms of a reaction mechanism from the most well-studied papain-like cysteine proteases. The question of what dominates the TEVp mechanism, nucleophile identity, or structural composition has never been previously addressed. Here, we use enhanced sampling multiscale modeling to uncover that TEVp combines the features of two worlds in such a way that potentially hampers its activity. We show that TEVp cysteine is strictly in the anionic form in a free enzyme similar to papain. Peptide binding shifts the equilibrium toward the nucleophile's protonated form, characteristic of chymotrypsin-like proteases, although the cysteinyl anion form is still present and interconversion is rapid. This way cysteine protonation generates enzyme states that are a diversion from the most effective course of action, with only 13.2% of Michaelis complex sub-states able to initiate the reaction. As a result, we propose an updated view on the reaction mechanism catalyzed by TEVp. We also demonstrate that AlphaFold is able to construct protease-substrate complexes with high accuracy. We propose that our findings open a way for its industrious use in enzymological tasks. Unique features of TEVp discovered in this work open a discussion on the evolutionary history and trade-offs of optimizing serine triad-associated folds to cysteine as a nucleophile.

Long-range allostery mediates the regulation of plasminogen activator inhibitor-1 by cell adhesion factor vitronectin

The serpin plasminogen activator inhibitor 1 (PAI-1) spontaneously undergoes a massive structural change from a metastable and active conformation, with a solvent-accessible reactive center loop (RCL), to a stable, inactive, or latent conformation, with the RCL inserted into the central β-sheet. Physiologically, conversion to the latent state is regulated by the binding of vitronectin, which hinders the latency transition rate approximately twofold. The molecular mechanisms leading to this rate change are unclear. Here, we investigated the effects of vitronectin on the PAI-1 latency transition using all-atom path sampling simulations in explicit solvent. In simulated latency transitions of free PAI-1, the RCL is quite mobile as is the gate, the region that impedes RCL access to the central β-sheet. This mobility allows the formation of a transient salt bridge that facilitates the transition; this finding rationalizes existing mutagenesis results. Vitronectin binding reduces RCL and gate mobility by allosterically rigidifying structural elements over 40 Å away from the binding site, thus blocking transition to the latent conformation. The effects of vitronectin are propagated by a network of dynamically correlated residues including a number of conserved sites that were previously identified as important for PAI-1 stability. Simulations also revealed a transient pocket populated only in the vitronectin-bound state, corresponding to a cryptic drug-binding site identified by crystallography. Overall, these results shed new light on PAI-1 latency transition regulation by vitronectin and illustrate the potential of path sampling simulations for understanding functional protein conformational changes and for facilitating drug discovery.

Deciphering the Alphabet of Disorder-Glu and Asp Act Differently on Local but Not Global Properties

Compared to folded proteins, the sequences of intrinsically disordered proteins (IDPs) are enriched in polar and charged amino acids. Glutamate is one of the most enriched amino acids in IDPs, while the chemically similar amino acid aspartate is less enriched. So far, the underlying functional differences between glutamates and aspartates in IDPs remain poorly understood. In this study, we examine the differential effects of aspartate and glutamates in IDPs by comparing the function and conformational ensemble of glutamate and aspartate variants of the disordered protein Dss1, using a range of assays, including interaction studies, nuclear magnetic resonance spectroscopy, small-angle X-ray scattering and molecular dynamics simulation. First, we analyze the sequences of the rapidly growing database of experimentally verified IDPs (DisProt) and show that glutamate enrichment is not caused by a taxonomy bias in IDPs. From analyses of local and global structural properties as well as cell growth and protein-protein interactions using a model acidic IDP from yeast and three Glu/Asp variants, we find that while the Glu/Asp variants support similar function and global dimensions, the variants differ in their binding affinities and population of local transient structural elements. We speculate that these local structural differences may play roles in functional diversity, where glutamates can support increased helicity, important for folding and binding, while aspartates support extended structures and form helical caps, as well as playing more relevant roles in, e.g., transactivation domains and ion-binding.

QM/MM Well-Tempered Metadynamics Study of the Mechanism of XBP1 mRNA Cleavage by Inositol Requiring Enzyme 1α RNase

A range of in silico methodologies were herein employed to study the unconventional XBP1 mRNA cleavage mechanism performed by the unfolded protein response (UPR) mediator Inositol Requiring Enzyme 1α (IRE1). Using Protein-RNA molecular docking along with a series of extensive restrained/unrestrained atomistic molecular dynamics (MD) simulations, the dynamical behavior of the system was evaluated and a reliable model of the IRE1/XBP1 mRNA complex was constructed. From a series of well-converged quantum mechanics molecular mechanics well-tempered metadynamics (QM/MM WT-MetaD) simulations using the Grimme dispersion interaction corrected semiempirical parametrization method 6 level of theory (PM6-D3) and the AMBER14SB-OL3 force field, the free energy profile of the cleavage mechanism was determined, along with intermediates and transition state structures. The results show two distinct reaction paths based on general acid-general base type mechanisms, with different activation energies that perfectly match observations from experimental mutagenesis data. The study brings unique atomistic insights into the cleavage mechanism of XBP1 mRNA by IRE1 and clarifies the roles of the catalytic residues H910 and Y892. Increased understanding of the details in UPR signaling can assist in the development of new therapeutic agents for its modulation.

Quantifying Energetic and Entropic Pathways in Molecular Systems

When examining dynamics occurring at nonzero temperatures, both energy and entropy must be taken into account to describe activated barrier crossing events. Furthermore, good reaction coordinates need to be constructed to describe different metastable states and the transition mechanisms between them. Here we use a physics-based machine learning method called state predictive information bottleneck (SPIB) to find nonlinear reaction coordinates for three systems of varying complexity. SPIB is able to correctly predict an entropic bottleneck for an analytical flat-energy double-well system and identify the entropy- and energy-dominated pathways for an analytical four-well system. Finally, for a simulation of benzoic acid permeation through a lipid bilayer, SPIB is able to discover the entropic and energetic barriers to the permeation process. Given these results, we thus establish that SPIB is a reasonable and robust method for finding the important entropy, energy, and enthalpy barriers in physical systems, which can then be used to enhance the understanding and sampling of different activated mechanisms.

Hybrid drug-screening strategy identifies potential SARS-CoV-2 cell-entry inhibitors targeting human transmembrane serine protease

The spread of coronavirus infectious disease (COVID-19) is associated with the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which has risked public health more than any other infectious disease. Researchers around the globe use multiple approaches to identify an effective approved drug (drug repurposing) that treats viral infections. Most of the drug repurposing approaches target spike protein or main protease. Here we use transmembrane serine protease 2 (TMPRSS2) as a target that can prevent the virus entry into the cell by interacting with the surface receptors. By hypothesizing that the TMPRSS2 binders may help prevent the virus entry into the cell, we performed a systematic drug screening over the current approved drug database. Furthermore, we screened the Enamine REAL fragments dataset against the TMPRSS2 and presented nine potential drug-like compounds that give us clues about which kinds of groups the pocket prefers to bind, aiding future structure-based drug design for COVID-19. Also, we employ molecular dynamics simulations, binding free energy calculations, and well-tempered metadynamics to validate the obtained candidate drug and fragment list. Our results suggested three potential FDA-approved drugs against human TMPRSS2 as a target. These findings may pave the way for more drugs to be exposed to TMPRSS2, and testing the efficacy of these drugs with biochemical experiments will help improve COVID-19 treatment.

Modelling Protein Plasticity: The Example of Frataxin and Its Variants

Frataxin (FXN) is a protein involved in storage and delivery of iron in the mitochondria. Single-point mutations in the FXN gene lead to reduced production of functional frataxin, with the consequent dyshomeostasis of iron. FXN variants are at the basis of neurological impairment (the Friedreich's ataxia) and several types of cancer. By using altruistic metadynamics in conjunction with the maximal constrained entropy principle, we estimate the change of free energy in the protein unfolding of frataxin and of some of its pathological mutants. The sampled configurations highlight differences between the wild-type and mutated sequences in the stability of the folded state. In partial agreement with thermodynamic experiments, where most of the analyzed variants are characterized by lower thermal stability compared to wild type, the D104G variant is found with a stability comparable to the wild-type sequence and a lower water-accessible surface area. These observations, obtained with the new approach we propose in our work, point to a functional switch, affected by single-point mutations, of frataxin from iron storage to iron release. The method is suitable to investigate wide structural changes in proteins in general, after a proper tuning of the chosen collective variable used to perform the transition.

Modeling the Structure and Interactions of Intrinsically Disordered Peptides with Multiple Replica, Metadynamics-Based Sampling Methods and Force-Field Combinations

Intrinsically disordered proteins play a key role in many biological processes, including the formation of biomolecular condensates within cells. A detailed characterization of their configurational ensemble and structure-function paradigm is crucial for understanding their biological activity and for exploiting them as building blocks in material sciences. In this work, we incorporate bias-exchange metadynamics and parallel-tempering well-tempered metadynamics with CHARMM36m and CHARMM22* to explore the structural and thermodynamic characteristics of a short archetypal disordered sequence derived from a DEAD-box protein. The conformational landscapes emerging from our simulations are largely congruent across methods and force fields. Nevertheless, differences in fine details emerge from varying combinations of force-fields and sampling methods. For this protein, our analysis identifies features that help to explain the low propensity of this sequence to undergo self-association in vitro, which are common to all force-field/sampling method combinations. Overall, our work demonstrates the importance of using multiple force-field and sampling method combinations for accurate structural and thermodynamic information in the study of disordered proteins.

SARS-CoV-2 Omicron spike glycoprotein receptor binding domain exhibits super-binder ability with ACE2 but not convalescent monoclonal antibody

SARS-CoV-2, the causative virus for COVID-19 has now super-mutated into the Omicron (Om) variant. On its spike (S) glycoprotein alone, more than 30 substitutions have been characterized with 15 within the receptor binding domain (RBD); It therefore calls to question the transmissibility and antibody escapability of Omicron. This study was setup to investigate the Omicron RBD's interaction with ACE2 (host receptor) and a SARS-CoV-2 neutralizing monoclonal antibody (mAb). In-silico mutagenesis was used to generate the Om-RBD in complex with ACE2 or mAb from the wildtype. HDOCK server was used to redock and score the mAbs in Om-RBD bound state relative to the wildtype. Stability of interaction between all complexes were investigated using all-atom molecular dynamics (MD). Analyses of trajectories showed that Om-RBD has evolved into an efficient ACE2 binder, via pi-pi (Om-RBD-Y501/ACE2-Y41) and salt-bridge (Om-RBD-K493/ACE2-Y41) interactions. Conversely, in binding mAb, it has become less efficient (Center of mass distance of RBD from mAb complex, wildtype ≈ 30 Å, Omicron ≈ 41 Å). Disruption of Om-RBD/mAb complex resulted from loose interaction between Om-RBD and the light chain complementarity-determining region residues. Omicron is expected to be better transmissible and less efficiently interacting with neutralizing convalescent mAbs with consequences on transmissibility provided other mutations within the S protein similarly promote cell fusion and viral entry.

Metadynamics Simulations to Study the Structural Ensembles and Binding Processes of Intrinsically Disordered Proteins

The structures of intrinsically disordered proteins (IDPs) are highly dynamic. It is hard to characterize the structures of these proteins experimentally. Molecular dynamics (MD) simulation is a powerful tool in the understanding of protein dynamic structures and function. This chapter describes the application of metadynamics-based enhanced sampling methods in the study of phosphorylation regulation on the structure of kinase-inducible domains (KID). The structural properties of free pKID and KID were obtained by parallel tempering metadynamics combined with well-tempered ensemble (PTMetaD WTE) method, and the binding free energy surfaces of pKID/KID and KIX were characterized by bias-exchanged metadynamics (BE-MetaD) simulations.

Markovian Weighted Ensemble Milestoning (M-WEM): Long-Time Kinetics from Short Trajectories

We introduce a rare-event sampling scheme, named Markovian Weighted Ensemble Milestoning (M-WEM), which inlays a weighted ensemble framework within a Markovian milestoning theory to efficiently calculate thermodynamic and kinetic properties of long-time-scale biomolecular processes from short atomistic molecular dynamics simulations. M-WEM is tested on the Müller-Brown potential model, the conformational switching in alanine dipeptide, and the millisecond time-scale protein-ligand unbinding in a trypsin-benzamidine complex. Not only can M-WEM predict the kinetics of these processes with quantitative accuracy but it also allows for a scheme to reconstruct a multidimensional free-energy landscape along additional degrees of freedom, which are not part of the milestoning progress coordinate. For the ligand-receptor system, the experimental residence time, association and dissociation kinetics, and binding free energy could be reproduced using M-WEM within a simulation time of a few hundreds of nanoseconds, which is a fraction of the computational cost of other currently available methods, and close to 4 orders of magnitude less than the experimental residence time. Due to the high accuracy and low computational cost, the M-WEM approach can find potential applications in kinetics and free-energy-based computational drug design.

CHARMM-GUI Drude prepper for molecular dynamics simulation using the classical Drude polarizable force field

Explicit treatment of electronic polarizability in empirical force fields (FFs) represents an extension over a traditional additive or pairwise FF and provides a more realistic model of the variations in electronic structure in condensed phase, macromolecular simulations. To facilitate utilization of the polarizable FF based on the classical Drude oscillator model, Drude Prepper has been developed in CHARMM-GUI. Drude Prepper ingests additive CHARMM protein structures file (PSF) and pre-equilibrated coordinates in CHARMM, PDB, or NAMD format, from which the molecular components of the system are identified. These include all residues and patches connecting those residues along with water, ions, and other solute molecules. This information is then used to construct the Drude FF-based PSF using molecular generation capabilities in CHARMM, followed by minimization and equilibration. In addition, inputs are generated for molecular dynamics (MD) simulations using CHARMM, GROMACS, NAMD, and OpenMM. Validation of the Drude Prepper protocol and inputs is performed through conversion and MD simulations of various heterogeneous systems that include proteins, nucleic acids, lipids, polysaccharides, and atomic ions using the aforementioned simulation packages. Stable simulations are obtained in all studied systems, including 5 μs simulation of ubiquitin, verifying the integrity of the generated Drude PSFs. In addition, the ability of the Drude FF to model variations in electronic structure is shown through dipole moment analysis in selected systems. The capabilities and availability of Drude Prepper in CHARMM-GUI is anticipated to greatly facilitate the application of the Drude FF to a range of condensed phase, macromolecular systems.

Probing the Suitability of Different Ca2+ Parameters for Long Simulations of Diisopropyl Fluorophosphatase

Organophosphate hydrolases are promising as potential biotherapeutic agents to treat poisoning with pesticides or nerve gases. However, these enzymes often need to be further engineered in order to become useful in practice. One example of such enhancement is the alteration of enantioselectivity of diisopropyl fluorophosphatase (DFPase). Molecular modeling techniques offer a unique opportunity to address this task rationally by providing a physical description of the substrate-binding process. However, DFPase is a metalloenzyme, and correct modeling of metal cations is a challenging task generally coming with a tradeoff between simulation speed and accuracy. Here, we probe several molecular mechanical parameter combinations for their ability to empower long simulations needed to achieve a quantitative description of substrate binding. We demonstrate that a combination of the Amber19sb force field with the recently developed 12-6 Ca2+ models allows us to both correctly model DFPase and obtain new insights into the DFP binding process.

Three Reasons Why Aspartic Acid and Glutamic Acid Sequences Have a Surprisingly Different Influence on Mineralization

Understanding the role of polymers rich in aspartic acid (Asp) and glutamic acid (Glu) is the key to gaining precise control over mineralization processes. Despite their chemical similarity, experiments revealed a surprisingly different influence of Asp and Glu sequences. We conducted molecular dynamics simulations of Asp and Glu peptides in the presence of calcium and chloride ions to elucidate the underlying phenomena. In line with experimental differences, in our simulations, we indeed find strong differences in the way the peptides interact with ions in solution. The investigated Asp pentapeptide tends to pull a lot of ions into its vicinity, and many structures with clusters of calcium and chloride ions on the surface of the peptide can be observed. Under the same conditions, comparatively fewer ions can be found in proximity of the investigated Glu pentapeptide, and the structures are characterized by single calcium ions bound to multiple carboxylate groups. Based on our simulation data, we identified three reasons contributing to these differences, leading to a new level of understanding additive-ion interactions.

How to Determine Accurate Conformational Ensembles by Metadynamics Metainference: A Chignolin Study Case

The reliability and usefulness of molecular dynamics simulations of equilibrium processes rests on their statistical precision and their capability to generate conformational ensembles in agreement with available experimental knowledge. Metadynamics Metainference (M&M), coupling molecular dynamics with the enhanced sampling ability of Metadynamics and with the ability to integrate experimental information of Metainference, can in principle achieve both goals. Here we show that three different Metadynamics setups provide converged estimate of the populations of the three-states populated by a model peptide. Errors are estimated correctly by block averaging, but higher precision is obtained by performing independent replicates. One effect of Metadynamics is that of dramatically decreasing the number of effective frames resulting from the simulations and this is relevant for M&M where the number of replicas should be large enough to capture the conformational heterogeneity behind the experimental data. Our simulations allow also us to propose that monitoring the relative error associated with conformational averaging can help to determine the minimum number of replicas to be simulated in the context of M&M simulations. Altogether our data provides useful indication on how to generate sound conformational ensemble in agreement with experimental data.

Transient Excursions to Membrane Core as Determinants of Influenza Virus Fusion Peptide Activity

Fusion of viral and host cell membranes is a critical step in the life cycle of enveloped viruses. In the case of influenza virus, it is mediated by subunit 2 of hemagglutinin (HA) glycoprotein whose N-terminal fragments insert into the target membrane and initiate lipid exchange. These isolated fragments, known as fusion peptides (HAfp), already possess own fusogenic activity towards liposomes. Although they have long been studied with the hope to uncover the details of HA-mediated fusion, their actual mechanism of action remains elusive. Here, we use extensive molecular dynamics simulations combined with experimental studies of three HAfp variants to fully characterize their free energy landscape and interaction with lipid bilayer. In addition to customary assumed peptides localization at lipid-water interface, we characterize membrane-spanning configurations, which turn out to be metastable for active HAfps and unstable for the fusion inactive W14A mutant. We show that, while the degree of membrane perturbation by surface peptide configurations is relatively low and does not show any mutation-related differences, the effect of deeply inserted configurations is significant and correlates with insertion depth of the N-terminal amino group which is the highest for the wild type HAfp. Finally, we demonstrate the feasibility of spontaneous peptide transition to intramembrane location and the critical role of strictly conserved tryptofan residue 14 in this process.

Distinct axial and lateral interactions within homologous filaments dictate the signaling specificity and order of the AIM2-ASC inflammasome

Inflammasomes are filamentous signaling platforms integral to innate immunity. Currently, little is known about how these structurally similar filaments recognize and distinguish one another. A cryo-EM structure of the AIM2PYD filament reveals that the architecture of the upstream filament is essentially identical to that of the adaptor ASCPYD filament. In silico simulations using Rosetta and molecular dynamics followed by biochemical and cellular experiments consistently demonstrate that individual filaments assemble bidirectionally. By contrast, the recognition between AIM2 and ASC requires at least one to be oligomeric and occurs in a head-to-tail manner. Using in silico mutagenesis as a guide, we also identify specific axial and lateral interfaces that dictate the recognition and distinction between AIM2 and ASC filaments. Together, the results here provide a robust framework for delineating the signaling specificity and order of inflammasomes.

Comparison of Different Chitosan Lipid Nanoparticles for Improved Ophthalmic Tetrandrine Delivery: Formulation, Characterization, Pharmacokinetic and Molecular Dynamics Simulation

In this study, three different chitosan, namely carboxymethyl chitosan (CMC), hydroxypropyl chitosan (HPC) and trimethyl chitosan (TMC) were used as cationic materials to prepare tetrandrine lipid nanoparticles (TET-LNPs) for the treatment of glaucoma. In vitro drug release and pre-corneal retention were used to select the optimal chitosan. In vitro drug release curves of three kinds of LNPs showed a sustained release and TMC-TET-LNPs were the slowest. Moreover, compared with CMC-TET-LNPs and HPC-TET-LNPs, TMC-TET-LNPs had longer corneal retention time. Afterwards, the characteristics of TMC-TET-LNPs were investigated. The ocular irritation study revealed no sign of irritation in rabbit eyes. The pharmacokinetic studies showed that the area under the curve of TMC-TET-LNPs was increased by 2.03 times than TET solution (p < 0.01). Furthermore, the drug biofilm interactions were evaluated by molecular dynamics (MD) simulation. In MD simulation, the strong hydrophobic group of TET interacted with the tail of POPC, making it hard to enter the hydrophobic region of the membrane, thereby restricting TET ocular bioavailability. The experiments and MD simulation indicated that TMC-TET-LNPs had great potential for ocular administration and MD simulation could predict transmembrane transport of drugs.