Identification of slow molecular order parameters for Markov model construction

Collective Variables and Facilitated Conformational Opening during Translocation of Human Mitochondrial RNA Polymerase (POLRMT) from Atomic Simulations

Collective variable (CV) identification is challenging in complex dynamical systems. To study the translocation of a single-subunit RNA polymerase (RNAP) during human mitochondrial transcription, we employed all-atom molecular dynamics (MD) as a vehicle to illustrate CV refinement in conformational samplings and dimension reduction analyses. RNAP translocation is an essential mechanical step of transcription elongation that dictates gene expression. The translocation generally follows from polymerization product release and proceeds to initial binding or preinsertion of incoming nucleotides. The human mitochondrial DNA-dependent RNAP (or POLRMT) plays a critical role in cellular metabolism and can be a key molecular off-target in the design of nucleotide analogue antiviral and antitumor drugs due to its structural similarities with many viral RNAPs or RNA-dependent RNA polymerases (RdRps). While POLRMT shares particularly high structural similarity with bacteriophage T7 RNAP, previous experimental studies and our current simulations suggest that POLRMT's mechanochemical coupling mechanisms may be distinct. In the current work, we modeled POLRMT elongation complexes and performed equilibrium MD simulations on the pre- and post-translocation models, with extensive samplings around two potential translocation paths (with or without coupling to the fingers subdomain conformational change). We then compared time-lagged independent component analysis (tICA) and the neural network implementation of the variational approach for Markov processes (VAMPnets) as dimensional reduction methods on selected atomic coordinate sets to best represent the sampled features from the MD simulations. Our results indicate that POLRMT translocation is likely coupled with NTP binding to enable fingers subdomain opening at post-translocation which would otherwise be nonstabilized, or the translocations may proceed futilely without the fingers opening for incoming NTP initial binding or incorporation. The time scale of the coupled translocation reaches over hundreds of microseconds, as predicted by the VAMPnets analyses. Such a time scale seems to match a last postcatalytic kinetic step suggested for the POLRMT elongation cycle by previous experimental measurements. Our MD simulation studies combining atomic coordinate refinements and dimension reduction analyses on top of extensive conformational samplings thus suggest a variation of Brownian ratcheting in POLRMT translocation, as if the Brownian motions of translocation are coupled with NTP binding, which captures transient fingers subdomain opening to couple the translocation with a sustained fingers opening.

Recognition of specific PIP2-subtype composition triggers the allosteric control mechanism for selective membrane targeting of cargo loading and release functions of the intracellular sterol transporter StarD4

We present a comprehensive, quantitative model of the allosteric molecular mechanisms of selective cholesterol (CHL) uptake and delivery by the StarD4 protein - an intracellular cholesterol trafficking protein that facilitates the crucial non-vesicular sterol transport between the plasma membrane and the endoplasmic reticulum. This sterol-specific transfer protein is essential for maintaining the healthy life of human cells. In its physiological function, StarD4 targets both sterol donor and acceptor membranes via interactions with anionic lipids. Experiments have illuminated the kinetics of this sterol transfer and shown it to be modulated by specific phosphatidylinositol phosphates (PIPs) on the target membrane, but the molecular mechanism of the recognition of the PIP2 subtype by StarD4, and how this affects the direction and kinetics of cholesterol transport remained unclear. By revealing a heretofore unrecognized allosteric mechanism that connects the sterol binding site to the part of the protein embedded in the membrane, we show here how StarD4 can respond with different actions to diverse organelle membranes based on their PIP2-subtype composition, in agreement with physiological and experimental evidence. The trajectories of extensive (millisecond range) molecular dynamics (MD) simulation of the StarD4-membrane interactions we calculated, were analyzed with advanced machine learning and information theory methods. Our findings outline how the specific molecular mechanism for recognizing PIP2-subtypes in membranes by StarD4 couples to the defined allosteric pathway that induces the CHL binding pocket to propagate the signal for either uptake or release of the sterol. The central role determined for allostery in these significant advances in the understanding of intracellular cholesterol trafficking by StarD4, aligns with experimentally determined properties of StarD4 function, and interprets them in experimentally testable atomistic terms that explain function-altering results of mutations.

Discovering cryptic pocket opening and binding of a stimulant derivative in a vestibular site of the 5-HT3A receptor

A diverse set of modulators, including stimulants and anesthetics, regulates ion channel function in our nervous system. However, structures of ligand-bound complexes can be difficult to capture by experimental methods, particularly when binding is dynamic. Here, we used computational methods and electrophysiology to identify a possible bound state of a modulatory stimulant derivative in a cryptic vestibular pocket of a mammalian serotonin-3 receptor. We first applied a molecular dynamics simulation-based goal-oriented adaptive sampling method to identify possible open-pocket conformations, followed by Boltzmann docking that combines traditional docking with Markov state modeling. Clustering and analysis of stability and accessibility of docked poses supported a preferred binding site; we further validated this site by mutagenesis and electrophysiology, suggesting a mechanism of potentiation by stabilizing intersubunit contacts. Given the pharmaceutical relevance of serotonin-3 receptors in emesis, psychiatric, and gastrointestinal diseases, characterizing relatively unexplored modulatory sites such as these could open valuable avenues to understanding conformational cycling and designing state-dependent drugs.

Mapping conformational landscape in protein folding: Benchmarking dimensionality reduction and clustering techniques on the Trp-Cage mini-protein

Quantitative characterization of protein conformational landscapes is a computationally challenging task due to their high dimensionality and inherent complexity. In this study, we systematically benchmark several widely used dimensionality reduction and clustering methods to analyze the conformational states of the Trp-Cage mini-protein, a model system with well-documented folding dynamics. Dimensionality reduction techniques, including Principal Component Analysis (PCA), Time-lagged Independent Component Analysis (TICA), and Variational Autoencoders (VAE), were employed to project the high-dimensional free energy landscape onto 2D spaces for visualization. Additionally, clustering methods such as K-means, hierarchical clustering, HDBSCAN, and Gaussian Mixture Models (GMM) were used to identify discrete conformational states directly in the high-dimensional space. Our findings reveal that density-based clustering approaches, particularly HDBSCAN, provide physically meaningful representations of free energy minima. While highlighting the strengths and limitations of each method, our study underscores that no single technique is universally optimal for capturing the complex folding pathways, emphasizing the necessity for careful selection and interpretation of computational methods in biomolecular simulations. These insights will contribute to refining the available tools for analyzing protein conformational landscapes, enabling a deeper understanding of folding mechanisms and intermediate states.

Improved Free-Energy Estimates for the Permeation of Bulky Antibiotic Molecules through Porin Channels Using Temperature-Accelerated Sliced Sampling

The estimation of accurate free energies for antibiotic permeation via the bacterial outer-membrane porins has proven to be challenging. Atomistic simulations of the process suffer from sampling issues that are typical of systems with complex and slow dynamics, even with the application of advanced sampling methods. Ultimately, the objective is to obtain accurate potential of mean force (PMF) for a large set of antibiotics and to predict permeation rates. Therefore, the computational expense becomes an important criterion as well. Simulation studies on the permeation process and similar complex processes have shown that both the sampling scheme employed and the procedure used for the generation of the initial states can critically affect the quality of the estimates obtained and the respective computational overhead. The temperature-accelerated sliced sampling method (TASS) has been shown to partly address the issues with efficient sampling of the important and slow degrees of freedom by enabling simultaneous biasing of a large number of collective variables. In this work, we investigate the effect of the procedure used for the generation of input conformations on the convergence of free-energy estimates obtained from TASS simulations. In particular, we compare the steered molecular dynamics (MD)-based procedure that has been used in previous TASS studies with the Monte Carlo pathway search method, which is used to obtain approximate permeation trajectories with minimum perturbation of the protein channel. We tested different input setups for enrofloxacin permeation through the porins OmpK35 and OmpE35. The best setup shows an improved agreement between independent PMFs in both cases at a much lower computational cost.

Deciphering the folding code of collagens

