Spectral Map: Embedding Slow Kinetics in Collective Variables

The dynamics of physical systems that require high-dimensional representation can often be captured in a few meaningful degrees of freedom called collective variables (CVs). However, identifying CVs is challenging and constitutes a fundamental problem in physical chemistry. This problem is even more pronounced when CVs need to provide information about slow kinetics related to rare transitions between long-lived metastable states. To address this issue, we propose an unsupervised deep-learning method called spectral map. Our method constructs slow CVs by maximizing the spectral gap between slow and fast eigenvalues of a transition matrix estimated by an anisotropic diffusion kernel. We demonstrate our method in several high-dimensional reversible folding processes.


TOC Graphic
2][3][4] However, this challenging task often requires relying on physical intuition or trial and error.For processes exhibiting many timescales, such as glass transitions 5 or crystallization, 6 CVs should capture the slow dynamics of rare transitions between long-lived metastable states.Failure to encode this can hinder the description of the underlying physical mechanisms.9][10][11][12][13][14][15][16][17] In many metastable systems, the long-term dynamics can often be approximately Markovian and modeled as diffusion along CVs in the presence of a free-energy landscape. 189][20][21] This results in the adiabatic elimination of fast variables, which are slaved to the slow variables or their statistics, leading to the effective slow dynamics of the system.However, selecting inappropriate CVs can result in non-Markovian dynamics with long memory effects. 22 this Letter, we propose spectral map, a straightforward method that can construct CVs that arise due to the adiabatic timescale separation in complex systems.Spectral map draws inspiration from recent developments in diffusion maps and parametric dimensionality reduction techniques. 1,15,23Our method is based on training a neural network to embed a high-dimensional system into a few slow CVs without supervision.Maximizing the spectral gap, we parametrize a low-dimensional representation corresponding to the slowest timescales relevant to the physical system.
In the following, we consider a high-dimensional system described by n configuration variables x = (x 1 , . . ., x n ) whose dynamics at temperature T is driven according to a potential energy function U (x), and sampled generally from an unknown equilibrium distribution.
However, if we represent the system using the microscopic coordinates, its dynamics proceeds according to a canonical equilibrium distribution given by the Boltzmann density p(x) = e −βU (x) /Z, where β = (k B T ) −1 is the inverse temperature and Z = dx e −βU (x) is the partition function of the system.
We simplify the high-dimensional configuration space by mapping it into a reduced space z = (z 1 , . . ., z d ) given by a set of d functions of the configuration variables, commonly referred to as CVs, where d ≪ n.Furthermore, we can encapsulate these functions into a parametrizable target mapping: where θ are adjustable parameters.The target mapping ensures that CVs represent the slowest variables in the system.By sampling the system in the CV space, its dynamics follows a marginal equilibrium density p(z) ∝ e −βF (z) , where F (z) is a free-energy landscape: up to an irrelevant constant, and δ(•) is the Dirac delta function.
To estimate the effective timescales characteristic of the system, we can model its dynamics as a Markov chain using kernel functions.However, as the marginal equilibrium density p(z) is unknown, we need a density-preserving kernel appropriate for data sampled from any underlying probability distribution.To achieve this, we employ an anisotropic diffusion kernel: 19,24,25 where ) is a kernel density estimate, and ε is a scale constant.Then, a Markov transition matrix can be constructed by normalizing eq 3: which describes a Markov chain in the CV space given by m kl = Pr {z t+1 = z l | z t = z k } that corresponds to a transition probability from z k to z l in an auxiliary time t.
As the Markov transition matrix is estimated in the reduced space, we consider the effective dynamics rather than the dynamics of the microscopic coordinates of the system.Thus, the Markov chain defined in the CV space is implicitly modeled as following an overdamped Langevin dynamics with the marginal equilibrium density p(z).7][28] However, by selecting CVs that arise from the timescale separation in the system as the reduced space, we can represent the dynamics of microscopic coordinates through the effective dynamics of slow CVs, which is approximately Markovian. 18,20nsequently, to ensure that the timescale separation is evident and that CVs closely follow the slow dynamics of the system, we perform a spectral decomposition of the Markov transition matrix in the CV space.This allows us to calculate its dominant eigenvalues . . .The difference between neighboring eigenvalues is called the spectral gap and measures the degree of the timescale separation between the slow and fast variables: where k > 0 indicates the number of metastable states in the CV space.As such, our objective is to reach the largest spectral gap between the slow and fast variables by adjusting the parameters of the target mapping ξ θ (eq 1).
Let us summarize the algorithm to calculate spectral map: 1.A dataset consisting of N configuration samples in a high-dimensional representation, , is obtained from a simulation; see Figure 1(a).
2. The target mapping ξ θ (x) represented as a neural network is trained through backpropagation by feeding the dataset X and mapping it to d CVs.The objective is  3.The trained target mapping is used to evaluate N samples {z k } N k=1 in the CV space and estimate the corresponding free-energy landscape F (z); see Figure 1(c).
As an initial application of our method, we consider the folding process of a ten-residue protein chignolin in solvent.A 100-µs unbiased molecular dynamics simulation at a temperature of 340 K of this process is obtained from ref 29.As a high-dimensional representation for chignolin, we use pairwise Euclidean distances between its Cα atoms, amounting to n = 45 configuration variables.The training set consists of 5000 samples extracted from the simulation every 2 ns.The training of the target mapping is carried out using k = 2 for 100 epochs with data batches consisting of 100 high-dimensional samples.Once the target mapping is trained, we evaluate all samples collected from the simulation (sampled every 200 ps) to construct the corresponding free-energy landscape.For additional technical details, we refer to the Supporting Information.
Our results are presented in Figure 2. By selecting k = 2 for the spectral gap (eq 5),  is carried out for k = 2 to 7 over 100 epochs, with data batches consisting of 100 highdimensional samples.Finally, the calculated spectral maps and free-energy landscapes are computed using the complete trajectory (sampled every 200 ps).For additional details, see the Supporting Information.
In this example, we aim to determine the optimal number of metastable states (k in eq 5) for maximizing the spectral gap.Our findings, presented in Figure 3, show spectral maps and corresponding free-energy landscapes for k = 2, 3, 4 metastable states.As can be seen in Figure 3(a), selecting k = 2 metastable states results in a distinct folded state and a loosely defined unfolded state.The spectral gap reaches its maximum value of σ = 0.91, indicating satisfactory timescale separation.For k = 3 and 4, the spectral gap decreases to σ = 0.67 and 0.61, respectively, as shown in Figure 3(b-c).We observe that the folded state remains unchanged for k = 3 and 4, while the unfolded state splits into several states due to its heterogeneity, resulting in worse timescale separation.
Additionally, in Figure 3, we can observe that as the number of metastable states (k) increases and the timescale separation worsens, the difficulty in separating the metastable states also results in underestimated free-energy barriers.Based on this, we propose a criterium to determine the optimal value of k.Specifically, we can calculate several spectral maps for increasing values of k and choose the one with the largest spectral gap, corresponding to the maximal separation between slow and fast eigenvalues.Moreover, such timescale separation induces the highest free-energy barriers between metastable states, indicating the superior quality of slow CVs.In the Supporting Information, we provide further details and demonstrate how the spectral gap decreases as the number of metastable states increases to k = 7.
For our final example, we examine the folding and unfolding processes of the BBA protein in solvent, which is simulated by a 200-µs molecular dynamics simulation at a temperature Our results regarding the BBA protein are presented in Figure 4. We focus on the CVs and corresponding free-energy landscape for k = 3 metastable states, as the maximal spectral gap for BBA is virtually the same for k = 2.For the number of metastable states k > 3, the separation between the slow and fast dynamics of BBA worsens (see the Supporting Information).The CVs found by spectral map uncover three distinct states, as shown in Figure 4(a-c).We find that the unfolded state, not the folded state, is the primary and  most populated free-energy basin, suggesting that the folded state is relatively less stable in solvent.This has also been observed in ref 29.Furthermore, we can see that spectral map also finds a β-hairpin state in addition to the folded and unfolded states of BBA.
Representative conformations of BBA in these three states are shown in Figure 4(d).
The capability of identifying the slow kinetics of the BBA protein by spectral map can also be demonstrated by comparing the learned slow CVs to physical descriptors used routinely for reversible folding processes, such as RMSD and the fraction of native contacts.As shown in the Supporting Information, these descriptors calculated in reference to the folded state of BBA are unable to distinguish between the unfolded and β-hairpin metastable states obtained by spectral map.
Interestingly, as can be seen in Figure 4(a), the minimum free-energy path found by spectral map corresponds to the transition directly between the folded and unfolded states.
In contrast, the transition between the folded and β-hairpin states is less likely to be sampled for BBA.This cannot be observed by examining the free-energy profile shown in Figure 4(b), where the transition from the folded to unfolded states proceeds through the β-hairpin state.
In contrast, the slow CVs found by spectral map reveal that the β-hairpin state is primarily reached from the unfolded state, indicating that this state is not an intermediate but a misfolded one.
In this Letter, we have introduced spectral map.Our technique identifies slow CVs from high-dimensional observations obtained from a molecular dynamics simulation.Spectral map employs a procedure to separate slow and fast eigenvalues of a Markov transition matrix that is defined in the CV space.The timescale separation is obtained by maximizing the spectral gap by training a deep neural network.Spectral map allows us to encode information about slow kinetics in the physical system in just a few CVs.Spectral map shares conceptual similarities with several recent techniques, such as spectral gap optimization of order parameters (SGOOP) 20 and the variational approach for Markov processes in a deep learning framework (VAMPnet). 30However, there are notable differences.For example, SGOOP constructs a linear combination of trial CVs and a Markov chain from counting transition between states.VAMPnet employs eigenvalues derived from a spectral decomposition of a lag-time dependent correlation matrix for its loss function.
In contrast, spectral map creates a Markov transition matrix from the anisotropic diffusion kernel, which does not require time series as input.This alternative approach could simplify the process of constructing slow CVs, eliminating the need to select a lag time in Markov state models that require time-dependent data to estimate the transition matrix.
For demonstration purposes, we have shown a method for identifying slow CVs solely from unbiased molecular dynamics simulations.However, only a simple adjustment is needed to expand spectral map and learn from enhanced sampling simulations.This adjustment involves a reweighting procedure for unbiasing Markov transition matrices with statistical weights from a biased simulation.Such a procedure is implemented in reweighted diffusion maps 1,15 and reweighted stochastic embedding techniques. 11,12,15We intend to further explore this approach in subsequent works.
Overall, spectral map shows promise in addressing the challenge of computing slow CV by effectively reducing the dimensionality of complex physical systems and has a large potential for further development.

