Non-Markovian modeling of protein folding

Memory-induced absolute negative mobility

Non-Markovian systems form a broad area of physics that remains greatly unexplored despite years of intensive investigations. The spotlight is on memory as a source of effects that are absent in their Markovian counterparts. In this work, we dive into this problem and analyze a driven Brownian particle moving in a spatially periodic potential and exposed to correlated thermal noise. We show that the absolute negative mobility effect, in which the net movement of the particle is in the direction opposite to the average force acting on it, may be induced by the memory of the setup. To explain the origin of this phenomenon, we resort to the recently developed effective mass approach to dynamics of non-Markovian systems.

A Gauss-Newton method for iterative optimization of memory kernels for generalized Langevin thermostats in coarse-grained molecular dynamics simulations

In molecular dynamics simulations, dynamically consistent coarse-grained (CG) models commonly use stochastic thermostats to model friction and fluctuations that are lost in a CG description. While Markovian, i.e., time-local, formulations of such thermostats allow for an accurate representation of diffusivities/long-time dynamics, a correct description of the dynamics on all time scales generally requires non-Markovian, i.e., non-time-local, thermostats. These thermostats typically take the form of a Generalized Langevin Equation (GLE) determined by a memory kernel. In this work, we use a Markovian embedded formulation of a position-independent GLE thermostat acting independently on each CG degree of freedom. Extracting the memory kernel of this CG model from atomistic reference data requires several approximations. Therefore, this task is best understood as an inverse problem. While our recently proposed approximate Newton scheme allows for the iterative optimization of memory kernels (IOMK), Markovian embedding remained potentially error-prone and computationally expensive. In this work, we present an IOMK-Gauss-Newton scheme (IOMK-GN) based on IOMK that allows for the direct parameterization of a Markovian embedded model.

Memory Corrections to Markovian Langevin Dynamics

Analysis of non-Markovian systems and memory-induced phenomena poses an everlasting challenge in the realm of physics. As a paradigmatic example, we consider a classical Brownian particle of mass M subjected to an external force and exposed to correlated thermal fluctuations. We show that the recently developed approach to this system, in which its non-Markovian dynamics given by the Generalized Langevin Equation is approximated by its memoryless counterpart but with the effective particle mass M∗ Collapse

Key Words

• colored noise
• correlated fluctuations
• memory
• non-Markovian

Collapse

MESH Headings Collapse

Grants

• 2022/45/B/ST3/02619 National Science Center

Collapse

Affiliation(s)

Mateusz Wiśniewski
Jerzy Łuczka
Jakub Spiechowicz Institute of Physics, University of Silesia, 41-500 Chorzów, Poland

Collapse

The role of memory-dependent friction and solvent viscosity in isomerization kinetics in viscogenic media

Molecular isomerization kinetics in liquid solvent depends on a complex interplay between the solvent friction acting on the molecule, internal dissipation effects (also known as internal friction), the viscosity of the solvent, and the dihedral free energy profile. Due to the absence of accurate techniques to directly evaluate isomerization friction, it has not been possible to explore these relationships in full. By combining extensive molecular dynamics simulations with friction memory-kernel extraction techniques we consider a variety of small, isomerising molecules under a range of different viscogenic conditions and directly evaluate the viscosity dependence of the friction acting on a rotating dihedral. We reveal that the influence of different viscogenic media on isomerization kinetics can be dramatically different, even when measured at the same viscosity. This is due to the dynamic solute-solvent coupling, mediated by time-dependent friction memory kernels. We also show that deviations from the linear dependence of isomerization rates on solvent viscosity, which are often simply attributed to internal friction effects, are due to the simultaneous violation of two fundamental relationships: the Stokes-Einstein relation and the overdamped Kramers prediction for the barrier-crossing rate, both of which require explicit knowledge of friction.

Accurate Memory Kernel Extraction from Discretized Time-Series Data

Memory effects emerge as a fundamental consequence of dimensionality reduction when low-dimensional observables are used to describe the dynamics of complex many-body systems. In the context of molecular dynamics (MD) data analysis, accounting for memory effects using the framework of the generalized Langevin equation (GLE) has proven efficient, accurate, and insightful, particularly when working with high-resolution time series data. However, in experimental systems, high-resolution data are often unavailable, raising questions about the impact of the data resolution on the estimated GLE parameters. This study demonstrates that direct memory extraction from time series data remains accurate when the discretization time is below the memory time. To obtain memory functions reliably, even when the discretization time exceeds the memory time, we introduce a Gaussian Process Optimization (GPO) scheme. This scheme minimizes the deviation of discretized two-point correlation functions between time series data and GLE simulations and is able to estimate accurate memory kernels as long as the discretization time stays below the longest time scale in the data, typically the barrier crossing time.

Mechanism of DNA origami folding elucidated by mesoscopic simulations

Many experimental and computational efforts have sought to understand DNA origami folding, but the time and length scales of this process pose significant challenges. Here, we present a mesoscopic model that uses a switchable force field to capture the behavior of single- and double-stranded DNA motifs and transitions between them, allowing us to simulate the folding of DNA origami up to several kilobases in size. Brownian dynamics simulations of small structures reveal a hierarchical folding process involving zipping into a partially folded precursor followed by crystallization into the final structure. We elucidate the effects of various design choices on folding order and kinetics. Larger structures are found to exhibit heterogeneous staple incorporation kinetics and frequent trapping in metastable states, as opposed to more accessible structures which exhibit first-order kinetics and virtually defect-free folding. This model opens an avenue to better understand and design DNA nanostructures for improved yield and folding performance.

Effective mass approach to memory in non-Markovian systems