Collagen proteins contain a characteristic structural motif called a triple helix. During the self-assembly of this motif, three polypeptides form a folding nucleus at the C-termini and then propagate towards the N-termini like a zip-chain. While polypeptides from human collagens contain up to a 1000 amino acids, those found in bacteria can contain up to 6000 amino acids. Additionally, the collagen polypeptides are also frequently interrupted by non-helical sequences that disrupt folding and reduce stability. Given the length of polypeptides and the disruptive interruptions, compensating mechanisms that stabilize against local unfolding during propagation and offset the entropic cost of folding are not fully understood. Here, we show that the information for the correct folding of collagen triple helices is encoded in their sequence as interchain electrostatic interactions, which likely act as molecular clamps that prevent local unfolding. In the case of humans, disrupting these electrostatic interactions is associated with severe to lethal diseases.

From closed to open: three dynamic states of membrane-bound cytochrome P450 3A4

Cytochrome P450 3A4 (CYP3A4) is a membrane bound monooxygenase. It metabolizes the largest proportion of all orally ingested drugs. Ligands can enter and exit the enzyme through flexible tunnels, which co-determine CYP3A4's ligand promiscuity. The flexibility can be represented by distinct conformational states of the enzyme. However, previous state definitions relied solely on crystal structures. We employed conventional molecular dynamics (cMD) simulations to sample the conformational space of CYP3A4. Five conformationally different crystal structures embedded in a membrane were simulated for 1 µs each. A Markov state model (MSM) coupled with spectral clustering (Robust Perron Cluster Analysis PCCA +) resulted in three distinct states: Two open conformations and an intermediate conformation. The tunnels inside CYP3A4 were calculated with CAVER3.0. Notably, we observed variations in bottleneck radii compared to those derived from crystallographic data. We want to point out the importance of simulations to characterize the dynamic behaviour. Moreover, we identified a mechanism, in which the membrane supports the opening of a tunnel. Therefore, CYP3A4 must be investigated in its membrane-bound state.

Deep generative modeling of temperature-dependent structural ensembles of proteins

Deep learning has revolutionized protein structure prediction, but capturing conformational ensembles and structural variability remains an open challenge. While molecular dynamics (MD) is the foundation method for simulating biomolecular dynamics, it is computationally expensive. Recently, deep learning models trained on MD have made progress in generating structural ensembles at reduced cost. However, they remain limited in modeling atomistic details and, crucially, incorporating the effect of environmental factors. Here, we present aSAM (atomistic structural autoencoder model), a latent diffusion model trained on MD to generate heavy atom protein ensembles. Unlike most methods, aSAM models atoms in a latent space, greatly facilitating accurate sampling of side chain and backbone torsion angle distributions. Additionally, we extended aSAM into the first reported transferable generator conditioned on temperature, named aSAMt. Trained on the large and open mdCATH dataset, aSAMt captures temperature-dependent ensemble properties and demonstrates generalization beyond training temperatures. By comparing aSAMt ensembles to long MD simulations of fast folding proteins, we find that high-temperature training enhances the ability of deep generators to explore energy landscapes. Finally, we also show that our MD-based aSAMt can already capture experimentally observed thermal behavior of proteins. Our work is a step towards generalizable ensemble generation to complement physics- based approaches.

Thermodynamic Interpolation: A Generative Approach to Molecular Thermodynamics and Kinetics

Using normalizing flows and reweighting, Boltzmann generators enable equilibrium sampling from a Boltzmann distribution, defined by an energy function and thermodynamic state. In this work, we introduce thermodynamic interpolation (TI), which allows for generating sampling statistics in a temperature-controllable way. We introduce TI flavors that work directly in the ambient configurational space, mapping between different thermodynamic states or through a latent, normally distributed reference state. Our ambient-space approach allows for the specification of arbitrary target temperatures, ensuring generalizability within the temperature range of the training set and demonstrating the potential for extrapolation beyond it. We validate the effectiveness of TI on model systems that exhibit metastability and nontrivial temperature dependencies. Finally, we demonstrate how to combine TI-based sampling to estimate free energy differences through various free energy perturbation methods and provide corresponding approximated kinetic rates, estimated through generator extended dynamic mode decomposition (gEDMD).

Structural Heterogeneity of Intermediate States Facilitates CRIPT Peptide Binding to the PDZ3 Domain: Insights from Molecular Dynamics and Markov State Model Analysis

Intrinsically disordered proteins (IDPs), characterized by a lack of defined tertiary structure, are ubiquitous and indispensable components of cellular machinery. These proteins participate in a diverse array of biological processes, often undergoing conformational transitions upon binding to their target, a phenomenon termed "folding-upon-binding." The finding raises the question of how to achieve rapid binding kinetics in the presence of intrinsic disorder, and the underlying molecular mechanism remains elusive. This study investigated the interaction between the C-terminal region of CRIPT and the third PDZ domain of PSD-95, a critical complex in neuronal development. Upon binding, the CRIPT peptide adopts a β-strand conformation, a process meticulously characterized through extensive molecular dynamics simulations totaling 67.7 μs. Our findings reveal a funnel-like binding landscape in which IDPs can adopt multiple conformations prior to binding, forming structurally heterogeneous intermediate complexes and leading to diverse binding pathways. The stabilization of these intermediate complexes necessitates a dynamic interplay of native and non-native interactions. Markov state model analysis underscores the important role of structural heterogeneity as it contributes to accelerated binding. These findings enrich the classical fly-casting mechanism and provide novel insights into the functional advantages conferred by intrinsic disorder.

Increasing the Efficiency of Ensemble Molecular Dynamics Simulations with Termination of Unproductive Trajectories Identified at Runtime

The application of molecular dynamics (MD) simulations to study increasingly larger and more complex systems is challenged by the required amounts of trajectory data needed to sample their conformational space appropriately. The analysis and interpretation phase of such massive data sets that have to be stored and fed to the various algorithms to reveal the dynamic behaviors of the systems and the underlying energetics in structural terms related to functional mechanisms are also a significant challenge. To develop computational means that can address these challenges, we are developing a software framework that can increase the efficiency of this process. We present one component of this framework that can reduce the size of the accumulating data set while maintaining the structural attributes, distribution, and relative probability ranking of the minima in the free energy map for the system. This framework component utilizes early termination of individual trajectories identified as unproductive in the sampling of conformational space. The criteria for termination are derived quantities such as collective variables (CVs) and secondary quantities calculated from the time series of CVs. They are computed and applied during the trajectory generation. The approach is illustrated with simulations of the FS peptide and evaluated from comparisons between the free energy surfaces calculated from ensembles of complete, unabridged simulations with those obtained from ensembles in which ∼5-50% of trajectories were terminated early. Our early termination approach can optimize computational efficiency while achieving a robust representation of conformational space.

Structural characterization of codon 129 polymorphism in prion peptide segments (PrP127-132) using the Markov State Models

The human prion protein gene (PRNP) consists of two common alleles that encode either methionine or valine residues at codon 129. Polymorphism at codon 129 of the prion protein (PRNP) gene is closely associated with genetic variations and susceptibility to specific variants of prion diseases. The presence of these different alleles, known as the PRNP codon 129 polymorphism, plays a significant role in disease susceptibility and progression. For instance, the prion fragment 127-132 (PrP127-132) has been implicated in the development of variant Creutzfeldt-Jakob disease (vCJD), due to the presence of methionine or valine at codon 129. This study aims to unravel the early structural changes brought by the presence of polymorphism at codon 129. Using molecular dynamics (MD) simulations, we present evidence highlighting a spectrum of structural transitions, uncovering the nuanced conformational heterogeneity governing the polymorphic behavior of the PrP127-132 chain. The Markov state model (MSM) analysis was able to predict several metastable states of these chains and established a kinetic network that describes transitions between these states. Additionally, the MSM analysis showed extra stability of the PrP-M129 polymorph due to less random-coiled motions, the formation of a salt bridge, and an increase in the number of native contacts. The pathogenicity of PrP-V129 can be attributed to enhanced random motion and the absence of a salt bridge.

Acceleration with Interpretability: A Surrogate Model-Based Collective Variable for Enhanced Sampling

