Kinetic distance and kinetic maps from molecular dynamics simulation

Exploring Conformational Landscapes Along Anharmonic Low-Frequency Vibrations

We aim to automatize the identification of collective variables to simplify and speed up enhanced sampling simulations of conformational dynamics in biomolecules. We focus on anharmonic low-frequency vibrations that exhibit fluctuations on time scales faster than conformational transitions but describe a path of least resistance toward structural change. A key challenge is that harmonic approximations are ill-suited to characterize these vibrations, which are observed at far-infrared frequencies and are easily excited by thermal collisions at room temperature. Here, we approached this problem with a frequency-selective anharmonic (FRESEAN) mode analysis that does not rely on harmonic approximations and successfully isolates anharmonic low-frequency vibrations from short molecular dynamics simulation trajectories. We applied FRESEAN mode analysis to simulations of alanine dipeptide, a common test system for enhanced sampling simulation protocols, and compared the performance of isolated low-frequency vibrations to conventional user-defined collective variables (here backbone dihedral angles) in enhanced sampling simulations. The comparison shows that enhanced sampling along anharmonic low-frequency vibrations not only reproduces known conformational dynamics but can even further improve the sampling of slow transitions compared to user-defined collective variables. Notably, free energy surfaces spanned by low-frequency anharmonic vibrational modes exhibit lower barriers associated with conformational transitions relative to representations in backbone dihedral space. We thus conclude that anharmonic low-frequency vibrations provide a promising path for highly effective and fully automated enhanced sampling simulations of conformational dynamics in biomolecules.

β-sheets mediate the conformational change and allosteric signal transmission between the AsLOV2 termini

Avena sativa phototropin 1 light-oxygen-voltage 2 domain (AsLOV2) is a model protein of Per-Arnt-Sim (PAS) superfamily, characterized by conformational changes in response to external environmental stimuli. This conformational change begins with the unfolding of the N-terminal A'α helix in the dark state followed by the unfolding of the C-terminal Jα helix. The light state is characterized by the unfolded termini and the subsequent modifications in hydrogen bond patterns. In this photoreceptor, β-sheets are identified as crucial components for mediating allosteric signal transmission between the two termini. Through combined experimental and computational investigations, the Hβ and Iβ strands are recognized as the most critical and influential β-sheets in AsLOV2's allosteric mechanism. To elucidate the role of these β-sheets, we introduced 13 distinct mutations (F490L, N492A, L493A, F494L, H495L, L496F, Q497A, R500A, F509L, Q513A, L514A, D515V, and T517V) and conducted comprehensive molecular dynamics simulations. In-depth hydrogen bond analyses emphasized the role of two hydrogen bonds, Asn482-Leu453 and Gln479-Val520, in the observed distinct behaviors of L493A, L496F, Q497A, and D515V mutants. This illustrates the role of β-sheets in the transmission of the allosteric signal upon the photoactivation of the light state.

Information Bottleneck Approach for Markov Model Construction

Markov state models (MSMs) have proven valuable in studying the dynamics of protein conformational changes via statistical analysis of molecular dynamics simulations. In MSMs, the complex configuration space is coarse-grained into conformational states, with dynamics modeled by a series of Markovian transitions among these states at discrete lag times. Constructing the Markovian model at a specific lag time necessitates defining states that circumvent significant internal energy barriers, enabling internal dynamics relaxation within the lag time. This process effectively coarse-grains time and space, integrating out rapid motions within metastable states. Thus, MSMs possess a multiresolution nature, where the granularity of states can be adjusted according to the time-resolution, offering flexibility in capturing system dynamics. This work introduces a continuous embedding approach for molecular conformations using the state predictive information bottleneck (SPIB), a framework that unifies dimensionality reduction and state space partitioning via a continuous, machine learned basis set. Without explicit optimization of the VAMP-based scores, SPIB demonstrates state-of-the-art performance in identifying slow dynamical processes and constructing predictive multiresolution Markovian models. Through applications to well-validated mini-proteins, SPIB showcases unique advantages compared to competing methods. It autonomously and self-consistently adjusts the number of metastable states based on a specified minimal time resolution, eliminating the need for manual tuning. While maintaining efficacy in dynamical properties, SPIB excels in accurately distinguishing metastable states and capturing numerous well-populated macrostates. This contrasts with existing VAMP-based methods, which often emphasize slow dynamics at the expense of incorporating numerous sparsely populated states. Furthermore, SPIB's ability to learn a low-dimensional continuous embedding of the underlying MSMs enhances the interpretation of dynamic pathways. With these benefits, we propose SPIB as an easy-to-implement methodology for end-to-end MSM construction.

Mechanistic insights into G-protein activation via phosphorylation mediated non-canonical pathway

Activation of heterotrimeric G-proteins (Gαβγ) downstream to receptor tyrosine kinases (RTKs) is a well-established crosstalk between the signaling pathways mediated by G-protein coupled receptors (GPCRs) and RTKs. While GPCR serves as a guanine exchange factor (GEF) in the canonical activation of Gα that facilitates the exchange of GDP for GTP, the mechanism through which RTK phosphorylations induce Gα activation remains unclear. Recent experimental studies revealed that the epidermal growth factor receptor (EGFR), a well-known RTK, phosphorylates the helical domain tyrosine residues Y154 and Y155 and accelerates the GDP release from the Gαi3, a subtype of Gα-protein. Using well-tempered metadynamics and extensive unbiased molecular dynamics simulations, we captured the GDP release event and identified the intermediates between bound and unbound states through Markov state models. In addition to weakened salt bridges at the domain interface, phosphorylations induced the unfolding of helix αF, which contributed to increased flexibility near the hinge region, facilitating a greater distance between domains in the phosphorylated Gαi3. Although the larger domain separation in the phosphorylated system provided an unobstructed path for the nucleotide, the accelerated release of GDP was attributed to increased fluctuations in several conserved regions like P-loop, switch 1, and switch 2. Overall, this study provides atomistic insights into the activation of G-proteins induced by RTK phosphorylations and identifies the specific structural motifs involved in the process. The knowledge gained from the study could establish a foundation for targeting non-canonical signaling pathways and developing therapeutic strategies against the ailments associated with dysregulated G-protein signaling.

Uncovering multiscale structure in the variability of larval zebrafish navigation

Animals chain movements into long-lived motor strategies, resulting in variability that ultimately reflects the interplay between internal states and environmental cues. To reveal structure in such variability, we build models that bridges across time scales that enable a quantitative comparison of behavioral phenotypes among individuals. Applied to larval zebrafish exposed to diverse sensory cues, we uncover a hierarchy of long-lived motor strategies, dominated by changes in orientation distinguishing cruising and wandering strategies. Environmental cues induce preferences along these modes at the population level: while fish cruise in the light, they wander in response to aversive (dark) stimuli or in search for prey. Our method enables us to encode the behavioral dynamics of each individual fish in the transitions among coarse-grained motor strategies. By doing so, we uncover a hierarchical structure to the phenotypic variability that corresponds to exploration-exploitation trade-offs. Within a wide range of sensory cues, a major source of variation among fish is driven by prior and immediate exposure to prey that induces exploitation phenotypes. However, a large degree of variability is unexplained by environmental cues, pointing to hidden states that override the sensory context to induce contrasting exploration-exploitation phenotypes. Altogether, our approach extracts the timescales of motor strategies deployed during navigation, exposing undiscovered structure among individuals and pointing to internal states tuned by prior experience.

Unbiasing Enhanced Sampling on a High-Dimensional Free Energy Surface with a Deep Generative Model