Recent pioneering experiments on non-Markovian dynamics done, e.g., for active matter have demonstrated that our theoretical understanding of this challenging yet hot topic is rather incomplete and there is a wealth of phenomena still awaiting discovery. It is related to the fact that typically for simplification the Markovian approximation is employed and as a consequence the memory is neglected. Therefore, methods allowing to study memory effects are extremely valuable. We demonstrate that a non-Markovian system described by the Generalized Langevin Equation (GLE) for a Brownian particle of mass M can be approximated by the memoryless Langevin equation in which the memory effects are correctly reproduced solely via the effective mass M^{*} of the Brownian particle which is determined only by the form of the memory kernel. Our work lays the foundation for an impactful approach which allows one to readily study memory-related corrections to Markovian dynamics.

Tutorial on how to build non-Markovian dynamic models from molecular dynamics simulations for studying protein conformational changes

Protein conformational changes play crucial roles in their biological functions. In recent years, the Markov State Model (MSM) constructed from extensive Molecular Dynamics (MD) simulations has emerged as a powerful tool for modeling complex protein conformational changes. In MSMs, dynamics are modeled as a sequence of Markovian transitions among metastable conformational states at discrete time intervals (called lag time). A major challenge for MSMs is that the lag time must be long enough to allow transitions among states to become memoryless (or Markovian). However, this lag time is constrained by the length of individual MD simulations available to track these transitions. To address this challenge, we have recently developed Generalized Master Equation (GME)-based approaches, encoding non-Markovian dynamics using a time-dependent memory kernel. In this Tutorial, we introduce the theory behind two recently developed GME-based non-Markovian dynamic models: the quasi-Markov State Model (qMSM) and the Integrative Generalized Master Equation (IGME). We subsequently outline the procedures for constructing these models and provide a step-by-step tutorial on applying qMSM and IGME to study two peptide systems: alanine dipeptide and villin headpiece. This Tutorial is available at https://github.com/xuhuihuang/GME_tutorials. The protocols detailed in this Tutorial aim to be accessible for non-experts interested in studying the biomolecular dynamics using these non-Markovian dynamic models.

Chemically Specific Systematic Coarse-Grained Polymer Model with Both Consistently Structural and Dynamical Properties

The coarse-grained (CG) model serves as a powerful tool for the simulation of polymer systems; its reliability depends on the accurate representation of both structural and dynamical properties. However, strong correlations between structural and dynamical properties on different scales and also a strong memory effect, enforced by chain connectivity between monomers in polymer systems, render developing a chemically specific systematic CG model a formidable task. In this study, we report a systematic CG approach that combines the iterative Boltzmann inversion (IBI) method and the generalized Langevin equation (GLE) dynamics. Structural properties are ensured by using conservative CG potentials derived from the IBI method. To retrieve the correct dynamical properties in the system, we demonstrate that using a combination of a Rouse-type delta function and a time-dependent short-time kernel in the GLE simulation is practically efficient. The former can be used to adjust the long-time diffusion dynamics, and the latter can be reconstructed from an iterative procedure according to the velocity autocorrelation function (ACF) from all-atomistic (AA) simulations. Taking the polystyrene as an example, we show that not only structural properties of radial distribution function, intramolecular bond, and angle distributions can be reproduced but also dynamical properties of mean-square displacement, velocity ACF, and force ACF resulted from our CG model have quantitative agreement with the reference AA model. In addition, reasonable agreements are observed in other collective properties between our GLE-CG model and the AA simulations as well.

Main role of fractal-like nature of conformational space in subdiffusion in proteins

Protein dynamics involves a myriad of mechanical movements happening at different time and space scales, which make it highly complex. One of the less understood features of protein dynamics is subdiffusivity, defined as sublinear dependence between displacement and time. Here, we use all-atoms molecular dynamics (MD) simulations to directly interrogate an already well-established theory and demonstrate that subdiffusivity arises from the fractal nature of the network of metastable conformations over which the dynamics, thought of as a diffusion process, takes place.

Accurate estimates of dynamical statistics using memory

Many chemical reactions and molecular processes occur on time scales that are significantly longer than those accessible by direct simulations. One successful approach to estimating dynamical statistics for such processes is to use many short time series of observations of the system to construct a Markov state model, which approximates the dynamics of the system as memoryless transitions between a set of discrete states. The dynamical Galerkin approximation (DGA) is a closely related framework for estimating dynamical statistics, such as committors and mean first passage times, by approximating solutions to their equations with a projection onto a basis. Because the projected dynamics are generally not memoryless, the Markov approximation can result in significant systematic errors. Inspired by quasi-Markov state models, which employ the generalized master equation to encode memory resulting from the projection, we reformulate DGA to account for memory and analyze its performance on two systems: a two-dimensional triple well and the AIB9 peptide. We demonstrate that our method is robust to the choice of basis and can decrease the time series length required to obtain accurate kinetics by an order of magnitude.

Log-periodic oscillations as real-time signatures of hierarchical dynamics in proteins

The time-dependent relaxation of a dynamical system may exhibit a power-law behavior that is superimposed by log-periodic oscillations. D. Sornette [Phys. Rep. 297, 239 (1998)] showed that this behavior can be explained by a discrete scale invariance of the system, which is associated with discrete and equidistant timescales on a logarithmic scale. Examples include such diverse fields as financial crashes, random diffusion, and quantum topological materials. Recent time-resolved experiments and molecular dynamics simulations suggest that discrete scale invariance may also apply to hierarchical dynamics in proteins, where several fast local conformational changes are a prerequisite for a slow global transition to occur. Employing entropy-based timescale analysis and Markov state modeling to a simple one-dimensional hierarchical model and biomolecular simulation data, it is found that hierarchical systems quite generally give rise to logarithmically spaced discrete timescales. By introducing a one-dimensional reaction coordinate that collectively accounts for the hierarchically coupled degrees of freedom, the free energy landscape exhibits a characteristic staircase shape with two metastable end states, which causes the log-periodic time evolution of the system. The period of the log-oscillations reflects the effective roughness of the energy landscape and can, in simple cases, be interpreted in terms of the barriers of the staircase landscape.

