Transition Manifolds of Complex Metastable Systems

We consider complex dynamical systems showing metastable behavior, but no local separation of fast and slow time scales. The article raises the question of whether such systems exhibit a low-dimensional manifold supporting its effective dynamics. For answering this question, we aim at finding nonlinear coordinates, called reaction coordinates, such that the projection of the dynamics onto these coordinates preserves the dominant time scales of the dynamics. We show that, based on a specific reducibility property, the existence of good low-dimensional reaction coordinates preserving the dominant time scales is guaranteed. Based on this theoretical framework, we develop and test a novel numerical approach for computing good reaction coordinates. The proposed algorithmic approach is fully local and thus not prone to the curse of dimension with respect to the state space of the dynamics. Hence, it is a promising method for data-based model reduction of complex dynamical systems such as molecular dynamics.


Introduction
With the advancement of computing power, we are able to simulate and analyze more and more complicated and high-dimensional models of dynamical systems, ranging from astronomical scales for the simulation of galaxies, over planetary and continental scales for climate and weather prediction, down to molecular and subatomistic scales via, e.g., molecular dynamics (MD) simulations aimed at gaining insight into complex biological processes. Particular aspects of such processes, however, can often be described by much simpler means than the full process, thus reducing the full dynamics to some essential behavior or effective dynamics in terms of some essential observables of the system. Extracting these observables and the related effective dynamics from a dynamical system, though, is one of the most challenging problems in computational modeling (Froyland et al. 2014).
One prominent example of dynamical reduction is arguably given by a variety of multiscale systems with explicit fast-slow time scale separation, mostly singularly perturbed systems, where either the fast component is considered in a quasi-stationary regime (i.e., the slow components are fixed and assumed not to change for the observation period), or the effective behavior of the fast components is injected into the slow processes, e.g., by averaging or homogenization (Pavliotis and Stuart 2008). Much of the recent attention has been directed to the case where the deduction of the slow (or fast) effective dynamics is not possible by purely analytic means, due to the lack of an analytic description of the system, or because the complexity of the system renders this task unfeasible (Froyland et al. 2014(Froyland et al. , 2016Coifman et al. 2008;Dsilva et al. 2016;Nadler et al. 2006;Singer et al. 2009;Crosskey and Maggioni 2017;Vanden-Eijnden 2007;Kevrekidis and Samaey 2009). However, all of these approaches still depend on some local form of time scale separation between the "fast" and the "slow" components of the dynamics.
The focus of this work is on specific multiscale systems without local dynamical slow-fast time scale separation, but for which a reduction to an effective dynamical behavior supported on some low-dimensional manifold is still possible. The dynamical property lying at the heart of our approach is that there is a time scale separation in the global kinetic behavior of the process, as opposed to the aforementioned slowfast behavior encoded in the local dynamics. Here, global kinetic behavior means that the multiple scales show up if we consider the Fokker-Planck equation associated with the dynamics, sayu = Lu, where the Fokker-Planck operator L will have several small eigenvalues, while the rest of its spectrum is significantly larger. Such dynamical systems exhibit metastable behavior, and the slow time scales are the time scales of statistical relaxation between the main metastable sets, while there is no time scale gap for the local dynamics within each of the metastable regions (Bovier et al. 2002;Schütte and Sarich 2014).
Global time scale separation induced by metastability has been analyzed for deterministic (Dellnitz and Junge 1999) and stochastic dynamical systems (Schütte et al. 1999;Huisinga et al. 2004) for more than a decade. A typical trajectory of a metastable dynamical system will spend most time within the metastable sets, while rare transitions between these sets happen as sudden "jumps" roughly along low-dimensional transition pathways that connect the metastable sets (Dellago and Bolhuis 2009;Noé et al. 2009;E and Vanden-Eijnden 2010). For an example, see Fig. 1.
The tool to describe the global kinetic behavior of a metastable system is the socalled transfer operator (the evolution operator of the Fokker-Planck equation), which acts on functions on the state space. The time scale separation we rely on here implies a spectral gap for this operator. This fact has been exploited to find low-dimensional representations of the global kinetics in the form of Markov chains whose (discrete) states represent the metastable sets, while the transition probabilities between the states approximate the jump statistics between the sets on long time scales. Under the name "Markov state models" (MSMs), this approach has led to a variety of methods (Bowman et al. 2014;Schütte and Sarich 2014) with broad application, e.g., in molecular dynamics, cf. Schütte et al. (1999), Pande et al. (2010), Schütte et al. (2011) and Chodera and Noé (2014). This reduction comes with a price: Since the relaxation kinetics is described just by jumps between the metastable sets in a (finite) discrete state space, any information about the transition process and its dynamical features is lost. A variety of approaches have been developed for complementing the MSM approach appropriately (Metzner et al. 2009), but a continuous (in time and space) low-dimensional effective description based on MSMs allowing to understand the transition mechanism is infeasible.
In another branch of the literature, again heavily influenced by molecular dynamics applications, model reduction techniques have been developed that assume the existence of a low-dimensional reaction coordinate or order parameter in order to construct an effective dynamics or kinetics: Examples are free energy-based techniques (Torrie and Valleau 1977;Laio and Gervasio 2008), trajectory-based sampling techniques (Faradjian and Elber 2004;Becker et al. 2012;Moroni et al. 2004;Pérez-Hernández et al. 2013), methods based on diffusive processes (Best and Hummer 2010;Zhang et al. 2016;Pavliotis and Stuart 2008), and many more that rely on the assumption that the reaction coordinates are known. The problem of actually constructing good reaction coordinates remains an area of ongoing research (Li and Ma 2014), to which this paper contributes. Typically, reaction coordinates are either postulated using system-specific expert knowledge (Camacho and Thirumalai 1993;Socci et al. 1996), an approximation to the dominant eigenfunctions of the transfer operator is sought (Schütte and Sarich 2014;Chodera and Noé 2014;Pérez-Hernández et al. 2013), or machine learning techniques are proposed (Ma and Dinner 2005). Froyland et al. (2014) show that these eigenfunctions are indeed optimal-in the sense of optimally representing the slow dynamics-but for high-dimensional systems computational reaction coordinate identification still is often infeasible. In the context of transition path theory (Vanden-Eijnden 2006), the committor function is known to be an ideal (Lu and Vanden-Eijnden 2014) reaction coordinate. In Pozun et al. (2012), the authors construct a level set of the committor using support vector machines, but the computation of reaction coordinates is infeasible for high-dimensional systems. The main problem in computing reaction coordinates for high-dimensional metastable systems results from the fact that all of these algorithms try to solve a global problem in the entire state space that cannot be decomposed easily into purely local computations.
In this article, we elaborate on the definition, existence and algorithmic identification of reaction coordinates for metastable systems: We define reaction coordinates as a small set of nonlinear coordinates on which a reduced system (Legoll and Lelièvre 2010;Zhang et al. 2016) can be defined having the same dominant time scales (in terms of transfer operator eigenvalues) as the original system. We then consider a lowdimensional state space on which the reduced dynamics is a Markov process. Thus, our approach utilizes concepts and transfer operator theory developed previously, but in our case the projected transfer operator is still infinite-dimensional, in stark contrast to its reduction to a stochastic matrix in the MSM approach.
The contribution of this paper is twofold: First, we develop a conceptual framework that identifies good reaction coordinates as the ones that parametrize a low-dimensional transition manifold M in the function space L 1 , which is the natural state space of the Fokker-Planck equationu = Lu associated with the dynamics. The property which defines M is that, on moderate time scales t fast < t t slow , the transition density functions of the dynamics concentrate around M. We provide evidence that such an M indeed exists due to metastability and the existence of transition pathways. Crucially, the dimension of M is often lower than the number of dominant eigenfunctions.
Second, we present an algorithm to construct approximate reaction coordinates. Our algorithm is data-driven and fully local, thus circumventing the main problem of previously proposed algorithms: In order to compute the value of the desired reaction coordinate ξ at a location x in the state space X, only the ability to simulate short trajectories initialized at x is needed. In particular, we assume no a priori knowledge of metastable sets, no global equilibration, and we do not need to resolve the slow time scales numerically. The algorithm is built on two pillars: 1. The simulation time scale t can be chosen a lot smaller than the dominant time scales t slow of the system, such that it is feasible to simulate many short trajectories of length t. 2. We utilize embedding techniques inspired by the seminal work of Whitney (1936) and the recent work Dellnitz et al. (2016) that allows one to take almost any mapping into a Euclidean space of more than twice the dimension of the manifold M and to obtain a one-to-one image of it.
These two pillars together with the low-dimensionality of M imply that we can represent the image of the reaction coordinate in a space with moderate (finite) dimension. Then, we can use established manifold learning techniques (Nadler et al. 2006;Coifman et al. 2008;Singer et al. 2009) to obtain a parametrization of the manifold in the embedding space and pull this parametrization back to the original state space, hence obtaining a reaction coordinate. The locality of the algorithm also implies that reaction coordinates are only computed in the region of state space where sampled points are available. This is a common issue with manifold learning algorithms; here, it manifests as the transition manifold being reliably learned only in regions we have good sampling coverage of. However, recently several methods have appeared in the literature that allow a fast exploration of the state space. These methods do not provide equilibrium sampling, but instead try to rapidly cover the essential part of the state space with sampling points. This can be achieved with enhanced sampling methods such as umbrella sampling (Kumar et al. 1992;Torrie and Valleau 1977), metadynamics (Laio and Gervasio 2008;Laio and Parrinello 2002), blue-moon sampling (Ciccotti et al. 2005), adaptive biasing force method (Darve et al. 2008) or temperature-accelerated molecular dynamics (Maragliano and Vanden-Eijnden 2006), as well as trajectory-based techniques such as milestoning (Faradjian and Elber 2004), transition interface sampling (Moroni et al. 2004) or forward flux sampling (Becker et al. 2012). Alternatively, several techniques such as the equation-free approach (Kevrekidis and Samaey 2009), the heterogeneous multiscale method (HMM) (E and Engquist 2003) and methods based on diffusion maps (Chiavazzo et al. 2016) have been developed to utilize short unbiased MD trajectories for extracting information that allows much larger time steps. This can be combined with reaction coordinate-based effective dynamics (Zhang et al. 2016;Zhang and Schuette 2017).
In principle, the method we present in this article may be combined with any enhanced sampling technique in order to generate sampling points that cover a large part of the state space. For simplicity, we will use long MD trajectories to generate our sampling points, but we do not require that the points are distributed according to an equilibrium distribution.
The paper is organized as follows: Sect. 2 introduces transfer operators, which describe the global kinetics of the stochastic process. Based on these transfer operators, we define metastability, i.e., the existence of dominant time scales. In Sect. 3, we describe the model reduction techniques Markov state modeling and coordinate projection that are designed to capture the dominant time scales of metastable systems. Furthermore, we characterize good reaction coordinates. In the first part of Sect. 4, we show that our dynamical assumption ensures the existence of good reaction coordinates, and then in the second part we describe our approach to compute them. Several numerical examples are given in Sect. 5. Concluding remarks and an outlook are provided in Sect. 6. using transfer operators associated with the system and their eigenfunctions. In this section, we will introduce different transfer operators needed for our considerations.

Transfer Operators
In what follows, P[ · | E] denotes probabilities conditioned on the event E and E[· | E] the expectation value. Furthermore, {X t } t≥0 is a stochastic process defined on a state space X ⊂ R n .
Definition 2.1 (Transition density function) Let A be any measurable set, then the transition density function p t : X × X → R ≥0 of a time-homogeneous stochastic process {X t } t≥0 is defined by That is, p t (x, y) is the conditional probability density of X t = y given that X 0 = x.
With the aid of the transition density function, we can now define transfer operators. Note, however, that the transition density is in general not known explicitly and needs to be estimated from simulation data. In what follows, we assume that there is a unique equilibrium density that is invariant under {X t } t≥0 , that is, it satisfies a.e. on X. Let μ denote the associated invariant measure dμ = dx.
Definition 2.2 (Transfer operators) Let p ∈ L 1 (X) be a probability density, 1 u = p/ ∈ L 1 μ (X) be a probability density with respect to the equilibrium density , and f ∈ L ∞ (X) be an observable of the system. For a given lag time t: (a) The Perron-Frobenius operator P t : L 1 (X) → L 1 (X) is defined by the unique linear extension of with respect to the equilibrium density is defined by the unique linear extension of to L 1 μ (X). 1 We denote by L q the space (equivalence class) of q-integrable functions with respect to the Lebesgue measure. L q ν denotes the same space of function, now integrable with respect to the measure ν.
All these are well-defined non-expanding operators on the respective spaces.
The equilibrium density satisfies P t = , that is, is an eigenfunction of P t with associated eigenvalue λ 0 = 1. The definition of T t relies on , and we have T t u = P t (u ). Instead of their natural domains from Definition 2.2, all our transfer operators are considered on the following Hilbert spaces: They are still well-defined nonexpansive operators on these spaces (Baxter and Rosenthal 1995;Schervish and Carlin 1992;Klus et al. 2017).
Furthermore, we will need the notion of reversibility for our considerations. Reversibility means that the process is statistically indistinguishable from its timereversed counterpart.

Definition 2.3 (Reversibility) A system is said to be
is satisfied for all x, y ∈ X.
In what follows, we will assume that the system is reversible. One prominent example for a class of SDEs satisfying uniqueness of the equilibrium density and reversibility is given by Here, V is called the potential, β is the non-dimensionalized inverse temperature, and W t is a standard Wiener process. The process generated by (2) is ergodic and thus admits a unique positive equilibrium density, given by (x) = exp(−βV (x))/Z , under mild growth conditions on the potential V . Note that the subsequent considerations hold for all stochastic processes that satisfy reversibility and ergodicity with respect to a unique positive invariant measure and are not limited to the class of dynamical systems given by (2). See Schütte and Sarich (2014) for a discussion of a variety of stochastic dynamical systems that have been considered in this context. As a result of the detailed balance condition, the Koopman operator K t and the Perron-Frobenius operator with respect to the equilibrium density T t become identical and we obtain P t f, g 1/μ = f, P t g 1/μ and T t f, g μ = f, T t g μ , i.e., all the transfer operators become self-adjoint on the respective Hilbert spaces from above. Here ·, · μ and ·, · 1/μ denote the natural scalar products on the weighted spaces L 2 μ and L 2 1/μ , respectively.

Spectral Decomposition
Due to the self-adjointness, the eigenvalues λ t i of P t and T t are real-valued and the eigenfunctions form an orthogonal basis with respect to ·, · 1/μ and ·, · μ , respectively. In what follows, we assume that the spectrum of T t is purely discrete given by (infinitely many) isolated eigenvalues. This assumption is made for the sake of simplicity. It is actually not required for the rest of our considerations; it would be sufficient to assume that the spectral radius R of the essential spectrum of T t is strictly smaller than 1, and some isolated eigenvalues of modulus larger than R exist. It has been shown that this condition is satisfied for a large class of metastable dynamical systems, see Schütte and Sarich (2014), Sec. 5.3 for details. For example, the process generated by (2) has purely discrete spectrum under mild growth and regularity assumptions on the potential V .
Under this condition, ergodicity implies that the dominant eigenvalue λ 0 is the only eigenvalue with absolute value 1 and we can thus order the eigenvalues so that The eigenfunction of T t corresponding to λ 0 = 1 is the constant function ϕ 0 = 1 X . Let ϕ i be the normalized eigenfunctions of T t , i.e., ϕ i , ϕ j μ = δ i j , then any function f ∈ L 2 μ (X) can be written in terms of the eigenfunctions as f = ∞ i=0 f, ϕ i μ ϕ i . Applying T t thus results in For more details, we refer to Klus et al. (2017) and references therein.

Implied Time Scales
The (implied) time scales on which the associated dominant eigenfunctions decay are given by If T t is a semigroup of operators, then there are κ i ≤ 0 with λ t i = exp(κ i t) such that t i = − κ −1 i holds. Assuming there is a spectral gap, the dominant time scales satisfy t 1 ≥ · · · ≥ t d t d+1 . These are the time scales of the slow dynamical processes, also called rare events, which are of primary interest in applications. The other fast processes are regarded as fluctuations around the relative equilibria (or metastable states) between which the relevant slow processes travel.

Projected Transfer Operators and Reaction Coordinates
The purpose of dimension reduction in molecular dynamics is to find a reduced dynamical model that captures the dominant time scales of the system correctly while keeping the model as simple as possible. In this section, we will introduce two different projections and the corresponding projected transfer operators. The goal is to find suitable projections onto the slow processes.

Galerkin Projections and Markov State Models
One frequently used approach to obtain a reduced model is Markov state modeling. The goal is to find a model that is as simple as possible and yet correctly reproduces the dominant time scales. Given a fixed t > 0, most authors (Noé et al. 2009;Pande et al. 2010 and it has been studied in detail under which condition this can be achieved (Djurdjevac et al. 2012;Sarich et al. 2010).
There are different ways of constructing an MSM, maybe the most intuitive one is also the simplest: Let the entries of T t be the transition rates between metastable sets. A typical molecular system with d dominant time scales will have d + 1 metastable sets C 1 , . . . , C d+1 (also called cores), and its dynamics is characterized by transitions between these sets and fluctuations inside the sets (see Fig. 1 for an illustration). Since the fluctuations are on faster time scales, we neglect them by setting (Schütte et al. 1999) where P μ denotes the probability measure conditioned to the initial condition X 0 being distributed according to μ. Thus, T t core,i j is the probability that the process in equilibrium jumps to the metastable set C j in time t, given that it started in the metastable set C i . Note that (5) can be equivalently rewritten as where 1 C i is the characteristic function of the set C i . Equation (6) readily suggests that T t core is a projection of the transfer operator T t , namely its Galerkin projection onto the space spanned by the characteristic functions 1 C 1 , . . . , 1 C d+1 (Schütte et al. 1999).
Definition 3.1 (Galerkin projection) Given a set of basis functions ψ 1 , . . . , ψ m ∈ L 2 μ (X), let V := span{ψ 1 , . . . , ψ m } and ψ := (ψ 1 , . . . , ψ m ) . The projection to V or, equivalently, to ψ, The residual projection is given by Equivalently, T t = ψ T t . We also denote the extension of T t to the whole L 2 μ (X), given by ψ T t ψ , by T t . Furthermore, we denote the matrix representation of T t with respect to the basis (ψ 0 , . . . , ψ d ) by T t as well. Either it will be clear from the context which of the objects T t is meant or it will not matter; e.g., the dominant spectrum is the same for all of them.
We see that T core is the matrix representation of the Galerkin projection with respect to the basis functions More general MSMs can be built by Galerkin projections of the transfer operator to spaces spanned by othernot necessarily piecewise constant-basis functions (Weber 2006;Schütte et al. 2011;Weber et al. 2017;Klus et al. 2016Klus et al. , 2017Pérez-Hernández et al. 2013;Noé and Nüske 2013). However, in some of these methods, one also often loses the interpretation of the entries of the matrix T t as probabilities.
Ultimately, the best MSM in terms of approximation quality in (4) is given by the Galerkin projection of T t onto the space spanned by its dominant eigenfunctions ϕ 0 , . . . , ϕ d . This space is invariant under T t since T t ϕ i = λ t i ϕ i and the dominant eigenvalues (and hence the time scales) are the same for the MSM and for T t . Due to the curse of dimensionality, however, the computation of the eigenfunctions ϕ i is in general infeasible for high-dimensional problems.
Remark 3.2 There are quantitative results assessing the error in (4) of the MSM in terms of the projection errors ⊥ ψ ϕ i L 2 μ , i = 0, . . . , d, cf. Schütte and Sarich (2014), Section 5.3. One can obtain a weaker, but similar result from our Lemma 3.5 in the next section.

Coordinate Projections and Effective Transfer Operators
While the MSMs from above successfully reproduce the dominant time scales of the original system, they often discard all other information about the system, such as the transition paths between metastable sets. Minimal coordinates that describe these transitions are called reaction coordinates, and reducing the dynamics onto these coordinates yields effective dynamics (Legoll and Lelièvre 2010;Zhang et al. 2016). The goal of the previous section-namely to retain the dominant time scales of the original dynamics in a reduced model-can now be reformulated for this lowerdimensional effective dynamics or, equivalently, for its (effective) transfer operator.
The so-called coarea formula (Federer 1969, Section 3.2), which can be considered as a nonlinear variant of Fubini's theorem, splits integrals over X into consecutive integrals over level sets of ξ and then over the range of ξ . For (7) where z = ξ(x) and σ z is the surface measure on L z . The coordinate projection, defined next, averages a given function along the level sets of a coordinate function ξ .
where μ z is a probability measure on L z with density (z) det(∇ξ ∇ξ) −1/2 with respect to σ z . Here, (z) is just the normalization constant so that μ z becomes a probability measure. The residual projection is given by To get a better feeling for the action of P ξ , note that , Or, in other words, μ z is the marginal of μ conditional to ξ(x) = z. Note, in particular, that P ξ f is itself a function on X, but it is constant on the level sets of ξ , and thus let us set P ξ f (ξ(x)) = P ξ f (x) for x ∈ L ξ(x) . It follows from the coarea formula (7) and (9) that Next, we state some properties of the coordinate projection.

Proposition 3.4
The coordinate projection has the following properties.
(a) P ξ is a linear projection, i.e., P 2 We use the coordinate projection to describe the dynamics-induced propagation of reduced distributions with respect to the variable ξ . To this end, we define the effective transfer operator T t ξ : We immediately obtain from the self-adjointness of T t (see Sect. 2) and Proposi- Thus, the spectrum of the effective transfer operator lies in the interval [−1, 1], too.
Returning to the purpose of these constructions, we call ξ a good reaction coordi- While the previously introduced Markov state model T t obtained by the Galerkin projection was approximating the dominant spectrum of the original transfer operator by a finite-dimensional operator (i.e., a matrix), the effective transfer operator still acts on an infinite-dimensional space. The reduction lies in the fact that T t operates on functions over X ⊆ R n , but the effective transfer operator T t ξ operates essentially on functions over ξ(X) ⊂ R k , although we embed those into X through the level sets of ξ .
As mentioned above, a Galerkin projection of the transfer operator onto its dominant eigenfunctions is a perfect MSM. In the same vein, we ask here how we can characterize a good reaction coordinate. We can make use of the following general result.
Lemma 3.5 Let H be a Hilbert space with scalar product ·, · and associated norm · ; let Q : H → H be some orthogonal projection on a linear subspace of H, with Q ⊥ = Id − Q. Let T : H → H be a self-adjoint non-expansive linear operator, and u with u = 1 its eigenvector, i.e., T u = λu for some λ ∈ R. If Q ⊥ u < ε, where ζ ≤ Q ⊥ u < ε since Q and T are non-expanding. Thus, u := Qu/ Qu satisfies T Q u = λu + ζ / Qu , and the orthogonality of Q gives Qu > √ 1 − ε 2 . Now, any orthogonal projection is self-adjoint, as shown in the proof of Proposition 3.4; hence, the operator QT Q is self-adjoint, too, and thus normal. From the theory of pseudospectra for normal operators (Trefethen and Embree 2005, Theorems 2.1, 2.2, and §4), we know that if With H = L 2 μ , Q = P ξ , and T = T t we immediately obtain the following result.
Corollary 3.6 implies that if the projection error of all dominant eigenfunctions is small, then ξ is a good reaction coordinate in the sense of (12). Very similar results are available for approximation of the eigenvalues of the infinitesimal generator of the Fokker-Planck equation associated with the transfer operator if the dynamical system under consideration is continuous in time (Zhang and Schuette 2017).
Under which conditions is the projection error small? Let us consider the case where there areφ i : . We then say that ϕ i is a function of ξ or that ξ parametrizes ϕ i . If ξ parametrizes ϕ i perfectly, the projection error obviously vanishes. Thus, trivially, by choosing ξ = ϕ = (ϕ 1 , . . . , ϕ d ) , we obtain a perfect reaction coordinate since withφ i (z) := z i we have ϕ i =φ i • ξ . However, the eigenfunctions are global objects, i.e., their computation is prohibitive in high dimensions. Since we are aiming at computing a reaction coordinate, we have to answer the question of whether there is a reaction coordinate ξ that can be evaluated based on local computations only, while it parametrizes the dominant eigenfunctions of T t well enough such that it leads to a small projection error. We will see next that this question can be answered by utilizing a common property of most metastable systems: The transitions between the metastable sets happen along the so-called reaction pathways, which imply the existence of transition manifolds in the space of transition densities. A suitable parametrization of this manifold results in a parametrization of the dominant eigenfunctions with a small error.

Identifying Good Reaction Coordinates
The goal is now to find a reaction coordinate ξ that is as low-dimensional as possible and results in a good projected transfer operator in the sense of (12). As we saw in the previous section, the condition P ⊥ ξ ϕ i L 2 μ ≈ 0 is sufficient. Thus, the idea to numerically seek ξ that parametrizes the dominant eigenfunctions of T t in the · L 2 μ -norm seems natural since this would lead to small projection error P ⊥ ξ ϕ i L 2 μ . In fact, eigenfunctions of transfer operators have been used before to compute reduced dynamics and reaction coordinates: In Froyland et al. (2014), methods to decompose multiscale systems into fast and slow processes and to project the dynamics onto these subprocesses based on eigenfunctions of the Koopman operator K t are proposed. In McGibbon et al. (2017), the dominant eigenfunctions of the transfer operator T t , which due to the assumed reversibility of the system is identical to K t , are shown to be good reaction coordinates. Also, committor functions (introduced in Appendix B), which are closely related to the dominant eigenfunctions, have been used as reaction coordinates in Du et al. (1998) and Lu and Vanden-Eijnden (2014).
However, we propose a fundamentally different path in defining and finding reaction coordinates, as working with dominant eigenfunctions has two major disadvantages: 1. The eigenproblem is global. Thus, if we wish to learn the value of an eigenfunction ϕ i at only one location x ∈ X, we need an approximation of the transfer operator T t that has to be accurate on all of X. The computational effort to construct such an approximation grows exponentially with dim(X); this is the curse of dimensionality. There have been attempts to mitigate this (Weber 2006;Junge and Koltai 2009;Weber 2012), but we aim to circumvent this problem entirely. Given two points x, y ∈ X, we will decide whether ξ(x) is close to ξ(y) or not by using only local computations around x and y (i.e., samples from the transition densities p t (x, ·) and p t (y, ·) for moderate t). 2. The number of dominant eigenfunctions (d + 1) equals the number of metastable states, and this number can be much larger than the dimension of the transition manifold. This fact is illustrated in Example 4.1.
Example 4.1 Let us consider a diffusion process of the form (2) with the circular multi-well potential shown in Fig. 2. Choosing a temperature that is not high enough for the central potential barrier to be overcome easily, transitions between the wells typically happen in the vicinity of a one-dimensional reaction pathway, the unit circle. The number of dominant eigenfunctions, however, corresponds to the number of wells. Nevertheless, projecting the system onto the unit circle would retain the dominant time scales of the system, cf. Sect. 5.

Parametrization of Dominant Eigenfunctions
If the (d + 1) dominant eigenfunctions do not depend fully on the phase space X, a lower-dimensional and ultimately easier to find reaction coordinate suffices for keeping the eigenvalue approximation error (12) small. It is easy to see that if there exists a function ξ : X → R k for some k so that the eigenfunctions ϕ are constant on the level sets of ξ , i.e., there exist functionsφ i : R k → R, i = 1, . . . , d, such that ϕ i =φ i • ξ , then the projection error P ⊥ ξ ϕ i L 2 μ is zero. A quantitative generalization of this is the statement that if the eigenfunctions ϕ i are almost constant on level sets of a ξ , then the projection error is small.

Lemma 4.2
Assume that there exists a function ξ : X → R k for some k and functions Proof Assuming (13) holds, there exists a function c i : Thus, we have For the projection error, we then obtain Remark 4.3 From the proof we see that the pointwise condition (13) can be replaced by the much weaker condition for all level sets L z of ξ .
From here on, we address the following two central questions:  Fig. 3 a, b The transition densities Q(x 1 ) and Q(x 2 ) are "similar" to Q(x * ) for some x * on the transition path (dashed line) that connects the metastable sets A and B. c The mapping Q can be thought of as mapping all points that are "similar" under Q to the same point in L 1 (X). The image of Q thus forms a r -dimensional manifold in L 1 (X) so-called reaction pathways, which are the low-dimensional dynamical backbone in the high-dimensional state space, connecting the metastable states via saddle points of the potential V (Freidlin and Wentzell 1998). From now on, we observe the system at an intermediate time scale t slow t t fast (where t slow and t fast are the implied time scales t d , t d+1 from Sect. 2.3) and thus assume that the process X t has already left the transition region (if it started there), equilibrated to a quasi-stationary distribution inside some metastable wells, but has not had enough time to equilibrate globally. At this time scale, starting in some x ∈ X, the transition density p t (x, ·) is observed to approximately depend only on progress along these reaction paths; see Fig. 3 for an illustration. This means that the density p t (x, ·) on the fiber perpendicular to the transition pathway is approximately the same as p t (x * , ·) for some x * on the transition pathway. As this pathway is low-dimensional, this means that the image Q(X) of the map is almost a low-dimensional manifold in L 1 (X).
The existence of this low-dimensional structure in the space of probability densities is exactly the assumption we need to ensure that the dominant eigenfunctions are lowdimensionally parametrizable, and thus that a low-dimensional reaction coordinate ξ exists. This assumption is made precise in Definition 4.4. To summarize, we will see that ξ is a good reaction coordinate if p t (x, ·) ≈ p t (y, ·) for ξ(x) = ξ(y).
Definition 4.4 ((ε, r )-reducibility and transition manifold) We call the process X t (ε, r )-reducible, if there exists a smooth closed r -dimensional manifold M ⊂ L 2 1/μ ⊂ L 1 (X) such that for t fast t t slow and all x ∈ X holds. We call M the transition manifold and the map Q : X → M, the mapping onto the transition manifold. We can set M = cl(Q(X)), where cl(Y) denotes the closure of the set Y. 3 Remark 4.5 While it is natural to motivate (ε, r )-reducibility by the existence of reaction pathways in phase space, it is not strictly necessary. There exist stochastic systems without low-dimensional reaction pathways whose densities still quickly converge to a transition manifold in L 1 . Future work includes the identification of necessary and sufficient conditions for the existence of transition manifolds (see the first point in Conclusion). We also further elaborate on the connection between reaction pathways and transition manifolds in Appendix B.
Remark 4.6 We recall from Sect. 2 that the Perron-Frobenius operator P t is also naturally defined on the space L 2 1/μ (Schervish and Carlin 1992). Further, with the Dirac distribution centered in x ∈ X, denoted by δ x , we formally have p t (x, ·) = P t δ x . Hence, the choice of norm in Definition 4.4 is natural. It should also be noted that since μ is a probability measure, the Hölder inequality yields f L 1 μ ≤ f L 2 μ . Using this we have which shows that if p t (x, ·) and p t (y, ·) are close in the L 2 1/μ norm, they are also close in the L 1 norm. We require the closeness of the respective p t (x, ·) in the L 2 1/μ norm for our theoretical considerations below, but otherwise we will think of them as functions in L 1 .
Note that we only need to evolve the system at hand for a moderate time t t slow , which has to be merely sufficiently large to damp out the fast fluctuations in the metastable states. This will be an important point later, allowing for numerical tractability.
Next, we show that (ε, r )-reducibility implies that dominant eigenfunctions are almost constant on the level sets of Q.
Lemma 4.7 If X t is (ε, r )-reducible, then for an eigenfunction ϕ i of T t with ϕ i L 2 μ = 1 and points x, y ∈ X with Q(x) = Q(y) we have Proof First note that for the transition densities p t (x, ·), p t (y, ·) it holds that (16) With this we can show the assertion: Applying (16), for some function e ∈ L 2 1/μ (X) with e L 2 1/μ ≤ 2ε, we get where in the last equation, we again used that due to reversibility K t = T t and that ϕ i is an eigenfunction. Thus, for the difference, we have Assuming that the eigenfunctions are normalized (which we do from now on), i.e., ϕ i L 2 μ = 1, and that ε is sufficiently small, Lemma 4.7 implies that the dominant eigenfunctions (i.e., |λ i | ≈ 1) are almost constant on the level sets of Q. This can now be used to show that the ϕ i are not fully dependent on X, but only on the level sets of Q (up to a small error), in a sense similar to Lemma 4.2.

Corollary 4.8 Let X t be (ε, r )-reducible. Then there exists a functionφ
Proof Fix x ∈ X, and let z = Q(x). Define the functionφ i bỹ

Since by Lemma 4.7 it holds that |ϕ
thus, our choice ofφ i gives

Embedding the Transition Manifold
In light of Corollary 4.8, one could say that Q is an "M-valued reaction coordinate." However, as we have no access to M so far, and a R k -valued reaction coordinate is more intuitive, we aim to obtain a more useful representation of the transition manifold through embedding it into a finite, possibly low-dimensional Euclidean space.
We will see that we are very free in the choice of the embedding mapping, even though the manifold M is not known explicitly. (We only assumed that it exists.) To achieve this, we will use an infinite-dimensional variant of the weak Whitney embedding theorem (Sauer et al. 1991;Whitney 1936), which, roughly speaking, states that "almost every bounded linear map from L 1 (X) to R 2r +1 will be one-to-one on M and its image." We first specify what we mean by "almost every" in the context of bounded linear maps, following the notions of Sauer et al. (1991).

Definition 4.9 (Prevalence) A Borel subset S of a normed linear space V is called
prevalent if there is a finite-dimensional subspace E of V such that for each v ∈ V, v + e belongs to S for (Lebesgue) almost every e in E.
As the infinite-dimensional embedding theorem from Hunt and Kaloshin (1999) is applicable not only to smooth manifolds, but to arbitrary subsets A ⊂ V of fractal dimension, it uses the concepts of box-covering dimension dim B (A) and thickness exponent τ (A) from fractal geometry. Intuitively, dim B (A) describes the exponent of the growth rate in the number of boxes of decreasing side length that are needed to cover A, and τ (A) describes how well A can be approximated using only finitedimensional linear subspaces of V. As these concepts coincide with the traditional measure of dimensionality in our setting, we will not go into detail here and point to Hunt and Kaloshin (1999) for a precise definition.
The general infinite-dimensional embedding theorem reads: Theorem 4.10 (Hunt and Kaloshin 1999, Theorem 3.9) Let V be a Banach space and A ⊂ V be a compact set with box-counting dimension d and thickness exponent τ . Let k > 2d be an integer, and let α be a real number with

Then for almost every (in the sense of prevalence) bounded linear function
where · 2 denotes the Euclidean 2-norm.
Note that (17) implies Hölder continuity of E −1 on E(A) and in particular that E is one-to-one on A and its image. Using that the box-counting dimension of a smooth r -dimensional manifold K is simply r and that the thickness exponent is bounded from above by the box-counting dimension, thus 0 ≤ τ (K) ≤ r , see Hunt and Kaloshin (1999), we get the following infinite-dimensional embedding theorem for smooth manifolds.
Corollary 4.11 Let V be a Banach space, let K ⊂ V be a smooth manifold of dimension r , and let k > 2r . Then almost every (in the sense of prevalence) bounded linear function E : V → R k is one-to-one on K and its image in R k .
Thus, since the transition manifold M is assumed to be a smooth r -dimensional manifold in L 1 (X), an arbitrarily chosen bounded linear map E : L 1 (X) → R 2r +1 can be assumed to be one-to-one on M and its image. In particular, E(M) is again an r -dimensional manifold (although not necessarily smooth). With this insight, we can now construct a reaction coordinate in Euclidean space: Corollary 4.12 Let X t be (ε, r )-reducible, and let E : L 1 (X) → R 2r +1 be one-to-one on M and its image. Define ξ : R n → R 2r +1 by Then there exists a functionφ i : R 2r +1 → R so that Proof As E is one-to-one on M and its image, it is invertible on E(M). Withφ i chosen as in the proof of Corollary 4.8, defineφ i : Then SinceM := E(M) is an r -dimensional manifold, ξ is effectively an r -dimensional reaction coordinate. Thus, if the right-hand side of (19) is small, the ϕ i are "almost parametrizable" by the r -dimensional reaction coordinate ξ . Using Lemma 4.2, we immediately see that this results in a small projection error P ⊥ ξ ϕ i , and due to Corollary 3.6 in a good transfer operator approximation; hence, ξ is a good reaction coordinate. The reaction coordinate ξ remains an "ideal" case, because we have no access to the map Q and hence to M, but only to Q(x) = p t (x, ·) ≈ Q(x). We show the construction of the ideal reaction coordinate ξ in Fig. 4.

Remark 4.13
The recent work of Dellnitz et al. (2016) uses similar embedding techniques to identify finite-dimensional objects in the state space of infinite-dimensional dynamical systems. They utilize the infinite-dimensional delay-embedding theorem of Robinson (2005), a generalization of the well-known Takens embedding theorem (Takens 1981), to compute finite-dimensional attractors of delay differential equations by established subdivision techniques (Dellnitz and Hohmann 1997).

Numerical Approximation of the Reaction Coordinate
Approximate Embedding of the Transition Manifold We now elaborate how to construct a good reaction coordinate ξ numerically. To use the central definition (18) in practice, two points have to be addressed: 1. How to choose the embedding E? 2. How to deal with the fact that we do not know Q?
For the choice of E, we restrict ourselves to linear maps of the form with arbitrarily chosen linearly independent functions η i ∈ L ∞ (X), where f, η i = f η i . In practice, we will choose the η i : X → R as linear functions themselves, i.e., η i (x) = a i x for some, usually randomly drawn, a i ∈ R n . Note that then η i / ∈ Fig. 5 How to make a good r -dimensional reaction coordinate out of E • Q? Given the analysis from the previous section, we would like to make the level sets L z of Q also the level sets of ξ (red line segment). Unfortunately, we have no access to these (Color figure online) L ∞ , but this is not a problem because we will embed the functions f = p t (x, ·), and p t (x, y) can be shown to decay exponentially as y 2 → ∞, cf. Bittracher et al. (2015), Theorem C.1. Thus, f, η i will exist. For linearly independent η i , these maps are still generic in the sense of the Whitney embedding theorem and thus still embed the transition manifold M. A natural choice for the approximation of the unknown map Q is the mapping to the transition probability density, The values on the right-hand side can in turn be approximated by a Monte Carlo quadrature, using only short-time trajectories of the original dynamics: where the (m) t (x) are independent realizations of X t with starting point X 0 = x, in practice realized by a stochastic integrator (e.g., Euler-Maruyama).
The Computationally Infeasible Reaction Coordinate ξ Note that E •Q is not yet an rdimensional reaction coordinate, since Q(X) is only approximately an r -dimensional manifold; more precisely, it lies in the ε-neighborhood of an r -dimensional submanifold M of L 1 . Hence, also E(Q(X)) is only approximately an r -dimensional manifold; see the magenta regions in Fig. 5.
The question now is how we can reduce E • Q to an r -dimensional good reaction coordinate. Since we know from above that ξ = E • Q is a good reaction coordinate, let us see what would be needed to construct it.
The property of ξ that we want is that it is constant along level sets L z of Q, i.e., ξ | L z = const (because this implies that it is a good reaction coordinate, cf. Corollary 4.12). Hence, if we could identifyM as an r -dimensional manifold in R 2r +1 , we would project E(Q(x)) along E(Q(L z )) ontoM-assuming thatM and E(Q(L z )) intersect in R 2r +1 -to obtain ξ(x) as the resulting point (see Fig. 5, where we would project along the red line on the right). Unfortunately, we have no access to Q (not to mention thatM and E(Q(L z )) need not intersect in R 2r +1 ) and hence to its level sets L z . Thus, this strategy seems infeasible.  (25). Our strategy will be to choose a specific r -dimensional reaction coordinate ξ and to show that in general it can be expected to be a good reaction coordinate.
Let us recall that, by assumption, the set Q(X) is contained in the ε-neighborhood of an unknown smooth r -dimensional manifold M ⊂ L 1 (X). Thus, a generic smooth map E : L 1 (X) → R 2r +1 will embed M into R 2r +1 , forming a diffeomorphism from M toM. Thus, E is going to map Q(X) to an O(ε)-neighborhood ofM. This means the r -dimensional manifold structure ofM should still be detectable and can be identified with standard manifold learning tools. We use the diffusion maps algorithm (see Sect. 4.4), which gives us a map : R 2r +1 → R r (the diffusion map). Then we define ξ as This is depicted on the right-hand side of Fig. 6, where the dashed red line shows the level setLˆz = {z ∈ R 2r +1 : (z) = (ẑ)}. Next, we consider the setLˆz := E −1 (Lˆz) ∩ Q(X). It holds thatLˆz = Q(x) ξ(x) = (ẑ) . Recall that E : M →M is one-to-one; thus,Lˆz intersects M in exactly one point. We define this one point as Q(x), and thus, Q is the projection onto M alongLˆz. We see that Q is well defined and that Q(x) = Q(y) ⇔ ξ(x) = ξ(y).
At this point we assume that E −1 is sufficiently well behaved in a neighborhood ofM and it does not "distort transversality" of intersections such that the diameter ofLˆz is O(ε) with a moderate constant in O(·). We will investigate a formal justification of this fact in a future work, here we assume it holds true, and we will see in the numerical experiments that the assumption is justified. This assumption implies that Q(x) − Q(y) L 2 1/μ = O(ε) for Q(x) = Q(y), i.e., for ξ(x) = ξ(y). Now, however, Lemma 4.7 implies that ϕ i is almost constant (up to an error O(ε)) on level sets of ξ , which, in turn, by Lemma 4.2 and Corollary 3.6 shows that ξ is a good reaction coordinate.

Identification ofM Through Manifold Learning
In this section, we describe how to identifyM numerically. The task is as follows: Given that we have computed E(Q(x i )) =ẑ i ∈ R 2r +1 for a number of sample points {x i } i=1 ⊂ X, we would like to identify the r -dimensional manifoldM, noting the points E(Q(x i )) are in a O(ε)-neighborhood ofM (see Sect. 4.3). Additionally, we would like an r -dimensional coordinate function : R 2r +1 → R r that parametrizesM (so that the level sets of are transversal toM). This is a default setting for which manifold learning algorithms can be applied. Many standard methods exist; we name multidimensional scaling (Kruskal 1964b, a), isomap (Tenenbaum et al. 2000), and diffusion maps ) as a few of the most prominent examples. Because of its favorable properties, we choose the diffusion maps algorithm here and summarize it briefly for our setting in what follows. For details, the reader is referred to Coifman and Lafon (2006), Nadler et al. (2006), Coifman et al. (2008 and Singer et al. (2009).
Given sample points {ẑ i } i=1 ⊂ R 2r +1 , diffusion maps proceed by constructing a similarity matrix W ∈ R × with where · 2 is the Euclidean norm in R 2r +1 , σ > 0 is a scale factor, and h : R → R + is a kernel function which is most commonly chosen as h(x) = exp(−x)1 x≤R with a suitably chosen cutoff R that sparsifies W and ensures that only local distances enter the construction. With D being the diagonal matrix containing the row sums of W , the similarity matrix is then normalized to giveW = D −1 W D −1 . Finally, the stochastic matrix P =D −1W is constructed, whereD is the diagonal matrix containing the row sums ofW . P is similar to the symmetric matrixD −1/2WD−1/2 ; thus, it has an orthonormal basis of eigenvectors {ψ i } −1 i=0 with real eigenvalues γ i . Since P is also stochastic, |γ i | ≤ 1. The diffusion map is then given by Using properties of the Laplacian eigenproblem onM, one can show that indeed parametrizes the r -dimensional manifoldM for suitably chosen σ .
Remark 4.14 The diffusion maps algorithm will only reliably identifyM based on the neighborhood relations between the embedded sample points z i , if the points cover all parts ofM sufficiently well. In particular, as p t (x, ·) and thus E • Q (x) vary strongly with x traversing the transition regions, a good coverage of those regions is required. For the various low-dimensional academic examples Sect. 5, this is ensured by choosing the x i to be a dense grid of points in X. For the high-dimensional example in Sect. 5.2, the evaluation points are generated as a subsample from a long equilibrated trajectory, essentially sampling μ. Both of these ad hoc methods are likely to be unapplicable in realistic high-dimensional systems with very long equilibration times. However, as we mentioned in Introduction, there exist multiple statistical and dynamical approaches to this common problem of quickly sampling the relevant parts of phase space, including the transition regions. Each of these sampling methods can be easily integrated into our proposed algorithm as a preprocessing step.
Fundamentally though, the central idea of our method does not depend crucially on the applicability of diffusion maps. Rather, the latter can be considered an optional post-processing step. Using the 2r + 1-dimensional reaction coordinate i.e., (26) without the manifold learning step, may in practice already represent a sufficient dimensionality reduction.
In addition, situations may occur where the a priori generation of evaluation points is not possible or desired. One of the final goals and the work currently in progress is the construction of an accelerated integration scheme that generates significant evaluation points and their reaction coordinate value "on the fly." This is related to the effective dynamics mentioned in fifth point in Conclusion. However, this also requires us to be able to evaluate the reaction coordinate at isolated points, independent of each other, and thus also necessitates the use of the above ξ instead of ξ .
2. Choose linearly independent functions η j ∈ L ∞ (X), j = 1, . . . , 2r + 1. The essential boundedness of the η j is not necessary, but |η j (x)| should not grow faster than a polynomial as x 2 → ∞. 3. In each point x i , start M simulations of length t and estimate E j Q(x i ) using (23) and (24), to obtain the pointẑ i ∈ R 2r +1 . We discuss the appropriate choice of M and t in Sect. 5.1. 4. Apply the diffusion maps technique from Sect. 4.4 for the point cloud {ẑ i } i=1 , and obtain : R 2r +1 → R r , a parametrization of the point in its r essential directions of variation. 5. By (27), we define the reaction coordinate as The numerical effort of this algorithm depends strongly on the third step. Given evaluation points, and a choice of M trajectories per point, the cost is mainly given by M · · c(t), where c(t) is the effort of a single numerical realization of the dynamics up to time t. The high-dimensional phase space only enters the algorithm as the domain of the observables η j . The cost of evaluating those typically very simple functions 4 at the M · end points of the trajectory is negligible. The cost of the method is thus essentially independent of n.
In order to demonstrate the efficacy of our method, we compute the reaction coordinates for three representative problems, namely a simple curved double-well potential, a multi-well potential defined on a circle, both in low and high dimensions, and two slightly different quadruple-well potentials stressing the difference between a oneand a two-dimensional reaction coordinate.

Curved Double-Well Potential
As a first verification, we consider a system with an analytically known reaction coordinate that is then used for comparison. Consider the two-dimensional drift-diffusion process (2) with potential and inverse temperature β = 0.5. This potential already served as a motivational example for the nature of reaction coordinates in Introduction and is shown in Fig. 1.
The system possesses two metastable sets around the minima (−1, 0) and (1, 0) , which are connected by the transition path {x ∈ R 2 | x 2 = 1 − x 2 1 }. The implied time scales, defined in (3), can be computed from the eigenvalues using a standard Ulam-type Galerkin discretization (Klus et al. 2016(Klus et al. , 2017  and are shown in Fig. 7a. 5 We observe a significant gap between t 1 and t 2 and thus identify t 1 as the last slow and t 2 as the first fast time scales. Choosing the lag time t = 2 then satisfies t slow > t > t fast . A visual inspection of a typical trajectory of length t starting in one of the two metastable sets as shown in Fig. 7b confirms that the respective set is sampled, yet a transition to the other set is a rare event. The low dimension of the system allows us to compute the reaction coordinate on a full regular grid over the phase space. We choose a 40×30 grid in the rectangular region [−2, 2] × [−1, 2] and denote the set of grid points by X. For this system, we expect a one-dimensional transition path and thus a one-dimensional reaction coordinate ξ . That is, r = 1 and 2r + 1 = 3. Thus, we choose three linear observables in our embedding function (21), e.g., whose coefficients were drawn uniformly from [−1, 1]. The expectation value in (23) is approximated by a Monte Carlo quadrature using M = 10 5 sample trajectories for each grid point, cf. (24). The parameter M was chosen such that the error in (24), commonly defined as the variance of the Monte Carlo sum, is sufficiently low. The resulting embedding of the grid points x into R 3 is shown in Fig. 8. The transition path seems to be already parametrized well by the individual components of E • Q.
For this example, the image of X under E • Q should form a compact neighborhood of the one-dimensional manifold E(M), as described in Sect. 4.3. The one-dimensional structure in E Q(X) is clearly visible, see Fig. 9a. To identify the one-dimensional coordinate along this set the diffusion map algorithm is used. Let 1 : E •Q)(X) → R  Legoll and Lelièvre (2010) show that the effective dynamics based on the reaction coordinate ξ * (x) = x 1 exp(−2x 2 ) accurately reproduces the long-time dynamics of the full process-although they do not use dominant eigenvalues of the transfer operator in their argumentation. It is easy to verify that the level sets of ξ * traverse the transition path orthogonally. Figure 10 shows the comparison of the level sets of ξ and ξ * . While the two reaction coordinates have different absolute values, their contour lines coincide well. As the projection operator P ξ only depends on the level sets of ξ , the projected transfer operators T t ξ and T t ξ * should be similar as well.
Projected Eigenvalue Error To conclude this example, we compute the dominant spectrum of the projected transfer operator and compare it to the spectrum of the full transfer operator. To discretize T t ξ , we use a simple Ulam-type discretization scheme based on a long equilibrated trajectory of the full dynamics. Recall from Fig. 10 Selected contour lines (black) of the newly identified reaction coordinate ξ and the reference reaction coordinate ξ * Sect. 3.2 that, although T t ξ formally acts as an operator on functions over X, it is constant along level sets of ξ and thus can be treated as an operator on functions over R r . For completeness, we state the rough outline of an algorithm that we used to approximate T t ξ . An introduction to Ulam-and other Galerkin-type discretization schemes for transfer operators can be found, e.g., in Klus et al. (2016).
1. Compute points X := { (kτ ) x 0 | k = 1, . . . , N }, a discrete trajectory with step size τ of the full phase space dynamics that adequately samples the invariant density . 2. Compute the reaction coordinate ξ on the points X. 3. Divide the neighborhood of ξ(X) into boxes or other suitable discretization elements {A 1 , . . . , A N } and sample the boxes from the trajectory, i.e., compute 4. Count the time-t-transitions within X between the boxes (where t is a multiple of τ ), i.e., compute the matrix 5. After row-normalization, the eigenvalues of T t ξ approximate the point spectrum of T t ξ .
Remark 5.1 Note that the equilibrated trajectory X is typically unavailable for more complex systems. In practice, one would replace steps 1 and 2 by directly computing a reduced trajectory Z = {z 1 , . . . , z N } ⊂ R r whose statistics approximate that of ξ X . The formulation of a reduced numerical integration scheme to realize this is a work currently in progress (see the fifth point in Conclusion).
For our example system, we compute X as a N = 10 6 step trajectory with step size τ = 10 −2 using the Euler-Maruyama scheme. However, to reduce the numerical We observe in Fig. 11 that the eigenvalues of T t ξ and T t are in excellent agreement. Not only the dominant eigenvalues λ 0 , λ 1 are approximated well (as predicted by Lemma 3.5), but also the further subdominant eigenvalues that are not covered by our theory. In particular, the reaction coordinate ξ provides a better approximation to the spectrum of T t than other, manually chosen reaction coordinates: Fig. 11 also shows the eigenvalues of the projected transfer operator associated with the reaction coordinates We see that these are consistently outperformed by the computed reaction coordinate ξ (although it appears that ζ 1 already is quite a good reaction coordinate).

Circular Potential
Let us now compute the reaction coordinates for the multi-well diffusion process described in Example 4.1. The corresponding k-well potential is defined as We use k = 7, for which the potential is shown in Fig. 2a. The potential as well as the dominant eigenvalues of the corresponding transfer operator clearly indicates the existence of seven metastable sets, yet a typical long-time trajectory, shown in Fig. 12a, suggests a one-dimensional transition path, the unit circle B 1 . We demonstrate that with our method, a reaction coordinate of minimal dimension can be computed. We again choose the inverse temperature β = 0.5 and perform the same analysis as in the previous subsection. For this system, a time scale gap between t 6 ≈ 1.53 and t 7 ≈ 0.05 can be found. We thus choose the intermediate time scale t = 0.1. Since we again expect a one-dimensional transition path, the three observables (29) are used for the embedding of M. We use the grid points of a 40 × 40 grid, denoted again by X, over the region [−2, 2] × [−2, 2] as our test points.
The individual components of the embedding E•Q are shown in Fig. 13. The embedded grid points, shown as the individual points in Fig. 14a, seem to concentrate around a one-dimensional circular manifold and thus reveal the one-dimensional nature of the reaction coordinate. Although slightly unintuitive, the diffusion maps algorithm now identifies two significant diffusion map components, as shown in Fig. 14a. The reason is that the circular manifold cannot be embedded into R 1 , so that a two-component coordinate is necessary to parametrize it. Figure 12b shows some contour lines (of equidistant values) of the two components of ξ . We see that ξ is almost constant on the seven metastable sets, but resolves the transition regions well.
Parametrization of the Dominant Eigenfunctions Next, we experimentally investigate how well the dominant eigenfunctions ϕ i of T t can be parametrized by the numerically computed reaction coordinate ξ . If the eigenfunctions are almost functions of ξ , then by Lemma 4.2 and Corollary 3.6 the reaction coordinate is suitable to reproduce all To this end, we compute the dominant eigenfunctions ϕ j , j = 0, . . . , d by the Ulam-type Galerkin method (as in the previous example) and plot ϕ j (x i ) against ξ(x i ). Note that due to the reasons discussed above, the range of ξ is a one-dimensional manifold in R 2 . If ϕ j can be parametrized byξ , we expect that The result is shown in Fig. 15. We clearly see the functional dependency of the first seven (i.e., the dominant) eigenfunctions on the reaction coordinate.

Circular Potential in Higher Dimensions
The identification of reaction coordinates is not limited to two dimensions. To show that our method can effectively find the reaction coordinates in high-dimensional systems, we extend the 7-well potential to ten dimensions by adding a quadratic term in x 3 , . . . , x 10 : V (x) = cos (7 arctan(x 2 , x 1 )) + 10 x 2 1 + x 2 2 − 1 2 + 10 10 j=3 x 2 j . Fig. 16 a, b Two components ξ 1 and ξ 2 on the sampling points X. The picture shows the first three dimensions of x, but is qualitatively the same when replacing x 3 by x j , j = 4, . . . , 10. c The values of ξ 1 and ξ 2 on X plotted against each other We expect the one-dimensional circle {x ∈ R 10 | x 2 1 + x 2 2 = 1, x j = 0, j = 3, . . . , 10} to be the transition path and accordingly choose a three-dimensional linear observable η(x) = A · x, A ∈ R 3×10 , where the coefficients A i j were again drawn uniformly from [−1, 1].
In ten dimensions, the computation of the reaction coordinate on all points of a regular grid is no longer possible due to the curse of dimensionality, and neither is the visualization of this grid. Instead, we compute ξ on 10 5 points sampled from the invariant measure and plot only the first three coordinates. Let this point cloud be called X.
Performing the standard procedure, i.e., embedding X into R 3 and identifying the one-dimensional core using diffusion maps, a two-component reaction coordinate is identified. Coloring the first three dimensions of X by ξ (Fig. 16a, b), we see that the expected reaction pathway is indeed parametrized. This pathway as well as the seven metastable states can also be recognized in a plot of the components of ξ(X) plotted against each other, indicating that the information about the dominant eigenfunctions, thus the long-time jump process, is indeed retained by ξ .

Two Quadruple-Well Potentials
Our theory is based on the existence of an r -dimensional transition manifold M in L 1 (X) around which the transition probability functions concentrate. In Appendix B, we argue that the existence of an r -dimensional transition path suffices to ensure the existence of M. Here we illustrate how the existence of the transition path is reflected in the embedding procedure.
For this we consider the "hilly" and "flat" quadruple-well potentials Fig. 17 Two quad-well potentials V 1 and V 2 possess qualitatively different transition regions Fig. 18 Embedding of the grid points for the a "hilly" and b "flat" four-well potential. A one-dimensional structure is only visible in a, i.e., in the presence of a one-dimensional transition path Both systems possess metastable sets around the four minima (±1, ±1), but V 1 confines its dynamics outside of the metastable sets onto a one-dimensional transition path, whereas V 2 does not impose such restrictions on the dynamics (see Fig. 17). For both potentials the time t = 1 lies inside the slow-fast time scale gap. Assuming a one-dimensional transition manifold (wrongfully for V 2 ), we use the three linear observables (29). A 40 × 40 grid on [−2, 2] × [−2, 2] is used as evaluation points for ξ . The embedding of these points by E • Q is shown in Fig. 18. We observe a onedimensional structure in the case of the "hilly" potential V 1 , whereas the embedding points of the "flat" potential V 2 lie on a seemingly two-dimensional manifold. As these embeddings are approximately one-to-one with the respective transition manifolds M, we conclude that in the case of V 1 the manifold M must be one-dimensional, whereas for V 2 it is two-dimensional.

Conclusion
Our main contributions in this paper are as follows: (a) We developed a mathematical framework to characterize good reaction coordinates for stochastic dynamical systems showing metastable behavior, but no local separation of fast and slow time scales. (b) We showed the existence of good low-dimensional reaction coordinates under certain dynamical assumptions on the system. (c) We proposed an algorithmic approach to numerically identify good reaction coordinates and the associated low-dimensional transition manifold based on local evaluation of short trajectories of the system only. Our numerical examples show how the procedure works, that it can be used in higher dimensions, and the examples give further evidence that the dynamical assumptions from (b) are valid in many realistic cases. The application of our approach to relevant biomolecular problems, e.g., in protein folding, is ongoing work.
Apart from the application to actual molecular systems, there are several open questions and challenges, which we will address in the future: • A rigorous mathematical justification for the dynamical assumption in Definition 4.4 in terms of the potential V and the noise intensity β −1 in (2) would be desirable. This seems to be a demanding task, as the interplay between potential landscape and the thermal forcing is non-trivial. For β −1 → 0 the problem can be handled by large deviation approaches; however, understanding increasing β −1 is challenging: The strength of noise increases, and additional transitions between metastable sets become more probable, as the barriers in the potential landscape become less significant, and thus the reaction coordinate may increase in dimension. • Also related to the previous point, the choice of the correct lag time t is crucial.
Choosing the time too small, the concentration of the transition densities near a low-dimensional manifold in L 1 may not have happened yet, but a too large lag time has severe consequences for the numerical expenses. If no expert knowledge of a proper lag time t is available, it has to be identified in a preprocessing step, for example using Markov state model techniques (Bowman et al. 2014). • As discussed in the last part of Sect. 4.3 and in Fig. 6, we need the embedding E not to distort transversality close to the transition manifold M too much, such that the realized reaction coordinate ξ is indeed a good one. Theoretical bounds shall be developed. This problem seems to be coupled with the problem of how to control the condition number of the embedding and its numerical realization.
• The dimension r of the reaction coordinate may not be known in advance; hence, we need an algorithmic strategy to identify this on the fly. Fortunately, once the sampling has been made, the evaluation of the embedding mapping E, and finding intrinsic coordinates on the set of data points embedded in R k , has a negligible numerical effort; hence, different embedding dimensions k can be probed via (21). Theorem 4.10 suggests that if the identified dimension of the reaction coordinate is smaller than k/2, then a reaction coordinate of sufficient dimension is found. • To benefit from the dimensionality reduction of the reaction coordinate ξ , the dynamics that generates the reduced transfer operator T t ξ has to be described in closed form. We are planning to employ techniques based on the Kramers-Moyal extension (Zhang et al. 2016) to again receive an SDE for a stochastic process on R r . • The embedding mapping E is evaluated by Monte Carlo quadrature (24). Although Monte Carlo quadrature is known to have a convergence rate independent of the underlying dimension n of X, there is still an impact of the dimension on the practical accuracy. This we shall investigate as well.

A Properties of P ξ
Proof of Proposition 3.4 (a) This property has been shown by Zhang et al. (2016) as well, and we include the short reasoning for completeness. The linearity of P ξ is obvious. The property P 2 ξ = P ξ follows from (8) by noting that P ξ f is constant on L z and that μ z is a probability measure for every z.
(b) From (10) we have for f, g ∈ L 2 μ (X) that where ( * ) follows from the linearity of P ξ and the fact that P ξ f | L ξ(x) = const; thus, P ξ g P ξ f (z) = P ξ f (z) P ξ g(z). Expression (30) is symmetric in f and g; hence, it follows that P ξ f, g μ = f, P ξ g μ .
(c) We first prove that P ξ is an orthogonal projection: and the claim follows.

B On the Existence of Reaction Coordinates
To motivate the existence of low-dimensional reaction coordinates, let us assume that the dynamics of consideration has d + 1 metastable regions C 0 , . . . , C d ⊂ X. Let C = i C i . For a selected lag time t > 0 we make the following two assumptions: 1) Fast local equilibration: If x is in (or close to) C i then we have where qs i is the quasi-stationary density of the metastable core C i : lim s→∞ P X s = y τ C i > s = qs i (y)dy with τ C i being the (random) exit time from the set C i . 2) Slow transitions: The typical transition time to reach C\C i when starting in C i is larger than t. In other words, t is such that if the process X s transitions from x to some C i , it did not transition through some other C j with high probability. These two assumptions essentially say that t is much larger than the fast time scales of the system, but smaller than the dominant time scales. It follows that, for any x ∈ X, where by assumption 2) the coefficients q i (x) are given by the committor functions q i (x) = P X t reaches C i before C\C i X 0 = x .
We say that P t δ x is an r -dimensional structure in L 1 (X) if there is a function ξ : X → R r that jointly parametrizes all the committor functions, i.e., q i =q i •ξ withq i : R r → R. If this is the case, then Q(x) = P t δ x ≈ d i=0q i (ξ(x)) qs i =: Q(x) Fig. 19 A potential energy landscape with four minima (black contours) and level sets of q 1 (colored contours). Left: The "hilly" landscape structure confines transition pathways to a narrow region close to the red square connecting the four minima. As a result, all committor level sets are orthogonal to this main transition path. Right: "Flat" landscape structure with more complicated level sets of the committors (Color figure online) and clearly dim(Q(X)) ≤ r since dim(ξ(X)) = r . Moreover, r ≤ d since we can explicitly construct ξ : X → R d as ξ = (q 1 , . . . , q d ). This obviously parametrizes q 1 , . . . , q d , and it also parametrizes q 0 since q 0 = 1 − d i=1 q i . However, r may also be smaller than d. As an example, consider the potential with 4 minima as shown in Fig. 19 on the left. At low temperatures, the "hilly" potential energy landscape confines all transitions between the minima C 0 , . . . , C 3 to a narrow region close to the red square connecting the four minima. Figure 19 shows the level sets of q 0 , the level sets of the other committors are given by the rotational symmetry of the problem. All four committors can be jointly parametrized by a single coordinate ξ which describes clockwise movement along the red square and is constant orthogonal to it. Therefore, r = 1. Figure 19 on the right shows the situation with a "flat" energy landscape. Transition paths between the minima are no longer confined to a onedimensional structure, and the committor level sets are more complicated. We can no longer parametrize all four committors with a single coordinate ξ , so r > 1. On the other hand, dim(X) = 2 so r = 2.
This structural difference of the potentials can also be seen when applying our algorithm to construct the reaction coordinate ξ , see Fig. 18 and Sect. 5.3.