Most enhanced sampling methods facilitate the exploration of molecular free energy landscapes by applying a bias potential along a reduced dimensional collective variable (CV) space. The success of these methods depends on the ability of the CVs to follow the relevant slow modes of the system. Intuitive CVs, such as distances or contacts, often prove inadequate, particularly in biological systems involving many coupled degrees of freedom. Machine learning algorithms, especially neural networks (NN), can automate the process of CV discovery by combining a large number of molecular descriptors and often outperform intuitive CVs in sampling efficiency. However, their lack of interpretability and high cost of evaluation during trajectory propagation make NN-CVs difficult to apply to large biomolecular processes. Here, we introduce a surrogate model approach using lasso regression to express the output of a neural network as a linear combination of an automatically chosen subset of the input descriptors. We demonstrate successful applications of our surrogate model CVs in the enhanced sampling simulation of the conformational landscape of alanine dipeptide and chignolin mini-protein. In addition to providing mechanistic insights due to their explainable nature, the surrogate model CVs showed a negligible loss in efficiency and accuracy, compared to the NN-CVs, in reconstructing the underlying free energy surface. Moreover, due to their simplified functional forms, these CVs are better at extrapolating to unseen regions of the conformational space, e.g., saddle points. Surrogate model CVs are also less expensive to evaluate compared to their NN counterparts, making them suitable for enhanced sampling simulation of large and complex biomolecular processes.

Computational Analysis of the Fully Activated Orexin Receptor 2 across Various Thermodynamic Ensembles with Surface Tension Monitoring and Markov State Modeling

In this study, we investigated the stability of the fully activated conformation of the orexin receptor 2 (OX2R) embedded in a pure POPC bilayer using MD simulations. Various thermodynamic ensembles (i.e., NPT, NVT, NVE, NPAT, μVT, and NPγT) were employed to explore the dynamical heterogeneity of the system in a comprehensive way. In addition, informational similarity metrics (e.g., Jensen-Shannon divergence) as well as Markov state modeling approaches were utilized to elucidate the receptor kinetics. Special attention was paid to assessing surface tension within the simulation box, particularly under NPγT conditions, where 21 nominal surface tension constants were evaluated. Our findings suggest that traditional thermodynamic ensembles such as NPT may not adequately control physical properties of the POPC membrane, impacting the plausibility of the OX2R model. In general, the performed study underscores the importance of employing the NPγT ensemble for computational investigations of membrane-embedded receptors, as it effectively maintains zero surface tension in the simulated system. These results offer valuable insights for future research aimed at understanding receptor dynamics and designing targeted therapeutics.

Hidden complexity of α7 nicotinic acetylcholine receptor desensitization revealed by MD simulations and Markov state modeling

The α7 nicotinic acetylcholine receptor is a pentameric ligand-gated ion channel that plays an important role in neuronal signaling throughout the nervous system. Its implication in neurological disorders and inflammation has spurred the development of numerous compounds that enhance channel activation. However, the therapeutic potential of these compounds has been limited by the characteristically fast desensitization of the α7 receptor. Using recent high-resolution structures from cryo-EM, and all-atom molecular dynamic simulations augmented by Markov state modeling, here we explore the mechanism of α7 receptor desensitization and its implication on allosteric modulation. The results provide a precise characterization of the desensitization gate and illuminate the mechanism of ion-pore opening/closing with an agonist bound. In addition, the simulations reveal the existence of a short-lived, open-channel intermediate between the activated and desensitized states that rationalizes the paradoxical pharmacology of the L247T mutant and may be relevant to type-II allosteric modulation. This analysis provides an interpretation of the signal transduction mechanism and its regulation in α7 receptors.

Classification and spatiotemporal correlation of dominant fluctuations in complex dynamical systems

The behaviors of many complex systems, from nanostructured materials to animal colonies, are governed by local events/rearrangements that, while involving a restricted number of interacting units, may generate collective cascade phenomena. Tracking such local events and understanding their emergence and propagation in the system is often challenging. Common strategies consist, for example, in monitoring over time parameters (descriptors) that are designed ad hoc to analyze certain systems. However, such approaches typically require prior knowledge of the system's physics and are poorly transferable. Here, we present a general, transferable, and agnostic analysis approach that can reveal precious information on the physics of a variety of complex dynamical systems starting solely from the trajectories of their constitutive units. Built on a bivariate combination of two abstract descriptors, Local Environments and Neighbors Shuffling and TimeSmooth Overlap of Atomic Position, such approach allows to (i) detect the emergence of local fluctuations in simulation or experimentally acquired trajectories of multibody dynamical systems, (ii) classify fluctuations into categories, and (iii) correlate them in space and time. We demonstrate how this method, based on the abstract concepts of local fluctuations and their spatiotemporal correlations, may reveal precious insights on the emergence and propagation of local and collective phenomena in a variety of complex systems ranging from the atomic- to the macroscopic-scale. This provides a general data-driven approach that we expect will be particularly helpful to study and understand the behavior of systems whose physics is unknown a priori, as well as to revisit a variety of physical phenomena under a new perspective.

Heterogeneous and multiple conformational transition pathways between pseudoknots of the SARS-CoV-2 frameshift element

Frameshifting is an essential mechanism employed by many viruses including coronaviruses to produce viral proteins from a compact RNA genome. It is facilitated by specific RNA folds in the frameshift element (FSE), which has emerged as an important therapeutic target. For SARS-CoV-2, a specific 3-stem pseudoknot has been identified to stimulate frameshifting. However, prior studies and our RNA-As-Graphs analysis coupled to chemical reactivity experiments revealed other folds, including a different pseudoknot. Although structural plasticity has been proposed to play a key role in frameshifting, paths between different FSE RNA folds have not been yet identified. Here, we capture atomic-level transition pathways between two key FSE pseudoknots by transition path sampling coupled to Markov State Modeling and our BOLAS free energy method. We reveal multiple transition paths within a heterogeneous, multihub conformational landscape. A shared folding mechanism involves RNA stem unpairing followed by a 5'-chain end release. Significantly, this pseudoknot transition critically tunes the tension through the RNA spacer region and places the viral RNA in the narrow ribosomal channel. Our work further explains the role of the alternative pseudoknot in ribosomal pausing and clarifies why the experimentally captured pseudoknot is preferred for frameshifting. Our capturing of this large-scale transition of RNA secondary and tertiary structure highlights the complex pathways of biomolecules and the inherent multifarious aspects that viruses developed to ensure virulence and survival. This enhanced understanding of viral frameshifting also provides insights to target key transitions for therapeutic applications. Our methods are generally applicable to other large-scale biomolecular transitions.

Markovian State Models uncover Casein Kinase 1 dynamics that govern circadian period

Circadian rhythms in mammals are tightly regulated through phosphorylation of Period (PER) proteins by Casein Kinase 1 (CK1, subtypes δ and ε). CK1 acts on at least two different regions of PER with opposing effects: phosphorylation of phosphodegron (pD) regions leads to PER degradation, while phosphorylation of the Familial Advanced Sleep Phase (FASP) region leads to PER stabilization. To investigate how substrate selectivity is encoded by the conformational dynamics of CK1, we performed a large set of independent molecular dynamics (MD) simulations of wildtype CK1 and the tau mutant (R178C) that biases kinase activity toward a pD. We used Markovian State Models (MSMs) to integrate the simulations into a single model of the conformational landscape of CK1 and used Gaussian accelerated molecular dynamics (GaMD) to build the first molecular model of CK1 and the unphosphorylated FASP motif. Together, these findings provide a mechanistic view of CK1, establishing how the activation loop acts as a key molecular switch to control substrate selectivity. We show that the tau mutant favors an alternative conformation of the activation loop and significantly accelerates the dynamics of CK1. This reshapes the binding cleft in a way that impairs FASP binding and would ultimately lead to PER destabilization and shorter circadian periods. Finally, we identified an allosteric pocket that could be targeted to bias this molecular switch. Our integrated approach offers a detailed model of CK1's conformational landscape and its relevance to normal, mutant, and druggable circadian timekeeping.