MoDAFold: a strategy for predicting the structure of missense mutant protein based on AlphaFold2 and molecular dynamics

Protein structure prediction is a longstanding issue crucial for identifying new drug targets and providing a mechanistic understanding of protein functions. To enhance the progress in this field, a spectrum of computational methodologies has been cultivated. AlphaFold2 has exhibited exceptional precision in predicting wild-type protein structures, with performance exceeding that of other methods. However, predicting the structures of missense mutant proteins using AlphaFold2 remains challenging due to the intricate and substantial structural alterations caused by minor sequence variations in the mutant proteins. Molecular dynamics (MD) has been validated for precisely capturing changes in amino acid interactions attributed to protein mutations. Therefore, for the first time, a strategy entitled 'MoDAFold' was proposed to improve the accuracy and reliability of missense mutant protein structure prediction by combining AlphaFold2 with MD. Multiple case studies have confirmed the superior performance of MoDAFold compared to other methods, particularly AlphaFold2.

Analysis of Velocity Autocorrelation Functions from Molecular Dynamics Simulations of a Small Peptide by the Generalized Langevin Equation with a Power-Law Kernel

Internal motions play an essential role in the biological functions of proteins and have been the subject of numerous theoretical and spectroscopic studies. Such complex environments are associated with anomalous diffusion where, in contrast to the classical Brownian motion, the relevant correlation functions have power law decays with time. In this work, we investigate the presence of long memory stochastic processes through the analysis of atomic velocity autocorrelation functions. Analytical expressions of the velocity autocorrelation function spectrum obtained through a Mori-Zwanzig projection approach were shown to be compatible with molecular dynamics simulations of a small helical peptide (8-polyalanine).

Investigation of Rare Protein Conformational Transitions via Dissipation-Corrected Targeted Molecular Dynamics

To sample rare events, dissipation-corrected targeted molecular dynamics (dcTMD) applies a constant velocity constraint along a one-dimensional reaction coordinate s, which drives an atomistic system from an initial state into a target state. Employing a cumulant approximation of Jarzynski's identity, the free energy ΔG(s) is calculated from the mean external work and dissipated work of the process. By calculating the friction coefficient Γ(s) from the dissipated work, in a second step, the equilibrium dynamics of the process can be studied by propagating a Langevin equation. While so far dcTMD has been mostly applied to study the unbinding of protein-ligand complexes, here its applicability to rare conformational transitions within a protein and the prediction of their kinetics are investigated. As this typically requires the introduction of multiple collective variables {xj} = x, a theoretical framework is outlined to calculate the associated free energy ΔG(x) and friction Γ(x) from dcTMD simulations along coordinate s. Adopting the α-β transition of alanine dipeptide as well as the open-closed transition of T4 lysozyme as representative examples, the virtues and shortcomings of dcTMD to predict protein conformational transitions and the related kinetics are studied.

Inferring free-energy barriers and kinetic rates from molecular dynamics via underdamped Langevin models

Rare events include many of the most interesting transformation processes in condensed matter, from phase transitions to biomolecular conformational changes to chemical reactions. Access to the corresponding mechanisms, free-energy landscapes and kinetic rates can in principle be obtained by different techniques after projecting the high-dimensional atomic dynamics on one (or a few) collective variable. Even though it is well-known that the projected dynamics approximately follows - in a statistical sense - the generalized, underdamped or overdamped Langevin equations (depending on the time resolution), to date it is nontrivial to parameterize such equations starting from a limited, practically accessible amount of non-ergodic trajectories. In this work we focus on Markovian, underdamped Langevin equations, that arise naturally when considering, e.g., numerous water-solution processes at sub-picosecond resolution. After contrasting the advantages and pitfalls of different numerical approaches, we present an efficient parametrization strategy based on a limited set of molecular dynamics data, including equilibrium trajectories confined to minima and few hundreds transition path sampling-like trajectories. Employing velocity autocorrelation or memory kernel information for learning the friction and likelihood maximization for learning the free-energy landscape, we demonstrate the possibility to reconstruct accurate barriers and rates both for a benchmark system and for the interaction of carbon nanoparticles in water.

Integrative generalized master equation: A method to study long-timescale biomolecular dynamics via the integrals of memory kernels

The generalized master equation (GME) provides a powerful approach to study biomolecular dynamics via non-Markovian dynamic models built from molecular dynamics (MD) simulations. Previously, we have implemented the GME, namely the quasi Markov State Model (qMSM), where we explicitly calculate the memory kernel and propagate dynamics using a discretized GME. qMSM can be constructed with much shorter MD trajectories than the MSM. However, since qMSM needs to explicitly compute the time-dependent memory kernels, it is heavily affected by the numerical fluctuations of simulation data when applied to study biomolecular conformational changes. This can lead to numerical instability of predicted long-time dynamics, greatly limiting the applicability of qMSM in complicated biomolecules. We present a new method, the Integrative GME (IGME), in which we analytically solve the GME under the condition when the memory kernels have decayed to zero. Our IGME overcomes the challenges of the qMSM by using the time integrations of memory kernels, thereby avoiding the numerical instability caused by explicit computation of time-dependent memory kernels. Using our solutions of the GME, we have developed a new approach to compute long-time dynamics based on MD simulations in a numerically stable, accurate and efficient way. To demonstrate its effectiveness, we have applied the IGME in three biomolecules: the alanine dipeptide, FIP35 WW-domain, and Taq RNA polymerase. In each system, the IGME achieves significantly smaller fluctuations for both memory kernels and long-time dynamics compared to the qMSM. We anticipate that the IGME can be widely applied to investigate biomolecular conformational changes.

