# Improving Estimation of the Koopman Operator with Kolmogorov-Smirnov
Indicator Functions

Van A. Ngo<EMAIL_ADDRESS>[  [  Yen Ting Lin  [<EMAIL_ADDRESS>Danny
Perez  [<EMAIL_ADDRESS>
###### Abstract

It has become common to perform kinetic analysis using approximate Koopman
operators that transforms high-dimensional time series of observables into
ranked dynamical modes. Key to a practical success of the approach is the
identification of a set of observables which form a good basis in which to
expand the slow relaxation modes. Good observables are, however, difficult to
identify a priori and sub-optimal choices can lead to significant
underestimations of characteristic timescales. Leveraging the representation
of slow dynamics in terms of Hidden Markov Model (HMM), we propose a simple
and computationally efficient clustering procedure to infer surrogate
observables that form a good basis for slow modes. We apply the approach to an
analytically solvable model system, as well as on three protein systems of
different complexities. We consistently demonstrate that the inferred
indicator functions can significantly improve the estimation of the leading
eigenvalues of the Koopman operators and correctly identify key states and
transition timescales of stochastic systems, even when good observables are
not known a priori.

ORNL]Advanced Computing for Life Sciences and Engineering, Computing and
Computational Sciences, National Center for Computational Sciences, Oak Ridge
National Lab, Oak Ridge, TN 37830 CNLS]Center for Nonlinear Studies (T-CNLS),
Theoretical Division, Los Alamos National Laboratory, NM 87545, USA
CCS]Information Sciences Group (CCS-3), Computer, Computational and
Statistical Sciences Division, Los Alamos National Laboratory, NM 87545, USA
T]Physics and Chemistry of Materials Group (T-1), Theoretical Division, Los
Alamos National Laboratory, NM 87545, USA

##  1 Introduction

Elucidating the kinetics describing rare structural or chemical reactions is
crucial to understand many biophysical and biochemical systems 1, 2, 3, 4, 5.
Even when long fully-resolved trajectories are available, e.g., via extensive
molecular dynamics (MD) simulations, extracting a reliable representation of
kinetics in terms of a handful of physical observables can be elusive due to
the high dimensionality and complexity of most application-relevant systems,
particularly in the absence of an intuitive reaction coordinate. A powerful
approach aimed at tackling such problems exploits the correspondence between
the spectral properties of the so-called Koopman operator and those of the
dynamics generator. It allows to efficiently reduce high dimensional
timeseries into a compact and tractable representation (See Sec. 2). Various
techniques and spectral analysis methods based on this mathematical formalism
include time-independent component analysis (TICA) for improving Markov State
Modeling (MSM)6, 7, 8, 9, the variational approach for Markov processes (VAMP)
incorporated in MSMbuilder10, 11 and PyEMMA12, as well as different variants
of the extended dynamical mode decomposition (EDMD) approach 13. These methods
underscore that analyzing the spectral properties of the Koopman operator is
powerful to understand and characterize the dynamics of complex systems.

In essence, an estimated Koopman operator based on EDMD (or VAMP or TICA) is
an approximation of the dynamics generator of measurable observables in an
infinite-dimensional Hilbert space14, 15, 16. This approximate Koopman
operator returns the expectation value at a time $t+\tau$ from a value of an
observable at time $t$17, 18. In practice, a finite set of timeseries of
observables is used to obtain an estimation of the Koopman operator in a
desirably small subspace relevant to slow kinetics. The eigenvalues and
eigenvectors of such approximate Koopman operator can then be used to describe
slow relaxation modes together with their corresponding characteristic
timescales. The accuracy of the approximation can however strongly depend on
the choice of observables.

In some cases, an intuitive reaction coordinate can be easily determined. For
instance, two torsional angles of the backbone of the alanine dipeptide as the
reaction coordinates are sufficient to describe the slow dynamics. In cases
where intuitive reaction coordinates cannot be determined, many possible
observables including root mean square displacement (RMSD), all possible
torsional angles, native contacts, and backbone distances can be considered
18. Depending on the target system’s complexity, these generic observables may
not form a good basis to extract transition timecales, requiring more
elaborate schemes to define a better observable 19. Identifying a compact yet
sufficiently-complete set of observables that is able to reliably approximate
the true relevant eigenfunctions of the Koopman operator remains an
outstanding challenge, although guidelines are gradually emerging 18, 20.

To understand this challenge better, we can define the properties of an
optimal set of observables. Mathematically, optimal observables should
effectively act as basis functions to approximate the eigenfunctions
corresponding to a “slow” subspace that describes rare transitions (e.g.,
protein folding/unfolding) or slowest relaxation modes of the exact Koopman
operator 21, which are of course a priori unknown. While general
eigenfunctions can be extremely complex, those that represent slow transition
between $M$ metastable states have simplified features: they can be shown to
be collectively approximated by linear combinations of $M$ indicator
functions, each of which takes non-zero values over one of the metastable sets
and zero otherwise 22. Even when a good basis to extract such indicator
functions is not available, it can be shown that these functions can be
inferred by representing the evolution of non-ideal observables in terms of a
Hidden Markov Model (HMM) with $M$ hidden states 23.

However, fitting an HMM model to high-dimensional data given an unknown $M$
number of hidden states is generally non-trivial, requiring iterative methods
such as Expectation Maximization 24, 25, 26 or Sequential Monte Carlo
methods27, all of which are considerably more complex than the methods based
on the traditional linear Koopman approaches. Note that the task of choosing
“right” observables can be alleviated by using non-linear optimization methods
such as neural networks 28, 29 or kernel methods 30, 31. These methods often
require very large amounts of data and careful regularization to avoid over-
fitting, and hence require more expertise and careful application/fine-tuning
than linear Koopman operator methods, which possess an appealing simplicity of
use and interpretation.

In this study, we propose a simple and scalable alternative to traditional
HMM-inference algorithms based on the two-sample Kolmogorov–Smirnov (KS)32, 33
test and agglomerative clustering 34, hereinafter referred to as _KS
clustering_. This KS clustering is used to identify good “surrogate”
observables that conceptually correspond to indicator functions over
hidden/metastable states. The key idea is that the statistics of the time
evolution of even imperfect observables should contain information that can be
used to distinguish the metastable states a system visits, hence allowing one
to infer a good basis for slow eigenfunctions from imperfect observables.

It is worth noting an important trade-off in this algorithm: while it produces
accurate characteristic timescales, the resulting surrogate observables are
not explicit functions of the degrees of freedom of the target dynamical
system, and so the eigenfunctions that are produced cannot directly be
interpreted mechanistically. The KS clustering however possess an advantage
over methods where HMM membership functions are defined in the observable
space, since it can implicitly construct surrogate functions that cannot be
explicitly expressed in the observable space.