Biased enhanced sampling methods that utilize collective variables (CVs) are powerful tools for sampling conformational ensembles. Due to their large intrinsic dimensions, efficiently generating conformational ensembles for complex systems requires enhanced sampling on high-dimensional free energy surfaces. While temperature-accelerated molecular dynamics (TAMD) can trivially adopt many CVs in a simulation, unbiasing the simulation to generate unbiased conformational ensembles requires accurate modeling of a high-dimensional CV probability distribution, which is challenging for traditional density estimation techniques. Here we propose an unbiasing method based on the score-based diffusion model, a deep generative learning method that excels in density estimation across complex data landscapes. We demonstrate that this unbiasing approach, tested on multiple TAMD simulations, significantly outperforms traditional unbiasing methods and can generate accurate unbiased conformational ensembles. With the proposed approach, TAMD can adopt CVs that focus on improving sampling efficiency and the proposed unbiasing method enables accurate evaluation of ensemble averages of important chemical features.

Physics-informed and data-driven discovery of governing equations for complex phenomena in heterogeneous media

Rapid evolution of sensor technology, advances in instrumentation, and progress in devising data-acquisition software and hardware are providing vast amounts of data for various complex phenomena that occur in heterogeneous media, ranging from those in atmospheric environment, to large-scale porous formations, and biological systems. The tremendous increase in the speed of scientific computing has also made it possible to emulate diverse multiscale and multiphysics phenomena that contain elements of stochasticity or heterogeneity, and to generate large volumes of numerical data for them. Thus, given a heterogeneous system with annealed or quenched disorder in which a complex phenomenon occurs, how should one analyze and model the system and phenomenon, explain the data, and make predictions for length and time scales much larger than those over which the data were collected? We divide such systems into three distinct classes. (i) Those for which the governing equations for the physical phenomena of interest, as well as data, are known, but solving the equations over large length scales and long times is very difficult. (ii) Those for which data are available, but the governing equations are only partially known, in the sense that they either contain various coefficients that must be evaluated based on the data, or that the number of degrees of freedom of the system is so large that deriving the complete equations is very difficult, if not impossible, as a result of which one must develop the governing equations with reduced dimensionality. (iii) In the third class are systems for which large amounts of data are available, but the governing equations for the phenomena of interest are not known. Several classes of physics-informed and data-driven approaches for analyzing and modeling of the three classes of systems have been emerging, which are based on machine learning, symbolic regression, the Koopman operator, the Mori-Zwanzig projection operator formulation, sparse identification of nonlinear dynamics, data assimilation combined with a neural network, and stochastic optimization and analysis. This perspective describes such methods and the latest developments in this highly important and rapidly expanding area and discusses possible future directions.

Surface-Accelerated String Method for Locating Minimum Free Energy Paths

We present a surface-accelerated string method (SASM) to efficiently optimize low-dimensional reaction pathways from the sampling performed with expensive quantum mechanical/molecular mechanical (QM/MM) Hamiltonians. The SASM accelerates the convergence of the path using the aggregate sampling obtained from the current and previous string iterations, whereas approaches like the string method in collective variables (SMCV) or the modified string method in collective variables (MSMCV) update the path only from the sampling obtained from the current iteration. Furthermore, the SASM decouples the number of images used to perform sampling from the number of synthetic images used to represent the path. The path is optimized on the current best estimate of the free energy surface obtained from all available sampling, and the proposed set of new simulations is not restricted to being located along the optimized path. Instead, the umbrella potential placement is chosen to extend the range of the free energy surface and improve the quality of the free energy estimates near the path. In this manner, the SASM is shown to improve the exploration for a minimum free energy pathway in regions where the free energy surface is relatively flat. Furthermore, it improves the quality of the free energy profile when the string is discretized with too few images. We compare the SASM, SMCV, and MSMCV using 3 QM/MM applications: a ribozyme methyltransferase reaction using 2 reaction coordinates, the 2'-O-transphosphorylation reaction of Hammerhead ribozyme using 3 reaction coordinates, and a tautomeric reaction in B-DNA using 5 reaction coordinates. We show that SASM converges the paths using roughly 3 times less sampling than the SMCV and MSMCV methods. All three algorithms have been implemented in the FE-ToolKit package made freely available.

Markov State Models: To Optimize or Not to Optimize

Markov state models (MSM) are a popular statistical method for analyzing the conformational dynamics of proteins including protein folding. With all statistical and machine learning (ML) models, choices must be made about the modeling pipeline that cannot be directly learned from the data. These choices, or hyperparameters, are often evaluated by expert judgment or, in the case of MSMs, by maximizing variational scores such as the VAMP-2 score. Modern ML and statistical pipelines often use automatic hyperparameter selection techniques ranging from the simple, choosing the best score from a random selection of hyperparameters, to the complex, optimization via, e.g., Bayesian optimization. In this work, we ask whether it is possible to automatically select MSM models this way by estimating and analyzing over 16,000,000 observations from over 280,000 estimated MSMs. We find that differences in hyperparameters can change the physical interpretation of the optimization objective, making automatic selection difficult. In addition, we find that enforcing conditions of equilibrium in the VAMP scores can result in inconsistent model selection. However, other parameters that specify the VAMP-2 score (lag time and number of relaxation processes scored) have only a negligible influence on model selection. We suggest that model observables and variational scores should be only a guide to model selection and that a full investigation of the MSM properties should be undertaken when selecting hyperparameters.

Microsecond Molecular Dynamics Simulations and Markov State Models of Mutation-Induced Allosteric Mechanisms for the Light-Oxygen-Voltage 2 Protein : Revealing Structural Basis of Signal Transmission Induced by Photoactivation of the Light Protein State

Avena Sativa phototropin 1 Light-oxygen-voltage 2 domain (AsLOV2) is the model protein of Per-Arnt-Sim (PAS) superfamily, characterized by conformational changes in response to external environmental stimuli. This conformational change is initiated by the unfolding of the N-terminal helix in the dark state followed by the unfolding of the C-terminal helix. The light state is defined by the unfolded termini and the subsequent modifications in hydrogen bond patterns. In this photoreceptor, β-sheets have been identified as crucial components for mediating allosteric signal transmission between the two termini. In this study, we combined microsecond all-atm molecular dynamics simulations and Markov state modeling of conformational states to quantify molecular basis of mutation-induced allostery in the AsLOV2 protein. Through a combination of computational investigations, we determine that the Hβ and Iβ strands are the most critical structural elements involved in the allosteric mechanism. To elucidate the role of these β-sheets, we introduced 13 distinct mutations (F490L, N492A, L493A, F494L, H495L, L496F, Q497A, R500A, F509L, Q513A, L514A, D515V, and T517V) and conducted comprehensive simulation analysis. The results highlighted the role of two hydrogen bond Asn482-Leu453 and Gln479-Val520 in the observed distinct behaviors of L493A, L496F, Q497A, and D515V mutants. The comprehensive atomistic-level analysis of the conformational landscapes revealed the critical functional role of β-sheet segments in the transmission of the allosteric signal upon the photoactivation of the light state.

Markov State Model of Solvent Features Reveals Water Dynamics in Protein-Peptide Binding

In this work, we investigate the role of solvent in the binding reaction of the p53 transactivation domain (TAD) peptide to its receptor MDM2. Previously, our group generated 831 μs of explicit-solvent aggregate molecular simulation trajectory data for the MDM2-p53 peptide binding reaction using large-scale distributed computing and subsequently built a Markov State Model (MSM) of the binding reaction (Zhou et al. 2017). Here, we perform a tICA analysis and construct an MSM with similar hyperparameters while using only solvent-based structural features. We find a remarkably similar landscape but accelerated implied timescales for the slowest motions. The solvent shells contributing most to the first tICA eigenvector are those centered on Lys24 and Thr18 of the p53 TAD peptide in the range of 3-6 Å. Important solvent shells were visualized to reveal solvation and desolvation transitions along the peptide-protein binding trajectories. Our results provide a solvent-centric view of the hydrophobic effect in action for a realistic peptide-protein binding scenario.

Free energy along drug-protein binding pathways interactively sampled in virtual reality

