Machine Learning of Slow Collective Variables and Enhanced Sampling via Spatial Techniques

I. Introduction

Complex systems in chemical physics often exhibit dynamics with multiple temporal scales, characterized by infrequent transitions between long-lived metastable states that occur on timescales orders of magnitude slower than fast molecular motions^[1][2][3]. This significant disparity in timescales is known as timescale separation. Understanding such physical processes depends on our ability to recognize patterns in molecular dynamics (MD) simulations. We typically simplify the dynamics by introducing a set of reaction coordinates, customarily referred to as order parameters or collective variables (CVs)^[4], which are meant to describe it on the macroscopic level. Then, we can estimate a free-energy landscape in CV space, which, to a large extent, is responsible for the thermodynamics and kinetics of physical processes^{[5][6][7][8][9][10][11][12][13]}.

However, the determination of CVs has proved challenging even for simpler systems^[14][15][16]. The most interesting properties of complex processes are often hidden in slow dynamics to which fast variables are adiabatically constrained. Therefore, CVs should describe transitions between states that occur when crossing free-energy barriers significantly higher than thermal energy ($\gg k_BT$). This picture is based on the transition state theory and Kramers’ theory for reaction dynamics, where the reactant and product states are separated by the energy barrier locating the transition state^[17]. Processes such as protein folding^[18][19], crystallization^[20], nucleation^[21], glass transitions^[22][23], aqueous systems^[24], catalysis^[25], or molecular recognition^{[26][27][28][29]} are only a few examples where these characteristics are present and that have frequently profited from such a reduced description.

Due to the rapid development of machine learning (ML) libraries^[30][31], using neural networks has become relatively straightforward and readily available for applications in chemical physics and MD. Interestingly, a fundamental challenge in ML is to develop simple and interpretable representations for complex data^{[32][33][34][35][36][37]}, which closely resembles the task of developing CVs for dynamical systems. Consequently, ML methods have been employed to extract meaningful information from simulations due to their ability to recognize statistical patterns^[38][39]. These techniques can be harnessed to devise algorithms for learning CVs hidden in data to explain the dynamics on the macroscopic scale. A variety of such data-driven methods has been developed at the intersection of statistical physics, MD, and ML. Many of them were recently reviewed^{[40][41][42][43][44][45][46][47][48][4][49][50][51][52]}.

Nonetheless, learning slow CVs remains a challenging task, presenting several difficulties. One key issue is that the quality of CVs is often significantly hampered by the inability to effectively capture longer timescales during standard simulations within a reasonable computing time. This is commonly known as the sampling problem in MD. As such, the construction of training datasets for ML techniques can be problematic as it cannot be known if every state is sufficiently sampled. Additionally, the scarcity of observations between states makes the representation of transition states in reduced space problematic. Enhanced sampling methods can partly alleviate the problem of poor statistics. However, they require correct data reweighting to obtain equilibrium characteristics, often tailored to a particular class of ML algorithms. These problems cause a circular dependency between sampling and learning that poses a major obstacle in developing these techniques.

In this brief review, we will focus on a specific type of ML methods for building slow CVs. Unlike other reviews that cover techniques using trajectories and their time-delayed versions as input to calculate kinetic quantities, such as correlation functions directly, our priority will be on unsupervised techniques that do not rely on temporal characteristics; instead, they estimate kinetics indirectly by analyzing the thermodynamic properties of MD data. These methods aim to learn the reduced space of CVs by capturing spatial characteristics of simulation data encoded in configuration or reduced space, such as the proximity between samples, density estimates, and weights derived from enhanced sampling simulations. We will explore various techniques, including spectral methods such as diffusion maps and their extensions and recently developed algorithms that leverage deep neural networks to learn slow CVs. Lastly, we will discuss potential avenues for future advancements in this field.

II. Background

A. Collective Variables