This article is organized as the follows. In Sections 2.1 and 2.2, we provide
theoretical background for the representation of the Koopman operator that can
be estimated from time series of observables obtained from stochastic systems.
Section 2.3 describes key concepts in HMM with multiple states that can be
included in indicator functions. Section 2.4 explains further why indicator
functions with the KS clustering algorithm work best in a reduced observable
space for defining hidden states. Section 2.5 illustrates the computation of
the Koopman operator with an HMM of two hidden states. In Section 3, we
demonstrate numerical results for the two-state HMM (Sec. 3.1) and apply the
KS clustering to three protein systems (Sec. 3.2). Finally, we discuss some
implications of the results in Section 4.

##  2 Theoretical background

###  2.1 Koopman representation of dynamical systems

Throughout this manuscript, we consider Markovian stochastic dynamics driving
the evolution of a thermal system whose microscopic state is denoted as
$\omega$ in a state space $\Omega$. For example, for a three-dimensional
$N$-atom molecular system with $\Omega=\mathbb{R}^{3N}$ evolving under
overdamped dynamics, the microscopic state can be fully characterized by
$\omega=(x_{1},\ldots,x_{N},y_{1}\ldots,y_{N},z_{1},\ldots,z_{N})$, which are
the Cartesian coordinates of all atoms. In the stochastic setting, an ensemble
of trajectories at time $t$ is characterized by $\rho(t,\omega)$ as a joint
probability density function in the continuous state space $\Omega$.

The dynamics is prescribed by an infinitesimal generator $\mathcal{L}$. To fix
ideas, we consider overdamped Langevin dynamics with Gaussian white noise,
$\mathcal{L}=-\sum_{i=1}^{3N}\left(\partial_{\omega_{i}}V\left(\omega\right)\right)\partial_{\omega_{i}}+2k_{B}T\partial_{\omega_{i}}^{2}$,
where $V(\omega)$ is the potential describing the interactions between the
atoms. For simplicity, we consider systems with detailed balance, which
guarantees reversibility 35. Note that the formalism itself is not specific to
Langevin dynamics, but is applicable to general reversible dynamics. The
infinitesimal generator $\mathcal{L}$ uniquely defines the evolution of the
probability density function,

$\frac{\partial}{\partial
t}\rho\left(t,\omega\right)=\mathcal{L}^{\dagger}\rho\left(t,\omega\right),$
(1)

where $\mathcal{L}^{\dagger}$ is the adjoint of $\mathcal{L}$. For reversible
dynamics, $\mathcal{L}=-\mathcal{L}^{\dagger}$. We assume that the stochastic
system is ergodic such that a unique stationary distribution,
$\rho_{\text{stat}}$, exists and satisfies
$\mathcal{L}^{\dagger}\rho_{\text{stat}}(\omega)=0$.

Let’s consider an observable $O$ that is a real-valued function of $\omega$
and define the so-called stochastic Koopman operator $\mathcal{K}_{t}$ 17 as

$\left(\mathcal{K}_{t}O\right)\left(\omega\right)\triangleq\mathbb{E}\left[O\left(\omega_{t}\right)|\omega_{0}=\omega\right],\forall\omega\in\Omega,$
(2)

where $\omega_{t}$ and $\omega_{0}$ denote the stochastic processes measured
at time $t$ and $0$, respectively. Using the semi-group notation, the finite-
time Koopman operator $\mathcal{K}_{t}=e^{t\mathcal{L}}$ maps the current
value of an observable $f$ to its expectation over the probability
distribution induced by the process $O_{t}$ at a later time $t$, given the
initial process $O_{0}$. A function $\phi$ is defined as a Koopman
eigenfunction if it satisfies $\left(\mathcal{K}_{t}\phi\right)=e^{\lambda
t}\phi$ or equivalently $\mathcal{L}\phi=\lambda\phi$ 14, 15, 36, 16. In this
study, we consider systems with point spectra, that is, systems with countable
$\lambda_{i}$, $i=1\ldots$, which can be ordered by their moduli, whose
corresponding eigenfunctions $\phi_{i}$ satisfy
$\mathcal{L}\phi_{i}=\lambda_{i}\phi_{i}$.

###  2.2 Data-driven estimation of the Koopman operator

The linearity of $\mathcal{K}_{t}$ and the correspondence between its
eigenvalues/eigenfunctions and those of the generator $\mathcal{L}$ discussed
in Sec. 2.1 can be leveraged to create powerful data-driven methods to
efficiently learn the characteristics of dynamics from timeseries of
observables. Methods such as TICA 6, VAMP 21, 17 and EDMD 13 provide linear
finite-dimensional approximations to the Koopman operator acting in the space
of selected observables. An estimated Koopman operator can be obtained by
minimizing the $L^{2}$-norm of the difference between the left- and right-hand
sides of Eq. (2) using a least-squares minimization over pairs of
configurations separated by a lag time $t$.

For brevity, we consider discrete-time Markov processes below, noting that it
is straightforward to generalize the analysis to continuous-time Markov
processes. In the discrete-time setting, an approximation to the Koopman
operator13, 37 given a sample path $\left\\{\omega_{t}\right\\}$,
$t=0,1,\ldots$, and a vector-valued observable $O$ is given by:

$K_{k}\equiv C(k)\cdot C^{-1}(0),$ (3)

where $C(k)$ with $k\in\mathbb{N}$ is the $k$-lag correlation of the
observable

$C(k):=\frac{1}{T}\sum_{s=0}^{T-1}O(\omega_{s+k})\otimes O(\omega_{s}),$ (4)

where the dyadic (outer) product is denoted by $\otimes$.

The eigenvalues of $K_{k}$ are approximations of $e^{\lambda_{i}k}$, and its
eigenvectors $\Phi_{i}$ can be used to approximate the eigenfunction
$\phi_{i}$ of the true Koopman eigenfunctions defined in Sec. 2.113, 37 as
$\phi_{i}(\omega_{s})\simeq O(\omega_{s})\cdot\Phi_{i}$, where $\cdot$ denotes
the inner product of the two vectors13. This simple and elegant procedure can
be shown to converge to the exact eigenvalues and eigenfunctions in the limit
of infinitely long timeseries of a set of observables $f$ which linearly span
a Koopman invariant subspace containing the corresponding set of Koopman
eigenfunctions $\left\\{\phi_{i}\right\\}_{i}$ 38, 37. Based on the Rayleigh-
Ritz method, Wu and Noé 17 showed that the data-driven estimation of the
Koopman operator is a variational problem when using methods such as TICA6, 7,
8, 9 and EDMD 13, with the approximate characteristic timescales approaching
the actual values from below in the infinite-data limit. When observables are
chosen as the first $M$ Koopman eigenfunctions,
$\left\\{\phi_{i}\right\\}_{i=1}^{M}$, the variational bound is tight, thus
recovering the optimal estimates of the Koopman eigenvalues and
eigenfunctions. In practice, selecting a finite number of observables is the
only option. This selection procedure is often system-specific, relying on
educated guesses or a priori information of an intuitive reaction coordinate.
The quality of the estimate can depend sensitively on a choice of observables,
as we now show.