Toward a structural identification of metastable molecular conformations

Interpreting high-dimensional data from molecular dynamics simulations is a persistent challenge. In this paper, we show that for a small peptide, deca-alanine, metastable states can be identified through a neural net based on structural information alone. While processing molecular dynamics data, dimensionality reduction is a necessary step that projects high-dimensional data onto a low-dimensional representation that, ideally, captures the conformational changes in the underlying data. Conventional methods make use of the temporal information contained in trajectories generated through integrating the equations of motion, which forgoes more efficient sampling schemes. We demonstrate that EncoderMap, an autoencoder architecture with an additional distance metric, can find a suitable low-dimensional representation to identify long-lived molecular conformations using exclusively structural information. For deca-alanine, which exhibits several helix-forming pathways, we show that this approach allows us to combine simulations with different biasing forces and yields representations comparable in quality to other established methods. Our results contribute to computational strategies for the rapid automatic exploration of the configuration space of peptides and proteins.

Anomalous diffusion, nonergodicity, non-Gaussianity, and aging of fractional Brownian motion with nonlinear clocks

How do nonlinear clocks in time and/or space affect the fundamental properties of a stochastic process? Specifically, how precisely may ergodic processes such as fractional Brownian motion (FBM) acquire predictable nonergodic and aging features being subjected to such conditions? We address these questions in the current study. To describe different types of non-Brownian motion of particles-including power-law anomalous, ultraslow or logarithmic, as well as superfast or exponential diffusion-we here develop and analyze a generalized stochastic process of scaled-fractional Brownian motion (SFBM). The time- and space-SFBM processes are, respectively, constructed based on FBM running with nonlinear time and space clocks. The fundamental statistical characteristics such as non-Gaussianity of particle displacements, nonergodicity, as well as aging are quantified for time- and space-SFBM by selecting different clocks. The latter parametrize power-law anomalous, ultraslow, and superfast diffusion. The results of our computer simulations are fully consistent with the analytical predictions for several functional forms of clocks. We thoroughly examine the behaviors of the probability-density function, the mean-squared displacement, the time-averaged mean-squared displacement, as well as the aging factor. Our results are applicable for rationalizing the impact of nonlinear time and space properties superimposed onto the FBM-type dynamics. SFBM offers a general framework for a universal and more precise model-based description of anomalous, nonergodic, non-Gaussian, and aging diffusion in single-molecule-tracking observations.

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.

A deep learning approach to the measurement of long-lived memory kernels from generalized Langevin dynamics

Memory effects are ubiquitous in a wide variety of complex physical phenomena, ranging from glassy dynamics and metamaterials to climate models. The Generalized Langevin Equation (GLE) provides a rigorous way to describe memory effects via the so-called memory kernel in an integro-differential equation. However, the memory kernel is often unknown, and accurately predicting or measuring it via, e.g., a numerical inverse Laplace transform remains a herculean task. Here, we describe a novel method using deep neural networks (DNNs) to measure memory kernels from dynamical data. As a proof-of-principle, we focus on the notoriously long-lived memory effects of glass-forming systems, which have proved a major challenge to existing methods. In particular, we learn the operator mapping dynamics to memory kernels from a training set generated with the Mode-Coupling Theory (MCT) of hard spheres. Our DNNs are remarkably robust against noise, in contrast to conventional techniques. Furthermore, we demonstrate that a network trained on data generated from analytic theory (hard-sphere MCT) generalizes well to data from simulations of a different system (Brownian Weeks-Chandler-Andersen particles). Finally, we train a network on a set of phenomenological kernels and demonstrate its effectiveness in generalizing to both unseen phenomenological examples and supercooled hard-sphere MCT data. We provide a general pipeline, KernelLearner, for training networks to extract memory kernels from any non-Markovian system described by a GLE. The success of our DNN method applied to noisy glassy systems suggests that deep learning can play an important role in the study of dynamical systems with memory.

Anisotropic Friction in a Ligand-Protein Complex

The effect of an externally applied directional force on molecular friction is so far poorly understood. Here, we study the force-driven dissociation of the ligand-protein complex biotin-streptavidin and identify anisotropic friction as a not yet described type of molecular friction. Using AFM-based stereographic single molecule force spectroscopy and targeted molecular dynamics simulations, we find that the rupture force and friction for biotin-streptavidin vary with the pulling angle. This observation holds true for friction extracted from Kramers' rate expression and by dissipation-corrected targeted molecular dynamics simulations based on Jarzynski's identity. We rule out ligand solvation and protein-internal friction as sources of the angle-dependent friction. Instead, we observe a heterogeneity in free energy barriers along an experimentally uncontrolled orientation parameter, which increases the rupture force variance and therefore the overall friction. We anticipate that anisotropic friction needs to be accounted for in a complete understanding of friction in biomolecular dynamics and anisotropic mechanical environments.

Selecting Features for Markov Modeling: A Case Study on HP35