We describe a two-step approach for combining interactive molecular dynamics in virtual reality (iMD-VR) with free energy (FE) calculation to explore the dynamics of biological processes at the molecular level. We refer to this combined approach as iMD-VR-FE. Stage one involves using a state-of-the-art 'human-in-the-loop' iMD-VR framework to generate a diverse range of protein-ligand unbinding pathways, benefitting from the sophistication of human spatial and chemical intuition. Stage two involves using the iMD-VR-sampled pathways as initial guesses for defining a path-based reaction coordinate from which we can obtain a corresponding free energy profile using FE methods. To investigate the performance of the method, we apply iMD-VR-FE to investigate the unbinding of a benzamidine ligand from a trypsin protein. The binding free energy calculated using iMD-VR-FE is similar for each pathway, indicating internal consistency. Moreover, the resulting free energy profiles can distinguish energetic differences between pathways corresponding to various protein-ligand conformations (e.g., helping to identify pathways that are more favourable) and enable identification of metastable states along the pathways. The two-step iMD-VR-FE approach offers an intuitive way for researchers to test hypotheses for candidate pathways in biomolecular systems, quickly obtaining both qualitative and quantitative insight.

Frequency-Selective Anharmonic Mode Analysis of Thermally Excited Vibrations in Proteins

Low-frequency molecular vibrations at far-infrared frequencies are thermally excited at room temperature. As a consequence, thermal fluctuations are not limited to the immediate vicinity of local minima on the potential energy surface, and anharmonic properties cannot be ignored. The latter is particularly relevant in molecules with multiple conformations, such as proteins and other biomolecules. However, existing theoretical and computational frameworks for the analysis of molecular vibrations have so far been limited by harmonic or quasi-harmonic approximations, which are ill-suited to describe anharmonic low-frequency vibrations. Here, we introduce a fully anharmonic analysis of molecular vibrations based on a time correlation formalism that eliminates the need for harmonic or quasi-harmonic approximations. We use molecular dynamics simulations of a small protein to demonstrate that this new approach, in contrast to harmonic and quasi-harmonic normal modes, correctly identifies the collective degrees of freedom associated with molecular vibrations at any given frequency. This allows us to unambiguously characterize the anharmonic character of low-frequency vibrations in the far-infrared spectrum.

Uncertainties in Markov State Models of Small Proteins

Markov state models are widely used to describe and analyze protein dynamics based on molecular dynamics simulations, specifically to extract functionally relevant characteristic time scales and motions. Particularly for larger biomolecules such as proteins, however, insufficient sampling is a notorious concern and often the source of large uncertainties that are difficult to quantify. Furthermore, there are several other sources of uncertainty, such as choice of the number of Markov states and lag time, choice and parameters of dimension reduction preprocessing step, and uncertainty due to the limited number of observed transitions; the latter is often estimated via a Bayesian approach. Here, we quantified and ranked all of these uncertainties for four small globular test proteins. We found that the largest uncertainty is due to insufficient sampling and initially increases with the total trajectory length T up to a critical tipping point, after which it decreases as 1 / T , thus providing guidelines for how much sampling is required for given accuracy. We also found that single long trajectories yielded better sampling accuracy than many shorter trajectories starting from the same structure. In comparison, the remaining sources of the above uncertainties are generally smaller by a factor of about 5, rendering them less of a concern but certainly not negligible. Importantly, the Bayes uncertainty, commonly used as the only uncertainty estimate, captures only a relatively small part of the true uncertainty, which is thus often drastically underestimated.

Toward a Benchmark for Markov State Models: The Folding of HP35

Adopting a 300 μs long MD trajectory of the folding of villin headpiece (HP35) by D. E. Shaw Research, we recently constructed a Markov state model (MSM) based on inter-residue contacts. The model reproduces the folding time and predicts that the native basin and unfolded region consist of metastable substates that are structurally well-characterized. Recognizing the need to establish well-defined benchmark problems, we study to what extent and in what sense this MSM can be employed as a reference model. Hence, we test the robustness of the MSM by comparing it to models that use alternative combinations of features, dimensionality reduction methods, and clustering schemes. The study suggests some main characteristics of the folding of HP35 that should be reproduced by other competitive models. Moreover, the discussion reveals which parts of the MSM workflow matter most for the considered problem and illustrates the promises and pitfalls of state-based models for the interpretation of biomolecular simulations.

An Efficient Path Classification Algorithm Based on Variational Autoencoder to Identify Metastable Path Channels for Complex Conformational Changes

Conformational changes (i.e., dynamic transitions between pairs of conformational states) play important roles in many chemical and biological processes. Constructing the Markov state model (MSM) from extensive molecular dynamics (MD) simulations is an effective approach to dissect the mechanism of conformational changes. When combined with transition path theory (TPT), MSM can be applied to elucidate the ensemble of kinetic pathways connecting pairs of conformational states. However, the application of TPT to analyze complex conformational changes often results in a vast number of kinetic pathways with comparable fluxes. This obstacle is particularly pronounced in heterogeneous self-assembly and aggregation processes. The large number of kinetic pathways makes it challenging to comprehend the molecular mechanisms underlying conformational changes of interest. To address this challenge, we have developed a path classification algorithm named latent-space path clustering (LPC) that efficiently lumps parallel kinetic pathways into distinct metastable path channels, making them easier to comprehend. In our algorithm, MD conformations are first projected onto a low-dimensional space containing a small set of collective variables (CVs) by time-structure-based independent component analysis (tICA) with kinetic mapping. Then, MSM and TPT are constructed to obtain the ensemble of pathways, and a deep learning architecture named the variational autoencoder (VAE) is used to learn the spatial distributions of kinetic pathways in the continuous CV space. Based on the trained VAE model, the TPT-generated ensemble of kinetic pathways can be embedded into a latent space, where the classification becomes clear. We show that LPC can efficiently and accurately identify the metastable path channels in three systems: a 2D potential, the aggregation of two hydrophobic particles in water, and the folding of the Fip35 WW domain. Using the 2D potential, we further demonstrate that our LPC algorithm outperforms the previous path-lumping algorithms by making substantially fewer incorrect assignments of individual pathways to four path channels. We expect that LPC can be widely applied to identify the dominant kinetic pathways underlying complex conformational changes.

Sequence controlled secondary structure is important for the site-selectivity of lanthipeptide cyclization

Lanthipeptides are ribosomally synthesized and post-translationally modified peptides that are generated from precursor peptides through a dehydration and cyclization process. ProcM, a class II lanthipeptide synthetase, demonstrates high substrate tolerance. It is enigmatic that a single enzyme can catalyze the cyclization process of many substrates with high fidelity. Previous studies suggested that the site-selectivity of lanthionine formation is determined by substrate sequence rather than by the enzyme. However, exactly how substrate sequence contributes to site-selective lanthipeptide biosynthesis is not clear. In this study, we performed molecular dynamic simulations for ProcA3.3 variants to explore how the predicted solution structure of the substrate without enzyme correlates to the final product formation. Our simulation results support a model in which the secondary structure of the core peptide is important for the final product's ring pattern for the substrates investigated. We also demonstrate that the dehydration step in the biosynthesis pathway does not influence the site-selectivity of ring formation. In addition, we performed simulation for ProcA1.1 and 2.8, which are well-suited candidates to investigate the connection between order of ring formation and solution structure. Simulation results indicate that in both cases, C-terminal ring formation is more likely which was supported by experimental results. Our findings indicate that the substrate sequence and its solution structure can be used to predict the site-selectivity and order of ring formation, and that secondary structure is a crucial factor influencing the site-selectivity. Taken together, these findings will facilitate our understanding of the lanthipeptide biosynthetic mechanism and accelerate bioengineering efforts for lanthipeptide-derived products.

Discovering Reaction Pathways, Slow Variables, and Committor Probabilities with Machine Learning

