Likelihood-based non-Markovian models from molecular dynamics.

SignificanceThe 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.

In different branches of science, the interpretation and mathematical modeling of both experimental and computational data requires the analysis of the system dynamics in terms of a reduced set of collective variables (CVs), or order parameters. Prominent examples include chemical reactions in solution, conformational changes in biomolecules, or phase transitions in condensed matter systems. A standard approach is to approximate the evolution of the CVs by an effective dynamics, namely, a closed equation in which the degrees of freedom beyond the CVs (forming the so-called environment or “bath”) do not appear explicitly. Such coarse-grained models not only provide a physical interpretation more accessible to understanding than the full system but also, from a numerical perspective, enable one to recover the desired dynamical properties with long but cheap simulations of the reduced system (while only shorter simulations of the large system are used to determine the effective dynamics).

The most widespread model for this task is the Langevin equation, which can be derived—in some particular cases—from the Hamiltonian dynamics of a small system interacting with a large environment. It describes the evolution of a Markov process, which requires that the decorrelation time of the environment be short compared to the characteristic times of the reduced system. However, many cases do not enter the validity range of this approximation, displaying memory effects (1–8). To go beyond the Markovian approximation, a popular class of processes is given by the generalized Langevin equation (GLE) (9–15)where x(t) is the value of the d-dimensional CV at time t; v(t) is its time derivative; M is an effective mass, is a mean force, usually deriving from a potential V identified with the free energy; K is a memory kernel; and R(t) is a (colored) noise.

This form of the GLE can be motivated from the dynamics of the original full system following the Mori–Zwanzig formalism (9, 16–18), even though it cannot be formally obtained as a controlled approximation of the exact coarse-grained dynamics, since a rigorous derivation generally results in a memory kernel that depends on the CVs (19, 20). Nevertheless, in practice, this simple form is the most widely used effective dynamics. While an analytical derivation of the memory kernel is possible only in a few cases (21), for more general systems, K can be estimated from a data-driven approach. In most cases, the goal is to extract the memory kernel from trajectories of the CV computed with all-atom simulations (4, 22–37).

As already mentioned, the solutions of Eq. are not Markov processes, except when K is the Dirac δ function and R is a white noise. Both for fitting the model and then for generating new trajectories of the effective dynamics, it is convenient to consider the subclass of models where an extended process (x, v, h) is Markovian, with h some hidden auxiliary variables (12, 38–45). Restricting further to the case where the evolution of the hidden variables and the coupling with the observed variables are linear, this leads to an equation of the formwhere are constant matrices and ξ and W are independent standard white noises. This gives a convenient class of models parametrized by the dimension d of h, the corresponding matrices, and the (rescaled) effective force . For equilibrium processes, the coefficients of Eq. are related by the so-called fluctuation–dissipation relation (40). Although we could enforce this condition, thereby reducing the number of parameters, we do not, since we also consider nonequilibrium systems in the following.

Integrating over the hidden variables, we recover Eq. with a memory kernel of the form of a finite Prony series (40, 41)where w and λ are (possibly complex) coefficients of the series derived from the matrices . In principle, on all finite time intervals, any kernel given as the sum of a Dirac function at zero and of a continuous function can be approximated arbitrarily accurately by a sum of the form Eq. . However, in practice d is relatively small, and a memory kernel with, e.g., an algebraic tail can only be approximated on a small time interval (38, 45).

The use of auxiliary variables in the form of Eq. has been abundantly used and studied as it allows efficient integration of GLE Eq. (40, 42, 46), even though other methods exist (28, 29, 33, 47, 48). The estimation of GLE parameters from simulations is an active field of research. The main method consists in a nonparametric estimation of the memory kernel via the Volterra integral equation (12, 25, 34, 37, 39, 46, 49, 50), but other methods have also been proposed (24, 32, 44, 45, 51). In the present work, we 1) introduce a parametric estimator of GLE coefficients, based on a maximum likelihood approach, and 2) show that it allows building faithful coarse-grained models of MD simulations in a cost-effective way (i.e., starting from a relatively small training data set), such that the dynamics is well reproduced.

Data-Driven Approach on Extended Dynamics