In statistical mechanics, we consider a system described by the microscopic coordinates $\mathbf{x} = (x_1, \ldots, x_n)$ whose dynamics at temperature $T$ evolves according to a potential energy function $U(\mathbf{x})$. This dynamics can be described by the following overdamped Langevin equation:

where $\beta = 1/k_BT$ is the inverse temperature, $k_B$ is the Boltzmann constant, and $d\mathbf{w}$ is the Brownian motion. The time-evolution of the system results in a canonical equilibrium distribution given by the Boltzmann density:

where $Z = \int d\mathbf{x} e^{-\beta U(\mathbf{x})}$ is the partition function of the system^[53]. We reduce the representation of the system by mapping it into reduced space defined by a set of $d$ functions of the microscopic coordinates, commonly referred to as CVs:

where $d \ll n$. The dynamics of the system in reduced space samples a marginal equilibrium density:

that is defined by weighting each slice through configuration space $\mathbf{x}$, denoted by the delta function $\delta(\cdot)$, with the Boltzmann factors $p(\mathbf{x}) \propto e^{-\beta U(\mathbf{x})}$. The marginal probability $p(\mathbf{z})$ contains information about the free-energy landscape:

Even for simple systems, the free-energy landscape contains many stable states that are separated by barriers much larger than thermal energy, leading to significant timescale disparities in the dynamics.

As summarized in a review by Peters^[54], a general requirement for optimal CVs is to preserve dynamical self-consistency: The dynamics projected onto the free-energy landscape should remain consistent with trajectories sampling configuration space. Taking this apart, we can list more specific characteristics that define optimal CVs:

1. CVs must accurately recognize metastability; that is, distinguish between long-lived metastable states and identify transition states^[55]. Accurate metastability recognition is often difficult to achieve due to the sampling problem^[5][10].
2. CVs need to model reduced dynamics as primarily corresponding to transitions on longer timescales, with the dynamics of fast variables being negligible^[56]. Slow and fast variables should be unmixed in such a way that they induce a significant separation of timescales. Moreover, CVs should not be degenerate; each should describe a different slow mode.
3. CVs need to, preferably, be Markovian for the ability to describe slow dynamics as evolution in the free-energy landscape with configuration-dependent diffusion coefficients^[57][58][59].
4. CVs must be applicable in CV-based enhanced sampling methods (i.e., smooth and differentiable), such as umbrella sampling^{[60][61][62][63]}, metadynamics^{[64][65][66][67]}, or variationally enhanced sampling^[68][69][70], to improve sampling in MD and drive it toward long-timescale processes.

B. Timescale Separation

To illustrate the problem of timescale separation, let us focus on the spectral theory of dynamical systems and reversible Markov processes^[71]. Consider the forward Fokker–Planck equation for time-propagation of a probability distribution $p(\mathbf{x},t)$: $\partial p/\partial t = -Lp$, where $L$ is the generator of a Markov process. Through a change of variables, the Fokker–Planck equation can be solved^[72]. The solution to this equation can be written in closed form as an eigenfunction expansion^[73]:

where $a_k$ are coefficients and $t$ is a time variable. For time $t \to \infty$, the solution of the forward equation converges to the equilibrium Boltzmann distribution $p(\mathbf{x})$.

Under general conditions, the generator of the diffusion process $L$ has a discrete eigenspectrum of eigenvalues $\mu_k$, and the corresponding eigenfunctions $\varphi_k(\mathbf{x})$. The zeroth eigenfunction is the equilibrium density $\varphi_0(\mathbf{x}) \propto e^{-\beta U(\mathbf{x})}$ with the eigenvalue $\mu_0 = 0$. The eigenvalues are non-negative and sorted in increasing order:

The dominant eigenvalues of the Markov generator decay exponentially and are linked to the slowest relaxation timescales in the system. Each eigenvalue can be matched with an effective timescale $t_k = 1/\mu_k$. In systems with timescale separation, only a few slow processes related to rare transitions between metastable states remain. As a result, the eigenspectrum of $L$ has a spectral gap, i.e., the largest difference between eigenvalues $\mu_{k+1}$ and $\mu_k$. This implies that the eigenvalues much lower than $\mu_{k+1}$ can be neglected as they correspond with rapid fluctuations within states and decay much faster, leading effectively to $k$ slow processes (Fig. 1).