A significant challenge faced by atomistic simulations is the difficulty, and often impossibility, to sample the transitions between metastable states of the free-energy landscape associated with slow molecular processes. Importance-sampling schemes represent an appealing option to accelerate the underlying dynamics by smoothing out the relevant free-energy barriers, but require the definition of suitable reaction-coordinate (RC) models expressed in terms of compact low-dimensional sets of collective variables (CVs). While most computational studies of slow molecular processes have traditionally relied on educated guesses based on human intuition to reduce the dimensionality of the problem at hand, a variety of machine-learning (ML) algorithms have recently emerged as powerful alternatives to discover meaningful CVs capable of capturing the dynamics of the slowest degrees of freedom. Considering a simple paradigmatic situation in which the long-time dynamics is dominated by the transition between two known metastable states, we compare two variational data-driven ML methods based on Siamese neural networks aimed at discovering a meaningful RC model─the slowest decorrelating CV of the molecular process, and the committor probability to first reach one of the two metastable states. One method is the state-free reversible variational approach for Markov processes networks (VAMPnets), or SRVs─the other, inspired by the transition path theory framework, is the variational committor-based neural networks, or VCNs. The relationship and the ability of these methodologies to discover the relevant descriptors of the slow molecular process of interest are illustrated with a series of simple model systems. We also show that both strategies are amenable to importance-sampling schemes through an appropriate reweighting algorithm that approximates the kinetic properties of the transition.

Identification of Catechins' Binding Sites in Monomeric Aβ42 through Ensemble Docking and MD Simulations

The assembly of the amyloid-β peptide (Aβ) into toxic oligomers and fibrils is associated with Alzheimer's disease and dementia. Therefore, disrupting amyloid assembly by direct targeting of the Aβ monomeric form with small molecules or antibodies is a promising therapeutic strategy. However, given the dynamic nature of Aβ, standard computational tools cannot be easily applied for high-throughput structure-based virtual screening in drug discovery projects. In the current study, we propose a computational pipeline-in the framework of the ensemble docking strategy-to identify catechins' binding sites in monomeric Aβ₄₂. It is shown that both hydrophobic aromatic interactions and hydrogen bonding are crucial for the binding of catechins to Aβ₄₂. Additionally, it has been found that all the studied ligands, especially EGCG, can act as potent inhibitors against amyloid aggregation by blocking the central hydrophobic region of Aβ. Our findings are evaluated and confirmed with multi-microsecond MD simulations. Finally, it is suggested that our proposed pipeline, with low computational cost in comparison with MD simulations, is a suitable approach for the virtual screening of ligand libraries against Aβ.

Girsanov Reweighting Enhanced Sampling Technique (GREST): On-the-Fly Data-Driven Discovery of and Enhanced Sampling in Slow Collective Variables

Molecular dynamics simulations of microscopic phenomena are limited by the short integration time steps which are required for numerical stability but which limit the practically achievable simulation time scales. Collective variable (CV) enhanced sampling techniques apply biases to predefined collective coordinates to promote barrier crossing, phase space exploration, and sampling of rare events. The efficacy of these techniques is contingent on the selection of good CVs correlated with the molecular motions governing the long-time dynamical evolution of the system. In this work, we introduce Girsanov Reweighting Enhanced Sampling Technique (GREST) as an adaptive sampling scheme that interleaves rounds of data-driven slow CV discovery and enhanced sampling along these coordinates. Since slow CVs are inherently dynamical quantities, a key ingredient in our approach is the use of both thermodynamic and dynamical Girsanov reweighting corrections for rigorous estimation of slow CVs from biased simulation data. We demonstrate our approach on a toy 1D 4-well potential, a simple biomolecular system alanine dipeptide, and the Trp-Leu-Ala-Leu-Leu (WLALL) pentapeptide. In each case GREST learns appropriate slow CVs and drives sampling of all thermally accessible metastable states starting from zero prior knowledge of the system. We make GREST accessible to the community via a publicly available open source Python package.

Mechanism of β-hairpin formation in AzoChignolin and Chignolin

AzoChignolin is a photoswitchable variant of the mini-protein Chignolin with an azobenzene (AMPP) replacing the central loop. AzoChignolin is unfolded with AMPP in the trans-isomer. Transition to the cis-isomer causes β-hairpin folding similar to Chignolin. The AzoChignolin system is excellently suited for comprehensive analysis of folding nucleation kinetics. Utilizing multiple long-time MD simulations of AzoChignolin and Chignolin in MeOH and water, we estimated Markov models to examine folding kinetics of both peptides. We show that while AzoChignolin mimics Chignolin's structure well, the folding kinetics are quite different. Not only folding times but also intermediate states differ, particularly Chignolin is able to fold in MeOH into an α-helical intermediate which is impossible to form in AzoChignolin. The Markov models demonstrate that AzoChignolin's kinetics are generally faster, specifically when comparing the two main microfolding processes of hydrophobic collapse and turn formation. Photoswitchable loops are used frequently to understand the kinetics of elementary protein folding nucleation. However, our results indicate that intermediates and folding kinetics may differ between natural loops and photoswitchable variants.

A variational approach to assess reaction coordinates for two-step crystallization

Molecule- and particle-based simulations provide the tools to test, in microscopic detail, the validity of classical nucleation theory. In this endeavor, determining nucleation mechanisms and rates for phase separation requires an appropriately defined reaction coordinate to describe the transformation of an out-of-equilibrium parent phase for which myriad options are available to the simulator. In this article, we describe the application of the variational approach to Markov processes to quantify the suitability of reaction coordinates to study crystallization from supersaturated colloid suspensions. Our analysis indicates that collective variables (CVs) that correlate with the number of particles in the condensed phase, the system potential energy, and approximate configurational entropy often feature as the most appropriate order parameters to quantitatively describe the crystallization process. We apply time-lagged independent component analysis to reduce high-dimensional reaction coordinates constructed from these CVs to build Markov State Models (MSMs), which indicate that two barriers separate a supersaturated fluid phase from crystals in the simulated environment. The MSMs provide consistent estimates for crystal nucleation rates, regardless of the dimensionality of the order parameter space adopted; however, the two-step mechanism is only consistently evident from spectral clustering of the MSMs in higher dimensions. As the method is general and easily transferable, the variational approach we adopt could provide a useful framework to study controls for crystal nucleation.

Chasing collective variables using temporal data-driven strategies

The convergence of free-energy calculations based on importance sampling depends heavily on the choice of collective variables (CVs), which in principle, should include the slow degrees of freedom of the biological processes to be investigated. Autoencoders (AEs), as emerging data-driven dimension reduction tools, have been utilised for discovering CVs. AEs, however, are often treated as black boxes, and what AEs actually encode during training, and whether the latent variables from encoders are suitable as CVs for further free-energy calculations remains unknown. In this contribution, we review AEs and their time-series-based variants, including time-lagged AEs (TAEs) and modified TAEs, as well as the closely related model variational approach for Markov processes networks (VAMPnets). We then show through numerical examples that AEs learn the high-variance modes instead of the slow modes. In stark contrast, time series-based models are able to capture the slow modes. Moreover, both modified TAEs with extensions from slow feature analysis and the state-free reversible VAMPnets (SRVs) can yield orthogonal multidimensional CVs. As an illustration, we employ SRVs to discover the CVs of the isomerizations of N-acetyl-N'-methylalanylamide and trialanine by iterative learning with trajectories from biased simulations. Last, through numerical experiments with anisotropic diffusion, we investigate the potential relationship of time-series-based models and committor probabilities.

Dynamic conformational states of apo, ATP and cabozantinib bound TAM kinases to differentiate active-inactive kinetic models