In statistics, a standard method to deal with hidden variables is the expectation-maximization (EM) algorithm, which belongs to the category of likelihood maximization algorithms (52, 53). It is frequently used to estimate parameters of time series models in the case of partial or noisy observations of the system, for either hidden Markov models (54) or state-space models (55). A first application in the context of GLE was proposed in ref. 38 to reconstruct the memory kernel in the absence of effective force and under more restrictive conditions than the method presented below.

The algorithm proceeds by alternating steps: in the E step, one determines the conditional probability law of the hidden variables given the observed ones at fixed parameters; in the M step, one optimizes the parameters to maximize the log-likelihood averaged with respect to these conditional laws. In the following we denote as Θ the whole set of parameters estimated after j iterations of the algorithm, which includes the mean force projected on some functional basis (which can be very large in general or reduced if prior knowledge on the system is available), the coefficients of the matrices of Eq. , and, for technical reasons discussed below, the mean value at time 0 of the hidden variables, .

EM Algorithm.

The available data, obtained from all-atom simulations, consists of a set of independent trajectories. For simplicity of the notation, we introduce the algorithm with only one trajectory for some time step and simulation time , the extension to the general case being straightforward. The statistical models we consider are Euler–Maruyama discretizations of Eq. with the same time step , for a fixed dimension d of auxiliary variables h. The state of the system at time will be denoted , and we write as a complete trajectory of the system. Hence, is the value of the known variables since from the choice of the Euler–Maruyama scheme, the velocity can be computed as . In the following we write as the probability density of a variable z; as the conditional probability density of z with respect to u; and, in both cases, to make explicit the value of the parameters if needed.

As the extended system is Markovian, we have for the probability density of a trajectoryand the form of Eq. and of the Euler–Maruyama scheme lead to a Gaussian transition kernel, characterized by its mean μ and variance Σ ().

E Step.

The first step is to compute the conditional law of the hidden variables given the observed variables at the current guess of the parameters, i.e., . Due to the Markovianity of the extended system, it is sufficient to compute the mean and variance of the Gaussian marginal laws for all . Taking advantage of the explicit form of the transition probability (), we apply an iterative predictor–corrector–smoother approach (also known as Kalman filter and Rauch–Tung–Striebel smoother) (56). Starting from the trajectory up to step k – 1, we determine the law of the hidden variable h conditioned on the past information only. We then use the expression of the transition probability to determine the current value of . These are the prediction and correction parts that are run forward on the trajectory, i.e., from k = 0 to k = N (arrow 1 in Fig. 1). The initial guess at k = 0 of uses the measured vector as the mean and an arbitrary variance (identity matrix). Such initial guess could be optimized, but we did not observe any influence on the final results. The second part of the E step, called the smoother part, computes and is run backward, i.e., from k = N to k = 0 (arrow 2 in Fig. 1), which finally gives the required probability law of conditioned to the full observed trajectory. Detailed formulas are presented in .

Fig. 1.

E step. (1) We first predict iteratively the history of h for the whole trajectory, using a predictor-corrector. (2) The values of are then smoothed iteratively backward from the end of the trajectory (see Data-Driven Approach on Extended Dynamics).

M Step.

For any set of parameters Θ, introduce the evidence lower bound after the jth iteration of the algorithm as the expectation with respect to of the log-likelihood of the full trajectory with parameter Θ, namely (see derivation in ),

The M step consists in setting to be the maximizer of this quantity. Notice that due to the particular form of Eq. , is an explicit function of Θ that can be easily optimized as described in .

Full Algorithm.

The algorithm is then run as follows. An initial random or informed guess is taken for the parameters. Such an informed guess could come from a previous execution of the algorithm with a different number of hidden dimensions. From parameters , a new set of parameters is computed through an iteration of E and M steps. Since maximizing the evidence lower bound increases the observed likelihood, the method is iterated until either a prescribed maximum number of EM steps or a convergence criterion is reached.

Assessing the Quality of a Given Model.