The spectral properties of reversible Markov processes can be related to the concept of metastability^[74]. Although this relation can be understood intuitively, Gaveau and Schulman^[75][76], drawing on the extensive work of Davies^[77][78][79], developed a spectral definition of metastability. They formally showed that dominant and nearly degenerate eigenvalues are related to metastable timescales. This concept relies on the presence of the spectral gap. If an eigenvalue is nearly degenerate, the equilibrium distribution separates into metastable states with infrequent transitions between them. Conversely, eigenvalue degeneracy exists if the equilibrium density breaks into metastable states separated by a free-energy barrier much larger than thermal energy. The eigenfunctions related to the dominant eigenvalues are linked to distributions that remain stable longer than transient processes. Furthermore, sign changes in these eigenfunctions indicate transitions between metastable states. The theory is summarized in a monograph by Bovier and Den Hollander^[80].

C. Enhanced Sampling

Acquiring an informative training dataset from unbiased MD trajectories is a crucial challenge. These trajectories need to spontaneously and repeatedly cross over all significant free-energy barriers in the system. However, the metastability leads to kinetic entrapment in a single state, making transitions between metastable states rare. To alleviate this issue, enhanced sampling methods can be used to improve sampling efficiency^{[6][7][8][9][10][11]}.

Enhanced sampling methods that require CVs to improve sampling are based on employing a nonphysical bias potential. To such methods, we can include umbrella sampling introduced by Torrie and Valeau^[60], adiabatic biasing force^[81], adiabatic free-energy dynamics^[82], metadynamics proposed by Laio and Parrinello^[64] and improved to the well-tempered variant by Barducci et al.^[65], mean-force dynamics^[83], or variationally enhanced sampling^[68]. Biasing the system can cause the probability distribution of collective variables (CVs) to significantly deviate from equilibrium, resulting in sampling according to a biased distribution:

where $V(\mathbf{z},t)$ is a time-dependent bias potential. To calculate equilibrium properties, such as free-energy landscapes, the bias must be reverted during postprocessing. This is customarily done by reweighting, where each sample is associated with a statistical weight to counter the effect of biasing. In general, the weights are given by the likelihood ratio between the equilibrium and the biased probability distributions (Eq. 8):

For methods using a quasi-stationary bias potential^[60][68][84] (e.g., umbrella sampling), or when the simulation is converged and the bias does not change significantly, the weights are given as:

In contrast, in metadynamics^[65], the bias potential changes over time and requires accounting for a time-dependent offset^[85]. Thus, the functional form of weights may vary depending on an enhanced sampling method and a reweighting algorithm^{[86][87][84][88]}. A summary of such methods was recently published by Kamenik et al^[89].

To efficiently sample and drive complex physical processes, high-quality CVs are required for biasing. However, learning CVs demands using exhaustively sampled data. This problem creates a challenging circular dependency, which is referred to as the “chicken-and-egg” problem^[45]. Advances in the determination of CVs help address this problem and contribute to the development and implementation of enhanced sampling methods.

III. Spatial Learning

Due to recent extensive advancements in data-driven temporal methods^{[90][91][92][69][93][94][95]}, there are numerous reviews summarizing this topic^{[45][46][49][50][51][52]}. In this work, however, we consider techniques that are “spatial,” i.e., algorithms for learning slow CVs that do not need to exploit temporal information in MD simulations. We can describe spatial techniques as those that rely on pairwise relations between samples in the dataset (usually through a distance metric) instead of counting transitions within a specified lag time. The development of such techniques can be traced back to the work of Shi and Malik^[96] on image segmentation and the classic Laplacian eigenmaps introduced by Belkin and Niyogi^{[97][98][99][100]}; and is closely related to graph spectral theory^[101] based on graphs, kernels, and random walks^{[102][103][104]}.