The dynamically active and inactive conformations of kinases play a crucial role in the activation of intracellular downstream signaling pathways. The all-atom molecular dynamics (MD) simulations at microsecond (µs) timescale and longer provide robust insights into the structural details of conformational alterations in kinases that contribute to their cellular metabolic activities and signaling pathways. Tyro3, Axl and Mer (TAM) receptor tyrosine kinases (RTKs) are overexpressed in several types of human cancers. Cabozantinib, a small molecule inhibitor constrains the activity of TAM kinases at nanomolar concentrations. The apo, complexes of ATP (active state) and cabozantinib (active and inactive states) with TAM RTKs were studied by 1 µs MD simulations followed by trajectory analyses. The dynamic mechanistic pathways intrinsic to the kinase activity and protein conformational landscape in the cabozantinib bound TAM kinases are revealed due to the alterations in the P-loop, α-helix and activation loop that result in breaking the regulatory (R) and catalytic (C) spines, while the active states of ATP bound TAM kinases are retained. The co-existence of dynamical states when bound to cabozantinib was observed and the long-lived kinetic transition states of distinct active and inactive structural models were deciphered from MD simulation trajectories that have not been revealed so far.Communicated by Ramaswamy H. Sarma.

Reweighted Manifold Learning of Collective Variables from Enhanced Sampling Simulations

Enhanced sampling methods are indispensable in computational chemistry and physics, where atomistic simulations cannot exhaustively sample the high-dimensional configuration space of dynamical systems due to the sampling problem. A class of such enhanced sampling methods works by identifying a few slow degrees of freedom, termed collective variables (CVs), and enhancing the sampling along these CVs. Selecting CVs to analyze and drive the sampling is not trivial and often relies on chemical intuition. Despite routinely circumventing this issue using manifold learning to estimate CVs directly from standard simulations, such methods cannot provide mappings to a low-dimensional manifold from enhanced sampling simulations, as the geometry and density of the learned manifold are biased. Here, we address this crucial issue and provide a general reweighting framework based on anisotropic diffusion maps for manifold learning that takes into account that the learning data set is sampled from a biased probability distribution. We consider manifold learning methods based on constructing a Markov chain describing transition probabilities between high-dimensional samples. We show that our framework reverts the biasing effect, yielding CVs that correctly describe the equilibrium density. This advancement enables the construction of low-dimensional CVs using manifold learning directly from the data generated by enhanced sampling simulations. We call our framework reweighted manifold learning. We show that it can be used in many manifold learning techniques on data from both standard and enhanced sampling simulations.

Temporally Coherent Backmapping of Molecular Trajectories From Coarse-Grained to Atomistic Resolution

Coarse-graining offers a means to extend the achievable time and length scales of molecular dynamics simulations beyond what is practically possible in the atomistic regime. Sampling molecular configurations of interest can be done efficiently using coarse-grained simulations, from which meaningful physicochemical information can be inferred if the corresponding all-atom configurations are reconstructed. However, this procedure of backmapping to reintroduce the lost atomistic detail into coarse-grain structures has proven a challenging task due to the many feasible atomistic configurations that can be associated with one coarse-grain structure. Existing backmapping methods are strictly frame-based, relying on either heuristics to replace coarse-grain particles with atomic fragments and subsequent relaxation or parametrized models to propose atomic coordinates separately and independently for each coarse-grain structure. These approaches neglect information from previous trajectory frames that is critical to ensuring temporal coherence of the backmapped trajectory, while also offering information potentially helpful to producing higher-fidelity atomic reconstructions. In this work, we present a deep learning-enabled data-driven approach for temporally coherent backmapping that explicitly incorporates information from preceding trajectory structures. Our method trains a conditional variational autoencoder to nondeterministically reconstruct atomistic detail conditioned on both the target coarse-grain configuration and the previously reconstructed atomistic configuration. We demonstrate our backmapping approach on two exemplar biomolecular systems: alanine dipeptide and the miniprotein chignolin. We show that our backmapped trajectories accurately recover the structural, thermodynamic, and kinetic properties of the atomistic trajectory data.

Simulation of ligand dissociation kinetics from the protein kinase PYK2

Early-stage drug discovery projects often focus on equilibrium binding affinity to the target alongside selectivity and other pharmaceutical properties. The kinetics of drug binding are ignored but can have significant influence on drug efficacy. Therefore, increasing attention has been paid on evaluating drug-binding kinetics early in a drug discovery process. Simulating drug-binding kinetics at the atomic level is challenging for the long time scale involved. Here, we used the transition-based reweighting analysis method (TRAM) with the Markov state model to study the dissociation of a ligand from the protein kinase PYK2. TRAM combines biased and unbiased simulations to reduce computational costs. This work used the umbrella sampling technique for the biased simulations. Although using the potential of mean force from umbrella sampling simulations with the transition-state theory over-estimated the dissociation rate by three orders of magnitude, TRAM gave a dissociation rate within an order of magnitude of the experimental value.

Conformational transitions in BTG1 antiproliferative protein and their modulation by disease mutants

B cell translocation gene 1 (BTG1) protein belongs to the BTG/transducer of ERBB2 (TOB) family of antiproliferative proteins whose members regulate various key cellular processes such as cell cycle progression, apoptosis, and differentiation. Somatic missense mutations in BTG1 are found in ∼70% of a particularly malignant and disseminated subtype of diffuse large B cell lymphoma (DLBCL). Antiproliferative activity of BTG1 has been linked to its ability to associate with transcriptional cofactors and various enzymes. However, molecular mechanisms underlying these functional interactions and how the disease-linked mutations in BTG1 affect these mechanisms are currently unknown. To start filling these knowledge gaps, here, using atomistic molecular dynamics (MD) simulations, we explored structural, dynamic, and kinetic characteristics of BTG1 protein, and studied how various DLBCL mutations affect these characteristics. We focused on the protein region formed by α2 and α4 helices, as this interface has been reported not only to serve as a binding hotspot for several cellular partners but also to harbor sites for the majority of known DLBCL mutations. Markov state modeling analysis of extensive MD simulations revealed that the α2-α4 interface in the wild-type (WT) BTG1 undergoes conformational transitions between closed and open metastable states. Importantly, we show that some of the mutations in this region that are observed in DLBCL, such as Q36H, F40C, Q45P, E50K (in α2), and A83T and A84E (in α4), either overstabilize one of these two metastable states or give rise to new conformations in which these helices are distorted (i.e., kinked or unfolded). Based on these results, we conclude that the rapid interconversion between the closed and open conformations of the α2-α4 interface is an essential component of the BTG1 functional dynamics that can prime the protein for functional associations with its binding partners. Disruption of the native dynamic equilibrium by DLBCL mutants leads to the ensemble of conformations in BTG1 that are unlikely structurally and/or kinetically to enable productive functional interactions with the binding proteins.

On Sampling Minimum Energy Path

Sampling the minimum energy path (MEP) between two minima of a system is often hindered by the presence of an energy barrier separating the two metastable states. As a consequence, direct sampling based on molecular dynamics or Markov Chain Monte Carlo methods becomes inefficient, the crossing of the energy barrier being associated to a rare event. Augmented sampling methods based on the definition of collective variables or reaction coordinates allow us to circumvent this limitation at the price of an arbitrary choice of the dimensionality reduction algorithm. We couple the statistical sampling techniques, namely, metadynamics and invertible neural networks, with autoencoders so as to gradually learn the MEP and the collective variable at the same time. Learning is achieved through a succession of two steps: statistical sampling of the most probable path between the two minima and redefinition of the collective variable from the updated data points. The prototypical Mueller potential with nearly orthogonal minima is employed to demonstrate the ability of such coupling to unravel a complex MEP.

Investigating the Bioactive Conformation of Angiotensin II Using Markov State Modeling Revisited with Web-Scale Clustering