###  2.3 Hidden Markov Chain as an effective model for describing systems with
multiple metastable states

In the following, we focus on the problem of characterizing the kinetics of
systems with metastable states, which are defined to have distinguishable
statistics, essential to the approximation of the Koopman operator discussed
in the previous sections. Let’s consider a system with $M$ such metastable
states and denote the $i^{\text{th}}$ metastable state, $i=1\ldots M$, by
$\Omega_{i}$ with $\Omega=\cup_{i=1}^{M}\Omega_{i}$, and
$\Omega_{i}\cap\Omega_{j}=\varnothing$ if $i\neq j$. Conceptually, a
metastable state is such that a typical trajectory would relax to quasi-
stationary distribution (QSD) within one such state much faster than it would
leave the state 39. Many systems in biology, chemistry, and materials science
exhibit strong metastability, which makes their study using direct simulation
methods such as molecular dynamics challenging due to the long waiting times
between state-to-state transitions. This setting implies the existence of a
slow subspace containing $M$ slow eigenvalues, well separated or statistically
distinguishable from the rest of the spectrum. As discussed above, indicator
functions of the form:

$\mathbf{1}_{\Omega_{i}}(\omega):=\left\\{\begin{array}[]{ll}1,&\text{if
}\omega\in\Omega_{i},\\\ 0,&\text{else},\end{array}\right.$ (5)

would form an excellent approximate basis for the global eigenfunctions of the
slow subspace of the Koopman operator $\mathcal{K}$.

Define the _local_ infinitesimal operator $\mathcal{L}_{i}^{\dagger}$ as the
operator describing the dynamics on $\Omega_{i}$, applying absorbing boundary
conditions on the boundary $\partial\Omega_{i}$ of $\Omega_{i}$, as in Ref.
39. Then, the largest eigenvalue and eigenfunction pair satisfies

$\mathcal{L}_{i}^{\dagger}\rho_{1}^{(i)}\left(\omega\right)=\lambda_{1}^{(i)}\rho_{1}^{(i)}\left(\omega\right),\quad
i=1\ldots M,$ (6)

where $\lambda_{1}^{(i)}<0$ implies decaying dynamics and
$\rho_{1}^{(i)}\left(\omega\right)=0$ if $\omega\notin\Omega_{i}$. The
eigenfunction $\rho_{1}^{(i)}(\omega)$ is referred to as the quasi-stationary
distribution (QSD) on state $\Omega_{i}$. Note that the $\lambda_{1}^{(i)}$
are now local quantities, in contrast to the discussion in the preceding
sections that focus on the eigenvalues of the global generator. Specifically,
the eigenvalue $\lambda_{1}^{(i)}$ quantifies the expected timescale
$-1/\lambda_{1}^{(i)}$, over which a system residing in state $\Omega_{i}$
would finally escape39. The metastability of the state can be quantified by
the ratio of the expected escaping time to the relaxation time to the QSD,
i.e., $(\lambda_{2}^{(i)}-\lambda_{1}^{(i)})/\lambda_{1}^{(i)}$. In the
following, we assume that all states are sufficiently metastable so that
$(\lambda_{2}^{(i)}-\lambda_{1}^{(i)})/\lambda_{1}^{(i)}\gg 1$, $\forall i$.

Let’s suppose we periodically observe the a trajectory of the system on some
timescale

$-1/(\lambda_{2}^{(i)}-\lambda_{1}^{(i)})\ll\tau\ll-1/\lambda_{1}^{(i)}$ (7)

We can set the timescale $\tau=1$ by choosing an appropriate unit of time.
With a probability $\sim 1-\exp(\lambda_{1}^{(i)})\lesssim 1$, we would
observe that the system remains in state $\Omega_{i}$; then, by construction,
the next observed configurations would be sampled from a distribution
approaching the QSD on state $\Omega_{i}$. Alternatively, with a probability
$\sim\exp(\lambda_{1}^{(i)})\ll 1$, we would observe that the system escaped
to another state $\Omega_{j}\neq\Omega_{i}$, where its configurations will be
sampled from the QSD on state $\Omega_{j}$. Crucially, when observed on the
timescale $\tau=1$ that is large compared to the internal relaxation time
within basins, the time series of any observable should be well approximated
by a sequence of i.i.d. random variables drawn from a distribution specific to
the QSD of the metastable state where the system is currently trapped. In
other words, such a time series will be well approximated by discrete-time
Hidden Markov Model (HMM) with $M$ hidden states, where the approximation
becomes increasingly good as each state becomes increasingly metastable and
statistically distinguishable (see Ref. 23 for an in-depth analysis). In this
representation, transitions between hidden states are described by a discrete-
time stochastic matrix $\mathbb{P}\left\\{S_{t+1}=i|S_{t}=j\right\\}$ with
$t=0,1,\ldots$ being the observation times and $S_{t}$ being a hidden state
realized in random processes. In HMM, we cannot directly measure $S_{t}$, but
rather, we measure an observable $O_{t}$ according to an observation model
$\mathbb{P}\left\\{O_{t}|S_{t}\right\\}$. In our case, the observation model
depends on the QSD measured on each state, i.e., $O_{t}=O(\omega_{t})$, where
$O$ is a prescribed deterministic function of the system’s state, and
$\omega_{t}$ is distributed according to the QSD of the states in which
$S_{t}$ resides.

This discussion highlights the important property of generic physical
observables: if hidden metastable states exist as described above, the
statistical distribution of generic observables should contain useful
information that can be used to implicitly reconstruct indicator functions
over hidden states. We now show how the statistics of these variations
following transitions between states dictate how well the corresponding slow
timescales can be estimated from the optimized Koopman operator.

###  2.4 Kolmogorov–Smirnov Clustering

The challenge in practice is to find a set of observables that can approximate
the indicator functions over metastable states as characterized in Sec. 2.3.
Our approach to this challenge is based on a simple characteristics of HMM: if
during two time intervals $[t_{1},t_{1}+\Delta t]$ and $[t_{2},t_{2}+\Delta
t]$ the dynamical system is in the same metastable state and delays between
observations obey Eq. (7), then the distributions of observables measured over
these time intervals should be statistically equivalent; and if the system is
in different metastable states during the two intervals, the corresponding
distributions of observables can presumably be statistically distinguishable
from one another. The “statistical distance” between distributions measured in
different time intervals can be quantified using the two-sample
Kolmogorov–Smirnov (KS) test 32. Namely, $D_{1,2}=\sup_{x}|F_{1}(x)-F_{2}(x)|$
(Figure 1a) is the KS statistic measuring the maximum difference between two
empirical cumulative distributions $F_{1}(x)$ and $F_{2}(x)$, which are
computed from the data collected in two time intervals $[t_{1},t_{1}+\Delta
t]$ and $[t_{2},t_{2}+\Delta t]$, respectively. A so-called distance matrix
contains all $D_{i,j}$ values, which are computed for all pairs of intervals.
Based on this distance matrix, indicator functions can be constructed (see
below).

In conventional applications, the KS statistics are used to reject the null
hypothesis that the two samples were drawn from the same underlying
distribution. In the present context, the KS statistics of the distributions
corresponding to every pair of intervals are instead used as a statistical
distance measure that allows for the clustering of all intervals into a number
of different groups. Specifically, each group will contain intervals that are
statistically similar to one another, while intervals that are very
statistically different will be assigned to different groups. This can be done
with any clustering methods that can operate from a user-provided pairwise
distance matrix; in the following, this was accomplished via hierarchical
agglomerative clustering40, which returns a hierarchy of clusters, i.e.,
clusters being merged in a bottom-up fashion until a preset number of cluster
or a critical inter-cluster distance threshold has been reached.

Figure 1: (Panel A) Illustration of the KS test. The test statistics
corresponds to the maximum difference between two cumulative distribution
functions measured. (Panel B) Schematic illustration of the proposed KS-
clustering approach.

Assuming that $M$ different groups are identified by such clustering, we can
build $M$ surrogate timeseries conceptually corresponding to indicator
functions over metastable states (rigorously, only $M-1$ such timeseries are
needed, since the $M^{\mathrm{th}}$ one can be expressed as a linear
combination of the $M-1$ others, up to an additive constant). Surrogate
indicator function $l_{i}$ with $1\leq i\leq M$ can be created by assigning a
value of $1$ to a given time interval when it was deemed a member of cluster
$i$ and a value of zero otherwise. This procedure can easily be generalized to
multiple observables using multi-dimensional generalizations of the KS test
41, or by defining the distance between two multi-dimensional distributions as
the maximal distance between any corresponding pair of one-dimensional
distributions, which we use in the following.

The overall algorithm is illustrated in Fig. 1. One first identifies a set of
base observables (1) which are processed using a conventional linear EDMD
procedure (2). The original timeseries are then compressed into a lower-
dimensional space via projection into the eigenvectors corresponding to the
slowest relaxation modes (3). In this reduced space, the KS clustering is
applied to identify indicator functions over the KS clusters (4), which are
then used to construct surrogate time-series (5). The projected descriptors
from step (3) are combined with the surrogate functions from step (5) and used
as input for a final EDMD analysis (6) to compute slow modes and corresponding
timescales (7).

###  2.5 A two-state HMM with Gaussian observation noise

To illustrate the arguments presented in Sec. 2.1-2.4, we consider an
analytically solvable HMM with two discrete state $S\in\left\\{1,2\right\\}$
and a single observable $O\in\mathbb{R}$ (note that this analysis can be
generalized to general $M$-state HMM). We use the standard notation that the
upper-case symbols with a subscript time stand for random processes, and
lower-case symbols stand for dummy variables or sample paths of random
processes. The Markov transition between the hidden states is characterized by
a Markov matrix $\bf{M}$ whose entries
$M_{ij}=\mathbb{P}\left\\{S_{t+1}=i|S_{t}=j\right\\}$:

$\bf{M}:=\begin{bmatrix}1-p_{+}&p_{-}\\\ p_{+}&1-p_{-}\end{bmatrix}.$ (8)

That is, with a probability $p_{+}$ (resp. $p_{-}$) the hidden state jumps
from 1 to 2 (resp. 2 to 1) in a single step. We note that the stationary
distribution of the hidden state
$\pi:=\left[\pi_{1},\pi_{2}\right]^{T}=\left[p_{-}/\left(p_{-}+p_{+}\right),p_{+}/\left(p_{-}+p_{+}\right)\right]^{T}$
satisfies $\mathbf{M}\cdot\pi=\pi$. As we are interested in systems whose
hidden states are metastable on the observation timescale, we have $p_{+}$,
$p_{-}\ll 1$. We consider a univariate Gaussian observation model, where the
observation $O_{t}$ at time $t$, which is a random variable (modeling the
quasi-stationary distribution), depends on only on the current hidden state
$S_{t}$:

$\displaystyle\rho\left(O_{t}=\omega|S_{t}=s\right)={}$
$\displaystyle\frac{1}{\sqrt{2\pi\sigma_{s}^{2}}}e^{-\frac{\left(\omega-\mu_{s}\right)^{2}}{2\sigma_{s}^{2}}}$
(9)

where $(\mu_{1},\sigma_{1})$ and $(\mu_{2},\sigma_{2})$ fully characterize the
observation model. We remark that the only timescale of the process is the
autocorrelation time of the hidden states, which is $-1/\log(1-p_{-}-p_{+})$.
Without loss of generality, we impose a zero-mean condition on the observable,
that is,
$\lim_{T\rightarrow\infty}\frac{1}{T}\sum_{i=0}^{T-1}\omega_{i}=\pi_{1}\mu_{1}+\pi_{2}\mu_{2}=0$.

Our goal is to analytically express the result of an EDMD procedure given an
infinitely long timeseries $O_{t}$’s. This corresponds to substituting the
analytical expression

$C(k):=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\sum_{s,s^{\prime}\in\left\\{1,2\right\\}}\,\mathbb{P}\left\\{\omega_{2}|s^{\prime}\right\\}\mathbb{P}\left\\{\omega_{1}|s\right\\}\,\left[\mathbf{M}^{k}\right]_{s^{\prime},s}\,\pi\left(s\right)\,\text{d}\omega_{1}\,\text{d}\omega_{2}.$
(10)

into Eq. 3. It is elementary to show that

$C(0)=\pi_{1}\left(\mu_{1}^{2}+\sigma_{1}^{2}\right)+\pi_{2}\left(\mu_{2}^{2}+\sigma_{2}^{2}\right)$
(11)

and for $k\in\mathbb{Z}_{\geq 1}$,

$C(k)=\left(1-p_{+}-p_{-}\right)^{k}\left(\mu_{1}-\mu_{2}\right)^{2}\pi_{1}\pi_{2}.$
(12)

Consequently, the estimated Koopman operator using this single observable with
$k$-lag is

$K_{k}=\gamma\left(1-p_{+}-p_{-}\right)^{k},$ (13)

where

$\gamma:=\frac{\left(\mu_{1}-\mu_{2}\right)^{2}\pi_{1}\pi_{2}}{\pi_{1}\left(\mu_{1}^{2}+\sigma_{1}^{2}\right)+\pi_{2}\left(\mu_{2}^{2}+\sigma_{2}^{2}\right)}.$
(14)

The only characteristic timescale of the process is then estimated as

$\tau_{k}=-\frac{k}{\log
K_{k}}=\frac{1}{-\log\left(1-p_{+}-p_{-}\right)+\frac{1}{k}\log\gamma}.$ (15)

For general $\sigma_{1}$, $\sigma_{2}\neq 0$ and a finite lag $k<\infty$,
$\gamma<1$ leads to an overestimation of the timescale $\tau_{k}$. This result
is consistent with the variational principle 17 which states that sub-optimal
observables result into overestimations of the characteristic timescales;
$-\log(\gamma)$ can be seen as a noise-to-signal metric which discounts the
timescale estimation by EDMD. While in principle the correct timescale is
recovered in the limit $k\rightarrow\infty$, approaching this limit could
require extremely long trajectories. In contrast, using zero-mean timeseries
of the indicator function over the hidden states as a basis results in an
unbiased estimation of $\tau_{k}=-1/\log\left(1-p_{+}-p_{-}\right)$ for any
value of $k$.

##  3 Numerical Results

Figure 2: Illustration of the KS-clustering approach using timeseries
generated with the Hidden Markov Models defined in Eqs. (8) and (9) with
$p_{+}=p_{-}=10^{-4}$ and $\mu_{1}=-0.5$, $\mu_{2}=0.5$. Algorithmic
parameters $\Delta t=500$ for the KS clustering procedure. Panels A-D:
$\sigma_{1}=\sigma_{2}=0.1$; Panels E-H: $\sigma_{1}=\sigma_{2}=0.5$. Panels A
and E are the time series generated by a standard kinetic Monte Carlo sampling
procedure, and the histogram of the time series are shown in panels B and F.
The red-colored indicator function based on the KS clustering is illustrated
in Panels C and G, as well as green-colored indicator functions (often called
discretized trajectories) based on a standard $k$-means clustering in Panels D
and H. 2, 42, 12

###  3.1 Two-state HMM with Gausssian noise

We now use the two-state HMM discussed above to numerically demonstrate that
the approximated indicator functions inferred from the KS clustering
significantly improve the timescale estimation. For each set of parameters
considered, a reference trajectory of $5\times 10^{5}$ steps is first
generated using a standard kinetic Monte Carlo procedure. In the following,
the transition probabilities between the hidden states was fixed at
$p_{+}=p_{-}=10^{-4}$; and the conditional means of the observation were set
to $\mu_{1}=-0.5$ and $\mu_{2}=0.5$, ensuring that the long-time observable
mean tends to 0. In the first parameter set, we consider
$\sigma_{1}=\sigma_{2}=0.1$, corresponding to a well-separated observation
model (see time series in Fig. 2A and empirical distribution in Fig. 2B). In
the second parameter set, we consider $\sigma_{1}=\sigma_{2}=0.5$; in this
case the stationary distribution is unimodal (see time series in Fig. 2E and
empirical distribution in Fig. 2F); this is representative of a situation
where the base observables are poor at distinguishing the two states (which is
a very common occurrence in practice).

In both parameter sets, the ground-truth characteristic timescale is
$-\log\left(1-p_{1}-p_{2}\right)\approx 5000$ discrete time steps. The
discounting factors $\gamma$ from Eq. (14) are $1/1.04$ and $1/2$,
respectively. These two cases are thus representative of good and bad
descriptors, where the convergence of EDMD with respect to the lag time would
be fast and slow as shown by the timescales (dashed red curve) in Fig. 3A and
B, respectively. Even for the “good” case (Fig. 3A), the slowest
characteristic timescale obtained from $\omega_{t}$ (raw data used for step
(2) in Fig. 1B) is severely underestimated for lag times less than 1000 steps.
The estimate is even worse in the “bad” case (Fig. 3B), where a 20%
underestimation still persists even for a lag of 10,000 steps.

We now demonstrate that the use of the KS indicator functions can greatly
improve the estimated timescales. To obtain the KS indicator functions (see
Sec. 2.4), we applied the procedure (steps 1-6) described in Fig. 1b. At step
(4), the data were split into time intervals of $\Delta t=500$, from which the
inter-interval KS distance matrix was computed. This matrix was then used as
an input to the hierarchical agglomerative clustering algorithm 40 with $M=2$
clusters. The KS indicator functions are plotted in Fig. 2C and G for both
cases. Note that since the dimension of the test system is one, steps (1-3)
can be skipped. As a point of comparison, we also considered a standard MSM
approach using k-mean clustering with 2 clusters; the corresponding
discretized trajectories are illustrated in Fig. 2D and H. For the “good”
case, Figures 2C and D show consistent indicator functions obtained from both
clustering algorithms. Fig. 3A shows that both of the clustering algorithms
yield accurate timescales at all lagtimes. When the distributions
corresponding to the two hidden states are not well separated (c.f. Fig. 2E
and F), Fig. 3B shows that the accuracy of the estimated timescale obtained
from the standard MSM with k-means clustering is basically identical to the
direct EDMD approach, while the KS-clustering approach accurately estimates
the slow timescale at all lagtimes. This results from the fact that the KS
clustering uses observable statistics over extended periods of time to
differentiate metastable states, while a k-means approach operates on each
sample separately, and hence cannot distinguish distributions that
significantly overlap with one another.

Figure 3: Timescales (Left axis) computed for the two-state Hidden Markov
Models. The dashed and solid black or red curves are the data obtained via
EDMD from $\omega_{t}$ and $\chi_{\text{ind}}$, respectively. Panel A:
$\sigma_{1}=\sigma_{2}=0.1$; Panel B: $\sigma_{1}=\sigma_{2}=0.5$. Label
$\omega_{t}$ indicates the estimated timescales obtained from the raw data via
EDMD without clustering. Label $\chi_{\text{ind}}$ indicates the timescales
obtained from the indicator functions computed via the KS clustering applied
to $\omega_{t}$. An ensemble of 100 time-series was used. Each timeseries has
$5\times 10^{5}$ steps. Figure 4: Timescale ($\tau_{k}$) computed for
different $\sigma_{1}=\sigma_{2}=\sigma$ and lagtime $k$ as functions of
$\Delta\mu=\mu_{2}-\mu_{1}=2\mu_{2}$, fixing $p_{1}=p_{2}=10^{-4}$. The dashed
green line is the ground true timescale, $\tau_{\text{truth}}=5000$ steps.
Here, we used the exact indicator functions ($\tilde{\chi}_{ind}$), which were
generated by the kinetic Monte Carlo procedure, for comparison with the
indicator function $\chi_{ind}$ generated by the KS clustering with number
$M=2$ of clusters used for all calculations here. Figure 5: Dependence of
timescales on $\Delta t$ and number of clusters ($M$), which is equal to the
number of indicator functions, $\chi_{\mathrm{ind}}$. Panels A-B:
$\sigma_{1}=\sigma_{2}=0.1$; Panels C-D: $\sigma_{1}=\sigma_{2}=0.5$; Panels A
and C: $M=2$; Panels B and D: $\Delta t=500$. Other parameters are the same as
in Figure 2. An ensemble of 10 timeseries was used to compute the confidence
intervals. Each timeseries has $5\times 10^{5}$ steps.

To further examine the effects of $\mu_{i}$ and $\sigma_{i}$ on the estimation
of timescales, we fixed $p_{+}=p_{-}=10^{-4}$ and varied
$\Delta\mu=\mu_{2}-\mu_{1}=2\mu_{2}$ and $\sigma_{1}=\sigma_{2}=\sigma$.
Equation (14) then simplifies to
$\gamma=1/\left[4(\sigma/\Delta\mu)^{2}+1\right]$. Figure 4 shows that the
results from direct EDMD on $\omega_{t}$ are underestimated in a way that is
in good agreement with the predictions of Eq. (15). In contrast, the use of
the indicator functions obtained from the KS clustering recovers reasonable
timescale even at short lagtimes and when the noise amplitude is very large
compared to $\Delta\mu$ in comparison with the exact indicator functions. The
very good agreement between the results obtained from inferred indicator
functions and actual hidden states (which are not normally accessible)
indicates that the fluctuations in predicted timescales result from the finite
trajectory lengths.

The KS clustering method requires selecting two tunable parameters: the length
of the bins $\Delta t$ used to partition the timeseries and the number of
target clusters $M$ used in the clustering algorithm. The impact of these two
parameters is shown in Fig. 5. Figures 5A and C show that $\Delta t$ around
500 steps produces accurate results, while larger values lead to an
overestimation of the characteristic timescale. This behavior can be
rationalized by considering the $\Delta t>>\tau$ limit. In this case, multiple
visits to different metastable states are averaged out within each interval.
This results into an underestimation of the transition probability and a
corresponding overestimation of the transition timescales. In order to be
accurate, on one hand, the algorithm requires that most time intervals contain
sections of trajectory remaining in the same hidden states, whose timescales
satisfy Eq. 7. On the other hand, $\Delta t$ should not be too small (e.g.,
$<50$ in our cases), so as to maximize the statistical power of the KS test
for accurately detecting transitions between metastable states. We therefore
recommend $\Delta t\simeq\tau_{M}/10$, where $\tau_{M}$ is the timescale of
the fastest mode in the target “slow” subspace. Of course, $\tau_{M}$ is not
known a priori, so $\Delta t$ can be estimated using the EDMD preprocessing
step discussed above. This estimate can be validated and readjusted after a
carrying out the whole procedure.

Figures 5B and D shows the effect of the number of clusters/hidden states $M$.
These results suggest that the estimated transition timescales can slightly be
overestimated when $M$ is larger than the number of hidden states (2 in this
case), although the dependence is rather weak. Consider the extreme limit
where $M$ is equal to the total number of time intervals, and hence each
interval will be assigned to its own cluster. In this case, it is possible for
EDMD to create a spurious linear combination of indicator functions whose
characteristic timescale approaches the trajectory length. We again recommend
using the EDMD preprocessing step discussed above to estimate the number of
slow characteristic timescales prior to the KS-clustering procedure.

###  3.2 Application to Protein Dynamics

Figure 6: Protein systems. A: Ribosomal Protein L9 (NTL9) (PDB: 2HBB), B:
Beta-beta-alpha Fold (BBA) (PDB: 1FME), and C: Trp-Cage (PDB: 2JOF). Lagtime =
100 steps or 200, 40 and 40 ns for NTL9, BBA and Trp-Cage, respectively;
$\Delta t=2100\text{ steps}$ or $4.2,0.84,\text{and }0.84\ \mu s$ for NTL9,
BBA and Trp-Cage, respectively. Middle Panels: The timeseries of the first
slowest components in the protein systems. Bottom Panels: the longest time
scale as a function of selected components.

To illustrate the performance of the KS clustering on more complex systems, we
consider three small proteins, NLT9, BBA, and Trp-Cage (Fig. 6A-C), each of
which was simulated for at least $200\text{\,}\mathrm{\SIUnitSymbolMicro
s}$19. We describe each trajectory by the time series of the torsional angles
of each amino acid in the proteins, indexed with $i$ 43, yielding 190, 152,
and 85 observables, respectively. We then applied EDMD to estimate the Koopman
operator for the three systems. The trajectories were then projected onto the
$M$ slowest eigenvectors of the Koopman operator. The resulting $M$ timeseries
were used as input features to the procedure (steps 1-6) shown in Figure 1B.
Note that when $M=1$, no surrogate indicator function is introduced and the
results are identical to that of a conventional EDMD approach.

Figure 6 shows the projection of the original trajectory into the slowest
eigenmode (middle row) of the Koopman operator inferred from the original data
(blue) and from the augmented data. While the data remains relatively noisy
when projected based on the torsional angles only (blue), the almost piecewise
character becomes evident when EDMD is augmented by the surrogate indicator
functions (black), which is a clear indication that improved eigenfunctions
are produced (see Sec. 2.5). This strongly suggests that the KS clustering
properly identifies the metastable states of the system. This improved
identification dramatically affects the estimated timescales compared to the
original EDMD (which corresponds to $M=1$), sometimes by up to an order of
magnitude (c.f. bottom row of Figures 6). The slowest characteristic timescale
is also observed to be relatively insensitive to a choice of $M$, beyond a
very large initial jump as soon as at least one surrogate indicator function
is added (i.e., when $M\geq 2)$.

As a more stringent test of the approach, we repeated the same procedure using
a single base observable, the RMSD with respect to the native state. Using
only the RMSD in an EDMD approach produces poor estimates of the timescales,
namely 0.25, 0.05, and $0.1\text{\,}\mathrm{\SIUnitSymbolMicro s}$ for NLT9,
BBA, and Trp-Cage, respectively. When combined with indicator functions
($M=2$, see SI) for NLT9, BBA, and Trp-Cage, the timescales are increased to
22.6, 2, and $3\text{\,}\mathrm{\SIUnitSymbolMicro s}$, respectively. These
estimated timescales from the RMSDs are close to the estimates from a much
richer set of observables 19 and timescales shown in Fig. 6. This suggests
that the KS-clustering is comparatively much less sensitive to the details of
the observable definition compared to the conventional EDMD approach, insofar
the observable distributions in different states is statistically
distinguishable.

##  4 Conclusion

EDMD-type methods are powerful kinetic analysis tools that are particularly
appealing due to their formal and practical simplicity. The results of this
type of analysis are however often sensitive to the ability of a given set of
observables to serve as a good basis for the slow subspace of the generator.
This would even be true of non-linear generalization of EDMD or even
conventional MSM approaches when the observable distributions corresponding to
different hidden states are not separable in the space spanned by the
observable. Rationally enriching the observable space is an outstanding
challenge that has yet to find a simple solution.

The main objective of this work was to propose a strategy to retain the
simplicity of linear EDMD-type methods, while improving the ability of the
method to extract accurate kinetic information from relatively poor
observables. The KS-clustering approach proposed here does not require the
distributions to be separable in the observable space, but only to be
sufficiently statistically distinct that those corresponding to different
hidden states can be identified using a two sample KS test. The surrogate
indicator functions obtained by clustering using the KS metric were shown
provide reliable estimates of known characteristic timescales, even when the
distribution of observables over the different metastable states strongly
overlap in the given observable space. It is however worth nothing that the
approximate eigenfunctions obtained by KS-clustering are not explicit
functions of the observables, which is certainly a limitation in terms of
their direct interpretability. However, the availability of timeseries
“labeled” in terms of putative hidden states can potentially be correlated
with known observables to interpret the nature of the transitions between
hidden states. When enriched with surrogate indicator functions over the
implicit metastable set, the simple linear EDMD procedure is shown to produce
very accurate timescale estimations at any lagtime, in contrast to
conventional linear approaches that can require very long lagtimes to produce
accurate estimates. The procedure is generic and scalable, and provides a
simple tool to improve linear approaches at low cost.

This work has been authored by employees of Triad National Security, LLC which
operates Los Alamos National Laboratory (LANL) under Contract No.
89233218CNA000001 with the U.S. Department of Energy/National Nuclear Security
Administration. The work has been supported by LDRD (Laboratory Directed
Research and Development) program at LANL under project 20190034ER (Massively-
Parallel Acceleration of the Dynamics of Complex Systems: a Data-Driven
Approach). V.A.N was partially supported by Director’s Postdoctoral
Fellowship, 20170692PRD4, for this work. V.A.N. is supported by Oak Ridge
National Laboratory, which is managed by UT-Battelle under Contract No. DE-
AC05-00OR22725 with the U.S. Department of Energy. This research used
resources of the Oak Ridge Leadership Computing Facility (OLCF). We also thank
DE Shaw for making their MD data available for this study.

Supporting information contains a PDF with one figure, and python scripts and
pandas pickles, which are made available at https://gitlab.com/ngoav/the-ks-
clustering.

## References

  * Noé and Fischer 2008 Noé, F.; Fischer, S. Transition Networks for Modeling the Kinetics of Conformational Change in Macromolecules. _Current Opinion in Structural Biology_ 2008, _18_ , 154–162
  * Pande et al. 2010 Pande, V. S.; Beauchamp, K.; Bowman, G. R. Everything You Wanted to Know about Markov State Models but Were Afraid to Ask. _Methods_ 2010, _52_ , 99–105
  * Prinz et al. 2011 Prinz, J.-H.; Wu, H.; Sarich, M.; Keller, B.; Senne, M.; Held, M.; Chodera, J. D.; Schütte, C.; Noé, F. Markov Models of Molecular Kinetics: Generation and Validation. _The Journal of Chemical Physics_ 2011, _134_ , 174105
  * Chodera and Noé 2014 Chodera, J. D.; Noé, F. Markov State Models of Biomolecular Conformational Dynamics. _Current Opinion in Structural Biology_ 2014, _25_ , 135–144
  * Wu et al. 2016 Wu, H.; Paul, F.; Wehmeyer, C.; Noé, F. Multiensemble Markov Models of Molecular Thermodynamics and Kinetics. _Proceedings of the National Academy of Sciences_ 2016, _113_ , E3221–E3230
  * Molgedey and Schuster 1994 Molgedey, L.; Schuster, H. G. Separation of a Mixture of Independent Signals Using Time Delayed Correlations. _Physical Review Letters_ 1994, _72_ , 3634–3637
  * Schwantes and Pande 2013 Schwantes, C. R.; Pande, V. S. Improvements in Markov State Model Construction Reveal Many Non-Native Interactions in the Folding of NTL9. _Journal of Chemical Theory and Computation_ 2013, _9_ , 2000–2009
  * Pérez-Hernández et al. 2013 Pérez-Hernández, G.; Paul, F.; Giorgino, T.; De Fabritiis, G.; Noé, F. Identification of Slow Molecular Order Parameters for Markov Model Construction. _The Journal of Chemical Physics_ 2013, _139_ , 015102
  * Naritomi and Fuchigami 2011 Naritomi, Y.; Fuchigami, S. Slow Dynamics in Protein Fluctuations Revealed by Time-Structure Based Independent Component Analysis: The Case of Domain Motions. _The Journal of Chemical Physics_ 2011, _134_ , 065101
  * Beauchamp et al. 2011 Beauchamp, K. A.; Bowman, G. R.; Lane, T. J.; Maibaum, L.; Haque, I. S.; Pande, V. S. MSMBuilder2: Modeling Conformational Dynamics on the Picosecond to Millisecond Scale. _Journal of Chemical Theory and Computation_ 2011, _7_ , 3412–3419
  * 11 Pande, V. S. MSMbuilder3.0
  * Scherer et al. 2015 Scherer, M. K.; Trendelkamp-Schroer, B.; Paul, F.; Pérez-Hernández, G.; Hoffmann, M.; Plattner, N.; Wehmeyer, C.; Prinz, J.-H.; Noé, F. PyEMMA 2: A Software Package for Estimation, Validation, and Analysis of Markov Models. _Journal of Chemical Theory and Computation_ 2015, _11_ , 5525–5542
  * Williams et al. 2015 Williams, M. O.; Kevrekidis, I. G.; Rowley, C. W. A data–driven approximation of the koopman operator: Extending dynamic mode decomposition. _Journal of Nonlinear Science_ 2015, _25_ , 1307–1346
  * Koopman 1931 Koopman, B. O. Hamiltonian systems and transformation in Hilbert space. _Proceedings of the National Academy of Sciences_ 1931, _17_ , 315–318
  * Koopman and Neumann 1932 Koopman, B. O.; Neumann, J. v. Dynamical systems of continuous spectra. _Proceedings of the National Academy of Sciences_ 1932, _18_ , 255–263
  * Mezic 2005 Mezic, I. Spectral Properties of Dynamical Systems, Model Reduction and Decompositions. _Nonlinear Dynamics_ 2005, _41_ , 309–325
  * Wu and Noé 2020 Wu, H.; Noé, F. Variational Approach for Learning Markov Processes from Time Series Data. _Journal of Nonlinear Science_ 2020, _30_ , 23–66
  * Scherer et al. 2019 Scherer, M. K.; Husic, B. E.; Hoffmann, M.; Paul, F.; Wu, H.; Noé, F. Variational Selection of Features for Molecular Kinetics. _The Journal of Chemical Physics_ 2019, _150_ , 194108
  * Lindorff-Larsen et al. 2011 Lindorff-Larsen, K.; Piana, S.; Dror, R. O.; Shaw, D. E. How Fast-Folding Proteins Fold. _Science_ 2011, _334_ , 517–520
  * Best et al. 2013 Best, R. B.; Hummer, G.; Eaton, W. A. Native Contacts Determine Protein Folding Mechanisms in Atomistic Simulations. _Proceedings of the National Academy of Sciences_ 2013, _110_ , 17874–17879
  * Wu et al. 2017 Wu, H.; Nüske, F.; Paul, F.; Klus, S.; Koltai, P.; Noé, F. Variational Koopman Models: Slow Collective Variables and Molecular Kinetics from Short off-Equilibrium Simulations. _The Journal of Chemical Physics_ 2017, _146_ , 154104
  * Röblitz and Weber 2013 Röblitz, S.; Weber, M. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. _Advances in Data Analysis and Classification_ 2013, _7_ , 147–179
  * Noé et al. 2013 Noé, F.; Wu, H.; Prinz, J.-H.; Plattner, N. Projected and Hidden Markov Models for Calculating Kinetics and Metastable States of Complex Molecules. _The Journal of Chemical Physics_ 2013, _139_ , 184114
  * Baum et al. 1970 Baum, L. E.; Petrie, T.; Soules, G.; Weiss, N. A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of Markov Chains. _The Annals of Mathematical Statistics_ 1970, _41_ , 164–171, Full publication date: Feb., 1970
  * Noé et al. 2013 Noé, F.; Wu, H.; Prinz, J.-H.; Plattner, N. Projected and hidden Markov models for calculating kinetics and metastable states of complex molecules. _The Journal of Chemical Physics_ 2013, _139_ , 184114
  * Rabiner 1989 Rabiner, L. A tutorial on hidden Markov models and selected applications in speech recognition. _Proceedings of the IEEE_ 1989, _77_ , 257–286
  * Smith 2013 Smith, A. _Sequential Monte Carlo methods in practice_ ; Springer Science & Business Media, 2013
  * Mardt et al. 2018 Mardt, A.; Pasquali, L.; Wu, H.; Noé, F. VAMPnets for deep learning of molecular kinetics. _Nature Communications_ 2018, _9_ , 5
  * Yeung et al. 2017 Yeung, E.; Kundu, S.; Hodas, N. O. Learning Deep Neural Network Representations for Koopman Operators of Nonlinear Dynamical Systems. _2019 American Control Conference (ACC)_ 2017, 4832–4839
  * Schwantes and Pande 2015 Schwantes, C. R.; Pande, V. S. Modeling Molecular Kinetics with tICA and the Kernel Trick. _Journal of Chemical Theory and Computation_ 2015, _11_ , 600–608, PMID: 26528090
  * Klus et al. 2019 Klus, S.; Husic, B. E.; Mollenhauer, M.; Noé, F. Kernel methods for detecting coherent structures in dynamical data. _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 2019, _29_ , 123112
  * Smirnov 1948 Smirnov, N. Table for Estimating the Goodness of Fit of Empirical Distributions. _The Annals of Mathematical Statistics_ 1948, _19_ , 279–281
  * Dimitrova et al. 2020 Dimitrova, D. S.; Kaishev, V. K.; Tan, S. Computing the Kolmogorov-Smirnov Distribution When the Underlying CDF Is Purely Discrete, Mixed, or Continuous. _Journal of Statistical Software_ 2020, _95_
  * Sokal and Rohlf 1962 Sokal, R. R.; Rohlf, F. J. THE COMPARISON OF DENDROGRAMS BY OBJECTIVE METHODS. _TAXON_ 1962, _11_ , 33–40
  * Klus et al. 2018 Klus, S.; Nüske, F.; Koltai, P.; Wu, H.; Kevrekidis, I.; Schütte, C.; Noé, F. Data-Driven Model Reduction and Transfer Operator Approximation. _Journal of Nonlinear Science_ 2018, _28_ , 985–1010
  * v. Neumann 1932 v. Neumann, J. Zur Operatorenmethode In Der Klassischen Mechanik. _Annals of Mathematics_ 1932, _33_ , 587–642
  * Lin et al. 2021 Lin, Y. T.; Tian, Y.; Anghel, M.; Livescu, D. Data-driven learning for the Mori-Zwanzig formalism: a generalization of the Koopman learning framework. 2021
  * Brunton et al. 2016 Brunton, S. L.; Brunton, B. W.; Proctor, J. L.; Kutz, J. N. Koopman Invariant Subspaces and Finite Linear Representations of Nonlinear Dynamical Systems for Control. _PLOS ONE_ 2016, _11_ , 1–19
  * Le Bris et al. 2012 Le Bris, C.; Lelievre, T.; Luskin, M.; Perez, D. A mathematical formalization of the parallel replica dynamics. _Monte Carlo Methods and Applications_ 2012,
  * Bouguettaya et al. 2015 Bouguettaya, A.; Yu, Q.; Liu, X.; Zhou, X.; Song, A. Efficient Agglomerative Hierarchical Clustering. _Expert Systems with Applications_ 2015, _42_ , 2785–2797
  * Fasano and Franceschini 1987 Fasano, G.; Franceschini, A. A multidimensional version of the Kolmogorov–Smirnov test. _Monthly Notices of the Royal Astronomical Society_ 1987, _225_ , 155–170
  * Lane et al. 2013 Lane, T. J.; Shukla, D.; Beauchamp, K. A.; Pande, V. S. To milliseconds and beyond: challenges in the simulation of protein folding. _Current Opinion in Structural Biology_ 2013, _23_ , 58–65, Folding and binding / Protein-nucleic acid interactions
  * McGibbon et al. 2015 McGibbon, R. T.; Beauchamp, K. A.; Harrigan, M. P.; Klein, C.; Swails, J. M.; Hernández, C. X.; Schwantes, C. R.; Wang, L.-P.; Lane, T. J.; Pande, V. S. MDTraj: A Modern Open Library for the Analysis of Molecular Dynamics Trajectories. _Biophysical Journal_ 2015, _109_ , 1528–1532