The primary difference between spatial and temporal techniques lies in how kinetics is estimated. Spatial techniques estimate kinetics indirectly by analyzing the thermodynamic characteristics of MD data, such as equilibrium probabilities, in contrast to temporal techniques. Additionally, in spatial techniques, we assume that MD data closely approximates overdamped Langevin dynamics (see Sec. II.A). For these reasons, we can refer to these methods as thermodynamics-informed learning.

A. Anisotropic Kernels

The core of most spatial learning methods involves establishing similarity between samples, typically through a distance metric and a kernel^[105]. For example, Laplacian eigenmaps construct a Gaussian kernel to model relations between $N$ samples in a dataset $X = \{\mathbf{x}_k\}_{k=1}^N$^{[97][98][99][100]}:

where $\varepsilon > 0$ is a scale parameter. This kernel is then used to define a Laplacian matrix and parametrize reduced space using its eigenvectors. However, methods that use a Gaussian kernel, such as Laplacian eigenmaps, cannot be used to compute slow CVs as their construction implicitly assumes that data is distributed uniformly. As the equilibrium density is often far from uniform, Laplacian eigenmaps have not seen many applications for analyzing trajectories. However, they are often used as a baseline for developing more advanced techniques.

Based on Laplacian eigenmaps, Coifman et al.^[106] developed the diffusion map algorithm that is especially suited for learning the reduced space of slow CVs. Diffusion maps use a density-preserving kernel for data sampled from any underlying probability distribution. For this, an anisotropic kernel is constructed on the dataset $X$^[107]:

where $\varepsilon$ is a scale parameter, $\rho(\mathbf{x}_k) = \sum_l G_\varepsilon(\mathbf{x}_k, \mathbf{x}_l)$ is a density estimate that allows us to include information about non-uniformly sampled data into the kernel, and $\alpha \in [0,1]$ is the anisotropic diffusion constant. Next, a Markov transition matrix is constructed by row-normalizing $K$:

to build a discrete Markov chain on the data:

that expresses a transition probability between $\mathbf{x}_k$ and $\mathbf{x}_l$. Note that this construction does not depend on the physical time. The local scale parameter $\varepsilon$ plays an important role in determining the quality of slow CVs, as it defines the scale within which the relation between two samples contributes to the Markov transition matrix.

Depending on the anisotropic diffusion constant $\alpha$, several kernel normalizations are available, which can change the long-time convergence of the Markov chain to a particular operator. This group of constructions is known as anisotropic diffusion maps^{[106][108][109][107][110]}. For example, with $\alpha = 1/2$, the Markov chain approaches the time asymptotics of the system by describing the dynamics by the Fokker–Planck anisotropic diffusion with the potential $U(\mathbf{x})$. As such, this normalization is commonly used to extract information from MD trajectories. Two other frequently considered values are $\alpha = 0$ and 1. The former results in the classical normalized graph Laplacian, while the latter yields the Laplace-Beltrami operator with a uniform probability density^{[106][108][109][107][110]}.

The advancements of the diffusion map algorithm and anisotropic Markovian kernels often involve using a kernel that captures more aspects of the data. For instance, self-tuning local kernels were introduced by Zelnik-Manor and Perona^[111]. Following this works by Rohrdanz et al.^[112] and Zhang et al.^[113][114] demonstrated that estimating the scale parameter as configuration-dependent $\varepsilon(\mathbf{x}_k)\varepsilon(\mathbf{x}_l)$, where each term can be calculated as the distance between $\mathbf{x}$ and its $n$-th nearest neighbor, improves the overall quality of slow CVs^[40]. A more general method for computing the local scale parameters was later proposed by Berry et al.^[115][116]