Molecular dynamics simulation is a powerful technique for studying the structure and dynamics of biomolecules in atomic-level detail by sampling their various conformations in real time. Because of the long timescales that need to be sampled to study biomolecular processes and the big and complex nature of the corresponding data, relevant analyses of important biophysical phenomena are challenging. Clustering and Markov state models (MSMs) are efficient computational techniques that can be used to extract dominant conformational states and to connect those with kinetic information. In this work, we perform Molecular Dynamics simulations to investigate the free energy landscape of Angiotensin II (AngII) in order to unravel its bioactive conformations using different clustering techniques and Markov state modeling. AngII is an octapeptide hormone, which binds to the AT1 transmembrane receptor, and plays a vital role in the regulation of blood pressure, conservation of total blood volume, and salt homeostasis. To mimic the water-membrane interface as AngII approaches the AT1 receptor and to compare our findings with available experimental results, the simulations were performed in water as well as in water-ethanol mixtures. Our results show that in the water-ethanol environment, AngII adopts more compact U-shaped (folded) conformations than in water, which resembles its structure when bound to the AT1 receptor. For clustering of the conformations, we validate the efficiency of an inverted-quantized k-means algorithm, as a fast approximate clustering technique for web-scale data (millions of points into thousands or millions of clusters) compared to k-means, on data from trajectories of molecular dynamics simulations with reasonable trade-offs between time and accuracy. Finally, we extract MSMs using various clustering techniques for the generation of microstates and macrostates and for the selection of the macrostate representatives.

Deeptime: a Python library for machine learning dynamical models from time series data

Generation and analysis of time-series data is relevant to many quantitative fields ranging from economics to fluid mechanics. In the physical sciences, structures such as metastable and coherent sets, slow relaxation processes, collective variables, dominant transition pathways or manifolds and channels of probability flow can be of great importance for understanding and characterizing the kinetic, thermodynamic and mechanistic properties of the system. Deeptime is a general purpose Python library offering various tools to estimate dynamical models based on time-series data including conventional linear learning methods, such as Markov state models (MSMs), Hidden Markov Models and Koopman models, as well as kernel and deep learning approaches such as VAMPnets and deep MSMs. The library is largely compatible with scikit-learn, having a range of Estimator classes for these different models, but in contrast to scikit-learn also provides deep Model classes, e.g. in the case of an MSM, which provide a multitude of analysis methods to compute interesting thermodynamic, kinetic and dynamical quantities, such as free energies, relaxation times and transition paths. The library is designed for ease of use but also easily maintainable and extensible code. In this paper we introduce the main features and structure of the deeptime software. Deeptime can be found under https://deeptime-ml.github.io/.

Distinct Binding Mechanisms for Allosteric Sodium Ion in Cannabinoid Receptors

The therapeutic potential of cannabinoid receptors is not fully explored due to psychoactive side effects and lack of selectivity associated with orthosteric ligands. Allosteric modulators have the potential to become selective therapeutics for cannabinoid receptors. Biochemical experiments have shown the effects of the allosteric Na⁺ binding on cannabinoid receptor activity. However, the Na⁺ coordination site and binding pathway are still unknown. Here, we perform molecular dynamic simulations to explore Na⁺ binding in the cannabinoid receptors, CB₁ and CB₂. Simulations reveal that Na⁺ binds to the primary binding site from different extracellular sites for CB₁ and CB₂. A distinct secondary Na⁺ coordination site is identified in CB₁ that is not present in CB₂. Furthermore, simulations also show that intracellular Na⁺ could bind to the Na⁺ binding site in CB₁. Constructed Markov state models show that the standard free energy of Na⁺ binding is similar to the previously calculated free energy for other class A GPCRs.

Mechanistic Origin of Partial Agonism of Tetrahydrocannabinol for Cannabinoid Receptors

Cannabinoid receptor 1 (CB₁) is a therapeutically relevant drug target for controlling pain, obesity, and other central nervous system disorders. However, full agonists and antagonists of CB₁ have been reported to cause serious side effects in patients. Therefore, partial agonists have emerged as a viable alternative as they can mitigate overstimulation and side effects. One of the key bottlenecks in the design of partial agonists, however, is the lack of understanding of the molecular mechanism of partial agonism itself. In this study, we examine two mechanistic hypotheses for the origin of partial agonism in cannabinoid receptors and predict the mechanistic basis of partial agonism exhibited by Δ⁹-Tetrahydrocannabinol (THC) against CB₁. In particular, we inspect whether partial agonism emerges from the ability of THC to bind in both agonist and antagonist-binding poses or from its ability to only partially activate the receptor. We used extensive molecular dynamics simulations and Markov state modeling to capture the THC binding in both antagonist and agonist-binding poses in the CB₁ receptor. Furthermore, we predict that binding of THC in the agonist-binding pose leads to rotation of toggle switch residues and causes partial outward movement of intracellular transmembrane helix 6 (TM6). Our simulations also suggest that the alkyl side chain of THC plays a crucial role in determining partial agonism by stabilizing the ligand in the agonist and antagonist-like poses within the pocket. Taken together, this study provides important insights into the mechanistic origin of the partial agonism of THC.

Collective variable discovery in the age of machine learning: reality, hype and everything in between

Understanding the kinetics and thermodynamics profile of biomolecules is necessary to understand their functional roles which has a major impact in mechanism driven drug discovery. Molecular dynamics simulation has been routinely used to understand conformational dynamics and molecular recognition in biomolecules. Statistical analysis of high-dimensional spatiotemporal data generated from molecular dynamics simulation requires identification of a few low-dimensional variables which can describe the essential dynamics of a system without significant loss of information. In physical chemistry, these low-dimensional variables are often called collective variables. Collective variables are used to generate reduced representations of free energy surfaces and calculate transition probabilities between different metastable basins. However the choice of collective variables is not trivial for complex systems. Collective variables range from geometric criteria such as distances and dihedral angles to abstract ones such as weighted linear combinations of multiple geometric variables. The advent of machine learning algorithms led to increasing use of abstract collective variables to represent biomolecular dynamics. In this review, I will highlight several nuances of commonly used collective variables ranging from geometric to abstract ones. Further, I will put forward some cases where machine learning based collective variables were used to describe simple systems which in principle could have been described by geometric ones. Finally, I will put forward my thoughts on artificial general intelligence and how it can be used to discover and predict collective variables from spatiotemporal data generated by molecular dynamics simulations.

Data driven collective variable discovery methods to capture conformational dynamics in biological macromolecules.

Molecular dynamics analysis of the structural properties of the transglutaminases of Kutzneria albida and Streptomyces mobaraensis

The microbial transglutaminase (TGase) from Streptomyces mobaraensis (MTGase) is widely used for industrial applications. However, in the last decades, TGases from other bacteria have been described. We focused our attention on TGase, from Kutzneria albida (KalbTGase), recently characterized as more selective than MTGase and proposed for applications in drug delivery. By comparison of the crystallographic structures, the volume of the catalytic site results smaller in KalbTGase. We compared KalbTGase and MTGase structural flexibility by molecular dynamics (MD) simulations at different conditions. KalbTGase is more rigid than MTGase at 300 K, but the catalytic site has a preserved conformation in both structures. Preliminary studies at higher temperatures suggest that KalbTGase acquires enhanced conformational flexibility far from the active site region. The volume of the catalytic active site pocket of KalbTGase at room temperature is smaller than that of MTGase, and decreases at 335 K, remaining stable after further temperature increase. On the contrary, in MTGase the pocket volume continues to decrease as the temperature increases. Overall, the results of our study suggest that at room temperature the enhanced specificity of KalbTGase could be related to a more closed catalytic pocket and lower flexibility than MTGase. Moreover, by preliminary results at higher temperature, KalbTGase structural flexibility suggests an adaptability to different substrates not recognized at room temperature. Lower adaptability of MTGase at higher temperature with a reduction of the catalytic pocket, instead, suggests a reduction of its activity.

Identifying the leading dynamics of ubiquitin: A comparison between the tICA and the LE4PD slow fluctuations in amino acids' position