Figure 1 :
Figure 1: Outline of spectral map.(a) Dataset X in a high-dimensional representation x = (x 1 , . . ., x n ) used to describe the system is taken as an input for the target mapping.(b) Target mapping z = ξ θ (x) is modeled as a neural network that embeds the system in its highdimensional representation to a low-dimensional map spanned by slow CVs z = (z 1 , . . ., z d ).An eigendecomposition of a Markov transition constructed from CV samples is performed (M ψ = λψ).The spectral gap σ is maximized based on the difference between neighboring eigenvalues {λ k } to separate the slow and fast timescales.(c) Trained neural network can be used to evaluate all available high-dimensional samples and calculate the corresponding free-energy landscape F (z).

Figure 2 :
Figure 2: Spectral map of chignolin folding in solvent calculated from a training dataset consisting of 5000 samples extracted every 2 ns from a 100-µs molecular dynamics simulation at a temperature of 340 K.The high-dimensional representation of chignolin is given by n = 45 pairwise Euclidean distances between Cα atoms.(a) Free-energy landscape estimated from the complete trajectory (sampled every 200 ps).The folded state is the main metastable state, with a less populated state representing the unfolded state.The states are separated by a free-energy barrier of around 12 kJ/mol.(b) Free-energy profile along the z 1 CV with z 2 integrated out.(c) The maximal separation between eigenvalues is obtained for the spectral gap of 0.8 at k = 2 (eq 5), with successive eigenvalues close to 0.