WeTICA: A directed search weighted ensemble based enhanced sampling method to estimate rare event kinetics in a reduced dimensional space

Estimating rare event kinetics from molecular dynamics simulations is a non-trivial task despite the great advances in enhanced sampling methods. Weighted Ensemble (WE) simulation, a special class of enhanced sampling techniques, offers a way to directly calculate kinetic rate constants from biased trajectories without the need to modify the underlying energy landscape using bias potentials. Conventional WE algorithms use different binning schemes to partition the collective variable (CV) space separating the two metastable states of interest. In this work, we have developed a new "binless" WE simulation algorithm to bypass the hurdles of optimizing binning procedures. Our proposed protocol (WeTICA) uses a low-dimensional CV space to drive the WE simulation toward the specified target state. We have applied this new algorithm to recover the unfolding kinetics of three proteins: (A) TC5b Trp-cage mutant, (B) TC10b Trp-cage mutant, and (C) Protein G, with unfolding times spanning the range between 3 and 40 μs using projections along predefined fixed Time-lagged Independent Component Analysis (TICA) eigenvectors as CVs. Calculated unfolding times converge to the reported values with good accuracy with more than one order of magnitude less cumulative WE simulation time than the unfolding time scales with or without a priori knowledge of the CVs that can capture unfolding. Our algorithm can be used with other linear CVs, not limited to TICA. Moreover, the new walker selection criteria for resampling employed in this algorithm can be used on more sophisticated nonlinear CV space for further improvements of binless WE methods.

Flow Matching for Optimal Reaction Coordinates of Biomolecular Systems

We present flow matching for reaction coordinates (FMRC), a novel deep learning algorithm designed to identify optimal reaction coordinates (RC) in biomolecular reversible dynamics. FMRC is based on the mathematical principles of lumpability and decomposability, which we reformulate into a conditional probability framework for efficient data-driven optimization using deep generative models. While FMRC does not explicitly learn the well-established transfer operator or its eigenfunctions, it can effectively encode the dynamics of leading eigenfunctions of the system transfer operator into its low-dimensional RC space. We further quantitatively compare its performance with several state-of-the-art algorithms by evaluating the quality of Markov state models (MSM) constructed in their respective RC spaces, demonstrating the superiority of FMRC in three increasingly complex biomolecular systems. In addition, we successfully demonstrated the efficacy of FMRC for bias deposition in the enhanced sampling of a simple model system. Finally, we discuss its potential applications in downstream applications such as enhanced sampling methods and MSM construction.

Molecular Mechanisms Underlying the Loop-Closing Dynamics of β-1,4 Galactosyltransferase 1

The β-1,4 galactosylation catalyzed by β-1,4 galactosyltransferases (β4Gal-Ts) is not only closely associated with diverse physiological and pathological processes in humans but also widely applied in the N-glycan modification of protein glycoengineering. The loop-closing process of β4Gal-Ts is an essential intermediate step intervening in the binding events of donor substrate (UDP-Gal/Mn2+) and acceptor substrate during its catalytic cycle, with a significant impact on the galactosylation activities. However, the molecular mechanisms in regulating loop-closing dynamics are not entirely clear. Here, we construct Markov state models (MSMs) based on approximately 20 μs of all-atom molecular dynamics simulations to explore the loop-closing dynamics for β-1,4 galactosyltransferase 1 (β4Gal-T1). Our MSM reveals five key metastable states of β4Gal-T1 upon substrate binding, indicating that the entire conformational transition occurs on a time scale of ∼10 μs. Moreover, a regulatory mechanism involving six conserved residues (R187, H190, F222, W310, I341, and D346) among β4Gal-Ts is validated to account for the loop-closing dynamics of the C-loop and W-loop by site-directed mutagenesis and enzymatic activity assays, exhibiting high consistency with our computational predictions. Overall, our research proposes detailed atomic-level insight into the loop-closing dynamics of the C-loop and W-loop on β4Gal-T1, contributing to a deeper understanding of catalytic mechanisms of β-1,4 galactosylation.

Adaptive sampling-based structural prediction reveals opening of a GABAA receptor through the αβ interface

γ-Aminobutyric acid type A (GABAA) receptors are ligand-gated ion channels in the central nervous system with largely inhibitory function. Despite being a target for drugs including general anesthetics and benzodiazepines, experimental structures have yet to capture an open state of classical synaptic α1β2γ2 GABAA receptors. Here, we use a goal-oriented adaptive sampling strategy in molecular dynamics simulations followed by Markov state modeling to capture an energetically stable putative open state of the receptor. The model conducts chloride ions with comparable conductance as in electrophysiology measurements. Relative to experimental structures, our open model is relatively expanded at both the cytoplasmic (-2') and central (9') gates, coordinated with distinctive rearrangements at the transmembrane αβ subunit interface. Consistent with previous experiments, targeted substitutions disrupting interactions at this interface slowed the open-to-desensitized transition rate. This work demonstrates the capacity of advanced simulation techniques to investigate a computationally and experimentally plausible functionally critical of a complex membrane protein yet to be resolved by experimental methods.

Exploring transition states of protein conformational changes via out-of-distribution detection in the hyperspherical latent space

Identifying transitional states is crucial for understanding protein conformational changes that underlie numerous biological processes. Markov state models (MSMs), built from Molecular Dynamics (MD) simulations, capture these dynamics through transitions among metastable conformational states, and have demonstrated success in studying protein conformational changes. However, MSMs face challenges in identifying transition states, as they partition MD conformations into discrete metastable states (or free energy minima), lacking description of transition states located at the free energy barriers. Here, we introduce Transition State identification via Dispersion and vAriational principle Regularized neural networks (TS-DAR), a deep learning framework inspired by out-of-distribution (OOD) detection in trustworthy artificial intelligence (AI). TS-DAR offers an end-to-end pipeline that can simultaneously detect all transition states between multiple free minima from MD simulations using the regularized hyperspherical embeddings in latent space. The key insight of TS-DAR lies in treating transition state structures as OOD data, recognizing that they are sparsely populated and exhibit a distributional shift from metastable states. We demonstrate the power of TS-DAR by applying it to a 2D potential, alanine dipeptide, and the translocation of a DNA motor protein on DNA, where it outperforms previous methods in identifying transition states.

The importance of sampling the dynamical modes: Reevaluating benchmarks for invariant and equivariant features of machine learning potentials for simulation of free energy landscapes

Machine learning interatomic potentials (MLIPs) are rapidly gaining interest for molecular modeling, as they provide a balance between quantum-mechanical level descriptions of atomic interactions and reasonable computational efficiency. However, questions remain regarding the stability of simulations using these potentials, as well as the extent to which the learned potential energy function can be extrapolated safely. Past studies have encountered challenges when MLIPs are applied to classical benchmark systems. In this work, we show that some of these challenges are related to the characteristics of the training datasets, particularly the inefficient exploration of the dynamical modes and the inclusion of rigid constraints. We demonstrate that long stability in simulations with MLIPs can be achieved by generating unconstrained datasets using unbiased classical simulations, provided that the important dynamical modes are correctly sampled. In addition, we emphasize that in order to achieve precise energy predictions, it is important to resort to enhanced sampling techniques for dataset generation, and we demonstrate that safe extrapolation of MLIPs depends on judicious choices related to the system's underlying free energy landscape and the symmetry features embedded within the machine learning models.

Augmenting Human Expertise in Weighted Ensemble Simulations through Deep Learning-Based Information Bottleneck

The weighted ensemble (WE) method stands out as a widely used segment-based sampling technique renowned for its rigorous treatment of kinetics. The WE framework typically involves initially mapping the configuration space onto a low-dimensional collective variable (CV) space and then partitioning it into bins. The efficacy of WE simulations heavily depends on the selection of CVs and binning schemes. The recently proposed state predictive information bottleneck (SPIB) method has emerged as a promising tool for automatically constructing CVs from data and guiding enhanced sampling through an iterative manner. In this work, we advance this data-driven pipeline by incorporating prior expert knowledge. Our hybrid approach combines SPIB-learned CVs to enhance sampling in explored regions with expert-based CVs to guide exploration in regions of interest, synergizing the strengths of both methods. Through benchmarking on alanine dipeptide and chignolin systems, we demonstrate that our hybrid approach effectively guides WE simulations to sample states of interest and reduces run-to-run variances. Moreover, our integration of the SPIB model also enhances the analysis and interpretation of WE simulation data by effectively identifying metastable states and pathways and offering direct visualization of dynamics.