Molecular Dynamics (MD) simulations of proteins implicitly contain the information connecting the atomistic molecular structure and proteins' biologically relevant motion, where large-scale fluctuations are deemed to guide folding and function. In the complex multiscale processes described by MD trajectories, it is difficult to identify, separate, and study those large-scale fluctuations. This problem can be formulated as the need to identify a small number of collective variables that guide the slow kinetic processes. The most promising method among the ones used to study the slow leading processes in proteins' dynamics is the time-structure based on time-lagged independent component analysis (tICA), which identifies the dominant components in a noisy signal. Recently, we developed an anisotropic Langevin approach for the dynamics of proteins, called the anisotropic Langevin Equation for Protein Dynamics or LE4PD-XYZ. This approach partitions the protein's MD dynamics into mostly uncorrelated, wavelength-dependent, diffusive modes. It associates with each mode a free-energy map, where one measures the spatial extension and the time evolution of the mode-dependent, slow dynamical fluctuations. Here, we compare the tICA modes' predictions with the collective LE4PD-XYZ modes. We observe that the two methods consistently identify the nature and extension of the slowest fluctuation processes. The tICA separates the leading processes in a smaller number of slow modes than the LE4PD does. The LE4PD provides time-dependent information at short times and a formal connection to the physics of the kinetic processes that are missing in the pure statistical analysis of tICA.

A quantitative paradigm for water-assisted proton transport through proteins and other confined spaces

Water-assisted proton transport through confined spaces influences many phenomena in biomolecular and nanomaterial systems. In such cases, the water molecules that fluctuate in the confined pathways provide the environment and the medium for the hydrated excess proton migration via Grotthuss shuttling. However, a definitive collective variable (CV) that accurately couples the hydration and the connectivity of the proton wire with the proton translocation has remained elusive. To address this important challenge-and thus to define a quantitative paradigm for facile proton transport in confined spaces-a CV is derived in this work from graph theory, which is verified to accurately describe water wire formation and breakage coupled to the proton translocation in carbon nanotubes and the Cl^-/H⁺ antiporter protein, ClC-ec1. Significant alterations in the conformations and thermodynamics of water wires are uncovered after introducing an excess proton into them. Large barriers in the proton translocation free-energy profiles are found when water wires are defined to be disconnected according to the new CV, even though the pertinent confined space is still reasonably well hydrated and-by the simple measure of the mere existence of a water structure-the proton transport would have been predicted to be facile via that oversimplified measure. In this paradigm, however, the simple presence of water is not sufficient for inferring proton translocation, since an excess proton itself is able to drive hydration, and additionally, the water molecules themselves must be adequately connected to facilitate any successful proton transport.

SGOOP-d: Estimating Kinetic Distances and Reaction Coordinate Dimensionality for Rare Event Systems from Biased/Unbiased Simulations

Understanding kinetics including reaction pathways and associated transition rates is an important yet difficult problem in numerous chemical and biological systems, especially in situations with multiple competing pathways. When these high-dimensional systems are projected on low-dimensional coordinates, which are often needed for enhanced sampling or for interpretation of simulations and experiments, one can end up losing the kinetic connectivity of the underlying high-dimensional landscape. Thus, in the low-dimensional projection, metastable states might appear closer or further than they actually are. To deal with this issue, in this work, we develop a formalism that learns a multidimensional yet minimally complex reaction coordinate (RC) for generic high-dimensional systems. When projected along this RC, all possible kinetically relevant pathways can be demarcated and the true high-dimensional connectivity is maintained. One of the defining attributes of our method lies in that it can work on long unbiased simulations as well as biased simulations often needed for rare event systems. We demonstrate the utility of the method by studying a range of model systems including conformational transitions in a small peptide Ace-Ala₃-Nme, where we show how two-dimensional and three-dimensional RCs found by our previously published spectral gap optimization method "SGOOP" [Tiwary, P. and Berne, B. J. Proc. Natl. Acad. Sci. 2016, 113, 2839] can capture the kinetics for 23 and all 28 out of the 28 dominant state-to-state transitions, respectively.

Determining Sequence-Dependent DNA Oligonucleotide Hybridization and Dehybridization Mechanisms Using Coarse-Grained Molecular Simulation, Markov State Models, and Infrared Spectroscopy

A robust understanding of the sequence-dependent thermodynamics of DNA hybridization has enabled rapid advances in DNA nanotechnology. A fundamental understanding of the sequence-dependent kinetics and mechanisms of hybridization and dehybridization remains comparatively underdeveloped. In this work, we establish new understanding of the sequence-dependent hybridization/dehybridization kinetics and mechanism within a family of self-complementary pairs of 10-mer DNA oligomers by integrating coarse-grained molecular simulation, machine learning of the slow dynamical modes, data-driven inference of long-time kinetic models, and experimental temperature-jump infrared spectroscopy. For a repetitive ATATATATAT sequence, we resolve a rugged dynamical landscape comprising multiple metastable states, numerous competing hybridization/dehybridization pathways, and a spectrum of dynamical relaxations. Introduction of a G:C pair at the terminus (GATATATATC) or center (ATATGCATAT) of the sequence reduces the ruggedness of the dynamics landscape by eliminating a number of metastable states and reducing the number of competing dynamical pathways. Only by introducing a G:C pair midway between the terminus and the center to maximally disrupt the repetitive nature of the sequence (ATGATATCAT) do we recover a canonical "all-or-nothing" two-state model of hybridization/dehybridization with no intermediate metastable states. Our results establish new understanding of the dynamical richness of sequence-dependent kinetics and mechanisms of DNA hybridization/dehybridization by furnishing quantitative and predictive kinetic models of the dynamical transition network between metastable states, present a molecular basis with which to understand experimental temperature jump data, and furnish foundational design rules by which to rationally engineer the kinetics and pathways of DNA association and dissociation for DNA nanotechnology applications.

CATBOSS: Cluster Analysis of Trajectories Based on Segment Splitting

Molecular dynamics (MD) simulations are an exceedingly and increasingly potent tool for molecular behavior prediction and analysis. However, the enormous wealth of data generated by these simulations can be difficult to process and render in a human-readable fashion. Cluster analysis is a commonly used way to partition data into structurally distinct states. We present a method that improves on the state of the art by taking advantage of the temporal information of MD trajectories to enable more accurate clustering at a lower memory cost. To date, cluster analysis of MD simulations has generally treated simulation snapshots as a mere collection of independent data points and attempted to separate them into different clusters based on structural similarity. This new method, cluster analysis of trajectories based on segment splitting (CATBOSS), applies density-peak-based clustering to classify trajectory segments learned by change detection. Applying the method to a synthetic toy model as well as four real-life data sets–trajectories of MD simulations of alanine dipeptide and valine dipeptide as well as two fast-folding proteins–we find CATBOSS to be robust and highly performant, yielding natural-looking cluster boundaries and greatly improving clustering resolution. As the classification of points into segments emphasizes density gaps in the data by grouping them close to the state means, CATBOSS applied to the valine dipeptide system is even able to account for a degree of freedom deliberately omitted from the input data set. We also demonstrate the potential utility of CATBOSS in distinguishing metastable states from transition segments as well as promising application to cases where there is little or no advance knowledge of intrinsic coordinates, making for a highly versatile analysis tool.

Markov state models of proton- and pore-dependent activation in a pentameric ligand-gated ion channel

Ligand-gated ion channels conduct currents in response to chemical stimuli, mediating electrochemical signaling in neurons and other excitable cells. For many channels, the details of gating remain unclear, partly due to limited structural data and simulation timescales. Here, we used enhanced sampling to simulate the pH-gated channel GLIC, and construct Markov state models (MSMs) of gating. Consistent with new functional recordings, we report in oocytes, our analysis revealed differential effects of protonation and mutation on free-energy wells. Clustering of closed- versus open-like states enabled estimation of open probabilities and transition rates, while higher-order clustering affirmed conformational trends in gating. Furthermore, our models uncovered state- and protonation-dependent symmetrization. This demonstrates the applicability of MSMs to map energetic and conformational transitions between ion-channel functional states, and how they reproduce shifts upon activation or mutation, with implications for modeling neuronal function and developing state-selective drugs.