In works by Dsilva et al.^[117][118] and Singer et al.^[119], it was proposed to use a heterogeneous Gaussian kernel to improve properties of the resulting CVs. Instead of using the Euclidean distance, this kernel introduces a Mahalanobis-like distance, which incorporates a covariance matrix. The implication of this is that the Mahalanobis kernel, by including the correlations in the dataset, can be used to remove the effect of observing the underlying space through a complex nonlinear function^{[119][117][118]}:

where the squared Mahalanobis distance is:

The local covariance matrix $\Sigma_k$ can be estimated as a sample covariance matrix at configuration $\mathbf{x}_k$ in its immediate neighborhood^{[120][119][117]} and $\dagger$ denotes a pseudo-inverse (as $\Sigma$ can be rank-deficient).

Subsequently, Berry and Sauer^[121] developed a generalization of diffusion maps to local kernels by introducing diffusion and drift terms in the distance metric, which should be additionally computed from the data^[122][123]. It was shown by Berry et al.^[124] that it is possible to improve anisotropic kernels by including Taken’s delay coordinates in datasets, especially when observations are scarce. Diffusion map was also embedded in a framework for coarse-graining and clustering^[125].

B. Reweighted Transitions

The concept of reweighting transition probabilities is crucial when using enhanced sampling algorithms to build the Markov transition matrix and, thus, CVs. A Markov chain constructed from biased data does not converge to the equilibrium density given by the Boltzmann distribution^[126][127]. This bias affects the Markov chain and leads to incorrect density and geometric relations between samples, which can result in reduced space that does not accurately represent the characteristics of the data. Reweighting pairwise probabilities counters the bias from the Markov matrix, yielding the unbiased Markov. While learning biased CVs can still be used to analyze, speed up, and drive the sampling of rare events^{[113][128][129]}, the necessity of a reweighting algorithm becomes apparent when we seek to restore the equilibrium properties of the system and compute slow CVs.

The initial approach to learning unbiased CVs from enhanced sampling simulations with the diffusion map algorithm was proposed by Ferguson et al.^[130], in which each configuration is weighted based on its importance in umbrella sampling simulations. A symmetric weighted Gaussian kernel was used by Zhang et al.^[114] to learn CVs from multiple metadynamics simulations. Building on the local kernels introduced by Berry and Sauer^[121], Banisch et al. and Trstanova et al. devised a general approach to reweighting transition probabilities based on target measure reweighting^[131][132]. This approach was later employed in works by Evans et al., where diffusion map with the Mahalanobis distance is constructed in $\mathbf{z}$ space^[131][132].

Zhang and Chen^[133] derived an alternative technique for reweighting, which Rydzewski et al.^[127] later generalized to multiple algorithms employing Markov transition kernels. They demonstrated that the anisotropic diffusion kernel as can be unbiased as:

where a transition reweighting factor $r_{kl} = w_k w_l$ incorporates importance weights from enhanced sampling simulations and $\rho$ are reweighted density estimates:

A detailed derivation with possible approximations is given by Rydzewski et al.^[127] As explained in Sec. II.C, the form of weight depends on the employed enhanced sampling and reweighting techniques^{[86][85][87][84][89]}.

Several approximate transition reweighting factors can be obtained depending on the scaling of the long-time asymptotics of the kernel with the constant $\alpha$^[127]. This kind of transition reweighting can be used for diffusion maps^[127] and deep learning^{[133][126][127]}. We refer to the review by Rydzewski et al.^[51] for a detailed discussion.

This idea was recently explored by Rydzewski^[134], who demonstrated that this form of transition reweighting in diffusion maps can be employed as a feature selection pipeline for further dimensionality reduction. This is done by leveraging the idea that the partial selection of variables should have a similar eigenspectrum to configuration space. This extension can provide an interpretable and explainable description by selecting physically important CVs for the given process^[134].

For a more general approach to dynamical transition reweighting, not limited to unbiasing transition probabilities in spatial techniques, see reviews by Chen and Chipot^[50], which discusses many reweighting methods for temporal techniques and Keller and Bolhuis^[135], where reweighting is examined from the perspective of Markov state models.

C. Eigendecomposition