Screening Fast-Mode Motion in Collective Variable Discovery for Biochemical Processes

Collective variables (CVs) describing slow degrees of freedom (DOFs) in biomolecular assemblies are crucial for analyzing molecular dynamics trajectories, creating Markov models and performing CV-based enhanced sampling simulations. While time-lagged independent component analysis (tICA) and its nonlinear successor, time-lagged autoencoder (tAE), are widely used, they often struggle to capture protein dynamics due to interference from random fluctuations along fast DOFs. To address this issue, we propose a novel approach integrating discrete wavelet transform (DWT) with dimensionality reduction techniques. DWT effectively separates fast and slow motion in protein simulation trajectories by decoupling high- and low-frequency signals. Based on the trajectory after filtering out high-frequency signals, which corresponds to fast motion, tICA and tAE can accurately extract CVs representing slow DOFs, providing reliable insights into protein dynamics. Our method demonstrates superior performance in identifying CVs that distinguish metastable states compared to standard tICA and tAE, as validated through analyses of conformational changes of alanine dipeptide and tripeptide and folding of CLN025. Moreover, we show that DWT can be used to improve the performance of a variety of CV-finding algorithms by combining it with Deep-tICA, a cutting-edge CV-finding algorithm, to extract CVs for enhanced-sampling calculations. Given its negligible computational cost and remarkable ability to screen fast motion, we propose DWT as a "free lunch" for CV extraction, applicable to a wide range of CV-finding algorithms.

Kinetics-Based State Definitions for Discrete Binding Conformations of T4 L99A in MD via Markov State Modeling

As a model system, the binding pocket of the L99A mutant of T4 lysozyme has been the subject of numerous computational free energy studies. However, previous studies have failed to fully sample and account for the observed changes in the binding pocket of T4 L99A upon binding of a congeneric ligand series, limiting the accuracy of results. In this work, we resolve the closed, intermediate, and open states for T4 L99A previously reported in experiment in MD and establish definitions for these states based on the dynamics of the system. From this analysis, we arrive at two primary conclusions. First, assignment of simulation trajectories into discrete states should not be done simply based on RMSD to crystal structures as this can result in misassignment of states. Second, the different metastable conformations studied here need to be carefully treated, as we estimate the time scales for conformational interconversion to be on the order of 102 to 103 ns─far longer than time scales for typical binding calculations. We conclude with a discussion on the need to develop enhanced sampling methods to generally account for significant changes in protein conformation due to relatively small ligand perturbations.

Structural insights into the high basal activity and inverse agonism of the orphan receptor GPR6 implicated in Parkinson's disease

GPR6 is an orphan G protein-coupled receptor with high constitutive activity found in D2-type dopamine receptor-expressing medium spiny neurons of the striatopallidal pathway, which is aberrantly hyperactivated in Parkinson's disease. Here, we solved crystal structures of GPR6 without the addition of a ligand (a pseudo-apo state) and in complex with two inverse agonists, including CVN424, which improved motor symptoms in patients with Parkinson's disease in clinical trials. In addition, we obtained a cryo-electron microscopy structure of the signaling complex between GPR6 and its cognate Gs heterotrimer. The pseudo-apo structure revealed a strong density in the orthosteric pocket of GPR6 corresponding to a lipid-like endogenous ligand. A combination of site-directed mutagenesis, native mass spectrometry, and computer modeling suggested potential mechanisms for high constitutive activity and inverse agonism in GPR6 and identified a series of lipids and ions bound to the receptor. The structures and results obtained in this study could guide the rational design of drugs that modulate GPR6 signaling.

Correlating enzymatic reactivity for different substrates using transferable data-driven collective variables

Machine learning (ML) is transforming the investigation of complex biological processes. In enzymatic catalysis, one significant challenge is identifying the reactive conformations (RC) of the enzyme:substrate complex where the substrate assumes a precise arrangement in the active site necessary to initiate a reaction. Traditional methods are hindered by the complexity of the multidimensional free energy landscape involved in the transition from nonreactive to reactive conformations. Here, we applied ML techniques to address this challenge, focusing on human pancreatic α-amylase, a crucial enzyme in type-II diabetes treatment. Using ML-based collective variables (CVs), we correlated the probability of being in a RC with the experimental catalytic activity of several malto-oligosaccharide substrates. Our findings demonstrate a remarkable transferability of these CVs across various compounds, significantly streamlining the modeling process and reducing both computational demand and manual intervention in setting up simulations for new substrates. This approach not only advances our understanding of enzymatic processes but also holds substantial potential for accelerating drug discovery by enabling rapid and accurate evaluation of drug efficacy across different generations of inhibitors.

Tutorial on Molecular Latent Space Simulators (LSSs): Spatially and Temporally Continuous Data-Driven Surrogate Dynamical Models of Molecular Systems

The inherently serial nature and requirement for short integration time steps in the numerical integration of molecular dynamics (MD) calculations place strong limitations on the accessible simulation time scales and statistical uncertainties in sampling slowly relaxing dynamical modes and rare events. Molecular latent space simulators (LSSs) are a data-driven approach to learning a surrogate dynamical model of the molecular system from modest MD training trajectories that can generate synthetic trajectories at a fraction of the computational cost. The training data may comprise single long trajectories or multiple short, discontinuous trajectories collected over, for example, distributed computing resources. Provided the training data provide sufficient sampling of the relevant thermodynamic states and dynamical transitions to robustly learn the underlying microscopic propagator, an LSS furnishes a global model of the dynamics capable of producing temporally and spatially continuous molecular trajectories. Trained LSS models have produced simulation trajectories at up to 6 orders of magnitude lower cost than standard MD to enable dense sampling of molecular phase space and large reduction of the statistical errors in structural, thermodynamic, and kinetic observables. The LSS employs three deep learning architectures to solve three independent learning problems over the training data: (i) an encoding of the high-dimensional MD into a low-dimensional slow latent space using state-free reversible VAMPnets (SRVs), (ii) a propagator of the microscopic dynamics within the low-dimensional latent space using mixture density networks (MDNs), and (iii) a generative decoding of the low-dimensional latent coordinates back to the original high-dimensional molecular configuration space using conditional Wasserstein generative adversarial networks (cWGANs) or denoising diffusion probability models (DDPMs). In this software tutorial, we introduce the mathematical and numerical background and theory of LSS and present example applications of a user-friendly Python package software implementation to alanine dipeptide and a 28-residue beta-beta-alpha (BBA) protein within simple Google Colab notebooks.

AMARO: All Heavy-Atom Transferable Neural Network Potentials of Protein Thermodynamics

All-atom molecular simulations offer detailed insights into macromolecular phenomena, but their substantial computational cost hinders the exploration of complex biological processes. We introduce Advanced Machine-learning Atomic Representation Omni-force-field (AMARO), a new neural network potential (NNP) that combines an O(3)-equivariant message-passing neural network architecture, TensorNet, with a coarse-graining map that excludes hydrogen atoms. AMARO demonstrates the feasibility of training coarser NNP, without prior energy terms, to run stable protein dynamics with scalability and generalization capabilities.

Improved Estimates of Folding Stabilities and Kinetics with Multiensemble Markov Models

