Deflation reveals dynamical structure in nondominant reaction coordinates

TL;DR

This paper introduces deflation into molecular kinetics analysis methods to better identify and model non-dominant reaction coordinates, enabling more targeted understanding of complex molecular processes.

Contribution

The paper presents a novel deflation technique integrated with variational methods to focus on non-slowest dynamical modes in molecular systems.

Findings

Deflation improves the informativeness of nondominant reaction coordinates.

dVAMP effectively isolates specific dynamical processes in protein folding.

Deflation can obscure undesired slow processes from kinetic models.

Abstract

The output of molecular dynamics simulations is high-dimensional, and the degrees of freedom among the atoms are related in intricate ways. Therefore, a variety of analysis frameworks have been introduced in order to distill complex motions into lower-dimensional representations that model the system dynamics. These dynamical models have been developed to optimally approximate the system's global kinetics. However, the separate aims of optimizing global kinetics and modeling a process of interest diverge when the process of interest is not the slowest process in the system. Here, we introduce deflation into state-of-the-art methods in molecular kinetics in order to preserve the use of variational optimization tools when the slowest dynamical mode is not the same as the one we seek to model and understand. First, we showcase deflation for a simple toy system and introduce the deflated variational approach to Markov processes (dVAMP). Using dVAMP, we show that nondominant reaction coordinates produced using deflation are more informative than their counterparts generated without deflation. Then, we examine a protein folding system in which the slowest dynamical mode is not folding. Following a dVAMP analysis, we show that deflation can be used to obscure this undesired slow process from a kinetic model, in this case a VAMPnet. The incorporation of deflation into current methods opens the door for enhanced sampling strategies and more flexible, targeted model building.

Equations

$$\mathbf{X}$$

$$\mathbf{Y}$$

$${\mathbf{C}}_{\tau }\mathbf{\Phi }={\mathbf{C}}_0\mathbf{\Phi \Lambda },$$

$${\mathbf{C}}_{\tau }$$

$${\mathbf{C}}_0$$

$${\mathbf{C}}_{00}$$

$${\mathbf{C}}_{0\tau }$$

$${\mathbf{C}}_{\tau \tau }$$

$$\tilde{\mathbf{X}}$$

$$\tilde{\mathbf{Y}}$$

$$\tilde{\mathbf{K}}=({\tilde{\mathbf{X}}}^{\top }\tilde{\mathbf{X}}{)}^{-1}{\tilde{\mathbf{X}}}^{\top }\tilde{\mathbf{Y}},$$

$$\tilde{\mathbf{K}}$$

$${\mathbf{U}}^{\top }{\mathbf{C}}_{00}\mathbf{U}$$

$${\mathbf{V}}^{\top }{\mathbf{C}}_{\tau \tau }\mathbf{V}$$

$${\xi }_c^{\prime }$$

$${\omega }_c^{\prime }$$

$$U(x)$$

$$P_2(n)$$

$$P_3(n)$$

$${\mathbf{X}}_{\text{new}}={\mathbf{X}}_{\text{original}}-{\xi }_c[{\mathbf{X}}_{\text{original}}^{\top }{\xi }_c({\xi }_c^{\top }{\xi }_c{)}^{-1}{]}^{\top },$$

$$\text{VAMP-}r\equiv \sum _{i=1}^k{\sigma }_i^r.$$

$${\varphi }_c$$

