{'fraction_non_alphanumeric': 0.03384698567180319, 'fraction_numerical': 0.010507344327295666, 'mean_word_length': 4.732617005889038, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identification of simple reaction coordinates from complex dynamics

Reaction coordinates are widely used throughout chemical physics to model and understand complex chemical transformations. We introduce a definition of the natural reaction coordinate, suitable for condensed phase and biomolecular systems, as a maximally predictive one-dimensional projection. We then show this criterion is uniquely satisfied by a dominant eigenfunction of an integral operator associated with the ensemble dynamics. We present a new sparse estimator for these eigenfunctions which can search through a large candidate pool of structural order parameters and build simple, interpretable approximations that employ only a small number of these order parameters. Example applications with a small molecule's rotational dynamics and simulations of bovine pancreatic trypsin inhibitor show that this approach can filter through statistical noise to identify simple reaction coordinates from complex dynamics.


I. INTRODUCTION
The reaction coordinate -a single collective variable that quantifies progress in a chemical reaction -is a ubiquitous concept in chemical kinetics. 1,2 Reaction coordinates are, for example, required for computing reaction rates using transition state theory, 1-3 computing kinetically meaningful free energy barriers, 4 and accelerating conformational sampling in many biomolecular simulation protocols. [5][6][7][8][9] Their most important use, however, is often in facilitating insight into chemical reaction mechanisms. [10][11][12] Implicit in the concept is the notion that the measurement of reaction coordinate is dynamically informative, and provides a proxy for the rate-limiting dynamical processes of the system. Reactions in soft matter and condensed phase systems, such as the folding of a protein or an enzyme-catalyzed chemical transformation take place in a high-dimensional phase space that may include many uninvolved solute and solvent coordinates. In this regime, identification of reaction coordinates is difficult. 12 Physical intuition may suffice to determine these critical degrees of freedom for low-dimensional systems, such as simple bimolecular gas-phase reactions. But in more complex processes involving tens of thousands or more atoms, rough energy landscapes, and/or solvent dynamics, methods to identify the reaction coordinate that rely merely on physical intuition or trial and error can be ad hoc and unsystematic. [13][14][15][16] We recognize that the identification of a system's reaction coordinate(s) is invaluable for physical interpretation of complex molecular systems, that researchers now have access to extremely large data sets of unbiased molecular dynamics (MD) simulations of biologically relevant macromolecules, and that the interpretation of these data is often a major bottleneck. 17 We therefore aim to develop a method to use these MD data sets to infer accurate and interpretable reaction coordinates. Our approach builds on time-structure independent components analysis (tICA). 18,19 But these tICA-derived reaction coordinates can be a black box; they are difficult to interpret physically because of their abstract construc-tion as linear combinations of a large number of structural features. In contrast, our new estimator explicitly adds a sparsity consideration into the formulation to filter through statistical noise and identify simple physical reaction coordinates from complex dynamics.
The structure of this paper is as follows: First, we define the natural reaction coordinate(s) based on a set of intuitive mathematical properties that these collective variables should satisfy. After introducing these properties, we discuss their relationship to other commonly used definitions of the reaction coordinate. Next, we prove that this definition is satisfied by the leading eigenfunctions of an integral operator governing the ensemble dynamics. * Finally, we introduce and demonstrate a practical new estimator which can approximate these reaction coordinates as extremely sparse, interpretable, regularized linear combinations of structural order parameters.

II. DEFINING THE NATURAL REACTION COORDINATE
Although (or perhaps because) the idea of the reaction coordinate is so widely used in chemical kinetics, the community has not always agreed on its precise meaning. A number of different definitions thereof have been proposed, including the minimum energy path or intrinsic reaction coordinate (MEP), [22][23][24][25] the minimum action path (MAP), [26][27][28][29][30] and the committor function. 31,32 In order to proceed in the face of this definitional ambiguity, we begin from first principles and propose a set of properties that a natural reaction coordinate should satisfy for any time-homogeneous, reversible, ergodic Markov process. This approach is geared towards conformational dynamics of soft matter systems, and we make none of the assumptions common in chemical kinetics about the existence of two metastable states, about the relative importance of entropic or enthalpic barriers, about low temperature, or about the number of pathways that are possible. This level of generality does come with a trade off; it makes it impossible to leverage quasiequilibrium approximations, and our algorithms will require equilibrium sampling. The mathematical properties which we specify require that the natural reaction coordinate (a) be a dimensionality reduction that (b) is defined only by the system's dynamics, and that (c) is the maximally predictive projection about the future evolution of the system. Below, we describe and define each of these criteria in detail.

A. A dimensionality reduction from Ω to R
The natural reaction coordinate should be a function which maps any point in the system's full phase space to a single real number. Notating the reaction coordinate as q, and phase space as Ω, † we may specify this as q : Ω → R.
The reason for this form is that it should be welldefined to calculate how "far along" the reaction coordinate any conformation is, or to speak about the mean value of the reaction coordinate for some equilibrium or non-equilibrium ensemble of conformations. Reaction coordinates taking this form include geometric or physical observables which could, in principle, be as simple as the distance between two specific atoms. Later on, we will show how the formulation also extends naturally to multiple orthogonal reaction coordinates.
On the other hand, path-based definitions of the reaction coordinate such as the MEP or MAP do not take this form. Instead of functions from Ω to R, a path through phase space is a function from R to Ω. These paths map an arc length to a phase space coordinate, and the value of the reaction coordinate is undefined for all conformations in Ω that are not on this path. For the minimum energy path, this issue was discussed by Natanson et al., 33 who showed that while a reaction coordinate of the form Ω → R could be defined by introducing a projection operator onto the MEP, there was considerable ambiguity in the choice of projection function. This † We use the phrase 'phase space' to refer to either a position, momenta phase space, or a position-only configuration space, depending on the underlying dynamics. For thermostatted Hamiltonian or Langevin dynamics, Ω = R 6N , where N is the number of atoms. For overdamped Langevin dynamics, also called Brownian or Smoluchowski dynamics, Ω = R 3N . In periodic boundary conditions, the position space is some 3N -dimensional torus, but the exact definition of Ω is not critical for our purposes.
ambiguity was present even for reactive systems containing only 3 atoms without roughness, and are exacerbated in high-dimensional and condensed phase systems. This is one factor which makes the Ω → R formulation more attractive than the R → Ω formulation.

B. Uniquely determined by the dynamics
The natural reaction coordinate should be uniquely defined by the equations of motion that govern the underlying dynamics in Ω, which include the system's Hamiltonian, boundary conditions, and integration scheme. We wish to define the natural reaction coordinate in a way that does not depend on particular "reaction" or "product" conformations or subsets of phase space.
Although it may appear intuitive to define the reaction coordinate in terms of two end points or two states, this definition has a number of formal and practical drawbacks. Subdividing phase space into non-overlapping reactant and product states, A ⊂ Ω, B ⊂ Ω, A ∩ B = ∅, is a useful device, but this is a construct imposed by the modeller, not the underlying Hamiltonian. All experimentally measurable observables, such as ensemble averages, single-molecule time series, or time-correlation functions of a spectroscopic quantity are independent of whether the modeller labels certain regions of phase space as A or B.
For systems containing a small number of atoms, it is often relatively obvious how these states should be determined: e.g. for a bond-forming reaction, one can simply measure whether the distance between the atoms is greater than a certain cutoff. And when the states are metastable, many quantities which might formally depend on the exact specification of the states' boundaries in fact have a very weak dependence thereon, as long as the perturbed state boundaries are still metastable. 34 But in high-dimensional systems where entropy plays a dominant role, and when confronted with significant roughness in the energy landscape on energy scales less than k B T , it can be very difficult in practice to identify these metastable states. Furthermore, many systems have more than two metastable states.
Consider protein folding dynamics, where A and B would generally be taken to be the protein's folded and unfolded states. A number of practical definitions of the folded or unfolded state, based on metrics including rootmean-square deviations to a crystal structure, numbers of native contacts, or radii of gyration are defensible. None, however, are obviously mandated. If the definition of the natural reaction coordinate depends on the exact linedrawing between folded and unfolded, each definition of the state boundaries would lead to a slightly different natural reaction coordinate, with no criteria to judge which is optimal.
In our view, a formal defintion of the natural reaction coordinate should be unique and independent of any partitioning of phase space into regions, and only a function of the system's underlying dynamics. As a dimensionality reduction, the natural reaction coordinate should teach us about the system's metastable states, not the other way around.

C. Maximally predictive projection
Finally, the key property that we use to define the natural reaction coordinate relates to its ability to optimally predict the dynamics. Of all possible one-dimensional measurements of the state of some high-dimensional dynamical system, the natural reaction coordinate should be the most informative about the future evolution of the system. This relates to the expectation, common in chemical kinetics, that the dynamics along the reaction coordinate are rate-limiting, and that all other degrees of freedom in the system equilibrate more rapidly. The maximally predictive single coordinate will measure progress with respect to the rate-limiting bottlenecks, as the orthogonal coordinates can more reliably be assumed to be at, or near, equilibrium.
We now formalize this notion mathematically. To begin, we define the following quantities: • The system has a unique equilibrium distribution over phase space, µ(x) : Ω → R. Note that ∀x, µ(x) > 0 and Ω dx µ(x) = 1.
• Initially, the state of an ensemble is described by a (generally non-equilibrium) probability distribution over phase space, p 0 (x) : Ω → R.
• We consider an ansatz reaction coordinate, q(x) : Ω → R, and an associated scalar, τ , which will be interpreted as a timescale of the dynamics along the ansatz reaction coordinate.
• The scalar projection of the initial distribution, p 0 , along the reaction coordinate is measured as q|p 0 = Ω dx q(x)p 0 (x).
• At some later time, t > 0, the system will have evolved from p 0 to a new distribution over phase space, p t (x) : Ω → R, according to the underlying equations of motion for the dynamics. Note that while p t is a probability distribution, it is not a random variable; it is produced deterministically from p 0 and the system's equations of motion. Now, consider the task of constructing an approximation to p t . This approximation,p t , is constrained to depend only on µ(x), q(x), τ , t, and q|p 0 . That is, given knowledge of the equilibrium distribution, the ansatz reaction coordinate, its timescale, and no other information about the current ensemble, p 0 , beyond its projection onto the ansatz reaction coordinate, our goal is to construct a prediction of the future ensemble at some later time t.
A basic dimensional analysis argument and the constraint that Ωp t = 1 is sufficient to establish that, assuming that q is measured in a system of units such that it has mean zero and unit variance in the equilibrium ensemble, the functional form ofp t given q must bẽ for some non-random function, f , that is independent of x. Later on, we will show that f is necessarily an exponential, f (t/τ ) = e −t/τ . For diffusion on a double well potential, a diagrammatic example of the family of predictions,p t , that can be made given a particular choice of q is shown in Fig. 1.
Even with full knowledge of the Hamiltonian and equations of motion, this prediction will not be exact because the one-dimensional measurement, q|p 0 , gives incomplete information about p 0 (x). We define the error in the prediction, E p0 [q], as the µ −1 -weighted mean squared error, Note that this error depends on the arbitrary initial distribution. To remove this dependency, we consider the worst-case error by maximizing over all possible initial distributions, The functional E[q] thus measures how well the measurement of an arbitrary collective variable can be used to predict the future state of the system. We define the natural reaction coordinate, q * , as the minimizer of E[q]. It is, in this sense, the collective variable which is maximally informative about the system's dynamics.

D. Alternative Definitions
The approach we have taken is not the only one possible. Note first the choice of error functional, Eq. (2). While it may not be initially intuitive, the µ −1 -weighting on the norm is the logical choice for a mean squared error. It is the µ −1 measure, combined with detailed balance, that ensures, for example, that minimizer, q * , is strictly independent of t (see Section III B). A different choice, like the Kullback-Leibler divergence of Wasserstein distance would be possible, 36 but lead to substantially different results. Additionally, observe that in contrast to many other formulations, 37-39 our approach is not based on the explicit construction of a one-dimensional Smoluchowski-like diffusion along the reaction coordinate.
Predictions,pt, made by the natural reaction coordinate, q * , for Smoluchowski diffusion on two-well potential, U (x) = (x − 1) 2 (x + 1) 2 with a uniform diffusion constant, D = 1. The upper left panel shows the stationary distribution, µ(x), and the lower left panel shows the natural reaction coordinate, q * , which changes sign between the two metastable states. The main panel shows the family of possible predictions,pt that can be made by Eq. (1) using this choice of q, indicating the variable partitioning of density between the two basins. For an arbitrary initial distribution, p0, this coordinate minimizes the worst-case predictive error about the future ensemble pt given only knowledge of the current ensemble's projection onto q. As discussed in Section III, q * was calculated from the second eigenfunction of the Smoluchowski operator, which was determined in this case using the FiPy PDE solver. 35 Next, we turn our discussion to an alternative reaction coordinate definition, the committor function. This quantity was first introduced by Onsager as the splitting probability for ion-pair recombination. 40 The committor is defined based on the prior identification of two nonoverlapping states, A ⊂ Ω, B ⊂ Ω, A ∩ B = ∅, which do not fully partition phase space, A ∪ B ⊂ Ω. Then, the committor, p A (x), is defined as the probability that a trajectory initiated from x would enter the set A before entering B. 31,32 In the context of protein folding, where A is taken to be the protein's folded state, the committor is often referred to as p-fold. 41,42 The committor, p A (x), takes a value of 1 for conformations inside A, and 0 for conformations inside B. The condition {x : p A (x) = 1/2} defines a transition state ensemble or separatrix -the set of conformations equally likely to commit to either state A or state B.
Using the concept of the ensemble of transition paths, which are defined as trajectory segments following the moment at which the system has exited the set A and up until the systems enters the set B, without re-entering A, Hummer proved an important result. 34 He showed that, for diffusive dynamics, the probability of being on a transition path given that the system is at x, P(TP AB |x), is determined by the committor alone, P(TP AB |x) = 2p A (x)(1 − p A (x)). This implies also that the separatrix can be identified as the set of conformations which are most likely to be on reaction paths.
A number of computational methods build approximations to the transition path ensemble, committor or iso-committor surfaces. These include transition path sampling (TPS), 31,32 transition interface sampling, 43 and the finite temperature string method. 44,45 Most of the existing algorithms that identify physical reaction coordinates from molecular simulations are based on committor analysis or TPS. [46][47][48][49][50][51][52][53] In the simplest version, one initializes a large collection of trajectories from isosurfaces of an ansatz reaction coordinate and measures which of the two basins, A or B, they commit to. If this coordinate is a good approximation to the committor, the measured splitting fraction will be narrowly peaked around the characteristic value. 54 Criteria based on this observation can then be used to screen an ansatz reaction coordinate, or optimize the parameters of a model for the reaction coordinate. 47 More efficient maximum likelihood method which fit a parametric model for the reaction coordinate from TPS data further refine this approach. 49,51 By design, these algorithms rest on the preidentification of the A and B states, and are not naturally suited to systems with more than two metastable states, although multiple-state extensions are available. 55 For the reasons discussed in Section II B, we thus dispute the claim that the committor should be taken as the perfect or exact reaction coordinate. 47,48,56,57 The authors' experiences with large-scale simulations of protein folding and activation on Folding@Home have shown that it can be difficult to locate and precisely define these metastable states. This suggests that, for an important class of problems, the metastable states should be constructed from the output of some model, as opposed to being treated as a modelling input. [58][59][60] These considerations motivate our formulation of the natural reaction coordinate in a manner independent of the choice to label certain regions of phase space as A or B.

III. A DOMINANT EIGENFUNCTION IS THE NATURAL REACTION COORDINATE
In this section, we demonstrate that the natural reaction coordinate, as defined by the minimizer of Eq. (4), is the second leading eigenfunction of an integral operator associated with a system's Markovian dynamics in Ω. For simplicity, we work here with a discrete-time Markov chain, {X 0 , X 1 , X 2 , . . .}, such as a typical all-atom MD simulation with a finite time step integrator, assuming only that the prediction interval t is greater than 1 step (typically on the order of 2 fs). Afterwards, we note why the same results apply if the underlying dynamics are a continuous-time Markov process, and discuss the natural generalization to multiple reaction coordinates.
We note that others have also used the term reaction coordinate for these eigenfunctions. 61,62

A. Preliminaries
The one-step dynamics of a system's Markovian evolution forward in time can be completely described in terms of a stochastic transition density kernel, where B (y) is the open -ball centered at y with infinitesimal measure dy. Essentially, this kernel measures the conditional probability of jumping from x to y in one step. ‡ Integrating over the initial ensemble, p t , gives a Chapman-Kolmogorov equation for the evolution of the ensemble to p t+1 , By assumption, we consider only ergodic and reversible Markov processes. Ergodicicity is the property that there do not exist two or more regions of Ω that are dynamically disconnected. That is, the integrated transition density is strictly positive, y∈A p(x, y) > 0 for all x and ‡ For example, if the underlying dynamics are overdamped Langevin on a potential energy function U (x) in units of kT with unit diffusion constant, simulated using an Euler-Maruyama integrator with a unit time step, the stochastic transition density kernel, p(x, y)dy, would be the probability density function of a Gaussian distribution with meanȳ = x − ∇U (x) and variance all non-empty subsets of Ω, A. The reversibility condition is that the Markov chain obeys a detailed balance equation with respect to its stationary measure, µ(x), For molecular dynamics, µ(x) is the equilibrium distribution associated with the thermodynamic ensemble that the system is sampling, such as the Boltzmann distribution at constant temperature, and reversibility can be interpreted as a type of generalized symmetry on the function p(x, y).
The form of our maximally predictive projection formulation suggests that the reaction coordinate acts like a perturbation to the equilibrium distribution. This suggests that we consider the equations for the time evolution of a new function, u t (x) ≡ p t (x)/µ(x), which measures the same information as p t (x), but encoded the excess or depletion of probability in an ensemble with respect to the stationary distribution. Applying the Chapman-Kolmogorov equation to the time evolution of u t , we have This equation is taken to define the action of the one step backward transfer operator, T (1), which is uniquely defined by the transition density kernel, The transfer operator has many properties -we refer the interested reader to the monograph of Schütte, Huisinga, and Deuflhard for mathematical details. 21 For our purposes, the most relevant properties are that T (1) is compact and self-adjoint, and thus has a complete, countable set of real eigenfunctions and eigenvalues, which we number in decreasing order by eigenvalue magnitude. Each ψ i can be assumed to be normalized such that they are orthonormal with respect to the µ-weighted inner product, Furthermore, the largest eigenvalue is λ 1 = 1, with associate eigenfunction ψ 1 (x) = 1, and the absolute values of the remaining eigenvalues lie within the unit interval, |λ i | < 1. 21 These properties imply that the action of T (1) on u t can be written as a spectral decomposition, By repeatedly applying the single-step T (1) operator, we can also build the multi-step T (t) operator. Because of the linearity of the operator and orthonormality of the eigenfunctions, each repeated application only pulls out another factor of the eigenvalue in the sum. The spectral decomposition of T (t) is thus B. Analysis of the error functional We now apply this spectral decomposition of the transfer operator to the analysis of the error functional from Section II C and prove why the natural reaction coordinate is equal to the second transfer operator eigenfunction, q * = ψ 2 .
First, observe that the form of the prediction about the future state of the system made using the reaction coordinate, Eq. (1), can also be written as some operator that maps p 0 →p t , or equivalently in u-notation as an approximate transfer operator, The approximate transfer operator,T (t), is rank 2; it has two non-zero eigenvalues, 1 and f (t/τ ), with associated eigenfunctions 1 and q(x) respectively.
Next, we rewrite the error the error functional, Eq. (2) in u-notation as well, where the maximum is understood to be taken over properly normalized u 0 , ||u 0 || µ = 1. Now, we prove by why ψ 2 = min q E[q]. Observe that for any q, there exists a function v(x) in the span of the first three eigenfunctions of T , v = a 1 ψ 1 + a 2 ψ 2 + a 3 ψ 3 , which is normalized, v|v µ = 1, and which is in the null space ofT , over all possible µ 0 , it also must be greater than the error incurred for this particular starting distribution, § To be concrete, set a 1 = 0, and choose a 2 and a 3 to satisfy a 2 ψ 2 |q µ = −a 3 ψ 3 |q µ and a 2 2 + a 2 3 = 1.
where the third line only includes a sum up to i = 3 because, by construction, v is in the span of the first three eigenfunctions. The final line follows because of the ordering of the eigenvalues and the normalization of v, implying 3 i=1 a 2 i = 1. Interpreting this inequality, we see that the worst-case prediction error for any ansatz reaction coordinate, q, is always greater than or equal to λ 2t 3 . Furthermore, for the particular choice q = ψ 2 and f (t/τ ) = λ t 2 , the equality is achieved, E[ψ 2 ] = λ 2t 3 . ¶ If we define τ ≡ −1/ ln λ 2 , we see also that f (t/τ ) can be written as f (t/τ ) = e −t/τ . Therefore, ψ 2 is the natural reaction coordinate, the minimizer of E[q].
The reader may recall that this argument is equivalent to the Eckart-Young Theorem on the optimal low-rank approximation of a matrix. 63

C. Continuous-time Markov processes
When the generating process is a continuous-time Markov process, T (1) has an infinitesimal generator, L, The set of eigenfunctions of L and T (t) are equivalent, so for these processes, ψ 2 can be defined in either manner.

D. Multiple reaction coordinates
One attractive property of this definition of the reaction coordinate is that it generalizes naturally to multiple orthogonal reaction coordinates ordered by timescale.
Recall that the maximally predictive projection criterion from Section II C assumed that the approximation, p t , was to be formed only from knowledge of the equilibrium distribution and the ansatz reaction coordinate. The multiple coordinate generalization follows from modifying this criteria to assume knowledge of µ and the first k − 1 eigenfunctions, µ and ψ 2 , . . . , ψ k−1 . Additionally, assume that the projection of the initial distribution onto each coordinate is available. Then, another application of the Eckart-Young Theorem shows that the maximally predictive remaining ansatz coordinate is ψ k . Multiple ¶ To demonstrate that E[ψ 2 ] = λ 2t 3 , note that for this choice of q and f (t/τ ),T (t) is equal to the sum of the first two terms in the spectral decomposition of T (t). The squared spectral norm of the difference between the two operators is the square of the largest eigenvalue of the difference operator. The first two eigenpairs having been subtracted out, the square of the largest remaining eigenvalue is λ 2t 3 . orthogonal natural reaction coordinates can thus be defined in a stepwise manner, and shown to be equal to the leading eigenfunctions, ψ 2 , . . . , ψ k . Systems contains k metastable states will have k − 1 eigenfunctions whose associated eigenvalues are close to one, separated from the remaining eigenvalues by a socalled spectral gap. 67

E. Two-dimensional example
In the left panel of Fig. 2, we show an example potential with two possible pathways between the dominant basins. The potential is given by the following expression, 38 For Smoluchowski dynamics at kT = 5 with a homogeneous diffusion constant, D = 1, the natural reaction coordinate, ψ 2 (x, y), is shown with solid contour lines in the right panel of Fig. 2. Although ψ 2 can be calculated without explicitly notating any two regions A and B as the reactant or product state, it provides a natural measure of progress of any conformation or ensemble between the two dominant metastable states in the upper left and lower right regions of the potential.

IV. THE TICA APPROXIMATOR
Markov state models (MSMs) and time-structure independent components analysis (tICA) are two widely used approximators for ψ 2 that can be parameterized directly from molecular dynamics trajectories. 18,19,67,68 Other popular estimators include diffusion maps and kernel tICA. 20,61,62,69,70 In the tICA method, the goal is to find the optimal variational approximation to ψ 2 using a linear combination of basis functions. These basis functions are generally structural order parameters that can be evaluated easily for each snapshot in a simulation, such as the distance between certain pairs of atoms or some nonlinear transformation thereof, torsion angles between quartets of atoms, or root-mean-squared deviations to certain landmark conformations.
Assume that there are m linearly-independent basis functions, where typical values of m are in the hundreds to thousands. Without loss of generality, we assume that the basis functions have been mean-subtracted, so that they have zero mean in the equilibrium ensemble. We label the collection of basis functions as {χ j } m j=1 . Because T is self-adjoint, it can be shown that the true eigenfunction, ψ 2 , satisfies a variational theorem, 71,72 Because inner products of the form q|T (t)•q µ can be interpreted as the value of the autocorrelation function of a mean-zero, unit variance observable at time t, 18,71,72 we see as well that ψ 2 , in addition to being the most predictive collective variable, as discussed above, is the most slowly decorrelating collective variable under the system's dynamics.
As in variational quantum chemistry methods, this quantity serves as a figure of merit for the optimization of a trial function. Expanding the ansatz as q = i a i χ i , the maximization is equivalent to the quadratic optimization problem The solution, a * , yielding the best approximation to ψ 2 in the span of the basis set is the generalized eigenvector associated with the largest generalized eigenvalue of the matrices C(t) and Σ. 73 The symmetric matrix, C(t), and positive-definite matrix, Σ, have elements given by, where the expectations are understood to be taken over the stochastic process. As discussed in detail by Schwantes and Pande 18 and Pérez-Hernández et al. 19 , the matrix elements can be estimated by empirical averages over the snapshots in molecular dynamics trajectories. The matrix C(t) is a collection of time-lagged correlations between the basis functions, and Σ is a covariance matrix of the basis functions. In Appendix A, we discuss the use of shrinkage estimators in approximating Σ from timeseries data.

V. A SPARSE APPROXIMATOR FOR THE DOMINANT EIGENFUNCTION
The tICA method has one obvious drawback: the solution, our approximate natural reaction coordinate, is a linear combination of all m basis functions, and the loadings are typically non-zero. This makes the solutions difficult to interpret in a mechanistic manner, because hundreds or thousands of different interatomic distances and/or torsion angles, for example, can be combined together into a single collective variable. Because an important property of reaction coordinates is their role in facilitating physical interpretation of the underlying molecular system, we consider it desirable to reduce the number of explicitly used variables.
These same interpretability issues arise with numerous methods in machine learning and statistics. For example, in multivariate linear regression, a response variable is modeled as the linear combination of input variables. Interpretable models, with only a small number of nonzero coefficients, can be obtained using variable selection methods such as the lasso. 74 The log-norm regularizer used in this work, log(1+|x|/ ) log(1+1/ ) , with = 10 −6 , as compared to the 1 norm. The log-norm is a closer approximation to the 0 norm, and is attractive computationally for this problem because it leads to a more efficient optimization algorithm than the 1.
In this section, we introduce a new sparse approximator for ψ 2 . The solution will share the same form as the tICA approximation, q = i a i χ i , except that the vast majority of the expansion coefficients, a i , will be zero.
One general approach for building sparsity-inducing estimators is to augment the objective function -in our case, Eq. (26) -with a regularization term that penalizes model complexity and steers the optimization towards solutions that fit the data well, but also remain simple. By scaling the strength of this term, the modeller can trade off between the two goals.
Arguably the most natural sparsity-inducing regularizer would be the 0 norm, a penalty proportional to the number of non-zero elements in the solution vector. Unfortunately, 0 -penalized problems generally require an NP-hard combinatorial search. For many problems, such as linear regression, the most common numericallytractable regularizers which lead to sparse solutions are based on the 1 norm, which is sometimes interpreted as a relaxation of 0 . 75,76 However, both the 0 and 1 versions of Eq. (26) are unsuitable. As discussed by Sriperumbudur, Torres, and Lanckriet 77 , the addition of either an 0 or 1 penalty to the Eq. (26) objective leads to the intractable problem of maximizing a non-concave objective function. They considered an alternative relaxation of the 0 penalty, Choosing a fixed > 0 yields a regularizer that is concave (see Fig. 3), which is a property that will allow the sparse tICA method with this choice regularizer to be optimized efficiently as a difference of convex programs. 78 Therefore, to define this sparse tICA algorithm, we adopt the following formulation: * * where ρ ≥ 0 is the regularization strength. At ρ = 0, the problem reduces to standard tICA. Larger values of ρ will induce sparsity in the solution vectors. Investigating sparse generalized eigenvalue problems, Sriperumbudur, Torres, and Lanckriet 77 showed that Algorithm 1 is a globally convergent method for solving Eq. (30). The algorithm is iterative, and refines an initial guess. Each iteration requires solving Eq. (31), a quadratically-constrained quadratic program (QCQP).

Algorithm 1 Sriperumbudur, Torres, and Lanckriet 77
Require: C ∈ S n , Σ ∈ S n ++ , ρ > 0, > 0 Choose τ > max(0, −λmin(C)), These QCQPs are convex. When the number of basis functions, m, is small (less than a few hundred), we have found that they can be solved quickly and with high accuracy by off-the-shelf convex optimization libraries. However, for sparse tICA, our interest is in searching for sparse linear combinations from libraries of many thousands of possible structural order parameters. In this regime, more efficient algorithms are necessary.

VI. AN ADMM SOLVER FOR THE QCQP SUBPROBLEM
We now derive a new, efficient solver for Eq. (31) using the alternating direction method of multipliers (ADMM). ADMM is a general method for constructing optimization algorithms for problems of the form minimize x,z f (x) + g(z) subject to Ax − Bz = c (32) * * At this point, we switch the notation slightly for clarity of presentation. x ∈ R m will be the vector of sparse tICA expansion coefficients being optimized, and we take the t-dependence of C(t) to be implicit, so we simply use the notation C.
where f (x) and g(z) are convex, but not necessarily smooth, functions. See Boyd et al. 79 for a comprehensive review. We take f (x) to be the original objective function from Eq. (31), where D(w) is matrix with the vector w along the diagonal, and g(z) to encode the constraint, where A = B = I n , and c = 0. The ADMM algorithm, in so-called scaled form, consists of the following iterations.
where is a scalar that acts like a step size parameter, and can be adjusted over the course of the optimization to maintain stability. By splitting the objective function into two parts, f and g, the algorithm can alternate taking steps that minimize over the variables x and z separately, with the u variable serving to pull these variables towards each other and enforce the constraint that x = z at convergence.
The advantage of this formulation is that, as we now show, both the x and the z optimization steps can be performed very efficiently.

A. ADMM x update
The x optimization, Eq. (35), can be rewritten as where v = z (k) − u (k) . This function is component-wise separable over the elements of x, f (x) = i f i (x i ). The minimization, Eq. (37), can thus be carried out as n separate scalar minimizations, Although this objective function is not differentiable, it is a simple application of subdifferential calculus to compute a simple closed-form expression for the minimizer. The explicit solution is where S, the soft-thresholding function, is defined as This simple form and component-wise separability means that the ADMM x update can be computed extremely rapidly.

B. ADMM z update
Because g(z) is a hard boundary function, the z update, Eq. (36), can be interpreted as the projection of a point a = x (k+1) + u (k+1) onto the constraint set, {z : z T Σz ≤ 1}, a hyper-ellipsoid. The problem can be rewritten as For the nontrivial case in which the point a lies outside the ellipsoid, a T Σa > 1, the solution, z * , is on the border of the ellipsoid, z * T Σz * = 1. By precomputing the eigendecomposition of Σ, this can be solved efficiently using Kiseliov's method which is detailed in Appendix B. 80 An open source implementation of the estimator is available in the MSMBuilder software package at http: //msmbuilder.org.

C. Further orthogonal reaction coordinates
Like tICA, our algorithm is not restricted to finding a single reaction coordinate, but can also identify sparse approximations to the other long-timescale eigenfunctions, ψ 3 , . . . , ψ k . Unlike in the tICA method, in which the full set of solutions can be computed simultaneously with a single call to a standard generalized eigensolver, each sparse reaction coordinate must be estimated with a separate calculation.
As with most iterative sparse principal components analysis methods, we obtain the remaining generalized eigenvectors by subtracting the influence of the solution from the matrix C, and then restarting optimization using the deflated matrix. The tradeoffs between methods for this deflation step have been discussed by Mackey. 81 Based on the recommendations therein, we have adopted Mackey's Schur complement deflation strategy.

D. Hyperparameter selection and implementation notes
In order to use sparse tICA in practice, a value of the regularization strength, ρ, must be chosen. When ρ = 0, sparse tICA reduces to the standard tICA algorithm, and An overcomplete set of 510 internal coordinates were measured from each frame, which included four dihedral angles (described by carbons 2-1-5-4, 2-1-5-6, 3-1-5-4, and 3-1-5-6) that described the inter-ring torsion angle. larger values of ρ will increase the sparsity. We recommend two possible methods of choosing ρ. First, with cross-validation, the modeller may split the data set into two or more portions, optimize the reaction coordinate at different values of ρ using one fraction of the data set, and check the value of the objective function on the left-out data set. For tICA and Markov state models, this approach was discussed McGibbon and Pande. 82 It is equally applicable to sparse tICA.
Alternatively, when the primary goal is to generate physically interpretable reaction coordinates, the modeller may choose the value of ρ to bring the number of non-zero loadings down to a pre-specified number that is amenable to interpretation. When employing this strategy, we recommend that modellers watch the value of the pseudoeiegenvalue (Rayleigh quotient), λ = x T Cx/x T Σx. It should decrease slightly with increasing ρ, but dramatic drops inλ may indicate overregularization.
The procedure also depends on > 0, which controls the shape of the regularizer. Lower values of lead to a tighter approximation of the 0 norm, but can also lead to numerical instabilities as the derivative of the regularizer near zero goes to infinity, as can be seen in Fig. 3. Empirically, we have found that = 10 −6 provides a suitable balance.
Finally, note that the scalar is required during the optimization as well. This parameter affects only the convergence rate of the solver, as opposed to the final solution, and can be dynamically adjusted over the course of the optimization using standard methods described by Boyd et al. 79 VII. EXAMPLES

A. Torsional reaction coordinate
We demonstrate our approach on molecular dynamics simulations of a simple 2-fluorobiphenyl derivative, shown in Fig. 4. This system is interesting as a toy example because chemical intuition suggests that the rotation of the rings with respect to one another will be hindered. We anticipate the dynamics of the aliphatic tails to be  Fig. 4 with increasing values of the regularization strength, ρ. The unregularized tICA results report a reaction coordinate which is a dense linear combination of all 510 input features. In contrast, with increasing values of the regularization strength, ρ, the sparse tICA algorithm filters out this noise to identify only the sines of the four dihedral angles that collectively characterize the inter-ring torsional reaction coordinate, with only a minor decrease in the psuedoeigenvalue,λ.
faster and uncoupled to the reaction coordinate. Can our algorithm recover this sparse reaction coordinate?
After parameterization with the generalized Amber forcefield, 83 we simulated the system in the gas phase for 250 ns at 290 K using a Langevin integrator with a friction coefficient of 1 ps −1 and timestep of 2 fs using OpenMM 6.3. 84 Snapshots from the simulation were saved every 20 ps. From each simulation snapshot, we recorded the values of an overcomplete set of 510 internal coordinates, which included the distances between all unique pairs of carbon atoms, measured in nanometers, the angles between pairs of bonded atoms, in radians, and the sine and cosine of the dihedral angles between all quartets of bonded atoms. After mean subtraction, these coordinates form our basis functions, χ i , for tICA and our sparse variant. Despite our chemical intuition, from an algorithmic perspective, finding the reaction coordinate for this system is something like finding a needle in a haystack.
In Fig. 5, we show the resulting dominant eigenvector as estimated by tICA and our new approach using increasing values of the regularizer, ρ. The pseudoeigenvalue,λ, is the Rayleigh quotient of the collective variable, related to its timescale byτ = −1/ ln λ i . In standard tICA, this value is maximized exclusively, whereas in sparse tICA, this objective is balanced against a penalty that favors zero coefficients. We see in Fig. 5 that the tICA solution, as expected, returns a collective variable that is a linear combination of all 510 input coordinates, with a nonzero component on each of the coordinates and significant noise.
In constrast, our sparse tICA algorithm suppresses this noise and identifies sparse collective variables that are formed from linear combinations of only a small number of the input degrees of freedom. This sparsity increases with larger values of the regularization strength, ρ, and only leads to a modest decrease in the approximated timescale associated with the coordinate. For ρ = 10 −3 and ρ = 10 −2 , only four input coordinates survive. Inspection of these coordinates shows that they are the sines of the four dihedral angles that cross between the rings (atoms 2-1-5-4, 2-1-5-6, 3-1-5-4, and 3-1-5-6 in Fig. 4). We interpret these results to show that sparse tICA has, without any prior chemical knowledge, filtered through a collection of structural order parameters, many of which are irrelevant in describing the slowest dynamical process of this molecule, and located the subset which can approximate the natural reaction coordinate.

B. Bovine pancreatic trypsin inhibitor (BPTI)
In this section, we apply the sparse tICA method to analyze the native state dynamics of the bovine pancreatic trypsin inhibitor (BPTI), a small 58-residue globular protein that has been extensively investigated by experimental and computational methods. We reanalyzed the one millisecond all-atom molecular dynamics simulation performed by D.E. Shaw Research at 300K with explicit solvent. 85 With its rigid disulfide bonds, the system remains folded over the course of the simulation, but samples a number of near-native states.
For each frame in the trajectory data set, sampled every 25 ns, we computed the value of an extensive set of 2880 structural order parameters from the backbone and side chain dihedral angles. For each of the 57 pro- The ARG 42 φ angle over the course of the 1 ms simulation of native state dynamics of BPTI performed by D.E. Shaw Research. 85 . Our sparse tICA identifies this as the reaction coordinate for a process that involves the opening and hydration of the protein's core.
tein backbone φ and ψ torsion angles, as well as the 46 χ 1 torsion angles, we computed 18 order parameters by evaluating the probability density function of the von Mises distribution at different values of its location parameter, evenly spaced around the unit circle at 20 • increments. A subset of these functions is shown in Fig. 6. These functions act like softened indicator functions that wrap appropriately on (−180 • , 180 • ). We hypothesized that this would be a suitable basis in which to expand the reaction coordinates for BPTI, because it is well suited for expressing a function representing flux between two regions on a Ramachandran plot. Each structural order parameters in our input basis set can thus be interpreted as roughly indicating whether a particular torsion angle is within one of 18 different ∼ 30 • windows.
Using these input features, we fit a sparse tICA model ������ ������ FIG. 8. The near-native (above) and ARG42-flipped (below) conformations of BPTI from the simulation trajectory. The first panel shows the near-native conformation sampled by the majority of the simulation with a ARG42 φ angle between −50 • and −150 • , with the expected four crystallographic waters. Nearly 800 µs into the simulation, the trajectory samples an alternate state in which the protein's core opens and hydrates and the crystallographic waters can exchange with the bulk. In this state, the ARG42 φ angle has flipped, putting its oxygen pointing into the now hydrated core.
with ρ = 0.005 and observed a surprising result. The first solution depends only on the φ dihedral angle of ARG 42. The timeseries of this angle over the course of the simulation is shown in Fig. 7, and we see that this degree of freedom makes a single dramatic flip over the course of the simulation. When we inspected conformations from this flipped state, we observed that the protein's core had opened and hydrated. While this large-scale structural change is obvious from visual inspection of the trajectory, fact that the φ angle of ARG 42 acts as a switch between these two states was unexpected. While many other degrees of freedom also change between these two states, such as the orientation of the upper disulfide linkage (visible in Fig. 8) these degrees of freedom also fluctuate within the near-native state. It is the rare inward flip of ARG 42 which we observe to draw in solvent to hydrate the protein's small core.

VIII. CONCLUSIONS
In this work, we have introduced a defintion of the natural reaction coordinate as a function that satisfies a set of simple mathematical properties: that it (a) is a dimensionality reduction that (b) is defined only by the system's dynamics, and that (c) is the maximally predictive projection about the future evolution of the system. The definition is particularly apt for soft-matter systems in which there may be more than two metastable states, or for systems in which identifying and structurally defining the metastable states is challenging. For any time-homogeneous, reversible, ergodic Markov chain such as thermostatted molecular dynamics, these properties are uniquely satisfied by a dominant eigenfunction of the transfer operator associated with the dynamics, ψ 2 . This eigenfunction is also the most slowly decorrelating collective variable in the system.
We developed a practical new estimator that builds upon the tICA method for estimating these eigenfunctions. Like tICA, this estimator is used to post-process molecular dynamics trajectories. Unlike the variational tICA method which constructs an approximation to these eigenfunctions using a linear combination of structural order parameters in which all of the coefficients are generally non-zero, our estimator finds sparse solutions. It is thus able both to filter through inevitable statistical noise and identify simple, interpretable strutural order parameters that approximate these natural reaction coordinates, without any prior knowledge of the system.
Application of this method to molecular dynamics simulations of a 2-fluorobiphenyl derivative and BPTI show that the approach can identify reaction coordinates for the slow dynamical processes in these data sets that are readily interpretable. In BPTI, we see that opening and hydration of the protein core is controlled by a flip of a single backbone φ angle at ARG 42.
We anticipate that this method will be useful for the analysis of today's large molecular dynamics data sets.
An implementation of this estimator is available in the MSMBuilder software package at http://msmbuilder. org/ under the GNU Lesser General Public License.

ACKNOWLEDGMENTS
The authors thank Brooke Husic and Thomas J. Lane for helpful discussions made during the preparation of this manuscript, Ariana Peck and Carlos X. Hernández for invaluable copy editing, and the National Institutes of Health under Nos NIH R01-GM62868 for funding. We graciously acknowledge D.E. Shaw Research for providing access to the BPTI trajectory data set.

Appendix A: Covariance matrix estimation
In this section, we discuss some issues related to the estimation of the covariance matrix, Σ, from timeseries data such as molecular dynamics simulations. If we consider a single trajectory of length N and collect the results of the evaluation of each of the zero-meaned m basis functions on each of the T snapshots into a matrix, χ ∈ R m×N , the standard estimator for Σ would be the sample covariance matrix, Covariance matrix estimation is a ubiquitous problem common to many fields of science and engineering, and a number of issues with this estimator are known. In particular, results from random matrix theory suggest that the eigenspectrum of the estimated covariance matrix,Ŝ, is over-dispersed with respect to the true value. That is, its large eigenvalues are too large, and its small eigenvalues are too small. For a fixed number of basis functions, m, the sample eigenvalues can be shown to converge to the true eigenvalues as N goes to infinity, 86 but when m is allowed to grow with N , keeping m/N fixed, results such as the Marčkenko-Pastur law suggest that the sample eigenvalues are not effective estimators, and do not converge to the true eigenvalues. 87 In the context of a weight matrix in a generalized eigenvalue problem, misestimation of the small eigenvalues of S is particularly problematic. The generalized eigenvalue problem requires that S be positive-definite -in the extreme case whenŜ is rank-deficient, the maximum value of Eq. (26) is not defined and we get the matrix equivalent of a division by zero.
The most popular class of stabilized covariance matrix estimators are called shrinkage estimators, and take the formΣ = (1 − γ)S + γ(Tr(S)/m)I, for some positive constant γ. The interpretation of this expression is that the shrunk covariance matrix is a convex combination of two estimators, the (low bias, but high variance) sample covariance matrix, and the (high bias, but low variance) estimator that assumes all basis functions have identical variances and zero covariance. An estimator of this form was first popularized by Ledoit and Wolf in the context of Markowitz portfolio selection. [88][89][90] Other shrinkage targets are possible beyond the scaled identity; we refer the reader to the excellent review by Schäfer and Strimmer. 91 The key insight of Ledoit and Wolf is that, under a Frobenius norm objective on the difference between the shrunk covariance matrix and the true covariance matrix, the asymptotically optimal value of the shrinkage constant, γ, can be estimated directly from S, without knowing the true covariance matrix. Thus, no extra tunable parameters need to be added to the algorithm, which is important for usability.
Further improvements to the Ledoit-Wolf (LW) estimator were made by Chen, Wiesel, and Hero III. 92 First, using the Rao-Blackwell theorem, 93 they produced a more accurate Rao-Blackwellized Ledoit-Wolf (RBLW) estimator for the optimal shrinkage constant that dominates the LW estimator. In addition, unlike the LW estimator, the RBLW estimator can be computed even more efficiently and essentially requires no significant computational work beyond the calculation of the sample covariance matrix, S. The expression for the RBLW-optimal shrinkage constant, γ, is where α, β, and U are given by We recommend this RBLW estimator for Σ for use with both tICA and sparse tICA.

Appendix C: Runtime performance
In addition to our ADMM-based solver, we implemented the sparse tICA algorithm using CVXPY and the off-the-shelf SCS solver to solve the QCQP. 95,96 In Fig. 9, we compare the runtime of these two approaches. For this comparison, we randomly generated the matrix Σ from a Wishart distribution with m degrees of freedom and an identity scale matrix, and initialized the ADMM solver from a vector, x, with elements drawn from the standard normal distribution. The error bars indicate standard deviations over 5 replicates. The timings were performed on a Mid 2014 Apple Macbook Pro laptop.
We see generally that our solver is roughly an order of magnitude faster on the QCQP than CVXPY with SCS. Our sparse tICA implementation, however, is also able to efficiently warm-start, because the vectors w and b also converge during the outer iteration of Algorithm 1. Because of this, we find that when we substitute in the off-the-shelf solver to Algorithm 1, the speedup achieved by our ADMM approach is even more substantial. For example, while converging the first sparse tICA solution with m = 500 using our ADMM implementation takes on the order of 0.1 seconds, the same optimization takes approximately 7 minutes using the off-the-shelf solver.  9. Comparison of the runtime of our specialized QCQP solver and a generic solver using CVXPY and SCS. 95, 96 We observe a speedup of approximately one order of magnitude. Efficient warm-starting of the QCQP in Algorithm 1 yields further improvements in runtime. Error bars indicate standard deviations over 5 replicates.