Markov State Models (MSMs) have been widely applied to understand protein folding mechanisms by predicting long time scale dynamics from ensembles of short molecular simulations. Most MSM estimators enforce detailed balance, assuming that trajectory data are sampled at an equilibrium. This is rarely the case for ab initio folding studies, however, and as a result, MSMs can severely underestimate protein folding stabilities from such data. To remedy this problem, we have developed an enhanced-sampling protocol in which (1) unbiased folding simulations are performed and sparse tICA is used to obtain features that best capture the slowest events in folding, (2) umbrella sampling along this reaction coordinate is performed to observe folding and unfolding transitions, and (3) the thermodynamics and kinetics of folding are estimated using multiensemble Markov models (MEMMs). Using this protocol, folding pathways, rates, and stabilities of a designed α-helical hairpin, Z34C, can be predicted in good agreement with experimental measurements. These results indicate that accurate simulation-based estimates of absolute folding stabilities are within reach, with implications for the computational design of folded miniproteins and peptidomimetics.

Calculating linear and nonlinear multi-ensemble slow collective variables for protein folding

Traditional molecular dynamics simulation of biomolecules suffers from the conformational sampling problem. It is often difficult to produce enough valid data for post analysis such as free energy calculation and transition path construction. To improve the sampling, one practical solution is putting an adaptive bias potential on some predefined collective variables. The quality of collective variables strongly affects the sampling ability of a molecule in the simulation. In the past, collective variables were built with the sampling data at a constant temperature. This is insufficient because of the same sampling problem. In this work, we apply the standard weighted histogram analysis method to calculate the multi-ensemble averages of pairs of time-lagged features for the construction of both linear and nonlinear slow collective variables. Compared to previous single-ensemble methods, the presented method produces averages with much smaller statistical uncertainties. The generated collective variables help a peptide and a miniprotein fold to their near-native states in a short simulation time period. By using the method, enhanced sampling simulations could be more effective and productive.

Molecular Modeling of Glycosylated Catalytic Domain of Human Carbonic Anhydrase IX

Glycans exhibit significant structural diversity due to the flexibility of glycosidic bonds linking their constituent monosaccharides and the formation of numerous hydrogen bonds. The present work searches a simulated ensemble of glycan chain conformations attached to the catalytic domain of N-glycosylated human carbonic anhydrase IX (HCA IX-c) to identify conformations pointed away or back-folded toward the protein surface guided by different amino acid residues. A series of classical molecular dynamics (MD) simulation studies for a total of 30 μs followed by accelerated MD simulations for a total of 2 μs have been performed using two different force fields to capture varying degrees of fluctuations of both glycan chain and HCA IX. From the underlying free energy profile and kinetics derived using hidden Markov state model, several stable glycan orientations are identified that extend away from the protein surface and convert among each other with rate constants of the order 107-1010 S-1. Most importantly, we have identified a rare glycan conformation which reaches close to a catalytically important amino acid residue, Glu-106. We further enlist the protein residues that couple such less frequent event of the glycan chain back-folding toward the surface of the protein.

Integration of Hydrogen-Deuterium Exchange Mass Spectrometry with Molecular Dynamics Simulations and Ensemble Reweighting Enables High Resolution Protein-Ligand Modeling

Hydrogen-Deuterium exchange mass spectrometry's (HDX-MS) utility in identifying and characterizing protein-small molecule interaction sites has been established. The regions that are seen to be protected from exchange upon ligand binding indicate regions that may be interacting with the ligand, giving a qualitative understanding of the ligand binding pocket. However, quantitatively deriving an accurate high-resolution structure of the protein-ligand complex from the HDX-MS data remains a challenge, often limiting its use in applications such as small molecule drug design. Recent efforts have focused on the development of methods to quantitatively model Hydrogen-Deuterium exchange (HDX) data from computationally modeled structures to garner atomic level insights from peptide-level resolution HDX-MS. One such method, HDX ensemble reweighting (HDXer), employs maximum entropy reweighting of simulated HDX data to experimental HDX-MS to model structural ensembles. In this study, we implement and validate a workflow which quantitatively leverages HDX-MS data to accurately model protein-small molecule ligand interactions. To that end, we employ a strategy combining computational protein-ligand docking, molecular dynamics simulations, HDXer, and dimensional reduction and clustering approaches to extract high-resolution drug binding poses that most accurately conform with HDX-MS data. We apply this workflow to model the interaction of ERK2 and FosA with small molecule compounds and inhibitors they are known to bind. In five out of six of the protein-ligand pairs tested, the HDX derived protein-ligand complexes result in a ligand root-mean-square deviation (RMSD) within 2.5 Å of the known crystal structure ligand.

Adaptive CVgen: Leveraging reinforcement learning for advanced sampling in protein folding and chemical reactions

Enhanced sampling techniques have traditionally encountered two significant challenges: identifying suitable reaction coordinates and addressing the exploration-exploitation dilemma, particularly the difficulty of escaping local energy minima. Here, we introduce Adaptive CVgen, a universal adaptive sampling framework designed to tackle these issues. Our approach utilizes a set of collective variables (CVs) to comprehensively cover the system's potential evolutionary phase space, generating diverse reaction coordinates to address the first challenge. Moreover, we integrate reinforcement learning strategies to dynamically adjust the generated reaction coordinates, thereby effectively balancing the exploration-exploitation dilemma. We apply this framework to sample the conformational space of six proteins transitioning from completely disordered states to folded states, as well as to model the chemical synthesis process of C60, achieving conformations that perfectly match the standard C60 structure. The results demonstrate Adaptive CVgen's effectiveness in exploring new conformations and escaping local minima, achieving both sampling efficiency and exploration accuracy. This framework holds potential for extending to various related challenges, including protein folding dynamics, drug targeting, and complex chemical reactions, thereby opening promising avenues for application in these fields.

On the humanization of VHHs: Prospective case studies, experimental and computational characterization of structural determinants for functionality

The humanization of camelid-derived variable domain heavy chain antibodies (VHHs) poses challenges including immunogenicity, stability, and potential reduction of affinity. Critical to this process are complementarity-determining regions (CDRs), Vernier and Hallmark residues, shaping the three-dimensional fold and influencing VHH structure and function. Additionally, the presence of non-canonical disulfide bonds further contributes to conformational stability and antigen binding. In this study, we systematically humanized two camelid-derived VHHs targeting the natural cytotoxicity receptor NKp30. Key structural positions in Vernier and Hallmark regions were exchanged with residues from the most similar human germline sequences. The resulting variants were characterized for binding affinities, yield, and purity. Structural binding modes were elucidated through crystal structure determination and AlphaFold2 predictions, providing insights into differences in binding affinity. Comparative structural and molecular dynamics characterizations of selected variants were performed to rationalize their functional properties and elucidate the role of specific sequence motifs in antigen binding. Furthermore, systematic analyses of next-generation sequencing (NGS) and Protein Data Bank (PDB) data was conducted, shedding light on the functional significance of Hallmark motifs and non-canonical disulfide bonds in VHHs in general. Overall, this study provides valuable insights into the structural determinants governing the functional properties of VHHs, offering a roadmap for their rational design, humanization, and optimization for therapeutic applications.

Drude2019IDPC polarizable force field reveals structure-function relationship of insulin

Intrinsically disordered proteins (IDPs) lack stable tertiary structures under physiological conditions, yet play key roles in biological processes and associated with human complex diseases. Their conformational characteristics and high content of charged residues make the use of polarizable force fields an advantageous for simulating IDPs. The Drude2019IDP polarizable force field, previously introduced, has demonstrated comprehensive enhancements and improvements in dipeptides, short peptides, and IDPs, achieving a balanced sampling between IDPs and structured proteins. However, the performance in simulating 5 dipeptides was found to be underestimate. Therefore, we individually performed reweighting and grid-based energy correction map (CMAP) optimization for these 5 dipeptides, resulting in the enhanced Drude2019IDPC force field. The performance of Drude2019IDPC was evaluated with 5 dipeptides, 5 disordered short peptides, and a representative IDP. The results demonstrated a marked improvement comparing with original Drude2019IDP. To further substantiate the capabilities of Drude2019IDPC, MD simulation and Markov state model (MSM) were applied to wild type and mutant for insulin, to elucidate the difference of conformational characteristics and transition path. The findings reveal that mutation can maintain the monomorphic characteristics, providing insights for engineered insulin development. These results indicate that Drude2019IDPC could be used to reveal the structure-function relationship for other proteins.

Dynamic Mechanism for Subtype Selectivity of Endocannabinoids