Markov state models represent a popular means to interpret molecular dynamics trajectories in terms of memoryless transitions between metastable conformational states. To provide a mechanistic understanding of the considered biomolecular process, these states should reflect structurally distinct conformations and ensure a time scale separation between fast intrastate and slow interstate dynamics. Adopting the folding of villin headpiece (HP35) as a well-established model problem, here we discuss the selection of suitable input coordinates or "features", such as backbone dihedral angles and interresidue distances. We show that dihedral angles account accurately for the structure of the native energy basin of HP35, while the unfolded region of the free energy landscape and the folding process are best described by tertiary contacts of the protein. To construct a contact-based model, we consider various ways to define and select contact distances and introduce a low-pass filtering of the feature trajectory as well as a correlation-based characterization of states. Relying on input data that faithfully account for the mechanistic origin of the studied process, the states of the resulting Markov model are clearly discriminated by the features, describe consistently the hierarchical structure of the free energy landscape, and─as a consequence─correctly reproduce the slow time scales of the process.

Memory Unlocks the Future of Biomolecular Dynamics: Transformative Tools to Uncover Physical Insights Accurately and Efficiently

Conformational changes underpin function and encode complex biomolecular mechanisms. Gaining atomic-level detail of how such changes occur has the potential to reveal these mechanisms and is of critical importance in identifying drug targets, facilitating rational drug design, and enabling bioengineering applications. While the past two decades have brought Markov state model techniques to the point where practitioners can regularly use them to glimpse the long-time dynamics of slow conformations in complex systems, many systems are still beyond their reach. In this Perspective, we discuss how including memory (i.e., non-Markovian effects) can reduce the computational cost to predict the long-time dynamics in these complex systems by orders of magnitude and with greater accuracy and resolution than state-of-the-art Markov state models. We illustrate how memory lies at the heart of successful and promising techniques, ranging from the Fokker-Planck and generalized Langevin equations to deep-learning recurrent neural networks and generalized master equations. We delineate how these techniques work, identify insights that they can offer in biomolecular systems, and discuss their advantages and disadvantages in practical settings. We show how generalized master equations can enable the investigation of, for example, the gate-opening process in RNA polymerase II and demonstrate how our recent advances tame the deleterious influence of statistical underconvergence of the molecular dynamics simulations used to parameterize these techniques. This represents a significant leap forward that will enable our memory-based techniques to interrogate systems that are currently beyond the reach of even the best Markov state models. We conclude by discussing some current challenges and future prospects for how exploiting memory will open the door to many exciting opportunities.

Building insightful, memory-enriched models to capture long-time biochemical processes from short-time simulations

The ability to predict and understand complex molecular motions occurring over diverse timescales ranging from picoseconds to seconds and even hours in biological systems remains one of the largest challenges to chemical theory. Markov state models (MSMs), which provide a memoryless description of the transitions between different states of a biochemical system, have provided numerous important physically transparent insights into biological function. However, constructing these models often necessitates performing extremely long molecular simulations to converge the rates. Here, we show that by incorporating memory via the time-convolutionless generalized master equation (TCL-GME) one can build a theoretically transparent and physically intuitive memory-enriched model of biochemical processes with up to a three order of magnitude reduction in the simulation data required while also providing a higher temporal resolution. We derive the conditions under which the TCL-GME provides a more efficient means to capture slow dynamics than MSMs and rigorously prove when the two provide equally valid and efficient descriptions of the slow configurational dynamics. We further introduce a simple averaging procedure that enables our TCL-GME approach to quickly converge and accurately predict long-time dynamics even when parameterized with noisy reference data arising from short trajectories. We illustrate the advantages of the TCL-GME using alanine dipeptide, the human argonaute complex, and FiP35 WW domain.

The effect of time resolution on the observed first passage times in diffusive dynamics

Single-molecule and single-particle tracking experiments are typically unable to resolve fine details of thermal motion at short timescales where trajectories are continuous. We show that, when a diffusive trajectory xt is sampled at finite time intervals δt, the resulting error in measuring the first passage time to a given domain can exceed the time resolution of the measurement by more than an order of magnitude. Such surprisingly large errors originate from the fact that the trajectory may enter and exit the domain while being unobserved, thereby lengthening the apparent first passage time by an amount that is larger than δt. Such systematic errors are particularly important in single-molecule studies of barrier crossing dynamics. We show that the correct first passage times, as well as other properties of the trajectories such as splitting probabilities, can be recovered via a stochastic algorithm that reintroduces unobserved first passage events probabilistically.

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

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

Diffusion in a two-dimensional energy landscape in the presence of dynamical correlations and validity of random walk model

Diffusion in a multidimensional energy surface with minima and barriers is a problem of importance in statistical mechanics and it also has wide applications, such as protein folding. To understand it in such a system, we carry out theory and simulations of a tagged particle moving on a two-dimensional periodic potential energy surface, both in the presence and absence of noise. Langevin dynamics simulations at multiple temperatures are carried out to obtain the diffusion coefficient of a solute particle. Friction is varied from zero to large values. Diffusive motion emerges in the limit of a long time, even in the absence of noise. Noise destroys the correlations and increases diffusion at small friction. Diffusion thus exhibits a nonmonotonic friction dependence at the intermediate value of the damping, ultimately converging to our theoretically predicted value. The latter is obtained using the well-established relationship between diffusion and random walk. An excellent agreement is obtained between theory and simulations in the high-friction limit but not so in the intermediate regime. We explain the deviation in the low- to intermediate-friction regime using the modified random walk theory. The rate of escape from one cell to another is obtained from the multidimensional rate theory of Langer. We find that enhanced dimensionality plays an important role. To quantify the effects of noise on the potential-imposed coherence on the trajectories, we calculate the Lyapunov exponent. At small friction values, the Lyapunov exponent mimics the friction dependence of the rate.

Compact and complete description of non-Markovian dynamics