Figure 3 :
Figure 3: Spectral map and free-energy landscapes of the folding of trp-cage in solvent calculated from a training dataset consisting of 10000 samples extracted every 2 ns from a 200-µs molecular dynamics simulation at a temperature of 290 K.In total, n = 190 pairwise Euclidean distances between Cα atoms are used as a high-dimensional representation.The spectral gaps for k = 2, 3, and 4 metastable states in (a), (b), and (c), respectively, indicate increasingly worsening separation between the timescales.In each spectral map, the unfolded state of trp-cage splits into additional metastable substates.The spectral gap values up to k = 7 and a comparison of the CVs found by spectral map to CVs often used for reversible folding are shown in the Supporting Information.

Figure 4 :
Figure 4: Spectral map and free-energy landscape of the folding of the BBA protein in solvent calculated from a training dataset consisting of 10000 samples extracted every 2 ns from a 200-µs molecular dynamics simulation at a temperature of 325 K.A high-dimensional representation given by n = 378 pairwise Euclidean distances between the Cα atoms of BBA.(a) Free-energy landscape showing metastable states spanned by CVs calculated by spectral map for k = 3 where the corresponding spectral gap reaches σ = 0.75.(b) Free-energy profile along the z 1 CV with z 2 integrated out.(c) Time series of z 1 showing changes between metastable states during the molecular dynamics simulation.(c) Representative conformations of the BBA protein corresponding to the folded state (F), the unfolded state (U), and the misfolded state in a β-hairpin structure (M).