Endocannabinoids are naturally occurring lipid-like molecules that bind to cannabinoid receptors (CB1 and CB2) and regulate many of human bodily functions via the endocannabinoid system. There is a tremendous interest in developing selective drugs that target the CB receptors. However, the biophysical mechanisms responsible for the subtype selectivity for endocannbinoids have not been established. Recent experimental structures of CB receptors show that endocannbinoids potentially bind via membrane using the lipid access channel in the transmembrane region of the receptors. Furthermore, the N-terminus of the receptor could move in and out of the binding pocket thereby modulating both the pocket volume and its residue composition. On the basis of these observations, we propose two hypothesis to explain the selectivity of the endocannabinoid, anandamide for CB1 receptor. First, the selectivity arises from distinct enthalpic ligand-protein interactions along the ligand binding pathway formed due to the movement of N-terminus and subsequent shifts in the binding pocket composition. Second, selectivity arises from the volumetric differences in the binding pocket allowing for differences in ligand conformational entropy. To quantitatively test these hypothesis, we perform extensive molecular dynamics simulations (∼0.9 milliseconds) along with Markov state modeling and deep learning-based VAMP-nets to provide an interpretable characterization of the anandamide binding process to cannabinoid receptors and explain its selectivity for CB1. Our findings reveal that the distinct N-terminus positions along lipid access channels between TM1 and TM7 lead to different binding mechanisms and interactions between anandamide and the binding pocket residues. To validate the critical stabilizing interactions along the binding pathway, relative free energy calculations of anandamide analogs are used. Moreover, the larger CB2 pocket volume increases the entropic effects of ligand binding by allowing higher ligand fluctuations but reduced stable interactions. Therefore, the opposing enthalpy and entropy effects between the receptors shape the endocannabinoid selectivity. Overall, the CB1 selectivity of anandamide is explained by the dominant enthalpy contributions due to ligand-protein interactions in stable binding poses. This study shed lights on potential selectivity mechanisms for endocannabinoids that would aid in the discovery of CB selective drugs.

Non-Markovian Dynamic Models Identify Non-Canonical KRAS-VHL Encounter Complex Conformations for Novel PROTAC Design

Targeted protein degradation (TPD) is emerging as a promising therapeutic approach for cancer and other diseases, with an increasing number of programs demonstrating its efficacy in human clinical trials. One notable method for TPD is Proteolysis Targeting Chimeras (PROTACs) that selectively degrade a protein of interest (POI) through E3-ligase induced ubiquitination followed by proteasomal degradation. PROTACs utilize a warhead-linker-ligand architecture to bring the POI (bound to the warhead) and the E3 ligase (bound to the ligand) into proximity. The resulting non-native protein-protein interactions (PPIs) formed between the POI and E3 ligase lead to the formation of a stable ternary complex, enhancing cooperativity for TPD. A significant challenge in PROTAC design is the screening of the linkers to induce favorable non-native PPIs between POI and E3 ligase. Here, we present a physics-based computational protocol to predict noncanonical and metastable PPI interfaces between an E3 ligase and a given POI, aiding in the design of linkers to stabilize the ternary complex and enhance degradation. Specifically, we build the non-Markovian dynamic model using the Integrative Generalized Master equation (IGME) method from ∼1.5 ms all-atom molecular dynamics simulations of linker-less encounter complex, to systematically explore the inherent PPIs between the oncogene homologue protein and the von Hippel-Lindau E3 ligase. Our protocol revealed six metastable states each containing a different PPI interface. We selected three of these metastable states containing promising PPIs for linker design. Our selection criterion included thermodynamic and kinetic stabilities of PPIs and the accessibility between the solvent-exposed sites on the warheads and E3 ligand. One selected PPIs closely matches a recent cocrystal PPI interface structure induced by an experimentally designed PROTAC with potent degradation efficacy. We anticipate that our protocol has significant potential for widespread application in predicting metastable POI-ligase interfaces that can enable rational design of PROTACs.

Molecular Dynamics Unveils Multiple-Site Binding of Inhibitors with Reduced Activity on the Surface of Dihydrofolate Reductase

The sensitivity to protein inhibitors is altered by modifications or protein mutations, as represented by drug resistance. The mode of stable drug binding to the protein pocket has been experimentally clarified. However, the nature of the binding of inhibitors with reduced sensitivity remains unclear at the atomic level. In this study, we analyzed the thermodynamics and kinetics of inhibitor binding to the surface of wild-type and mutant dihydrofolate reductase (DHFR) using molecular dynamics simulations combined with Markov state modeling. A strong inhibitor of methotrexate (MTX) showed a preference for the active site of wild-type DHFR with minimal binding to unrelated (secondary) sites. Deletion of a side-chain fragment in MTX largely destabilized the active site-bound state, with clear evidence of binding to secondary sites. Similarly, the F31V mutation in DHFR diminished the specificity of MTX binding to the active site. These results reveal the presence of multiple-bound states whose stabilities are comparable to or higher than those of the unbound state, suggesting that a reduction in the binding affinity for the active site significantly elevates the fractions of these states. This study presents a theoretical model that more accurately interprets the altered drug sensitivity than the traditional two-state model.

Symmetry-adapted Markov state models of closing, opening, and desensitizing in α 7 nicotinic acetylcholine receptors

α7 nicotinic acetylcholine receptors (nAChRs) are homopentameric ligand-gated ion channels with critical roles in the nervous system. Recent studies have resolved and functionally annotated closed, open, and desensitized states of these receptors, providing insight into ion permeation and lipid binding. However, the process by which α7 nAChRs transition between states remains unclear. To understand gating and lipid modulation, we generated two ensembles of molecular dynamics simulations of apo α7 nAChRs, with or without cholesterol. Using symmetry-adapted Markov state modeling, we developed a five-state gating model. Free energies recapitulated functional behavior, with the closed state dominating in absence of agonist. Open-to-nonconducting transition rates corresponded to experimental open durations. Cholesterol relatively stabilized the desensitized state, and reduced open-desensitized barriers. These results establish plausible asymmetric transition pathways between states, define lipid modulation effects on the α7 nAChR conformational cycle, and provide an ensemble of structural models applicable to rational design of lipidic pharmaceuticals.

Enhanced Sampling of Biomolecular Slow Conformational Transitions Using Adaptive Sampling and Machine Learning

Biomolecular simulations often suffer from the "time scale problem", hindering the study of rare events occurring over extended time scales. Enhanced sampling techniques aim to alleviate this issue by accelerating conformational transitions, yet they typically necessitate well-defined collective variables (CVs), posing a significant challenge. Machine learning offers promising solutions but typically requires rich training data encompassing the entire free energy surface (FES). In this work, we introduce an automated iterative pipeline designed to mitigate these limitations. Our protocol first utilizes a CV-free count-based adaptive sampling method to generate a data set rich in rare events. From this data set, slow modes are identified using Koopman-reweighted time-lagged independent component analysis (KTICA), which are subsequently leveraged by on-the-fly probability enhanced sampling (OPES) to efficiently explore the FES. The effectiveness of our pipeline is demonstrated and further compared with the common Markov State Model (MSM) approach on two model systems with increasing complexity: alanine dipeptide (Ala2) and deca-alanine (Ala10), underscoring its applicability across diverse biomolecular simulations.

A Mode Evolution Metric to Extract Reaction Coordinates for Biomolecular Conformational Transitions

The complex, multidimensional energy landscape of biomolecules makes the extraction of suitable, nonintuitive collective variables (CVs) that describe their conformational transitions challenging. At present, dimensionality reduction approaches and machine learning (ML) schemes are employed to obtain CVs from molecular dynamics (MD)/Monte Carlo (MC) trajectories or structural databanks for biomolecules. However, minimum sampling conditions to generate reliable CVs that accurately describe the underlying energy landscape remain unclear. Here, we address this issue by developing a Mode evolution Metric (MeM) to extract CVs that can pinpoint new states and describe local transitions in the vicinity of a reference minimum from nonequilibrated MD/MC trajectories. We present a general mathematical formulation of MeM for both statistical dimensionality reduction and machine learning approaches. Application of MeM to MC trajectories of model potential energy landscapes and MD trajectories of solvated alanine dipeptide reveals that the principal components which locate new states in the vicinity of a reference minimum emerge well before the trajectories locally equilibrate between the associated states. Finally, we demonstrate a possible application of MeM in designing efficient biased sampling schemes to construct accurate energy landscape slices that link transitions between states. MeM can help speed up the search for new minima around a biomolecular conformational state and enable the accurate estimation of thermodynamics for states lying on the energy landscape and the description of associated transitions.