Generalized master equations provide a theoretically rigorous framework to capture the dynamics of processes ranging from energy harvesting in plants and photovoltaic devices to qubit decoherence in quantum technologies and even protein folding. At their center is the concept of memory. The explicit time-nonlocal description of memory is both protracted and elaborate. When physical intuition is at a premium, one would desire a more compact, yet complete, description. Here, we demonstrate how and when the time-convolutionless formalism constitutes such a description. In particular, by focusing on the dissipative dynamics of the spin-boson and Frenkel exciton models, we show how to: easily construct the time-local generator from reference reduced dynamics, elucidate the dependence of its existence on the system parameters and the choice of reduced observables, identify the physical origin of its apparent divergences, and offer analysis tools to diagnose their severity and circumvent their deleterious effects. We demonstrate that, when applicable, the time-local approach requires as little information as the more commonly used time-nonlocal scheme, with the important advantages of providing a more compact description, greater algorithmic simplicity, and physical interpretability. We conclude by introducing the discrete-time analog and a straightforward protocol to employ it in cases where the reference dynamics have limited resolution. The insights we present here offer the potential for extending the reach of dynamical methods, reducing both their cost and conceptual complexity.

Pair-Reaction Dynamics in Water: Competition of Memory, Potential Shape, and Inertial Effects

When described by a one-dimensional reaction coordinate, pair-reaction rates in a solvent depend, in addition to the potential barrier height and the friction coefficient, on the potential shape, the effective mass, and the friction relaxation spectrum, but a rate theory that accurately accounts for all of these effects does not exist. After a review of classical reaction-rate theories, we show how to extract all parameters of the generalized Langevin equation (GLE) and, in particular, the friction memory function from molecular dynamics (MD) simulations of two prototypical pair reactions in water, the dissociation of NaCl and of two methane molecules. The memory exhibits multiple time scales and, for NaCl, pronounced oscillatory components. Simulations of the GLE by Markovian embedding techniques accurately reproduce the pair-reaction kinetics from MD simulations without any fitting parameters, which confirms the accuracy of the approximative form of the GLE and of the parameter extraction techniques. By modification of the GLE parameters, we investigate the relative importance of memory, mass, and potential shape effects. Neglect of memory slows down NaCl and methane dissociation by roughly a factor of 2; neglect of mass accelerates reactions by a similar factor, and the harmonic approximation of the potential shape gives rise to slight acceleration. This partial error cancellation explains why Kramers' theory, which neglects memory effects and treats the potential shape in harmonic approximation, describes reaction rates better than more sophisticated theories. In essence, all three effects, friction memory, inertia, and the potential shape nonharmonicity, are important to quantitatively describe pair-reaction kinetics in water.

Explicit analytical form for memory kernel in the generalized Langevin equation for end-to-end vector of Rouse chains

The generalized Langevin equation (GLE) provides an attractive theoretical framework for investigating the dynamics of conformational fluctuations of polymeric systems. While the memory kernel is a central function in the GLE, explicit analytical forms for this function have been challenging to obtain, even for the simple models of polymer dynamics. Here, we achieve an explicit analytical expression for the memory kernel in the GLE for the end-to-end vector of Rouse chains in the overdamped limit. Our derivation takes advantage of the finding that the dynamics of the end-to-end vector of Rouse chains with both free ends are equivalent to those of Rouse chains with one free end and the other fixed. For the latter model, we first show that the equations of motion of the Rouse modes as well as their statistical properties can be obtained under the boundary conditions where the free end is held fixed temporarily. We then analytically solve the terms associated with intrachain interactions in the GLE. By formally comparing these terms with the GLE based on the Rouse modes, we obtain an explicit expression for the memory kernel, along with analytical forms for the potential field and the random colored noise force. Our analytical memory kernel is confirmed by numerical calculations in the Laplace space and is shown to yield asymptotic behaviors that are consistent with previous studies. Finally, we utilize our analytical result to simulate the cyclization dynamics of Rouse chains and discuss the scaling of the cyclization time with chain length.

Assessing transition rates as functions of environmental variables

We present a method to estimate the transition rates of molecular systems under different environmental conditions that cause the formation or the breaking of bonds and require the sampling of the Grand Canonical Ensemble. For this purpose, we model the molecular system in terms of probable "scenarios," governed by different potential energy functions, which are separately sampled by classical MD simulations. Reweighting the canonical distribution of each scenario according to specific environmental variables, we estimate the grand canonical distribution, then use the Square Root Approximation method to discretize the Fokker-Planck operator into a rate matrix and the robust Perron Cluster Cluster Analysis method to coarse-grain the kinetic model. This permits efficiently estimating the transition rates of conformational states as functions of environmental variables, for example, the local pH at a cell membrane. In this work, we formalize the theoretical framework of the procedure, and we present a numerical experiment comparing the results with those provided by a constant-pH method based on non-equilibrium Molecular Dynamics Monte Carlo simulations. The method is relevant for the development of new drug design strategies that take into account how the cellular environment influences biochemical processes.

Generalized Langevin dynamics simulation with non-stationary memory kernels: How to make noise

We present a numerical method to produce stochastic dynamics according to the generalized Langevin equation with a non-stationary memory kernel. This type of dynamics occurs when a microscopic system with an explicitly time-dependent Liouvillian is coarse-grained by means of a projection operator formalism. We show how to replace the deterministic fluctuating force in the generalized Langevin equation by a stochastic process, such that the distributions of the observables are reproduced up to moments of a given order. Thus, in combination with a method to extract the memory kernel from simulation data of the underlying microscopic model, the method introduced here allows us to construct and simulate a coarse-grained model for a driven process.

Microscopic origin of diffusive dynamics in the context of transition path time distributions for protein folding and unfolding