The number of hidden dimensions d is an important parameter of the algorithm. It can be chosen using a model validation approach, classically by dividing the set of trajectories between a training and a validation set. However, here we simply compute the optimal parameters for several values of d and compare the predictions of the corresponding models for a number of observable properties, such as the memory kernel, velocity-autocorrelation functions (VACF), or first passage times (FPTs). Similarly, the quality of the model depends on the time step used for the coarse-grained dynamics. This choice depends among other things on the numerical scheme for the propagator. For a given underlying dynamics of the full system, the most accurate choice for the coarse-grained one is to use the same time step , but as a compromise for the amount of data, one can also use (i.e., using only every m step), with m a small integer.

Efficient Sampling of New Trajectories.

Once the model has been optimized by the EM algorithm, it can be used to generate new trajectories in the CV space. Due to their limited computational cost compared to MD trajectories, such synthetic data grant easier access to well-converged average properties, in the form of static and dynamic observables. As an example, in Results the mean FPTs (as well as their probability densities) of a Lennard–Jones (LJ) dimer in a bath are estimated based on the GLE model and compared with the corresponding ones extracted from expensive MD simulations.

Results

We first present the result of the algorithm on a simple yet nontrivial test case with a one-dimensional (1D) system following the extended dynamics of Eq. , with five hidden dimensions and a quadratic potential well , using 20 trajectories of steps and a time step of . The effective force is fitted as a linear function of x. Fig. 2 compares the result of our algorithm to the true memory function that can be computed from Eq. and the one obtained by the Volterra method (Materials and Methods). It demonstrates that the present EM method is able to reproduce the true memory kernel. Furthermore, the parametric structure of the fitted model enforces the decay to zero of the memory kernel, whereas the Volterra method is unstable for times longer than a few tens. Fig. 2 finally shows that the method accurately reproduces the VACF.

Fig. 2.

Equilibrium 1D case. (A) Memory kernel divided by the mass M: the true kernel used to generate the reference trajectories (dark blue line) is compared with the predictions of the Volterra method (cyan solid line, with shaded area indicating uncertainties computed from a bootstrap analysis) and of the present EM method (dash-dotted yellow line). (B) Velocity autocorrelation function, from the reference trajectories (dark blue line) and from new trajectories sampled using the fitted EM model (dash-dotted yellow line).

The algorithm also applies to multidimensional and nonequilibrium systems. This is illustrated in Fig. 3 for a 2D system with two different thermal noises along each axis, with temperatures T  =  1 and T  =  5 and a quadratic potential whose principal axes are not aligned with the x and y axes, leading to nonequilibrium conditions. This setup is inspired by a similar Markovian model used to describe nonequilibrium experiments on cold atoms (57). We run 20 trajectories of steps with a time step of . The effective 2D force is fitted as a linear combination of x and y. The corresponding quadratic potential, illustrated in Fig. 3, is in good agreement with the one used to generate the trajectories. Fig. 3 then shows that the algorithm correctly estimates the memory kernel (in the present case, a simple one with a single hidden dimension for each visible dimension). In particular, the presence of a strong Markovian component is captured by the algorithm but missed by the Volterra method. Finally, the dynamics of the system is well reproduced, as demonstrated for the VACF in Fig. 3.

Fig. 3.

Nonequilibrium 2D case. (A) Locus of for the original potential and the one estimated by the EM algorithm. (B) xx component of the reconstructed memory kernel divided by the mass M: the true kernel used to generate the reference trajectories (dark blue line) is compared with the predictions of the Volterra method (cyan solid line, with shaded area indicating uncertainties computed from a bootstrap analysis) and of the present EM method (dash-dotted yellow line); the left peak represents the Dirac function of Eq. . (C) Velocity autocorrelation function (for the x component of the velocity), from the reference trajectories (dark blue line) and from new trajectories sampled using the fitted EM model (dash-dotted yellow line).

The present approach, based on the reproduction of the dynamics of trajectories rather than the memory kernel or the VACF, conveniently provides other observables beyond these two. Indeed, by generating new trajectories corresponding to the fitted GLE model, one has in principle access to all properties that can be computed from the time evolution of the CVs. As an illustration, Fig. 4 shows for the same nonequilibrium 2D case the stationary probability distribution and the average velocity as a function of the position, estimated using either the initial trajectories used to fit the GLE model (Fig. 4) or the same number of trajectories generated with the latter (Fig. 4). Despite the relatively small number (only 20) of original trajectories used to fit the model and to compute the properties, those computed from the fitted GLE model are in very good agreement with the original ones.