Inhibitor binding and disruption of coupled motions in MmpL3 protein: Unraveling the mechanism of trehalose monomycolate transport

Mycobacterial membrane protein Large 3 (MmpL3) of Mycobacterium tuberculosis (Mtb) is crucial for the translocation of trehalose monomycolate (TMM) across the inner bacterial cell membrane, making it a promising target for anti-tuberculosis (TB) drug development. While several structural, microbiological, and in vitro studies have provided significant insights, the precise mechanisms underlying TMM transport by MmpL3 and its inhibition remain incompletely understood at the atomic level. In this study, molecular dynamic (MD) simulations for the apo form and seven inhibitor-bound forms of Mtb MmpL3 were carried out to obtain a thorough comprehension of the protein's dynamics and function. MD simulations revealed that the seven inhibitors in this work stably bind to the central channel of the transmembrane domain and primarily forming hydrogen bonds with ASP251, ASP640, or both residues. Through dynamical cross-correlation matrix and principal component analysis analyses, several types of coupled motions between different domains were observed in the apo state, and distinct conformational states were identified using Markov state model analysis. These coupled motions and varied conformational states likely contribute to the transport of TMM. However, simulations of inhibitor-bound MmpL3 showed an enlargement of the proton channel, potentially disrupting coupled motions. This indicates that inhibitors may impair MmpL3's transport function by directly blocking the proton channel, thereby hindering coordinated domain movements and indirectly affecting TMM translocation.

Cyclopamine modulates smoothened receptor activity in a binding position dependent manner

Cyclopamine, a natural alkaloid, can act as an agonist when it binds to the Cysteine-Rich Domain (CRD) of Smoothened receptor and as an antagonist when it binds to the Transmembrane Domain (TMD). To study the effect of cyclopamine binding to each site experimentally, mutations in the other site are required. Hence, simulations are critical for understanding the WT activity due to binding at different sites. Using multi-milliseconds long aggregate MD simulations combined with Markov state models and machine learning, we explore the dynamic behavior of cyclopamine's interactions with different domains of WT SMO. A higher population of the active state at equilibrium, a lower free energy barrier of ~2 kcal/mol, and expansion of hydrophobic tunnel to facilitate cholesterol transport agrees with cyclopamine's agonistic behavior when bound to CRD. A higher population of the inactive state at equilibrium, a higher free energy barrier of ~4 kcal/mol and restricted hydrophobic tunnel shows cyclopamine's antagonistic behavior when bound to TMD. With cyclopamine bound to both sites, there is a slightly larger inactive population at equilibrium and an increased free energy barrier (~3.5 kcal/mol) exhibiting an overall weak antagonistic effect. These findings show cyclopamine's domain-specific modulation of SMO regulates Hedgehog signaling and cholesterol transport.

Coarse-Graining Conformational Dynamics with Multidimensional Generalized Langevin Equation: How, When, and Why

A data-driven ab initio generalized Langevin equation (AIGLE) approach is developed to learn and simulate high-dimensional, heterogeneous, coarse-grained (CG) conformational dynamics. Constrained by the fluctuation-dissipation theorem, the approach can build CG models in dynamical consistency (DC) with all-atom molecular dynamics. We also propose practical criteria for AIGLE to enforce long-term DC. Case studies of a toy polymer, with 20 CG sites, and the alanine dipeptide, with two dihedral angles, elucidate why one should adopt AIGLE or its Markovian limit for modeling CG conformational dynamics in practice.

The impact of molecular configuration on the bond breaking rates of hydrocarbons: a computational study

The dissociation of hydrocarbon bonds plays a pivotal role in their utilization, whether through fuel combustion or the thermo-cracking of large hydrocarbons in petroleum refinement. Previous studies have primarily focused on the effects of temperature, pressure, and chemical environment on hydrocarbon reactions. However, the influence of molecular configuration on bond breaking rates has not been thoroughly explored. In this study, we propose an approach to compute bond dissociation rates, and apply it to the reactive molecular dynamics simulation (ReaxFF) trajectories of three molecules: n-tridecane, n-pentane, and 1,3-propanediol. Our results reveal that the bond dissociation rate depends not only on the bond position in the chain, but also on the molecular configuration. Stretched configurations exhibit higher dissociation rates, particularly favoring the breaking of central bonds. Conversely, when the molecule is coiled, resulting in a reduced size, terminal bonds exhibit higher dissociation rates. This research contributes to a deeper understanding of molecular dissociation properties in the oxidation of hydrocarbons, and provides a way to explore the bond breaking properties of other molecules.

Optimized Collective Variable for Collapse Transition in Linear Hydrophobic Polymers: Importance of Hydration Water and End-to-End Distance

Choosing an appropriate collective variable (CV) for any biomolecular process is a challenging task. Researchers are developing methods to solve this issue using a variety of methodologies, most recently using machine learning (ML) methods. In this work, we investigate the mechanism of collapse transition across various lengths of polymer systems through adaptively sampled multiple short trajectories utilizing the Time Lagged Independent Component Analysis (TICA) framework. From TICA analysis, it is revealed that the radius of gyration (Rg) and end-to-end distance serve as good order parameters (OPs) for these systems describing overall energy landscapes. Markov state model (MSM) and mean first passage time (MFPT) analysis suggest that hydration water (Nw) plays a determining role in dictating the time scale and barrier for the collapsed transition for the C40 system. P-fold analysis on identifying transition state ensembles (TSE) identified by committor analysis also strengthens the role of Nw in such a transition. TICA, MSM, and committor analyses on the collapse transition for C45 reveal similarities with C40 systems in different aspects. Furthermore, we propose a pipeline integrating XGBoost regression along with an interpretable ML model, Shapley Additive exPlanation (SHAP) to precisely elucidate the contribution of each OP locally at the TSE. Through this approach, we observe that the collapse transition is primarily driven by Nw for both polymer systems. A carefully designed protocol for the collapsed transition of C60 systems indirectly reiterates the above result. Overall, our results suggest that while the end-to-end distance should be considered for better resolution of metastable states in the landscape, Nw is the crucial coordinate to be used in enhanced sampling for the exploration of actual collapse transitions for linear hydrophobic polymer systems. The Python code for analyzing the contribution of different OPs in the TSE using an ML-aided protocol is available on GitHub (https://github.com/saikat-ai/linear_polymer_project).

Electric field modulated configuration and orientation of aqueous molecule chains

Understanding how external electric fields (EFs) impact the properties of aqueous molecules is crucial for various applications in chemistry, biology, and engineering. In this paper, we present a study utilizing molecular dynamics simulation to explore how direct-current (DC) and alternative-current (AC) EFs affect hydrophobic (n-triacontane) and hydrophilic (PEG-10) oligomer chains. Through a machine learning approach, we extract a 2-dimensional free energy (FE) landscape of these molecules, revealing that electric fields modulate the FE landscape to favor stretched configurations and enhance the alignment of the chain with the electric field. Our observations indicate that DC EFs have a more prominent impact on modulation compared to AC EFs and that EFs have a stronger effect on hydrophobic chains than on hydrophilic oligomers. We analyze the orientation of water dipole moments and hydrogen bonds, finding that EFs align water molecules and induce more directional hydrogen bond networks, forming 1D water structures. This favors the stretched configuration and alignment of the studied oligomers simultaneously, as it minimizes the disruption of 1D structures. This research deepens our understanding of the mechanisms by which electric fields modulate molecular properties and could guide the broader application of EFs to control other aqueous molecules, such as proteins or biomolecules.