Experimentally measured transition path time distributions are usually analyzed theoretically in terms of a diffusion equation over a free energy barrier. It is though well understood that the free energy profile separating the folded and unfolded states of a protein is characterized as a transition through many stable micro-states which exist between the folded and unfolded states. Why is it then justified to model the transition path dynamics in terms of a diffusion equation, namely the Smoluchowski equation (SE)? In principle, van Kampen has shown that a nearest neighbor Markov chain of thermal jumps between neighboring microstates will lead in a continuum limit to the SE, such that the friction coefficient is proportional to the mean residence time in each micro-state. However, the practical question of how many microstates are needed to justify modeling the transition path dynamics in terms of an SE has not been addressed. This is a central topic of this paper where we compare numerical results for transition paths based on the diffusion equation on the one hand and the nearest neighbor Markov jump model on the other. Comparison of the transition path time distributions shows that one needs at least a few dozen microstates to obtain reasonable agreement between the two approaches. Using the Markov nearest neighbor model one also obtains good agreement with the experimentally measured transition path time distributions for a DNA hairpin and PrP protein. As found previously when using the diffusion equation, the Markov chain model used here also reproduces the experimentally measured long time tail and confirms that the transition path barrier height is ∼3k_BT. This study indicates that in the future, when attempting to model experimentally measured transition path time distributions, one should perhaps prefer a nearest neighbor Markov model which is well defined also for rough energy landscapes. Such studies can also shed light on the minimal number of microstates needed to unravel the experimental data.

Barrier-crossing times for different non-Markovian friction in well and barrier: A numerical study

We introduce a generalized Langevin model system for different non-Markovian effects in the well and barrier regions of a potential, and use it to numerically study the barrier-crossing time. In the appropriate limits, our model interpolates between the theoretical barrier-crossing-time predictions by Grote and Hynes (GH), as well as by Pollak et al., which for a single barrier memory time can differ by several orders of magnitude. Our model furthermore allows one to test an analytic rate theory for space-inhomogeneous memory, which disagrees with our numerical results in the long well-memory regime. In this regime, we find that short barrier memory decreases the barrier-crossing time as compared to long barrier memory. This is in contrast with the short well-memory regime, where both our numerical results and the GH theory predict an acceleration of the barrier crossing time with increasing barrier memory time. Both effects, the "Markovian-barrier acceleration" and GH "non-Markovian-barrier acceleration," can be understood from a committor analysis. Our model combines finite relaxation times of orthogonal degrees of freedom with a space-inhomogeneous coupling to such degrees and represents a step towards more realistic modeling of reaction coordinates.

Free Energy Landscapes, Diffusion Coefficients, and Kinetic Rates from Transition Paths

We address the problem of constructing accurate mathematical models of the dynamics of complex systems projected on a collective variable. To this aim we introduce a conceptually simple yet effective algorithm for estimating the parameters of Langevin and Fokker-Planck equations from a set of short, possibly out-of-equilibrium molecular dynamics trajectories, obtained for instance from transition path sampling or as relaxation from high free-energy configurations. The approach maximizes the model likelihood based on any explicit expression of the short-time propagator, hence it can be applied to different evolution equations. We demonstrate the numerical efficiency and robustness of the algorithm on model systems, and we apply it to reconstruct the projected dynamics of pairs of C₆₀ and C₂₄₀ fullerene molecules in explicit water. Our methodology allows reconstructing the accurate thermodynamics and kinetics of activated processes, namely free energy landscapes, diffusion coefficients, and kinetic rates. Compared to existing enhanced sampling methods, we directly exploit short unbiased trajectories at a competitive computational cost.

Effects of Hydrodynamic Backflow on the Transmission Coefficient of a Barrier-Crossing Brownian Particle

The slow power law decay of the velocity autocorrelation function of a particle moving stochastically in a condensed-phase fluid is widely attributed to the momentum that fluid molecules displaced by the particle transfer back to it during the course of its motion. The forces created by this backflow effect are known as Basset forces, and they have been found in recent analytical work and numerical simulations to be implicated in a number of interesting dynamical phenomena, including boosted particle mobility in tilted washboard potentials. Motivated by these findings, the present paper is an investigation of the role of backflow in thermally activated barrier crossing, the governing process in essentially all condensed-phase chemical reactions. More specifically, it is an exact analytical calculation, carried out within the framework of the reactive-flux formalism, of the transmission coefficient κ(t) of a Brownian particle that crosses an inverted parabola under the influence of a colored noise process originating in the Basset force and a Markovian time-local friction. The calculation establishes that κ(t) is significantly enhanced over its backflow-free limit.

Position-dependent memory kernel in generalized Langevin equations: Theory and numerical estimation

Generalized Langevin equations with non-linear forces and position-dependent linear friction memory kernels, such as commonly used to describe the effective dynamics of coarse-grained variables in molecular dynamics, are rigorously derived within the Mori–Zwanzig formalism. A fluctuation–dissipation theorem relating the properties of the noise to the memory kernel is shown. The derivation also yields Volterra-type equations for the kernel, which can be used for a numerical parametrization of the model from all-atom simulations.

Generalized Langevin equation with a nonlinear potential of mean force and nonlinear memory friction from a hybrid projection scheme