Fig. 4.

Nonequilibrium 2D case: beyond the kernel and the VACF. Stationary probability distribution (colors) and average velocity (arrows) as function of the position for (A) the original dynamics and (B) the GLE model estimated by the EM algorithm. The two ellipses are the same as in Fig. 3 and represent the locus of for the original potential (blue line) and the one estimated by the EM algorithm (dashed yellow line).

As a final illustration, we apply our algorithm to a more realistic 3D system composed of 512 LJ particles at reduced temperature and reduced density . Two of the LJ particles are singled out to form a dimer (24), the others constituting the solvent. The CV of interest is the distance r between the two particles forming the dimer. LJ parameters for all interactions are taken as ϵ  =  1 and σ  =  1 (in LJ units), except between the two particles forming the dimer, with and σ  =  1. The size of the cubic simulation box is , with periodic boundary conditions in all directions. The dynamics is integrated with a time step of (in LJ units) in the microcanonical ensemble (NVE; constant number of particles N, volume V and energy E) using the LAMMPS simulation package (58). We run 20 trajectories with length of 105 time steps, and CV values are extracted every 2 steps.

We fit GLE models defined by Eq. with the EM algorithm for a number of hidden dimensions ranging from three to six. In all cases, the effective force determined from the MD trajectory is used as the single function of the above-mentioned functional basis, so that fitting this part reduces to determining a single prefactor. Our aim is to test the ability of these models to reproduce, in the statistical sense, the properties of the original simulations. In order to check the importance of the hidden variables, we also provide an analysis for a Markovian model, fitted using a maximum likelihood algorithm with 0 hidden dimensions (corresponding to the M step of the above EM algorithm). For each fitted GLE model, we generate 75 new trajectories of length 105 time steps using Eq. , to compute the observable properties and compare them with those obtained from the original set of MD trajectories.

We first compare the stationary distribution for the various GLE models in Fig. 5, which shows the free energy as a function of the r coordinate computed from the histogram of each set of new trajectories (i.e., not the one corresponding to the fitted effective force). The good agreement with the MD free energy profile demonstrates 1) that the coefficient multiplying the model free energy profile of each model is fitted precisely and 2) that the numerical integration of the GLE models is performed accurately. Notice that free energy beyond is affected by the size of the periodic box. The free energy displays two potential wells at and , corresponding to the contact pair (CP, i.e., the dimer) and the solvent shared pair (SSP, with solvent atoms belonging to the solvation shells of both solutes), whose dynamics is investigated below.

Fig. 5.

LJ fluid: two solutes in an explicit solvent. In all panels, results are shown for the reference MD trajectory (dark blue line) and for trajectories generated by the estimated Markov model (dashed cyan lines) and EM models with three to six hidden dimensions (from dark green to purple). Unless specified, all quantities are in LJ units. (A) Free energy (in units of the thermal energy, k) estimated from histograms of the distance r between the two solutes. The various curves are shifted by physically irrelevant constants for clarity. (B) Memory kernel and (C) velocity autocorrelation function. Inset in B shows a zoom on intermediate times. (D) Mean FPT to reach r = 2.0 starting from the distance r. (E) Distribution of the FPT for trajectories starting at and ending at r = 2.0 (Inset shows a zoom on the tail of the distributions, on semilogarithmic scale).

We then consider dynamical observables in Fig. 5, which shows the memory kernel estimated from the Volterra method (4) for MD as well as GLE trajectories, and Fig. 5, which illustrates the VACF (the velocities being computed numerically from positions both in the MD and GLE trajectories). In both cases, increasing the number of hidden dimensions increases the fidelity of the model with respect to the original data. The latter are correctly reproduced for five and six hidden dimensions. The plot also shows the poor quality of the Markovian model, which confirms the necessity of introducing some hidden variables (37).

Finally, we study the transition kinetics between the CP and SSP states, as a stringent test requiring accurate reproduction of both thermodynamic and dynamic properties of the system. Fig. 5 represents the mean FPT to reach the SSP state starting from smaller r distances, whereas Fig. 5 represents the FPT distribution for trajectories starting from the CP state and reaching the SSP state. Clearly, a sufficient number of hidden dimensions (in this case five to six) allows us to quantitatively reproduce the detailed transition statistics. This demonstrates again both the importance of memory effects and the ability of the present algorithm to reconstruct an accurate GLE model.