In learning algorithms that use a few eigenvectors of the Markov transition matrix to span $\mathbf{z}$ space, a mapping into $\mathbf{z}$ space is obtained by solving an eigendecomposition problem:

where $\lambda_k$ and $\psi_k$ are the eigenvalues and corresponding eigenvectors of the Markov transition matrix $M$, respectively. As explained in Sec. II.B, as a result of the existence of the spectral gap between neighboring eigenvalues $\lambda_k$, slow CVs can be approximated by the following truncated mapping:

where $d$ is the dimension of $\mathbf{z}$ space. The eigenvalues of the Markov transition matrix $M$ are (sorted in non-ascending order):

where the eigenvalue $\lambda_0$ corresponds to the equilibrium distribution of the Markov chain given by the eigenvector $\psi_0$. The dominant eigenvalues related to the slowest relaxation timescales in the system^[74] and the fast eigenvalues have a negligible contribution to slow CVs. In the case of anisotropic diffusion maps, the eigenvalues $\lambda_k$ are related to the eigenvalues of the Fokker–Planck generator $\mu_k$ by the relation $\lambda_k = e^{-\mu_k}$.

Several techniques use the mapping provided by diffusion maps as an initial guess to improve slow CVs iteratively. For instance, the eigenvectors of the Markov transition matrix $M$ can serve as a basis to approximate kinetic quantities such as the transfer operator. This approach was exploited in works by Boninsegna et al.^[136], Noe and Clementi^[137][138], and more recently by Thiede et al. using a Galerkin approximation^[139].

Algorithms that use an eigendecomposition to construct $z$ space require an out-of-sample extension to map samples outside of the dataset. Specifically, for diffusion maps the Nyström extension^[140][141], Laplacian pyramids^[117], and geometric harmonics^[142][143] interpolators were used. A detailed analysis of out-of-sample algorithms was published by Bengio et al^[144].

D. Reweighted Stochastic Embedding

Reweighted stochastic embedding (RSE) is a recent framework for the parametric learning of slow CVs, introduced by Rydzewski et al.^[127], which employs algorithms to construct unbiased Markov transition matrices with transition reweighting (Sec. III.2), allowing for the estimation of CVs from data collected in enhanced sampling simulations. Building on the work of van Maaten, Hinton, and Roweis^{[145][146][147][148]}, RSE optimizes a loss function to learn the mapping to reduced space. Specifically, it projects samples into $\mathbf{z}$ space using a neural network, while ensuring that the statistical distance between transition matrices estimated in both configuration space and $\mathbf{z}$ space is minimized (Fig. 2).

The first technique of this framework is stochastic kinetic embedding (StKE), which was proposed by Zhang and Chen^[133]. StKE combines modeling a slow manifold with parametric dimensionality reduction, building upon the reweighted anisotropic diffusion kernel. As such, StKE can learn slow CVs from biased data sampled in enhanced sampling simulations. In addition, it uses an iterative procedure incorporating temperate-accelerated MD^[149][150] to alleviate the circular dependency^[133], allowing the use of this algorithm on the fly in atomistic simulations^{[133][49][127][151]}. Subsequently, Rydzewski and Valsson introduced a RSE technique called multiscale reweighted stochastic embedding (MRSE)^[126] that shares similarities with StKE^[127][51]. The main difference between StKE and MRSE is, as in many methods discussed in this review, boils down to using other kernels to estimate Markov transition matrices. In MRSE, the process of constructing unbiased transition probabilities from enhanced sampling simulations involves adaptively estimating a kernel in $\mathbf{x}$ space based on information theory principles. In contrast, StKE employs a fixed anisotropic diffusion kernel, as used in diffusion maps (see Sec. III.2). This topic is discussed in detail in the review by Rydzewski et al^[51].