$$({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}({\mathbf{X}}^{\top }\mathbf{Y})({\mathbf{Y}}^{\top }\mathbf{Y}{)}^{-1}({\mathbf{Y}}^{\top }\mathbf{X}),$$

$$\mathbf{Y}{\mathbf{C}}_{\tau \tau }^{-\frac{1}{2}}$$

$$\mathbf{Y}$$

$${\mathbf{K}}_f$$

$${\mathbf{K}}_b$$

$$({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}({\mathbf{X}}^{\top }\mathbf{Y})({\mathbf{Y}}^{\top }\mathbf{Y}{)}^{-1}({\mathbf{Y}}^{\top }\mathbf{X})$$

$$({\mathbf{Y}}^{\top }\mathbf{Y}{)}^{-1}({\mathbf{Y}}^{\top }\mathbf{X})({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}({\mathbf{X}}^{\top }\mathbf{Y})$$

Full text

I Introduction

Molecular dynamics (MD) simulations provide a detailed view of a system that describes its atomistic dynamics: how a protein folds or misfolds, how a ligand binds, or how conformational changes occur. However, interpretable dynamical motifs are encoded in collective variables that must be computed from the simulation output. A wealth of research has been dedicated to extracting thermodynamic and kinetic information from time-series simulation data, often under the assumption that a Markovian lag time can be defined that is long enough for transitions to be memoryless but short enough to resolve the dynamics of interest (see Ref. Husic and Pande, 2018 for a review).

When models are Markovian, the future state of the system depends only on its present state, and not on its history. Thus, the system dynamics can be completely described by time-lagged pairs of data points. The true Markovian dynamics are described by some unknown dynamical operator whose eigen- or singular functions can be approximated in a data-driven way using the time-lagged representation Schütte and Sarich (2015).

Recently, two variational principles were developed for computing the optimal approximations to the dominant eigen- or singular functions of the unknown operator, which represent the slow dynamical modes of the system. First, Noé and Nüske (2013) introduced the variational approach to conformational dynamics (VAC) for reversible dynamics; more recently, Wu and Noé (2017) derived the variational approach to Markov processes (VAMP) for irreversible dynamics, which is a generalization of the VAC.

When analyzing time-lagged data, suitable feature representations must be chosen for the data, which often involve transformations of the original spatial coordinates such that the features are rotation- and translation-invariant. For proteins, common features include backbone dihedral angles and pairwise contact distances. From the time-lagged feature representation, VAMP can be used in two ways. First, VAMP yields the optimal linear model from those features. Second, VAMP can be used to compare the approximation quality of nonlinear models of the features, such as Markov state models (MSMs) Husic and Pande (2018), Koopman models Wu et al. (2017), or VAMPnets Mardt et al. (2018). A method for comparing feature choices directly has also recently been presented Scherer et al. (2019).

The introduction of variational methods into kinetic model construction has shifted the modeling community away from heuristic parameter choices and towards objective model rankings. However, the VAC and the VAMP optimize the global kinetics, which are characterized by the slowest dynamical modes of the system. It is a common criticism of these methods that the slow modes are not always the same as the modes of interest Husic and Pande (2018); Scherer et al. (2019). For example, the optimization method introduced by Banushkina and Krivov (2015) identifies a rare event as the optimal reaction coordinate instead of the desired mode, and the method derived by McGibbon, Husic, and Pande (2017) applied to a protein folding dataset finds an irrelevant distance contraction instead of the folding process.

Thus, applying a variational method when modeling a time-series dataset will optimize the approximation of the slowest modes, regardless of whether they are interesting. While in some cases it can be straightforward to modify the kinetic model a posteriori in order to remove the undesired modes after model construction, it becomes messier if not impossible to remove the same slow mode from every candidate model when multiple models are ranked and their modes are nonequivalent.

Furthermore, constraints of commonly used kinetic models require that the slow modes they identify are mutually orthogonal, which places restrictions on the modeling of the process of interest if it is not the slowest process. Finally, researchers sometimes find that the reaction coordinates of desired modes feature structures also present in slower reaction coordinates of undesired modes, thus obscuring the nondominant process.

Inspired by these limitations, in this work we consider the use of deflation to address two problems at the frontier of molecular kinetics analyses. First, we investigate whether the incorporation of deflation into the variational approach can increase the interpretability of nondominant reaction coordinates. Deflated reaction coordinates are no longer required to be orthogonal to previously computed coordinates; thus, if previously computed reaction coordinates are undesired, deflation may prevent the contamination of interesting modes with uninteresting ones. Second, we examine whether the deflation of a specific process can be used construct a new kinetic model that optimizes the description of only the desired slow dynamics. In this case, deflation is used to create a new basis, from which multiple models can be constructed and ranked with the guarantee that they are not optimizing for the accuracy of the uninteresting slow mode. From our results, we suggest several next steps for sampling and model construction.

II Theory

II.1 Time-lagged independent component analysis

The bio- and chemical physics communities often analyze time-series data in which the ordering of data points is not arbitrary but instead is indexed by time. A suite of algorithms developed for MD simulations involves the analysis of time-lagged data using feature transformations such that the dynamics are Markovian—i.e., independent of the system history for a defined lag time $\tau$. Since the dynamics do not depend on system history, the entire dataset can be represented by $T-\tau$ time-lagged pairs, for $T$ total simulation time.

Equivalently, time-lagged data can be represented by two matrices,

[TABLE]

where $\mathbf{z}_{t}$ is the column vector of mean-free feature values at time $t$, where $\mathbf{X}$ represents the first $T-\tau$ data points and $\mathbf{Y}$ represents their values after a time shift of $\tau$. For MD analyses, we assume that the length of $\mathbf{z}_{t}$ is smaller than $T$; i.e., $\mathbf{X}$ and $\mathbf{Y}$ have more rows than columns. For example, $\mathbf{z}_{t}$ might contain the (mean-free) inter-residue distances or side chain torsional angles of a protein.

Time-lagged independent component analysis (TICA)Naritomi and Fuchigami (2011); Schwantes and Pande (2013); Pérez-Hernández et al. (2013) leverages the matrix representations (1) and (2) to solve the generalized eigenvalue problem,

[TABLE]

for symmetric correlation and covariance matrices $\mathbf{C}_{\tau}$ and $\mathbf{C}_{0}$, where,

[TABLE]

respectively.

The TICA solution (3) gives a set of weights stored in the columns of $\mathbf{\Phi}$ (i.e., the eigenvector solutions) which can be used to transform the data from its original space to a new space. Each solution is a reaction coordinateMcGibbon, Husic, and Pande (2017): it characterizes a dynamical process within the data. We will see later that TICA is a special case of canonical correlation analysis (CCA).

It is clear from definitions (4) and (5) that a strict reversibility requirement is imposed upon the TICA formulation. In particular, (4) assumes that viewing the data in the backwards direction in time would be equally valid if it were to have occurred in the forward direction, which introduces bias that may not be appropriate for systems with rare events. This reversibility requirement has both mathematical and physical implications: the eigenvalues of the solution will be real and the dynamical system modeled is assumed to be in thermodynamic equilibrium. The two interpretations are connected: the processes in the simulated system, represented by the weights (eigenvector solutions), relax to equilibrium when the dynamics are reversible.

TICA was introduced into MD analysis in 2011 as a method to identify slow kinetic modes in proteins Naritomi and Fuchigami (2011). Later, TICA was used as a preprocessing step in the construction of MSMs Schwantes and Pande (2013). Most importantly for our purposes, TICA is equivalent to the linear VAC: namely, it has been shown to produce the optimal linear approximations for the slow modes of reversible systems Pérez-Hernández et al. (2013). In other words, by transforming our time-lagged features using TICA solutions, we obtain the best possible linear description of the system kinetics using those features—as long as we are comfortable assuming reversible dynamics.

II.2 Variational approach to Markov processes

However, what if we are not comfortable assuming reversible dynamics, but still want to obtain linearly optimal approximations of the slow modes in our features? In 2017, Wu and Noé (2017) solved this problem: instead of seeking the eigenvector solutions of TICA, which requires a reversibility assumption, the optimal solutions instead come from the singular vectors of the so-called Koopman matrix.

The Koopman matrix $\mathbf{K}$ minimizes the error in the regression problem $\mathbf{Y}=\mathbf{XK}$. For time-lagged data, $\mathbf{K}$ thus approximates a dynamical propagator: the right hand side approximates the propagation of the data after the lag time $\tau$.111Here, we adopt the statistics convention by writing our linear equation $\mathbf{Y}=\mathbf{X}\mathbf{K}$. In the propagator context, it is more common to write $\mathbf{Y}=\mathbf{K}\mathbf{X}$, and thus the equations for $\mathbf{K}$ in this work are equal to $\mathbf{K}^{\top}$ with the latter convention. Approximating the singular functions of the propagator by analyzing $\mathbf{K}$ is the linear VAMP Wu and Noé (2017); Paul et al. (2019).

To compute $\mathbf{K}$ we first calculate three matrices from our mean-free, time-lagged dataset (without any reversibility requirements):

[TABLE]

Then, we perform a whitening transformation on the matrices $\mathbf{X}$ and $\mathbf{Y}$,

[TABLE]

in order to remove their internal covariances. Instead of determining $\mathbf{K}$ directly, we instead calculate $\tilde{\mathbf{K}}$, which propagates $\tilde{\mathbf{X}}\rightarrow\tilde{\mathbf{Y}}$ in the whitened space.

The solution to the regression problem $\tilde{\mathbf{Y}}=\tilde{\mathbf{X}}\tilde{\mathbf{K}}$ is given by the ordinary least squares estimatorWehmeyer and Noé (2018),

[TABLE]

which, by using substitutions from (6) through (10), we can show yields222We assume in  (9)-(12) that the matrices $\mathbf{X}^{\top}\mathbf{X}$, etc., are full rank.,

[TABLE]

VAMP proceeds by performing the singular value decomposition (SVD) $\tilde{\mathbf{K}}=\mathbf{U^{\prime}SV^{\prime}}^{\top}$, and then transforming the resulting matrices of left ($\mathbf{U^{\prime}}$) and right ($\mathbf{V^{\prime}}$) singular values into the non-whitened space. This procedure is summarized in Alg. 1.

The outputs $\mathbf{U}$ and $\mathbf{V}$ store the reaction coordinates in their columns—the weights for the transformations of $\tilde{\mathbf{X}}$ and $\tilde{\mathbf{Y}}$, respectively. By the definition of the SVD, we have the following orthogonality conditions:

[TABLE]

When $\mathbf{C}_{00}=\mathbf{C}_{\tau\tau}$ and $\mathbf{C}_{0\tau}=\mathbf{C}_{0\tau}^{\top}$, the SVD of $\tilde{\mathbf{K}}$ is equal to is eigendecomposition. In this case, the data is reversible, and VAMP is identical to TICA Naritomi and Fuchigami (2011); Schwantes and Pande (2013); Pérez-Hernández et al. (2013) and the linear VAC Noé and Nüske (2013).

II.3 Iterative CCA with NIPALS

CCA is a well-known statistical algorithm Hotelling (1936); Wold (2006) that finds transformations of $\mathbf{X}$ and $\mathbf{Y}$ such that their representations in the new latent space (i.e., defined by the transformations) are maximally correlated. Thus, VAMP is a version of CCA where the data is time-lagged Wu and Noé (2017); Noé (2018). For VAMP, we saw in Alg. 1 that a single SVD is performed to obtain all the weights simultaneously.

CCA can be performed with a single SVD as in VAMP, or alternatively by using iterative deflation. This produces different weights due to the deflation process. The iterative deflation algorithm developed by Wold Wold (2006) proceeds as follows:

The nonlinear iterative least squares (NIPALS) method in lines 1–9 of Alg. 2 is equivalent to the power method for computing dominant singular vectors Wold (1973); Wegelin et al. (2000). In lines 11–12, the data is deflated such that new dominant singular vectors can be calculated in the next iteration. In the molecular kinetics context, we can call Alg. 2 “dVAMP”, i.e., VAMP with deflation. The first reaction coordinates (weights) calculated from dVAMP are identical to the first VAMP coordinates (see Appendices A and B). The crucial difference between VAMP and dVAMP is that the nondominant reaction coordinates (i.e., all but the first coordinate) are different due to deflation.

Deflation is common to statistical algorithms such as partial least squares (PLS) and principal component analysis (PCA)Wold (1975); Wegelin et al. (2000), both of which can be performed using Alg. 2 with changes only to lines 4 and 7.444PCA is equivalent to PLS when $\mathbf{X}=\mathbf{Y}$. Deflation can be viewed as subtracting the rank-one approximations of $\mathbf{X}$ and $\mathbf{Y}$ at step $c$ from $\mathbf{X}$ and $\mathbf{Y}$ at the previous update Wegelin et al. (2000). This is intended to minimize the effect of a previously computed weight on the subsequent weights by removing the presence of the associated singular value and vectorWhite (1958). Using deflation to obtain the weights destroys the orthogonality conditions in (13) and (14). Furthermore, there is no persistent matrix of which the singular vectors are computed (see also Appendix C).

By applying the deflated weights to the non-deflated data,

[TABLE]

dVAMP produces transformations that differ from VAMP.555 Note that if the deflated weights are applied to the the deflated data, the transformations will be the same in VAMP and dVAMP. The nonequivalence of VAMP and dVAMP transformations occurs when applying the deflated weights to nondeflated data.

We expect the weights identified using Alg. 2 to be minimally contaminated by previously identified coordinates.

To demonstrate the application of dVAMP in elucidating nondominant reaction coordinates, we consider a low-dimensional toy system. Then, for an atomistic protein folding system, we explore the use of deflation by performing a dVAMP analysis in order to deflate a specific slow dynamical mode from the feature data before applying a kinetic model. This strategy allows selected modes to be obscured from the model optimization procedure, because the associated timescale no longer appears slow.

III Experiments

III.1 Asymmetric double wells

Due to deflation, we expect that reaction coordinates obtained with dVAMP will minimally contaminate subsequently calculated coordinates. To investigate this, we design a dataset with dependencies among its degrees of freedom in order to mimic a real system, but with few enough dimensions that we can assess the results relative to the original data.

First, we perform three independent simulations along the following asymmetric double well potential,

[TABLE]

which is visualized with a red line in the left plot of Fig. 1A. The first simulation produces the histogram of particle positions $P_{1}(n)$ for $n$ up to $50000$ simulation steps. The histogram is illustrated on the same plot, and represents the first dimension of our toy dataset.

The other two independent simulations yield particle positions $P^{\prime}_{2}(n)$ and $P^{\prime}_{3}(n)$. We then generate the second and third dimensions of our toy dataset from $P^{\prime}_{2}$ and $P^{\prime}_{3}$ so that they are dependent upon the first dimension $P_{1}$:

[TABLE]

where $\mathcal{N}$ is the normal distribution. The histograms of particle positions $P_{2}(n)$ and $P_{3}(n)$ are shown in the center and right plots of Fig. 1A.

Now we wish to compare the results of the standard VAMP algorithm with the results of dVAMP. To do this, we input the mean-free three-dimensional coordinates into Algs. 1 and 2 and compute the VAMP and dVAMP reaction coordinates. Figures 1B and C illustrate the toy system results. The first reaction coordinate is identical to the first deflated reaction coordinate because no deflation has been performed. The singular values of the coordinates, provided in Fig. 1B, are identical.

The difference between VAMP and dVAMP becomes apparent in the nondominant reaction coordinates. The second VAMP coordinate is roughly Gaussian; in contrast, the second dVAMP coordinate reveals meaningful structures within the original data. Similarly, the third dVAMP transformation presents a more familiar structure with respect to the original data, whereas the third VAMP transformation is still “contaminated” by a Gaussian near the center of the coordinate. The Spearman rank correlation coefficients Spearman (1904) between pairs of corresponding reaction coordinates are also provided in Fig. 1C, and the coefficients’ distances from $1$ are consistent with qualitative assessments of their differences in structure.

III.2 Folding of villin

III.2.1 Modeling with TICA

Now we are interested in exploring how deflation can be used to tackle setbacks to kinetic modeling in a higher-dimensional protein example. To do this, we analyze an MD simulation of the $35$-residue villin headpiece (PDB ID: 2f4k) performed by Lindorff-Larsen et al. (2011). The folded state is depicted in the purple structures in Fig. 2A. The simulation dataset consists of a single $125$ $\mu$s trajectory with $34$ instances each of folding and unfolding. For our analysis, we subsample the original dataset by a factor of $10$ to produce $2$ ns timesteps.

Before constructing a kinetic model, we convert the Cartesian coordinates output from the simulation into features; in this case, a set of pairwise contact distances, where the distance is defined by the closest heavy atoms in the pair. For our analysis, we consider the pairwise distances among a subset of twelve residues at three-residue intervals starting with the first (LEU 1, GLU 4, …, LEU 34) for a total of $66$ distances.

Now, we want to use TICA to build a kinetic model for our data. To perform TICA, we can use VAMP (Alg. 1) with the reversibility conditions in Sec. II.1. This is achieved by letting $\mathbf{X}\leftarrow\mathbf{X^{\prime}}^{\frown}\mathbf{Y^{\prime}}$ and $\mathbf{Y}\leftarrow\mathbf{Y^{\prime}}^{\frown}\mathbf{X^{\prime}}$, where $\mathbf{X}^{\prime}$ and $\mathbf{Y}^{\prime}$ are the original time-lagged data and ⌢ signifies concatenation along the row axis. We use a lag time of $\tau=50$ ns, or $25$ time steps. The TICA results are shown in a two-dimensional histogram in Fig. 2B. The first TICA coordinate separates the majority of configuration space from a rare misfolded state (green in Fig. 2A). The second TICA coordinate separates the denatured state from the folded state (blue and purple in Fig. 2A, respectively).

In order to demonstrate the application of dVAMP and subsequent analysis with deflation, we will now assert that we are not interested in the rare misfolded state due to its poor sampling, and that it is therefore an artifact that is disruptive to our kinetic model. If we wanted to stop at TICA, we could ignore the first TIC. However, we may be interested in any of the following further analyses:

1. (i)

creating an MSM Husic and Pande (2018), especially when using kinetic mapping Noé and Clementi (2015) (in which the slowest modes are most important to the MSM state decomposition); 2. (ii)

evaluating the approximation quality of that MSM using the VAC Noé and Nüske (2013); McGibbon and Pande (2015); 3. (iii)

macrostating our TICA model or MSM using an eigenvector-based approach such as Perron cluster cluster analysis (PCCA) Deuflhard et al. (2000); 4. (iv)

constructing a VAMPnet Mardt et al. (2018) to obtain a kinetic model (similar to performing (i)-(iii) in sequence);

or the application of a variety of other tools.

What (i)-(iv) have in common is that if an undesired slow process is present in the model, then it will correspondingly affect analysis in the following ways:

1. (i)

the MSM state decomposition will be influenced by the artifact, 2. (ii)

MSM variational scores will be a function of how accurately the artifact kinetics are modeled, 3. (iii)

divisive clustering methods such as PCCA will divide macrostates according to the artifact, and 4. (iv)

the VAMPnet optimization will involve optimally modeling the artifact.

III.2.2 Removing undesired processes with deflation

The state of the art in Markovian kinetic modeling has not yet provided a solution to these problems other than the acquiescence that if we are approximating system kinetics then we must work with whatever the global kinetics are. However, the use of deflation provides a new option. In Sec. III.1, we saw how the use of dVAMP can reveal structure in nondominant reaction coordinates.

Now, we will selectively use deflation to prevent our model from optimizing the kinetics of the precise process we do not want. Rewriting line 11 in Alg. 2 for this purpose, we have,

[TABLE]

where $c$ is the index of the coordinate that is undesired, and $\mathbf{\xi}_{c}$ is the score of that coordinate, which was calculated in the course of Alg. 2.666This is different from $\mathbf{\xi}^{\prime}_{c}$ in Eqn. 15, which is not calculated in the course of Alg. 2 on data that has been deflated $c$ times, but instead is calculated on the original data afterward. This removes the effect of the $c$th coordinate on the dynamics by setting the eigen- or singular value of the associated eigen- or singular vector of the matrix product $\tilde{\mathbf{K}}$ (12) to zero.

For our purposes, this means that deflation has rendered the corresponding dynamical process fast: it appears as noise in the new basis. Because the VAC, VAMP, and related methods seek to optimally approximate the slowest processes in the system, deflation effectively “removes” this process from the optimization procedure of the global kinetic model by lowering its timescale in the new feature basis.

To achieve this for villin, we perform dVAMP as in Sec. III.1 through the first dVAMP coordinate in order to obtain the scores for the undesired (slowest) process. Then, we apply (17) with $c=1$.777 Note that if the undesired process is not the slowest process, i.e., $c>1$, we can deflate it exclusively using (17) once for $c$ only. However, we must perform dVAMP through the $c$th coordinates to obtain the necessary weights. Recall that our data $\mathbf{X}_{\text{original}}$ are the $66$ pairwise contact distances, and so $\mathbf{X}_{\text{new}}$ are corresponding deflated distances.888To calculate all of the distances, we perform the same transformation for $\mathbf{Y}$ using $\omega_{c}$ in order to obtain the last $\tau$ data points, which are not present in $\mathbf{X}$ due to the time lag. The deflated distances form a basis in which the eigenvalue of an undesired process is set to zero. They can now be input into any model ranking framework (e.g., MSMs or VAMPnets with variable hyperparameters) such that the resulting kinetic models will not be optimizing the artifact kinetics alongside the interesting ones.

III.2.3 Constructing VAMPnets for villin folding

We now investigate the difference in modeling results between the original and deflated bases using VAMPnets Mardt et al. (2018). VAMP (Alg. 1) calculates an SVD in line 1 which produces a diagonal $k\times k$ matrix $\mathbf{S}$ of positive singular values $\sigma_{i},i\leq k$. The sum of the highest $k$ singular values raised to the $r$th power was proven by Wu and Noé (2017) to be a variational score—the VAMP-$r$ score—for the approximation quality of the kinetic model:

[TABLE]

It follows from the existence of a scalar variational bound that machine learning algorithms can be designed to optimize it McGibbon and Pande (2015). Using neural networks, VAMPnets achieve exactly this by constructing few-state Koopman models directly from features by optimizing a VAMP score Mardt et al. (2018). Here, we build VAMPnets for both our original (mean-free) contact distances and our new, deflated basis. We chose VAMPnets because, due to the neural network architecture, there is no opportunity to manually modify a model while it is being built (e.g., by manually removing an undesired TIC).

For our VAMPnets, we use a lag time of $300$ ns and choose two metastable states. The training and validation sets are split to be approximately equal in size ($51$% and $49$% of the data, respectively), and $10$ epochs are performed when training each model. The remaining hyperparameters are the defaults at the time of modeling (batch size, 1000; network depth, 6; layer width, 100; learning rate, 0.0001). The VAMPnet output is a list of state assignment probabilities for each frame in the trajectory; binary (“hard”) memberships are determined by assigning each frame to the state with the greater probability.

Our VAMPnet results are illustrated in Fig. 3. The top row of histograms summarizes the VAMPnet constructed on the original contact distances. The upper left plot shows distributions for each (hard) VAMPnet state across the values of the first TIC from the original TICA analysis. The upper right plot shows distributions for the same states across root-mean-square deviation (RMSD) values to the folded structure of villin. It is clear from these plots that the two-state VAMPnet separates our artifact (4.6% of frames) from the rest of the ensemble.

The second row of histograms shows the same analysis for the VAMPnet constructed on the deflated contact distances with the same hyperparameters. The lower left plot shows that both VAMPnet states contain structures from the artifact region. The lower right plot shows that the new two-state VAMPnet generally isolates the folded structure (26.8%) from the denatured ensemble. Thus, by using deflated contact distances, we obtain a VAMPnet that captures folding and unfolding using the same parameters with which we were unable to capture folding and unfolding before the basis deflation.

We can see from the lower left histogram in Fig. 3 that the “artifact” is present in the folded state, and further analysis shows that the majority of structures assigned to the folded state with RMSDs greater than $7.5$ Å to the folded structure are “artifact” structures (as defined by their location along the first TIC). We do not claim that this analysis—a two-state VAMPnet with deflated distances—is the best one for villin (for example, previous research suggests a three-state model is better for describing the same dataset Beauchamp et al. (2012)). Instead, we have demonstrated that dVAMP can be used to prevent a kinetic model from optimizing the accuracy of a (subjectively) undesired process a priori.

IV Discussion

We have examined incorporation of deflation into state-of-the-art methods in molecular kinetics. It is already known that VAMP is equivalent to CCA without deflation. In this study we explore the use of dVAMP—VAMP with deflation—for the development of reaction coordinates. Then, we demonstrate the selective use of deflation to remove a specific slow process from a kinetic model; in this case, a VAMPnet.

In the toy example in Sec. III.1, we see that nondominant reaction coordinates obtained from dVAMP qualitatively provide more information about the system than their nondeflated VAMP counterparts. In addition to increased interpretability, these coordinates might be used for enhanced sampling. For example, MD simulations can be started from the top of a barrier revealed only in a deflated reaction coordinate.

In the protein folding example, we see that deflation can obscure a long-timescale process from a kinetic model in order to facilitate further analysis that is desired to be independent of that process. In MD simulations, the dominant processes may not be the same as the processes of interest due to low sampling or force field artifacts Banushkina and Krivov (2015); McGibbon, Husic, and Pande (2017). Thus, deflation presents a systematic way of removing the effects of these undesired modes on model building. Using deflated data, MSMs, Koopman models, or VAMPnets can be constructed such that the eigen- or singular vectors of the resulting model capture only those dynamics which are under investigation. Furthermore, multiple models can be compared with the guarantee none of them are optimizing the accuracy of the undesired slow mode.

Frameworks such as MSMs and Koopman models require that the dynamical modes are orthogonal. Thus, if dynamical modes of interest are faster than an undesired process, they are constrained by this orthogonality condition. By deflating the basis, we gain the flexibility to optimally model the dynamical modes of interest without requiring that they are orthogonal to something uninteresting.

This work opens the door not only for improved analysis but also for related methods developments. Often we are interested in some separate observable that changes in time with the data and its relationship to the features. The implementation of a deflation method in the context of explaining the relationships between a feature and an observable is a valuable pursuit and has already been explored in the time-independent context Hub and De Groot (2009); Krivobokova et al. (2012). The adaptation of such methods to time-series dynamics would be of great utility to the field and the advance presented here may serve as a starting point.

Several open source Python software packages facilitated this study. The double well potential in Sec. III.1 was simulated using PyEMMA Wehmeyer et al. (2018). The 2D histogram in Fig. 2 was generated using Corner Foreman-Mackey (2016), contact distances for villin were calculated with MDTraj McGibbon et al. (2015), VAMPnets were calculated using the deeptime package Mardt et al. (2018). An example Jupyter notebook Kluyver et al. (2016) is provided to reproduce the analysis in Sec. III.1 which also uses NumPy Van Der Walt, Colbert, and Varoquaux (2011), SciPy Jones et al. (01), Scikit-learn Pedregosa et al. (2011), Matplotlib Hunter (2007), and Seaborn Waskom et al. (2017). Protein trajectories were visualized with VMD Humphrey, Dalke, and Schulten (1996).

Acknowledgements

BEH is extremely grateful to Fabian Paul, Muneeb Sultan, Moritz Hoffmann, Andreas Mardt, Luca Pasquali, Stefan Klus, Tim Hempel, and Chris Wehmeyer for helpful discussions and feedback. We thank D. E. Shaw Research for providing the protein simulation dataset. BEH and FN acknowledge funding from the European Commission (ERC CoG 772230 “ScaleCell”) and the Deutsche Forschungsgemeinschaft (SFB1114/A04).

Appendix A Intuition for the relationship between VAMP and canonical CCA

We can demonstrate the equivalence of VAMP and canonical CCA (“dVAMP” for time-lagged data) by considering the calculation of the first component, which is identical in both formulations.

By making substitutions from lines 4–8 of Alg. 2, we can write,

[TABLE]

Further substituting the definition of the Moore-Penrose pseudoinverse Penrose (1955), $\mathbf{X}^{+}\equiv(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}$, we obtain,

[TABLE]

which has been shown to be equivalent to the power method for identifying the dominant eigenvector of a matrix Wegelin et al. (2000). Thus we see that $\phi_{c}$ is the dominant eigenvector of the expression,

[TABLE]

for $\mathbf{X}$ at a given iteration (i.e., possibly deflated).

To obtain an intuition for the connection of this expression (21) to the VAMP solution in Alg. 1, we first substitute (9), (10), and (12) into $\tilde{\mathbf{Y}}=\tilde{\mathbf{X}}\tilde{\mathbf{K}}$,

[TABLE]

which simplifies to,

[TABLE]

Thus the matrix that converts $\mathbf{X}$ to $\mathbf{Y}$ (i.e., propagates the time-series if $\mathbf{X}$ and $\mathbf{Y}$ represent time-lagged data) is,

[TABLE]

where the subscript $f$ indicates that, in the time-lagged interpretation, the matrix propagates the data forward in time. It is easy to see that for the backward direction,

[TABLE]

Now, we can see that (21) is equal to $\mathbf{K}_{f}\mathbf{K}_{b}$, as well as the analogous solution for the $\mathbf{Y}$ weights:

[TABLE]

CCA with deflation via NIPALS obtains the dominant eigenvector of the left-hand sides, whereas VAMP calculates the dominant eigenvector of the right-hand sides. The dominant eigenvalues of the expressions in (26) and (27) are the squares of the dominant singular values obtained from VAMP and dVAMP.

In general, the equivalences above apply when $\mathbf{X}$ and $\mathbf{Y}$ are the same for both algorithms, and only the dominant component is considered.

Appendix B CCA with iterative deflation via the SVD

Algorithm 2 can be rewritten with an SVD in place of the inner NIPALS loop as follows:

However, the implementation of this algorithm contains multiple inverse square roots, which will be numerically unstable with deflation, leading to errors that are propagated with each new component. Regularization can be employed to avoid ill-conditioned covariance matrices, but then the solution should no longer be expected to match that of the algorithm with NIPALS. Indeed, it is generally a good idea to regularize covariance matrices, including in the NIPALS algorithm and TICA. The interested reader is referred to Appendix B in Ref. McGibbon, Husic, and Pande, 2017.

Appendix C Connection to Koopman operator approximation

The matrix $\mathbf{K}$ is a finite-dimensional linear approximation to the continuous integral Koopman operator $\mathcal{K}$ and can be called the Koopman matrix Wu et al. (2017); Wu and Noé (2017); Paul et al. (2019); Koopman (1931); Mezić (2005); Klus et al. (2018). In the time-lagged context, the operator propagates a view of the system at time $t$ to a possibly different view of the system at time $t+\tau$. While we have focused on obtaining reaction coordinates, it is useful to note that the weight (singular) vectors obtained in VAMP (Alg. 1) are approximations to the singular functions of the Koopman operator. The first weight vector calculation in dVAMP (Alg. 2) approximates the dominant singular functions because it is equivalent to the single SVD calculation in VAMP (Alg. 1; see Appendices A and  B).

However, if a relationship between the subsequent weight vectors calculated in dVAMP and $\mathcal{K}$ exists, it is not straightforward. In a sense, a new operator is approximated with each deflation for a matrix built from new, deflated datasets $\mathbf{X}_{c}$ and $\mathbf{Y}_{c}$, in which the presence of the previously calculated $c$ weight vectors have been removed. Thus, each new set of weight vectors approximates the dominant singular functions of some new operator $\mathcal{K}_{c}$, which describes the deflated dynamics.