TYROSINE KINASES: COMPLEX MOLECULAR SYSTEMS CHALLENGING COMPUTATIONAL METHODOLOGIES

Classical molecular dynamics (MD) simulations based on atomic models play an increasingly important role in a wide range of applications in physics, biology, and chemistry. Nonetheless, generating genuine knowledge about biological systems using MD simulations remains challenging. Protein tyrosine kinases are important cellular signaling enzymes that regulate cell growth, proliferation, metabolism, differentiation, and migration. Due to the large conformational changes and long timescales involved in their function, these kinases present particularly challenging problems to modern computational and theoretical frameworks aimed at elucidating the dynamics of complex biomolecular systems. Markov state models have achieved limited success in tackling the broader conformational ensemble and biased methods are often employed to examine specific long timescale events. Recent advances in machine learning continue to push the limitations of current methodologies and provide notable improvements when integrated with the existing frameworks. A broad perspective is drawn from a critical review of recent studies.

Unsupervised Learning Methods for Molecular Simulation Data

Unsupervised learning is becoming an essential tool to analyze the increasingly large amounts of data produced by atomistic and molecular simulations, in material science, solid state physics, biophysics, and biochemistry. In this Review, we provide a comprehensive overview of the methods of unsupervised learning that have been most commonly used to investigate simulation data and indicate likely directions for further developments in the field. In particular, we discuss feature representation of molecular systems and present state-of-the-art algorithms of dimensionality reduction, density estimation, and clustering, and kinetic models. We divide our discussion into self-contained sections, each discussing a specific method. In each section, we briefly touch upon the mathematical and algorithmic foundations of the method, highlight its strengths and limitations, and describe the specific ways in which it has been used-or can be used-to analyze molecular simulation data.

Folding and misfolding of potassium channel monomers during assembly and tetramerization

The dynamics and folding of potassium channel pore domain monomers are connected to the kinetics of tetramer assembly. In all-atom molecular dynamics simulations of Kv1.2 and KcsA channels, monomers adopt multiple nonnative conformations while the three helices remain folded. Consistent with this picture, NMR studies also find the monomers to be dynamic and structurally heterogeneous. However, a KcsA construct with a disulfide bridge engineered between the two transmembrane helices has an NMR spectrum with well-dispersed peaks, suggesting that the monomer can be locked into a native-like conformation that is similar to that observed in the folded tetramer. During tetramerization, fluoresence resonance energy transfer (FRET) data indicate that monomers rapidly oligomerize upon insertion into liposomes, likely forming a protein-dense region. Folding within this region occurs along separate fast and slow routes, with τ_fold ∼40 and 1,500 s, respectively. In contrast, constructs bearing the disulfide bond mainly fold via the faster pathway, suggesting that maintaining the transmembrane helices in their native orientation reduces misfolding. Interestingly, folding is concentration independent despite the tetrameric nature of the channel, indicating that the rate-limiting step is unimolecular and occurs after monomer association in the protein-dense region. We propose that the rapid formation of protein-dense regions may help with the assembly of multimeric membrane proteins by bringing together the nascent components prior to assembly. Finally, despite its name, the addition of KcsA's C-terminal "tetramerization" domain does not hasten the kinetics of tetramerization.

Tetranucleosome Interactions Drive Chromatin Folding

The multiscale organizational structure of chromatin in eukaryotic cells is instrumental to DNA transcription, replication, and repair. At mesoscopic length scales, nucleosomes pack in a manner that serves to regulate gene expression through condensation and expansion of the genome. The particular structures that arise and their respective thermodynamic stabilities, however, have yet to be fully resolved. In this study, we combine molecular modeling using the 1CPN mesoscale model of chromatin with nonlinear manifold learning to identify and characterize the structure and free energy of metastable states of short chromatin segments comprising between 4- and 16-nucleosomes. Our results reveal the formation of two previously characterized tetranucleosomal conformations, the "α-tetrahedron" and the "β-rhombus", which have been suggested to play an important role in the accessibility of DNA and, respectively, induce local chromatin compaction or elongation. The spontaneous formation of these motifs is potentially responsible for the slow nucleosome dynamics observed in experimental studies. Increases of the nucleosome repeat length are accompanied by more pronounced structural irregularity and flexibility and, ultimately, a dynamic liquid-like behavior that allows for frequent structural reorganization. Our findings indicate that tetranucleosome motifs are intrinsically stable structural states, driven by local internucleosomal interactions, and support a mechanistic picture of chromatin packing, dynamics, and accessibility that is strongly influenced by emergent local mesoscale structure.

What Markov State Models Can and Cannot Do: Correlation versus Path-Based Observables in Protein-Folding Models

Markov state models (MSMs) have been widely applied to study the kinetics and pathways of protein conformational dynamics based on statistical analysis of molecular dynamics (MD) simulations. These MSMs coarse-grain both configuration space and time in ways that limit what kinds of observables they can reproduce with high fidelity over different spatial and temporal resolutions. Despite their popularity, there is still limited understanding of which biophysical observables can be computed from these MSMs in a robust and unbiased manner, and which suffer from the space-time coarse-graining intrinsic in the MSM model. Most theoretical arguments and practical validity tests for MSMs rely on long-time equilibrium kinetics, such as the slowest relaxation time scales and experimentally observable time-correlation functions. Here, we perform an extensive assessment of the ability of well-validated protein folding MSMs to accurately reproduce path-based observable such as mean first-passage times (MFPTs) and transition path mechanisms compared to a direct trajectory analysis. We also assess a recently proposed class of history-augmented MSMs (haMSMs) that exploit additional information not accounted for in standard MSMs. We conclude with some practical guidance on the use of MSMs to study various problems in conformational dynamics of biomolecules. In brief, MSMs can accurately reproduce correlation functions slower than the lag time, but path-based observables can only be reliably reproduced if the lifetimes of states exceed the lag time, which is a much stricter requirement. Even in the presence of short-lived states, we find that haMSMs reproduce path-based observables more reliably.

A kinetic ensemble of the Alzheimer's Aβ peptide

The conformational and thermodynamic properties of disordered proteins are commonly described in terms of structural ensembles and free energy landscapes. To provide information on the transition rates between the different states populated by these proteins, it would be desirable to generalize this description to kinetic ensembles. Approaches based on the theory of stochastic processes can be particularly suitable for this purpose. Here, we develop a Markov state model and apply it to determine a kinetic ensemble of Aβ42, a disordered peptide associated with Alzheimer's disease. Through the Google Compute Engine, we generated 315-µs all-atom molecular dynamics trajectories. Using a probabilistic-based definition of conformational states in a neural network approach, we found that Aβ42 is characterized by inter-state transitions on the microsecond timescale, exhibiting only fully unfolded or short-lived, partially folded states. Our results illustrate how kinetic ensembles provide effective information about the structure, thermodynamics and kinetics of disordered proteins.

Markov State Models to Elucidate Ligand Binding Mechanism

Molecular dynamics simulations can now routinely access the microsecond timescale, making feasible direct sampling of ligand association events. While Markov State Model (MSM) approaches offer a useful framework for analyzing such trajectory data to gain insight into binding mechanisms, accurate modeling of ligand association pathways and kinetics must be done carefully. We describe methods and good practices for constructing MSMs of ligand binding from unbiased trajectory data and discuss how to use time-lagged independent component analysis (tICA) to build informative models, using as an example recent simulation work to model the binding of phenylalanine to the regulatory ACT domain dimer of phenylalanine hydroxylase. We describe a variety of methods for estimating association rates from MSMs and discuss how to distinguish between conformational selection and induced-fit mechanisms using MSMs. In addition, we review some examples of MSMs constructed to elucidate the mechanisms by which p53 transactivation domain (TAD) and related peptides bind the oncoprotein MDM2.