RSE employs building transition matrices in both $\mathbf{x}$ and $\mathbf{z}$ spaces (Fig. 2). As with many neutral network-based techniques for learning CVs, $\mathbf{x}$ space can comprise variables other than the microscopic coordinates, which are called features or descriptors. The transition matrix $M$ constructed in $\mathbf{x}$ space remains constant throughout learning, while the matrix $Q$ in $\mathbf{z}$ space is adjusted depending on a neural network that performs dimensionality reduction, i.e., $f_w(\mathbf{x}) = \mathbf{z}$. Most generally, in RSE, a weighted Gaussian mixture is used to construct the transition matrix in $\mathbf{x}$ space^[126]:

where the sum goes over scale parameters. In $\mathbf{z}$ space, the transition matrix can be given, for example, by a $t$-distribution kernel^[126]:

RSE minimizes the Kullback–Leibler divergence^[152] to learn CVs, which can be interpreted as a “distance” between probability distributions. Thus, after the training converges, the transition probabilities in both spaces should be approximately equal. More details about these algorithms can be found in Refs.^[127][126].

E. Spectral Map

The first technique devised to maximize timescale separation to find CVs in complex systems was proposed by Tiwary and Berne^[153] and subsequently expanded^{[154][155][156][157]}. Their technique, called spectral gap optimization of order parameters (SGOOP), is based on constructing a transition matrix using the principles of the maximum caliber framework^[158]. As opposed to the techniques reviewed here, SGOOP explicitly uses time information to construct slow CVs.

A recent unsupervised statistical learning technique for learning slow CVs that is also based on maximizing timescale separation is spectral map, developed by Rydzewski^[159]. It is modeled using an overdamped Langevin diffusion in $\mathbf{z}$ space^[160]. Spectral map proceeds by mapping the dynamics into $\mathbf{z}$ space using a neural network and constructing a Markov transition matrix by row-normalizing the anisotropic diffusion kernel (Eq. 12), however, from data in $\mathbf{z}$ space:

where $\mathbf{z} = f_w(\mathbf{x})$ given by the neural network with learnable parameters $w$. The transition matrix is then spectrally decomposed to estimate the degree of timescale separation from the spectral gap in its eigenspectrum. The spectral gap is used as a score function for the neural network and maximized during learning:

where $\lambda_k$ are the eigenvalues of $Q$ sorted in decreasing order and $m$ is the number of metastable states in the system. As the spectral gap is maximized, $\mathbf{z}$ space is adjusted accordingly by improving the parameters of the neural network. At the end, $\mathbf{z}$ space corresponds to slow CVs. A simplified diagram spectral map and comparison to RSE is given in Fig. 2.

Rydzewski and Gökdemir^[161] showed that maximizing timescale separation in spectral map results in the dynamics in $\mathbf{z}$ space with Markovian characteristics. In their work, it was shown that it is possible to construct a high-quality Markov state model based on the learned slow CVs and estimate kinetics accurately. In another work, Rydzewski showed that the framework can be easily extended for learning the transition state ensembles^[160] (Fig. 3), which is demanding for complex systems due to the scarcity of transitions between states^{[162][163][164]}.

Using the transition state ensemble to count transition paths^[165], Rydzewski^[160] showed that a slow CV learned by spectral map closely approaches a Markovian limit for overdamped Langevin dynamics^[59]. Moreover, it was illustrated that spectral maps can estimate the quality of the reduced representations with commonly used physical descriptors by comparing their spectral gaps. It was demonstrated that spectral map can be used to construct interpretable reaction coordinates for protein folding with a linear model instead of a deep neural network, and they are slower than the fraction of native contacts or end-to-end distance^[160].

F. Enhanced Sampling via Neural Networks

After the training procedure, a neural network representing CVs can be used for the purposes of ehnanced sampling. To bias such a neural network, a biasing force must be applied in CV space. This force is equal to the negative derivative of the biasing potential with respect to the CVs, which can be estimated using the chain rule:

where the second term on the right-hand side is automatically computed through backpropagation. By accumulating the biasing potential in CV space, the neural network can be used to push the system out of local minima. Such CVs, in the form of a neural network, can be integrated into several advanced MD simulation codes, such as PLUMED^{[166][167][168][169][170]}.