Conclusions

In this work we addressed the construction of reduced mathematical models of the dynamics of complex molecular systems. Projecting the phase-space trajectories on a reduced set of CVs leads to a powerful framework for the prediction of thermodynamic and kinetic properties of experimental interest. However, the key problems in this context consist of the identification of a suitable dynamical equation and its parametrization. We developed an approach combining GLEs, their numerically efficient representation via Markovian equations including hidden variables, and a powerful machine-learning algorithm borrowed from the field of statistical modeling and data science. Starting from non-Markovian trajectories (e.g., projected all-atom molecular dynamics trajectories in condensed-matter applications), we maximize the likelihood of an extended Markovian model employing the EM algorithm. The advantage of obtaining an explicit parametrization allows for inexpensive sampling of synthetic trajectories, which can be used for the direct computation of quantitative observables (beyond the standard memory kernel and VACF) such as stationary currents in nonequilibrium situations or the distribution of FPTs between metastable states, generally hard to access through atomistic simulations.

Several features distinguish our approach from others existing in the literature. First, the model we optimize includes an explicit parametrization of both the friction and the noise, ensuring consistency between the analysis of the MD trajectories and the generation of new projected trajectories. Second, our method is based on a maximum likelihood procedure, which is well justified from a mathematical perspective. In particular, instead of estimating a nonparametric kernel which is then parametrized (as e.g., in Volterra-based approaches), the parametric model is directly fitted on the data; this should limit the accumulation of errors. Third, we do not enforce equilibrium conditions (such as the fluctuation dissipation theorem) on the model, so that the present approach offers the possibility to investigate nonequilibrium systems. Finally, the present approach readily applies to multidimensional CVs and corresponding matrix memory kernels.

The maximum likelihood approach offers a versatile strategy to implement various extended Markovian models, which could be extended in particular to position-dependent GLEs and higher-order discretization schemes. Overall, the present work provides an efficient way to generate reduced dynamical models for multidimensional CVs, with the same memory kernels as the underlying complex system, enabling the accurate sampling of the long-time dynamics of the latter at a dramatically reduced computational cost.

Materials and Methods

Estimate of the (Potential of) Mean Force.

In the first two examples, the coefficients of the quadratic potentials in the EM method follow from those of the corresponding forces, which are the ones determined numerically along with the parameters related to the memory (). Potentials of mean force in the Volterra method result from quadratic fits of the logarithm of histograms of the position. We obtain the results for the memory kernel with the Volterra approach using the memtools package [https://github.com/jandaldrop/memtools (4)] in the 1D case and the multidimensional version of ref. 39 (see our implementation at https://github.com/HadrienNU/VolterraBasis) in the 2D case.

MD Simulation Details for the LJ Dimer.

The dynamics is integrated with a time step of (in LJ units) in the NVE ensemble with the velocity Verlet algorithm using the LAMMPS simulation package (58). We run 20 trajectories of 105 time steps, and CV values are extracted every 2 steps.

EM Convergence.

Initial values of the parameters are taken randomly. For all examples, we stop the EM iterations if the difference of log-likelihood between two EM steps is less than or if the number of EM steps exceeds 2,000.

Density and Average Velocity for the Nonequilibrium 2D Case.

The density is estimated by kernel density estimation using the positions along the trajectories. The average velocities are estimated conditionally on the positions using kernel regression. The same Gaussian kernel is used in both cases, with a bandwidth of 1. For Fig. 4, 20 new trajectories of steps with a time step of were sampled from the fitted GLE model and compared to the original 20 trajectories of Fig. 4.

Mean FPT Estimation.

The FPT is estimated for molecular dynamics starting by restraining the initial position with a parabolic potential as a function of r using PLUMED (59). Two thousand trajectories are generated from different restrained positions. Gaussian kernel estimates (with bandwidth of ) of the mean FPT as well as the FPT density are then obtained conditioned on the realized starting position. The FPT and mean FPT from the fitted models are computed using 1,500 trajectories per initial value of the distance, again employing kernel estimates.