We introduce a hybrid projection scheme that combines linear Mori projection and conditional Zwanzig projection techniques and use it to derive a generalized Langevin equation (GLE) for a general interacting many-body system. The resulting GLE includes (i) explicitly the potential of mean force (PMF) that describes the equilibrium distribution of the system in the chosen space of reaction coordinates, (ii) a random force term that explicitly depends on the initial state of the system, and (iii) a memory friction contribution that splits into two parts: a part that is linear in the past reaction-coordinate velocity and a part that is in general nonlinear in the past reaction coordinates but does not depend on velocities. Our hybrid scheme thus combines all desirable properties of the Zwanzig and Mori projection schemes. The nonlinear memory friction contribution is shown to be related to correlations between the reaction-coordinate velocity and the random force. We present a numerical method to compute all parameters of our GLE, in particular the nonlinear memory friction function and the random force distribution, from a trajectory in reaction coordinate space. We apply our method on the dihedral-angle dynamics of a butane molecule in water obtained from atomistic molecular dynamics simulations. For this example, we demonstrate that nonlinear memory friction is present and that the random force exhibits significant non-Gaussian corrections. We also present the derivation of the GLE for multidimensional reaction coordinates that are general functions of all positions in the phase-space of the underlying many-body system; this corresponds to a systematic coarse-graining procedure that preserves not only the correct equilibrium behavior but also the correct dynamics of the coarse-grained system.

Non-Markovian rate theory on a multidimensional reaction surface: Complex interplay between enhanced configuration space and memory

A theory of barrier crossing rate on a multidimensional reaction energy surface is presented. The theory is a generalization of the earlier theoretical schemes to higher dimensions, with the inclusion of non-Markovian friction along both the reactive and the nonreactive coordinates. The theory additionally includes the bilinear coupling between the reactive and the nonreactive modes at the Hamiltonian level. Under suitable conditions, we recover the rate expressions of Langer and Hynes and establish a connection with the rate treatment of Pollak. Within the phenomenology of generalized Langevin equation description, our formulation provides an improvement over the existing ones because we explicitly include both the non-Markovian effects along the reaction coordinate and the bilinear coupling at the Hamiltonian level. At intermediate-to-large friction, an increase in dimensionality by itself tends to reduce the rate, while the inclusion of the memory effects increases the rate. The theory predicts an increase in rate when off-diagonal friction terms are included. We present a model calculation to study isomerization of a stilbene-like molecule using the prescription of Hochstrasser and co-workers on a two-dimensional reaction energy surface, employing Zwanzig–Bixon hydrodynamic theory of frequency-dependent friction. The calculated rate shows a departure from the predictions of Langer’s theory and also from the two-dimensional transition state theory.

Likelihood-based non-Markovian models from molecular dynamics

The analysis of complex systems with many degrees of freedom generally involves the definition of low-dimensional collective variables more amenable to physical understanding. Their dynamics can be modeled by generalized Langevin equations, whose coefficients have to be estimated from simulations of the initial high-dimensional system. These equations feature a memory kernel describing the mutual influence of the low-dimensional variables and their environment. We introduce and implement an approach where the generalized Langevin equation is designed to maximize the statistical likelihood of the observed data. This provides an efficient way to generate reduced models to study dynamical properties of complex processes such as chemical reactions in solution, conformational changes in biomolecules, or phase transitions in condensed matter systems.

We introduce a method to accurately and efficiently estimate the effective dynamics of collective variables in molecular simulations. Such reduced dynamics play an essential role in the study of a broad class of processes, ranging from chemical reactions in solution to conformational changes in biomolecules or phase transitions in condensed matter systems. The standard Markovian approximation often breaks down due to the lack of a proper separation of time scales, and memory effects must be taken into account. Using a parametrization based on hidden auxiliary variables, we obtain a generalized Langevin equation by maximizing the statistical likelihood of the observed trajectories. Both the memory kernel and random noise are correctly recovered by this procedure. This data-driven approach provides a reduced dynamical model for multidimensional collective variables, enabling the accurate sampling of their long-time dynamical properties at a computational cost drastically reduced with respect to all-atom numerical simulations. The present strategy, based on the reproduction of the dynamics of trajectories rather than the memory kernel or the velocity-autocorrelation function, conveniently provides other observables beyond these two, including, e.g., stationary currents in nonequilibrium situations or the distribution of first passage times between metastable states.

Particle dynamics in viscoelastic media: Effects of non-thermal white noise on barrier crossing rates

The growing interest in the dynamics of self-driven particle motion has brought increased attention to the effects of non-thermal noise on condensed phase diffusion. Thanks to data recently collected by Ferrer et al. on activated dynamics in the presence of memory [Phys. Rev. Lett. 126, 108001 (2021)], some of these effects can now be characterized quantitatively. In the present paper, the data collected by Ferrer et al. are used to calculate the extent to which non-thermal white noise alters the time taken by single micron-sized silica particles in a viscoelastic medium to cross the barrier separating the two wells of an optically created bistable potential. The calculation-based on a generalized version of Kramers's flux-over-population approach-indicates that the added noise causes the barrier crossing rate (compared to the noise-free case) to first increase as a function of the noise strength and then to plateau to a constant value. The precise degree of rate enhancement may depend on how the data from the experiments conducted by Ferrer et al. are used in the flux-over-population approach. As claimed by Ferrer et al., this approach predicts barrier crossing times for the original silica-fluid system that agree almost perfectly with their experimental counterparts. However, this near-perfect agreement between theory and experiment is only achieved if the theoretical crossing times are obtained from the most probable values of a crossing time distribution constructed from the distributions of various parameters in Kramers's rate expression. If the mean values of these parameters are used in the expression instead, as would be commonly done, the theoretical crossing times are found to be as much as 1.5 times higher than the experimental values. However, these times turn out to be consistent with an alternative model of viscoelastic barrier crossing based on a mean first passage time formalism, which also uses mean parameter values in its rate expression. The rate enhancements predicted for barrier crossing under non-thermal noise are based on these mean parameter values and are open to experimental verification.