We want to underline that there might be more requirements for slow CVs represented by a neural network (not only limited to spatial techniques). An often overlooked issue that can harm the convergence of biasing methods is that the neural network may learn a function where $\nabla_\mathbf{x} f(\mathbf{x}) \approx 0$ in basins. According to Darve et al.^[81], biasing a CV can be imagined in terms of an object that is pulled or pushed, where the CV has a “mass” attached to it that is related to the inverse of $\nabla_\mathbf{x} f(\mathbf{x})$. Consequently, applying bias to neural networks with $\nabla_\mathbf{x} f(\mathbf{x}) \approx 0$ in energy minima might be inefficient due to the large mass and lead to numerical stability issues in MD simulations.

IV. Summary

Overall, we think further research in spatial techniques will follow by carefully incorporating more thermodynamical information into ML. Due to rapid developments in physics-informed algorithms, we expect that the primary effort will be directed toward solving the problem of constructing interpretable and explainable reaction coordinates for complex systems in chemical physics.

To address this issue, we can examine the theoretical progress in modeling slow dynamics in the context of timescale separation in CV space^{[171][172][173][174]}. By investigating slow dynamics using overdamped Langevin dynamics in a free-energy landscape with configuration-dependent diffusion coefficients, we can propose a Markovian interpretation of the physical process. The diffusion tensors, which depend on the coordinates, are important for reduced dynamics and can impact the free-energy landscape by altering transition states and barrier height^{[175][176][177][178][179]}. To account for this in spatial techniques, we can incorporate information about them in anisotropic kernels and transition matrices. Additionally, analyzing spatial techniques from the perspective of spectral graph theory^[101], especially the long-term behavior of Markov chains, the asymptotic rate of convergence to equilibrium, and mixing rates^[180][181], can lead to improvements.

For spatial techniques to learn from enhanced sampling simulations, slow CVs should be computed using unbiased Markov matrices through a transition reweighting algorithm, such as those presented in the review, to capture equilibrium information accurately. It would be interesting to explore the relationship between the reweighting of Markov transition matrices and dynamical path reweighting, for example, based on the Girsanov theorem^{[182][183][184]}. To improve sampling and drive it toward complex physical processes, spatial techniques can be extended with a general iterative learning-sampling framework where rounds of learning slow CVs (including reweighting) are followed by biasing using an enhanced sampling technique. Such iterative approaches have already been implemented using ML to learn from MD simulations^{[92][185][186][133][187][93][188][169][189]}.

Finally, we underline that apart from spatial techniques, many others can be used to study complex processes in the fields of chemical physics and MD^{[190][191][192][193][194][195][196][197][198][199][200][201][202][203][204][205][206][207][208][209][210][211][212]}. A detailed introduction to such ML methods can be found in recent reviews^{[40][41][42][43][44][45][46][47][48][49][50][51][52]}. We think, however, that recent results in spatial techniques for learning slow CVs are worthy of further development and could provide a valuable alternative to temporal techniques for understanding the physics of complex systems.

Statements and Declarations

Acknowledgements

The research was supported by the National Science Center in Poland (Sonata 2021/43/D/ST4/00920, “Statistical Learning of Slow Collective Variables from Atomistic Simulations”). J. R. acknowledges funding the Ministry of Science and Higher Education in Poland. The authors acknowledge insightful feedback from Haochuan Chen, Luke Evans, Luigi Bonati, and Omar Valsson.

Conflicts of Interest

The authors have no conflicts to disclose.

Author Contributions

Tuğçe Gökdemir: Conceptualization (equal); Supervision (equal); Writing - original draft (equal); Writing - review & editing (equal). Jakub Rydzewski: Conceptualization (equal); Supervision (equal); Writing - original draft (equal); Writing - review & editing (equal).

Data Availability

Data sharing is not applicable to this article as no new data were created or analyzed in this study.