# Exact Topology of Dynamic Probability Surface of an Activated Process by
Persistent Homology

Farid Manuchehrfar  Huiyu Li  Wei Tian  Ao Ma<EMAIL_ADDRESS>Jie Liang
<EMAIL_ADDRESS>Center for Bioinformatics and Quantiative Biology and
Department of Bioengneering, University of Illinois at Chicago, Chicago, IL
60607.

###### Abstract

To gain insight into reaction mechanism of activated processes, we introduce
an exact approach for quantifying the topology of high-dimensional probability
surfaces of the underlying dynamic processes. Instead of Morse indexes, we
study the homology groups of a sequence of superlevel sets of the probability
surface over high-dimensional configuration spaces using persistent homology.
For alanine-dipeptide isomerization, a prototype of activated processes, we
identify locations of probability peaks and connecting-ridges, along with
measures of their global prominence. Instead of a saddle-point, the transition
state ensemble (TSE) of conformations are at the most prominent probability
peak after reactants/products, when proper reaction coordinates are included.
Intuition-based models, even those exhibiting a double-well, fail to capture
the dynamics of the activated process. Peak occurrence, prominence, and
locations can be distorted upon subspace projection. While principal component
analysis account for conformational variance, it inflates the complexity of
the surface topology and destroy dynamic properties of the topological
features. In contrast, TSE emerges naturally as the most prominent peak beyond
the reactant/product basins, when projected to a subspace of minimum dimension
containing the reaction coordinates. Our approach is general and can be
applied to investigate the topology of high-dimensional probability surfaces
of other activated process.

###### keywords:

surface topology, landscape analysis, persistent homology, active process,
transition state, energy flow

These authors contributed equally These authors contributed equally These
authors contributed equally

##  1 Introduction

Activated processes are ubiquitous in molecular systems, ranging from chemical
reactions of small molecules to dynamic conformational changes and enzymatic
reactions of proteins. In proteins, all functionally important processes are
activated processes, which provide well-defined rates essential for proteins
to carry out their roles in the cellular context, as proper timing is required
for proper function.

The prevalent picture describing an activated process is that of a transition
between two meta-stable basins on the free energy landscape separated by a
barrier, whose height is large compared to thermal energy 1. The slow time
scale of activation arises from the fact that the molecular system can rarely
accumulate sufficient energy in the relevant degrees of freedom (DoFs) to
surpass the transition barrier. This simple and elegant picture originates
from reaction rate theories, such as the well-known transition state theory
and Kramer’s theory 1, 2, 3, 4, 5, 6, 7 developed in studies of the dynamics
of chemical reactions of small molecules.

A key concept in reaction rate theories is that of reaction coordinates: a few
special coordinates exists that can fully determine the progress of a reaction
process 8, 9, 10. A requirement for reaction coordinates is that they must
accurately locate the transition barrier. Accordingly, the numerous degree of
freedoms (DoFs) in a complex molecular system (e.g., a protein molecule, a
system of solute and solvent) can be divided into reaction coordinates and
heat bath. Reaction coordinates play central roles as they determine both the
mechanism and the rate of activation. For example, to modify the activity of
an enzyme, one should modify residues involved in the reaction coordinates of
the enzyme activities 11, 12, as this will modify both the reaction pathway
and the barrier height for activation. In contrast, modifying residues that
belong to the heat bath will not alter the enzymatic activity, as the role of
the heat bath is to provide energy to the reaction coordinates to cross the
activation barrier during rare fluctuations, which is largely a non-specific
process.

Given such significance, it is important to develop a rigorous and
quantitative criterion for determining the correct reaction coordinates. This
task was accomplished with the development of the procedure of committor test,
which is characterized by the committor value $p_{B}$ 8, 9, 13, 14: the
probability that a dynamic trajectory initiated from a given configuration to
reach the product basin before visiting the reactant basin. By definition, the
reactant and product states have committor values of 0 and 1, respectively,
whereas the optimal transition state coincides with $p_{B}=0.5$. The committor
value $p_{B}$ therefore provides a rigorous parameterization of the reaction
process. Thus, the intuitive, albeit qualitative, notion of reaction
coordinates translates into a rigorous definition of the few coordinates that
are sufficient for determining the committor value of any given configuration.
Du et. al first adopted this rigorous definition of reaction coordinates in
the context of protein folding 9; Chandler and co-workers established its
usage as a standard practice in the general context of activated processes 8.

While this rigorous criterion has been well accepted, identifying the correct
reaction coordinates turns out to be rather difficult, even for systems of
modest complexity. One example is the $C_{7eq}\rightarrow C_{7az}$
isomerization reaction of the alanine dipeptide in vacuum, a prototype for
studying biomolecular conformational transitions. Alanine dipeptide is the
smallest molecule that satisfies the criterion that distinguishes complex
molecules from small molecules: the non-reaction coordinates in the system
constitute a large enough thermal bath to provide the reaction coordinates
with adequate energy to cross the activation barrier. As a result, the
$C_{7eq}\rightarrow C_{7az}$ transition displays features of activated
dynamics that are unique to complex molecules but absent in small molecules.
It was first found by Bolhuis et. al 15 that the conventional Ramachandran
torsional angles $\phi$ and $\psi$, while sufficient for distinguishing the
two stable basins, are inadequate for locating the transition state. Instead,
another torsional angle $\theta_{1}$ was found to be an essential reaction
coordinate––a rather counter-intuitive finding 16, 17. The counter-intuitive
nature of reaction coordinates turned out to be more often the norm than
exception in complex systems 15, 18, 19, posing a formidable challenge, as
sans intuition, there appears no guidance in sight.

The challenge in identifying reaction coordinates had motivated efforts in
developing rigorous methods for their identification in complex systems since
the early 2000’s 8, 9, 15, 18, 20, 21. Beyond unreliable intuition and trial-
and-error, the first systematic method was that of machine-learning, in which
a neural network was used to automatically identify the optimal reaction
coordinates from a prepared pool of candidates 21. This method was used to
successfully identify the key solvent coordinate that controls the
isomerization dynamics of an alanine dipeptide in solution, which had defied
prior intuition-based trial-and-error efforts. The success of this machine
learning based approach lead a series of further developments along similar
lines 10, 18, 22, 23, 24, 25, 26, 27, 28, 19.

However, a major deficiency of machine learning-based methods is that they
cannot answer the real question concerning reaction coordinates – why some
coordinates are more important for activation than the others? Instead, these
methods can only inform us empirically which coordinates appear to be
important based on well-defined criteria. Recently, Ma and co-workers
developed a rigorous theory for mapping out the flow of potential energy
through individual coordinates 29, 17. It was found that the reaction
coordinates are the coordinates that carry high energy flows during the
activation process. This result suggested an appealing physical picture:
energy flows from fast coordinates into slow coordinates during activation, so
that adequate energy can accumulate in the slow coordinates, enabling them to
cross the activation barrier on path of these slow coordinates. This physical
picture also suggested that reaction coordinates are preferred channels of
energy flows and are encoded in the protein structure. Through analysis of
energy flow, one can obtained a prioritized list of coordinates that likely
play most significant roles in the activation process.

The most celebrated concept in reaction rate theories is that of the
transition state (TS), which is the dynamic bottleneck of an activated
process. Conventional thoughts are that they are located at a critical point
with Mores index of 1 on the high-dimensional potential energy landscape of
the molecule. If all the reaction coordinates of an activated process are
known, the TS will be an index-1 critical point on the free energy surface of
the reaction coordinates, the sole direction with negative Hessian pointing
along the ideal one-dimensional reaction coordinate valid at the top of the
activation barrier. Based on this picture, the surface of the multi-
dimensional probability distribution of reaction coordinates constructed from
an ensemble of reactive trajectories should be highly structured, with the
important dynamic states (e.g. TS) manifesting as critical topological
features.

To gain insights into reaction mechanism, a common practice is to study
features of the free energy surface. However, all relevant information of an
activated process is contained in the reactive trajectories. Instead of the
free energy surface, one can construct the dynamic probability surface of the
transition state. For certain systems, this can be achieved using the
transition path sampling (TPS) method. Once a large ensemble of reactive
trajectories are generated, an ensemble of configurations that the system
sampled during the transition process can be harvested. From this ensemble,
one can generate a dynamic probability landscape or transition state surface,
which is usually high-dimensional.

The focus of this work is the analysis of the exact topology of the high-
dimensional dynamic probability landscape, and the establishment of its
relationship with the transition state of the active process. Instead of Morse
index, we characterize the high-dimensional transition state surface by its
topological structures in homology groups. We analyze topological changes in
the superlevel set of the probability surface at different probability levels.
Using the technique of persistent homology, we identify the locations of
probability peaks in the high-dimensional configuration space and ridges
connecting them, along with measures of their global prominence.

We apply this approach to study the active process of $C_{7eq}\rightarrow
C_{7ax}$ isomerization of the alanine dipeptide, a well-characterized model
system for studying protein conformational changes 15, 29, 17. After the exact
topological structures of the dynamic probability surface are constructed, we
identify the location of the ensemble of transition state conformations.
Instead of a saddle point with a Morse index of 1, the transition state
ensemble (TSE) is found to be located at the top of the most prominent peak
after those of the reactant and product. In addition, the dynamically
important topological structures are retained when the surface is projected
onto the 2-dimensional plane of ($\phi-\theta_{1}$), which are known to be the
reaction coordinates 17. In contrast, when projected to the intuition-derived
($\phi-\theta$)-plane, the topological features of the probability surface no
longer contain dynamic information on the transition state. Furthermore, we
find that PCA dimension reduction distorts surface topologies, such that the
transition state ensemble cannot be recovered from PCA-derived topological
features: Instead of simplification, PCA destroys the dynamic properties of
topological features of the original transition state surface.

Overall, we introduce a novel approach for quantifying the exact topology of
high-dimensional surfaces. With this approach, we have characterized the
precise topology of the transition state surface of the active process of
alanine dipeptide isomerization. We have also established that the TSE are
located at the most prominent probability peak beyond those of the reactant
and the product, when the subspace of projection contain the proper reaction
coordinates. The new language of homology group and the technique of
persistent homology on superlevel sets of high-dimensional probability
surfaces introduced here are general and can be applied to investigations of
the topology of high-dimensional surfaces over configuration space encountered
in other problems of activated process.

##  2 Theory, Models and Methods

We briefly discuss how the dynamic probability surface can be constructed for
the molecular system of our interests. We then discuss the problem of
understanding 118

the topological structures of the probability surface over the relevant
configuration space. Our focus will be the introduction of the method of
homology groups for analyzing the exact topological structures of the dynamic
probability surface. This approach is based on the homological structures of
the sequence of the superlevel sets of a probability surface, and differs from
previous efforts that are based on critical points, Morse theory, and Euler
characteristics.

###  2.1 Constructing dynamic probability surface: transition state ensemble
and path sampling

We use the transition path sampling (TPS) method 8 to generate a sufficiently
large ensemble of reactive trajectories for the active process we study. To
ensure the transition process is fully covered without bias, the duration of
the trajectories is much longer than that of the transition process. Along
each trajectories, we further harvested a number of system configuration with
a specific time intervals to generate an ensemble of the relevant
configurations that the system samples during the transition process. From
this ensemble, we construct a dynamic probability landscape over a
$d$-dimensional subspace, namely, the joint probability distribution over the
$d$-coordinates.

###  2.2 Configuration space and probability surface.

We now discuss the general problem of analyzing probability surface over a
configuration space of arbitrary dimension. We first introduce a few relevant
concepts, and then review current approaches. This is followed by an
exposition of the basics of homology groups and persistent homology. We focus
on developing these concepts in the setting of cubics and superlevel sets of
the probability surface. We will also discuss the key concept of filtration
and its construction, and how it can be used to uncover the exact high-
dimensional topological features.

#### Configuration space.

We begin our discussion in a general setting. We use $d$ number of features
that describe the configurations of a molecule. These can be bond lengths,
bond angles, and torsional angles describing the structure of the molecule.
For alanine dipeptide, there are a total of $60$ features that fully describes
the configuration of the molecule. For a molecule of finite size, the
configuration space $\mathbb{M}\subset{}^{d}$ is compact. It is likely that
$\mathbb{M}$ lies in a subspace of the Euclidean space d, as there are
coupling between different degrees of freedom.

#### Probability surface.

Each configuration
${\mbox{\boldmath$x$}}=(x_{1},\,x_{2},\cdots,x_{d})\in\mathbb{M}$ has a
probability $f({\mbox{\boldmath$x$}})\in[0,\,1]$ associated with it. Namely,
we have a function $f:\mathbb{M}\rightarrow{}_{[0,\,1]}$ that assigns the
probability value $p({\mbox{\boldmath$x$}})$ to each configuration $x$.

#### Superlevel sets and sublevel sets.

For a value $0\leq a\leq 1$, we can identify all $x$ whose probability values
$f({\mbox{\boldmath$x$}})\geq a$, which is called the superlevel set
$\mathbb{M}_{f\geq a}$:

$\mathbb{M}_{f\geq
a}\equiv\\{{\mbox{\boldmath$x$}}\in\mathbb{M}|f({\mbox{\boldmath$x$}})\geq
a\\}=f^{-1}([a,\,1)).$

Similarly, we can have the the sublevel set $\mathbb{M}_{f\leq a}$:

$\mathbb{M}_{f\leq
a}\equiv\\{{\mbox{\boldmath$x$}}\in\mathbb{M}|f({\mbox{\boldmath$x$}})\leq
a\\}=f^{-1}((0,\,a]).$

###  2.3 Topology of probability surface: Critical points and Morse indices.

It is of great importance to understand the topological structures of the
probability surface $f({\mbox{\boldmath$x$}})$. The approach of analyzing the
critical points is well-practiced for exploring the topological structures of
high-dimensional surfaces.

#### Critical Points.

The critical points of $f(\mathbb{M})$ are where all first derivatives of $f$
vanishes:

$\frac{\partial f({\mbox{\boldmath$x$}})}{\partial x_{1}}=0,\,\frac{\partial
f({\mbox{\boldmath$x$}})}{\partial x_{2}}=0,\,\cdots,\,\frac{\partial
f({\mbox{\boldmath$x$}})}{\partial x_{d}}=0.$

Critical points are coordinates independent and can be further classified into
different types by the secondary derivatives. We can organize the secondary
derivatives into a $d\times d$ Hessian matrix
${\mbox{\boldmath$H$}}f({\mbox{\boldmath$x$}})$, whose entries are:

$({\mbox{\boldmath$H$}}f)_{i,\,j}=\frac{\partial^{2}f}{\partial x_{i}\partial
x_{j}}.$

At non-degenerative critical points, where the Hessian matrices are non-
singular, the Hessian will have a mixture of positive and negative
eigenvalues. The number of negative eigenvalues $\eta$ is the Morse index of
the critical point.

#### Topology of surface by critical points and Morse theory.

At a critical point, the topology of sublevel sets changes. For a critical
point $x$ with $f({\mbox{\boldmath$x$}})=a$, consider the sublevel sets
$\mathbb{M}_{f<a}$ and the sublevel set $\mathbb{M}_{f<(a+\epsilon)}$ slightly
above it by a small amount of $\epsilon>0$. The topology of $\mathbb{M}_{f<a}$
and $\mathbb{M}_{f<(a+\epsilon)}$ are different, as one cannot be deformed
into another. However, once an $\eta$-dimensional handle is attached to
$\mathbb{M}_{f<a}$ at the critical point $x$, it has the same homotopy type as
$\mathbb{M}_{f<(a+\epsilon)}$ 30, namely, these two can be deformed into one
another. In fact, the homotopy class of $\mathbb{M}$ can be characterized by
toplological changes at critical points. The problem of determining the
topology of the surface $f({\mbox{\boldmath$x$}})$ on $\mathbb{M}$ then
becomes the problem of determining the critical points of $f(x)$ on
$\mathbb{M}$ and their Morse indices.

#### Euler characteristics.

A special case of the celebrated Morse inequalities relates the number of
critical points of different indices to the Euler characteristics, which is

$\chi(\mathbb{M})=\sum_{k=1}^{d}(-1)^{k}m_{k},$

where $m_{k}$ is the number of critical points of index $k$. The Euler
characteristics provides information on the number of holes of various
dimensions, and can also be written as:

$\chi(\mathbb{M})=\sum_{k=1}^{d}(-1)^{k}\beta_{k},$

where $\beta_{k}$ is the $k$-th Betti number, which counts the number of
$k$-dimensional holes. For example, in 3, $\beta_{0}$ counts the number of
connected components, $\beta_{1}$ counts the number of tunnels, and
$\beta_{3}$ counts the number of voids. Voids and tunnels have been
extensively examined in the studies of protein structures and functions 31,
32, 33, 34, 35, 36, 37, 38, 39, 40.

Persistent homology has also been applied in analysis of structures of
chemical compounds and proteins, and for correlation with their biophysical
properties as well as data analysis 41, 42, 43.

#### Prior work on characterizing topological surface of molecular system.

The topology of surfaces over configuration space has been the subject of many
investigations 44, 45, 46, 47, 48 . For example, thermodynamic phase
transitions have been explored from the viewpoint of topological changes of
submanifolds of the configuration space. Changes in Euler characteristics have
been found to signal the occurrence of phase transitions for certain systems
45, 46.

However, it is difficult in practice to understand the global surface topology
using critical points and Morse indices. There are only a few model systems
where all critical points are known analytically 49. Numerical computation
faces numerous challenges. First, it is difficult to identify all critical
points. Methods based on Newton-Ralphson and other techniques require initial
guesses and do not guarantee identifications of all critical points.
Furthermore, many initial guesses fall into the same large basins of
attraction and yield no new information. Second, as the probability surface
can only be constructed approximated from points sampled in the high-
dimensional configuration space, the degree of sampling may not be
sufficiently detailed to capture the original topology of the configuration
space. Third, sampling has to be sufficiently detailed to accurately measure
the first and second derivatives along each coordinate direction at all
locations. Fourth, as derivatives reflect local properties of the probability
surface, there can be numerous critical points that may be trivial and of
little importance if the probability surface is rugged. It is difficulty to
distinguish those reflecting the global features of the surface from those
reflecting local dimples or even those that are due to noise in sampling.
While there have been efforts in visualizations of potential surface of
molecules and other heuristics (see 47 and references within), to our
knowledge, it is not yet possible to characterize all critical points and the
homotopy type of the probability or potential surface of a molecule in three
dimensional space.

###  2.4 Topology of probability surface: Homology group and persistent
homology

#### Background and overview.

Instead of analyzing the critical points, we adopt a different approach. We
are interested in global features such as the occurrence of different peaks,
and how they are connected. This can be achieved by examining globally the
structures of holes of various dimension and how such structures changes at
different sublevel or superlevel sets of the probability landscape. We adopt
an approach based on the theory of homology group and persistent homology 50,
51, 52. Homology group studies holes in topological spaces and is a classic
topic in algebraic topology 53, 54. Persistent homology computes these holes
and measures their scales at different spatial resolutions 51, 52. Compared to
homotopy, homology groups are more amenable to computation. Our approach is
only feasible due to recent progress in computational topology and topological
data analysis 51, 55, 52, 50, 56.

To provide an intuitive picture illustrating this approach, we envision a sea
level on top of the probability landscape over the configuration space (Fig
1). We are interested in how mountain peaks emerge from the sea when the sea
level is lowered gradually, and how independent mountain peaks become
connected by land-ridges when the sea level is lowered further. These are
related to $0$-dimensional holes, which are connected components.

Figure 1: Sea levels of the probability landscape $f({\mbox{\boldmath$x$}})$
on the configuration space $\mathbb{M}$. The superlevel sets
$\mathbb{M}_{i}=\mathbb{M}_{f\geq a_{i}}=f^{-1}(\geq a_{i})$ at different sea
level (white regions) can have different topology, with different number of
components shown in white.

#### Complex and chain.

We first discuss how to represent the $d$-dimensional configuration space
$\mathbb{M}$ 57. In this study, we use cubic complexes 57, 56. A
$d$-dimensional cubic complex $K$ is constructed from a union of points, line
segments, squares, cubes, and their $k$-dimensional counterparts glue together
properly, where $k\leq d$ and all have unit length (except points, which have
no lengths). We call each of these a $k$-cell or a $k$-cube (see Fig 2a for a
3-cell). While the topology of $\mathbb{M}$ is invariant whether it is
represented by cubic complexes or other complexes such as simplicial
complexes, the nature of grid representation of the molecular configurations
makes this choice convenient 56.

We can build up our cubic complexe $K$ from cubes to represent the
configuration spaces. Consider a set of $k$-cells, we can sum them up. We call
the total summation of a set of $k$ cells a $k$-chain. Fig 2b shows an example
where two 3-cells are summed up to form a 3-chain. Fig 2c shows how nine
2-cells are summed up to form a 2-chain. Here the binary operation of summing
over two $k$-cells is orientation sensitive: two $k$-cells of the same
underlying space but opposite orientation cancel out each other when summed
up. If a set is equipped with a binary operation satisfying certain
requirements, it is called a group mathematically. The set of $k$-chains from
the $K$ complex with our binary operation of summation therefore form a chain
group $\mathsf{C}_{k}(K)$.

Figure 2: An illustration of cubic complex. a). A 3-cubic cell, with the
orientation of its 2-faces shown. b). Two 3-cubes are summed to form a
3-chain. The internal square is contributed twice from the two cubes. As each
surface is oriented (i.e., counter-clock-wise by the outward surface normal),
these two squares have opposite orientations and cancel each other when
summed. c). An example of a 2-chain formed by 2-cells and its boundary.

#### Boundaries.

We now set boundaries. The boundary of an individual $k$-cell is the set of
its $(k-1)$-dimensional faces, which by definition forms a $(k-1)$-chain. The
boundary of a 3-cubic cell is shown in Fig 2a, which is the set of the six
oriented squares. The boundary of a $k$-chain is the sum of the boundaries of
its element $k$-cells. Because of the nature of our sum operation, internal
structures cancel out. Consider the boundary of the two neighboring three-
dimensional cubes in Fig 2b. The interfacial square is contributed twice, once
from each cube, but with opposite orientation as both are counter clock-wise
around their outwards normals. When these two cubes are glued together, these
two boundary squares are summed up, and they cancel each other out. The
overall outcome of this summation is indeed the outer boundary of the union of
the two neighboring cubes. Fig 2c shows a 2-chain and its boundary. This holds
true in other dimensions as well, namely, a $(k-1)$-dimensional face from two
neighboring $k$-cells have opposite orientations and cancels each other out
upon summation.

With this summation, we obtain the boundary of a $k$-chain from $K$ by
applying the boundary operator $\partial_{k}$:

$\partial_{k}:\mathsf{C}_{k}(K)\rightarrow\mathsf{C}_{k-1}(K).$ (1)

#### Cycles.

There are certain $k$-chains that have no boundaries. They are called
$k$-dimensional cycles or $k$-cycles. With the binary operation of summation
discussed earlier, the set of $k$-cycles from $K$ form the cycle group
$\mathsf{Z}_{k}(K)$:

$\mathsf{Z}_{k}(K)\equiv\\{c\in\mathsf{C}_{k}(K)|\,\partial_{k}c=\emptyset\\}.$

As an example, we consider the three-dimensional cube again (Fig 2a). We take
its six surface squares that fully enclose the solid cube. These square form a
2-chain. As a whole, this 2-chain itself does not have boundaries, as the six
squares are glued together along the borders and there are no openings.

#### Kernel and image.

Analogous to the null space or kernel in linear algebra, the cycle group
$\mathsf{Z}_{k}(K)$ is the kernel of the operator $\partial_{k}$, as each of
its member $k$-chain has no boundary by definition and $\partial_{k}$ will
send it to null:

$\mathsf{Z}_{k}(K)=\ker\mathsf{C}_{k}(K)\equiv\\{c\in\mathsf{C}_{k}(K)|\partial_{k}c=\emptyset\\}.$

We now move one dimension up and consider boundaries of $(k+1)$-chains in $K$.
Boundaries are one dimension lower and therefore the boundaries of
$(k+1)$-chains are $k$-chains. They are called the $k$-boundaries of $K$ and
form the $k$-boundary group $\mathsf{B}_{k}(K)$. As each $k$-boundary is
obtained when $\partial_{k+1}$ is applied to a $(k+1)$-chain, collectively
they are the image of $\partial_{k+1}$:

$\mathsf{B}_{k}(K)=\text{im}\,\partial_{k+1}(\mathsf{C}_{k+1}(K))\equiv\\{c\in\mathsf{C}_{k}(K)|\;\partial_{k+1}(c^{\prime})=c,\,c^{\prime}\in\mathsf{C}_{k+1}(K)\\}.$

It turns out that all $k$-boundaries themselves have no $(k-1)$-boundaries.
This is due to a fundamental property of the boundary homomorphisms in
topology, which states that for any $k\geq 1$ 53,

$\partial_{k-1}\circ\,\partial_{k}=\emptyset.$ (2)

From our previous example in Fig 2a, the 2-chain of the six squares that
enclose the solid cube is the boundary of a 3-chain (a lone 3-cell in this
case). They themselves do not have any opening, and hence the boundary of this
2-chain is $\emptyset$. A consequence of this general property is that we have
$\mathsf{B}_{k}(K)\subseteq\mathsf{Z}_{k}(K)\subseteq\mathsf{C}_{k}(K)$.

Figure 3: An illustration of homology classes of $k$-cycles in a
$(k+1)$-manifold $\mathbb{M}_{k+1}$. a). The two $k$-cycles in brown and dark
brown enclose a $(k+1)$-body ($k$-chains in light brown). They contain no
holes and are boundary cycles. They belong to the same equivalence class of
$[\emptyset]$. The $k$-cycles in green/light green, purple/light purple, and
blue/light blue each contain $k$-holes $h_{\alpha},\,h_{\beta}$ and
$h_{\gamma}$, respectively, and are part of equivalence classes of
$[h_{\alpha}],\,[h_{\beta}]$ and $[h_{\gamma}]$, respectively. b). The
$k$-cycle $l_{1}$ contains a hole. The boundary cycle $l_{2}$ contains no hole
but a $(k+1)$ body. Note that $l_{1}$ and $l_{2}$ share a common piece of
boundary but in opposite orientations. c). When $l_{1}$ and $l_{2}$ are summed
up, we obtain the $k$-cycle $l_{3}$, which contains the same hole as $l_{1}$.
Both $l_{1}$ in b) and $l_{3}$ in c) belong to the same equivalence class of
$k$-cycles containing this hole.

#### Homology group and Betti number.

There are two types of $k$-cycles: Those enclose $(k+1)$-bodies and those
enclose $(k+1)$-holes (Fig 3a). The former are boundaries of the enclosed
bodies of $(k+1)$-chains and can be collapsed into a point (Fig 3a, $k$-cycles
in dark brown/brown enclosing $(k+1)$-chains shaded in light brown ). The
latter are not boundaries of $(k+1)$ bodies and cannot be collapsed into a
point (Fig 3a, other cycles) . We will first distinguish these two types of
cycles. Furthermore, among cycles that do not enclose a body, they may be so
due to different reasons, as the holes they contain may be different (Fig 3a).
We will distinguish these different situations as well.

We consider all cycles containing the same hole essentially the same and group
them into one equivalence class of cycles. As an illustration, green/light
green $k$-cycles in Fig 3a form an equivalence class $[h_{\alpha}]$ as they
all contain hole $h_{\alpha}$. So do the purple/light purple $k$-cycles (class
$[h_{\beta}]$ containing hole $h_{\beta}$), and the blue/light blue $k$-cycles
(class $[h_{\gamma}]$ containing hole $h_{\gamma}$). A special equivalence
class are cycles containing the $\emptyset$ hole or no hole (Fig 3a,
brown/dark brown cycles, class $[\emptyset]$). We call cycles in each
equivalence class homologous to each other. If they encircle different holes,
they belong to different equivalence classes.

We elaborate on this. Among all elements of $\mathsf{Z}_{k}(K)$, which are
$k$-cycles, we identify all $k$-boundaries, which contain $(k+1)$ bodies and
are elements of $\mathsf{B}_{k}(K)$. Because of Eqn (2), they have no
$(k+1)$-holes. We put them into a class denoted as $[\emptyset]$ as they
contain no holes (or $\emptyset$-hole) (Fig 3a, brown/dark brown cycles
$\in[\emptyset]$). For the remaining $k$-cycles of $\mathsf{Z}_{k}(K)$, they
are not in $[\emptyset]$ but may contain different holes. We identify those
contain hole $h_{a}$, and put them into the equivalence class denoted as
$[h_{a}]$ (Fig 3a, green/light green cycles $\in[h_{a}]$). Remaining
$k$-cycles that contain hole $h_{b}$ are put into the class $[h_{b}]$, and so
on (Fig 3a, purple/light purple cycles $\in[h_{b}]$, and blue/light blue
cycles $\in[h_{c}]$). Each element of the set
$\\{[\emptyset],\,[h_{a}],\,[h_{b}],\,\cdots\\}$ is an equivalence class.

As these equivalence classes themselves form a set, and the outcome of the
binary operation of summation on elements of $\mathsf{Z}_{k}$ is preserved,
this set form a new group. This new group is called a quotient group, as it is
obtained from $\mathsf{Z}_{k}(K)$ after factoring out the boundaries
$\mathsf{B}_{k}(K)$. The  $k$-th homology group $\mathsf{H}_{k}(K)$ is this
quotient group:

$\mathsf{H}_{k}(K)\equiv\mathsf{Z}_{k}(K)/\mathsf{B}_{k}(K).$

The elements of $\mathsf{H}_{k}(K)$ are equivalence classes of homologous
cycles representing the holes (or lack of) they encloses.

Two $k$-cycles are homologous to each other if they contain the same hole, or
equivalently, if one can be obtained from another by adding a $k$-boundary
(Fig 3b and Fig 3c). To illustrate this, note that the cycle labeled $l_{3}$
in Fig 3c can be obtained by adding the $k$-boundary of a $(k+1)$-body labeled
$l_{2}$ to the $k$-cycle labeled $l_{1}$ enclosing a hole (Fig3b). Due the
cancellation nature of our summation, the commonly shared piece of the
boundaries is cancelled out and we have the larger $k$-cycle $l_{3}$ in Fig 3c
enclosong the same hole. Upto the difference of a boundary of a solid
$(k+1)$-body, these two $k$-cycles are the same and belong to the same
equivalence class. It is not difficult to see that we can repeat this
operations of adding certain $k$-boundaries and convert any homologous
$k$-cyles between one another.

The number of the equivalent classes, or the number of independent
$k$-dimensional holes, is counted by the dimension of the homology group. It
is called the $k$-th Betti number $\beta_{k}(K)$:

$\beta_{k}(K)=\dim(\mathsf{H}_{k}(K))$

Figure 4: The probability landscape $f(\mathbb{M})$ and the topology of its
superlevel set $\mathbb{M}_{f\geq}$. a) The landscape and a sea level. The
superlevel sets $\mathbb{M}_{f\geq}$ are the regions of the domain
$\mathbb{M}\subset\mathbb{{}^{d}}$ (shown as a plane) whose landscape value is
above the sea. b) At $f({\mbox{\boldmath$x$}})=1$, all is below the sea level
and $\mathbb{M}_{f\geq 1}=\emptyset$. At
$f({\mbox{\boldmath$x$}})=b_{1},\,b_{2}$ and $d_{2}$, the topology of
$\mathbb{M}_{f\geq}$ (shown in white) changes. At
$f({\mbox{\boldmath$x$}})=0$, all is above sea level and we have
$\mathbb{M}_{f\geq}=\mathbb{M}$. c) The persistent diagram of the birth and
death value of $f({\mbox{\boldmath$x$}})$ for the $0$-th homology group
representing the two peaks. The sublevel sets below the sea $\mathbb{M}_{f<}$
are shown in blue.

#### Filtration.

We now examine the topological structures of holes in the probability
landscape on the configuration space, when we restrict to configurations all
with probabilities above certain value. By gradually adjusting this value, we
will be able to trace out the details of topological changes. For an
illustration, envision a sea on top of the probability landscape (Fig 4). At
the level of $f({\mbox{\boldmath$x$}})=1$, it covers the whole landscape. The
domain of the part of the landscape above the sea level is $\emptyset$. We
gradually lower the sea level to value $b_{1}$, when the first peak emerges
from the sea (birth of the first peak). At this time, we have the superlevel
set $\mathbb{M}_{f\geq b_{1}}$, which are the set of points
$\\{{\mbox{\boldmath$x$}}\in\mathbb{M}|f({\mbox{\boldmath$x$}})\geq b_{1}\\}$.
They form the white region(s) in Fig 4. We further lower the sea level to
$b_{2}$ when another peak emerges (birth of the second peak), at which time we
have the superlevel set $\mathbb{M}_{f\geq b_{2}}$. Suppose we continue this
process until sea level reaches $d_{2}$ where the two peaks are merged
together (death place of the second peak) by a land ridge that has just
emerged above the sea level. At this sea level, we have $\mathbb{M}_{f\geq
d_{2}}$. At each of these levels, the topology of the superlevel set changes,
namely, one component, two components, and then one component again. These
changes are captured by the changing homology groups and the Betti numbers.

We now generalize. We have a descending sequence of probability values
corresponding to the lowering sea level:

$1=a_{0}>a_{1}>a_{2}>\cdots>a_{n}=0,$

and the corresponding superlevel sets, or the domains of the part of the
landscape above the sea level, which are subspaces of $\mathbb{M}$:

$\emptyset=\mathbb{M}_{0}\subset\mathbb{M}_{1}\subset\mathbb{M}_{2}\cdots\subset\mathbb{M}_{n}=\mathbb{M}.$

Recall we have the full configuration space $\mathbb{M}$ represented by a
cubic complex $K$. Each superlevel set $\mathbb{M}_{i}$ is represented by a
subcomplex $K_{i}\subset K$, which can be derived from the original full
complex $K$. We then have the corresponding sequence of subcomplexes:

$\emptyset=K_{0}\subset K_{1}\subset K_{2}\cdots\subset K_{n}=K.$

This sequence of subcomlexes is called a filtration.

We are interested in how the topology of $\mathbb{M}_{f\geq a_{i}}$ evolves at
different $a_{i},i=0,\cdots,n$. This is represented by the corresponding
sequence of homology groups connected by linear maps:

$0=\mathsf{H}_{k}(K_{0})\rightarrow\mathsf{H}_{k}(K_{1})\rightarrow\cdots\rightarrow\mathsf{H}_{k}(K_{n})=\mathsf{H}_{k}(K).$

#### Persistence and Persistent diagram.

As we move from $K_{i-1}$ to $K_{i}$, we may gain a new equivalence class,
(e.g., a new peak for $0$-th homology as in our example), or we may lose one
(e.g. when a peak is merged with another one). We say that an equivalence
class of a $k$-cycle $[\alpha_{i}]$ is born at $a_{i}$ if its equivalence
class is present in $K_{i}$ but absent in $K_{i-1}$ for any value of
$a_{i-1}<a_{i}$. The class dies at $a_{i}$ if it is present in $K_{i-1}$ for
any value of $a_{i-1}<a_{i}$ but not at $a_{i}$. We record the location and
the value of $a_{i}$, namely, the corresponding $k$-cube and its probability
value whose inclusion lead to the birth and death events.

The prominence of the topological feature of a $k$-cycle is encoded in its
life-time or persistence. Denote the birth value and the death value of class
$[\alpha_{i}]$ as $b_{i}$ and $d_{i}$, respectively. The persistence of class
$[\alpha_{i}]$ is then $b_{i}-d_{i}$.

In the example shown in Fig 4, the equivalence class of $0$-cycles
(components) associated with the first peak is born at
$f({\mbox{\boldmath$x$}})=b_{1}$. The equivalent class associated with the
second peak is born at $f({\mbox{\boldmath$x$}})=b_{2}$. At
$f({\mbox{\boldmath$x$}})=d_{2}$, these two components merge together. We say
that the second peak dies at $d_{2}$, and its persistence is $b_{2}-d_{2}$.
The first peak dies at $f({\mbox{\boldmath$x$}})=0$, and its persistence is
$b_{1}-0=b_{1}$.

We record the birth and death events of homology classes in a two-dimensional
plot, which is called the persistent diagram. Each homology class is
represented by a point in this diagram, where the birth value $b_{i}$ and the
death value $d_{i}$ are taken as its coordinates ($b_{i},\,d_{i}$). Fig 4c
shows the persistent diagram of our illustrative example.

In general, we have the $k$-th persistent diagram $\mathcal{P}_{k}(f)$ of
$k$-cycles for our probability function $f:\mathbb{M}\rightarrow{}_{[0,1]}$.
It is the set of points such that each point $(x_{1},x_{2})$ represents a
distinct topological feature of $k$-cycle, which is present in
$\mathsf{H}_{k}(\mathbb{M}_{f\geq a})=\mathsf{H}_{k}(f^{-1}([a,1))$ for
$a\in[x_{1},x_{2})$.

#### Computation.

The key to study homology groups in high-dimensional space is the construction
of the $K$-complex to represent the configuration space $\mathbb{M}$. In this
study, we use the cubic algorithm of 56 with modifications 58, so it can be
applied to higher dimensions. We consider only the $0$-th persistent homology
groups, which records the birth and death of probability peaks. The locations
$x$ where birth and death events occur, namely, the corresponding $k$-cubes
are also computed.

##  3 Results

###  3.1 Model system and computation

#### Ensemble of reactive trajectories and conformations.

The isomerization of alanine dipeptide in vacuum provides a tractable system
for understanding the process of activation in details, and has been well
studied as a model for understanding protein conformational changes 21, 29,
59.

Figure 5: The isomerization reaction of alanine dipeptide. (a) Conformations
from the reactant and product basins before and after the isomerization. (b)
The six reaction coordinates of the isomerization process of alanine peptide
examined in this study.

Using transition path sampling 8, we harvested 6 million reactive
trajectories. Each trajectory is of 2.5 ps duration. We further collect
conformations every 50 steps at 1 fs/step along each trajectory. Altogether,
we have a total of $1.5\times 10^{11}$ conformations. All simulations are
conducted using the molecular dynamics software suite GROMACS4.5.4 60.
Amber$94$ force field was used to facilitate the comparison with previous
results. The simulation was performed with constant total energy $36$ KJ/mol,
such that the averaged temperature is $300$K for the transition path ensemble.
Note that the transition portion of each reactive trajectory is around $0.2$
ps, thus the majority of our $2.5$ ps trajectories are within the two stable
basins. Here the reactant basin is defined in radian as
$(\phi,\,\psi)\in[(-3.49,\,-0.96)\times(-1.57,\,3.32)]$ and the product basin
is defined as $(\phi,\,\psi)\in[(0.87,\,1.74)\times(-1.39,\,0)]$.

#### Constructing dynamic probability surface of transition state.

We then construct the dynamic probability surface of the isomeriztion reaction
from the sampled $1.5\times 10^{11}$ conformations. With a balanced
consideration of the available MD simulation trajectories and the
dimensionality, we construct a 5-dimensional configuration space for this
study. Based on previous analysis using the energy flow theory 17, we selected
the top 5 coordinates ($\phi,\theta_{1},\psi,\alpha,\beta$) that contribute
most to the activation dynamics. The original 60 dimensional space is then
projected onto this 5 dimensional space, where each dimension is divided into
15 bins. This leads to $15^{5}=759,375$ 5-dimensional hypercubes.

#### Computing topological structures of the dynamic probability landscape.

We then carry out persistent homology analysis. Computations are conducted on
a machine with a 20-core Xeon E5-2670CPU of 2.5 GHz, with a cache size of 20
MB and memory of 128 GB Ram. The computing time for finding the significant
peaks and ridges connecting them is $\approx$ 30 seconds.

#### Committor test for conformations at selected locations of the
configuration space.

We carry out committor test for configurations of the dipeptide identified by
persistent homology. The committor value of a configuration is defined as the
probability that a dynamic trajectory initialed from this configuration, with
initial momenta drawn from the Boltzmann distribution, reaches the product
basin before the reactant basin. A configuration with committor value
$p_{B}=0.5$ is regarded as a member of the transition state ensemble.

In a committor test 8, 9, we need to generate an ensemble of tentative
transition state conformations from locations in configurtaion space where
probability peaks and ridges identified by persistent homology are located.
These conformations all share the same target values for the selected
coordinates (e.g., $\phi$, $\theta_{1}$) that correspond to the location of
the selected peak, with the other coordinates sampling the equilibrium
distribution. For this, we add harmonic restraint potentials on the selected
coordinates to the system potential energy function. The minima of the
harmonic restraints are at the target values. Equilibrium MD simulations are
then carried out. Conformations harvested from such simulations are filtered
to generate an ensemble of conformations that all share the same target values
for the selected coordinates. The restraint potential is used to enrich
conformations that satisfy this criterion.

###  3.2 Topology and dynamic properties of 5-d dynamic probability surfaces

#### Dynamic probability surface on configuration space of
($\phi,\,\theta_{1},\,\psi,\alpha,\beta$).

We examine the topological structures of this 5-d dynamic probability surface.
There are four peaks, each located in a 5-d cube (Fig. 6a, red dots for birth
locations of peaks and blue dots for their ridges or death locations, and
values are listed in SI Table S1). The most prominent peak with the largest
persistence is peak $b_{1}$ (see persistent diagram in Fig. 6b), which
corresponds to the product basin. The second most prominent $b_{2}$ is the
reactant basin. The third prominent peak at $b_{3}$ (Fig 9a) has roughly the
same probability as the reactant basin $b_{2}$, but a shorter persistence,
namely, it does not stand out from the surrounding landscape as much. It
subsequently merges with the peak at the product basin $b_{1}$. As expected,
peaks have all negative eigenvalues for their Hessian matrices, and ridges
have one positive and four negative eigenvalues (see Supplementary Info).

We then take conformations from the four peaks and the three ridges (SI Table
S1), and carry out committor tests. We find that the most prominent peak
$b_{3}$ beyond the reactant and product basins fully captures the transition
state ensemble ($p_{B}$ centered at 0.5, Fig 6c): Trajectories initiated from
configurations at this location have equal probability towards the reactant or
the product basin.

Figure 6: The 5-d dynamic probability surface on the ($\phi-\theta_{1}$)
plane, its topological structures, and the committor values. (a) The 5-d
dynamic probability surface $p(\phi,\psi,\theta_{1},\alpha,\beta)$ shown on
the ($\phi-\theta_{1}$) plane. Red and blue dots are locations of probability
peaks and ridges (see also SI Table S1). (b) The persistent diagram recording
the birth and death probabilities $p(b_{i})$ and $p(d_{i})$ of the peaks in
$y$ and $x$, respectively. (c-e) Distributions of committor values $p_{B}$ for
trajectories from locations of $b_{3}$, $d_{2}$, and $d_{3}$, respectively.
(c) The transition state ensemble is located at $b_{3}$.

In contrast, all committor values $p_{B}$ for locations of $b_{2}$, $d_{2}$,
and $d_{4}$ are found to be 0. Reaction trajectories starting from
conformations at these locations fall back to the reactant basin. The
committor values for peak $b_{1}$ are all 1.0: Trajectories from this location
all go to the product basin. The committor values for conformations at the
ridges $d_{2}$ and $d_{3}$ follow one-sided distributions (Fig 6d-6e,
respectively). Only a negligible amount of conformations have $p_{B}=0.5$.

These results demonstrate that the dynamic properties of the transition state
ensemble is capture by the topological features of this dynamic probability
surface. Furthermore, the transition state ensemble is located at the most
prominent peak outside of the reactant and product basins, instead of a saddle
point.

#### Dynamic probability surface on configuration space
($\phi,\,\psi,\alpha,\beta,\,\theta_{2}$).

To examine the importance of proper choice of the coordinates, we construct
another 5-d probability surface by omitting the coordinate $\theta_{1}$ and
replacing it with $\theta_{2}$, which is on the other end of the molecule in a
position symmetric to $\theta_{1}$. Fig. 7a shows the projection of the 5-d
surface in $-\ln p({\mbox{\boldmath$x$}})$ on the ($\phi-\theta_{2}$) plane.

There are four significant peaks (Fig. 7a, red and blue dots), which are also
shown in the persistent diagram (Figure 7b), each located in a 5-d cube (see
SI Table S1). Peaks $b_{1}$ and $b_{2}$ are the most and second most prominent
peaks, corresponding to the product and the reactant basins, respectively. In
this projection, peaks are separated only in $\phi$ and they have almost
identical $\theta_{2}$ values. This differs from the 5-d surface containing
$\theta_{1}$ (Fig 6).

We then sample conformations from the location of each of the four peaks and
the three ridges (SI Table S1) and perform committor tests. The committor
values for peak $b_{1}$ are all 1.0, where trajectories starting here all go
to the product basin. All committor values for conformations at $d_{2}$,
$b_{4}$, and $d_{4}$ are 0.0. Reaction trajectories from these locations all
go to the reactant basin. The committor values for $d_{2}$ and $d_{3}$ follow
one-sided distributions (Fig. 7d-7e). Only a tiny amount of the conformations
have $p_{B}=0.5$, indicating that the transition state ensemble are located
elsewhere.

The committor values for $b_{3}$ has a flat distribution. While there are
conformations with $p_{B}$ values close to 0.5, their frequency is similar to
any other $p_{B}$ values. This indicates that the 5-d cube where $b_{3}$ is
located contains some transition state conformations as its $\phi$ value is
correct, but also a mixture of other conformations with diverse dynamic
properties. Overall, without the reaction coordinate $\theta_{1}$, this 5-d
dynamic probability surface does not describe the activated process
adequately, and cannot be used to identify the transition state ensemble.

Figure 7: A different 5-d dynamic probability surface with $\theta_{2}$
replacing $\theta_{1}$ on the ($\phi-\theta_{2}$) plane, its topological
structure represented in persistent diagram, and distributions of committor
values. (a) The 5-d dynamic probability surface
$p(\phi,\psi,\alpha,\beta,\theta_{2})$ projected onto the ($\phi-\theta_{2}$)
plane. Red and blue dots are locations of probability peaks and ridges. (b)
The persistent diagram recording the birth and death probabilities $p(b_{i})$
and $p(d_{i})$ of the peaks in $y$ and $x$, respectively. (c-e) Distribution
of the committor values $p_{B}$ for trajectories from locations of $b_{3}$,
$d_{2}$, and $d_{3}$.

###  3.3 Topology and dynamic properties of projected 2-d dynamic probability
surfaces.

As it is difficult to directly study the topology of a high-dimensional
surface, a common practice is to project the surface to a lower dimensional
subspace and analyze the topology of the projected surface instead. The caveat
of this practice is that the original topological features may be lost, new
features that are artifacts may arise due to the marginalization of the
probability distributions (See Fig. 9 for details).

To assess how well the dynamic properties is retained in a subspace, we
project the 5-d surface on ($\phi$, $\theta_{1}$, $\psi$, $\alpha$, $\beta$)
to 2-d planes. We then analyze the topological structures of the projected
surfaces, and assess the dynamic behavior of the identified topological
features. We carried out this analysis using the 2-d planes of
($\phi-\theta_{1}$) and ($\phi-\psi$). These choices are based on prior
knowledge of the role of $\phi$ and $\theta_{1}$ in alanine-dipeptide, and the
widely practiced ($\phi-\psi$) based Ramachandran plot.

#### Projecting 5-d dynamic probability surface to the ($\phi-\theta_{1}$)
plane.

After projection, there are only three significant probability peaks (Fig. 8a,
Fig 9a, red and blue dots, and SI Table 2), instead of four peaks found in the
full $5$-d probability surface of Fig. 7. The most and the next prominent
peaks $b_{1}$ and $b_{2}$ shown in the persistent diagram of Fig. 8b are the
product and reactant basins, respectively.

It is informative to compare peak locations on the 5-d surface to that on the
projected surface (SI Table S1 and SI Table S2). The $\phi$ coordinate for the
product basin $b_{1}$ is altered from 1.25 to 0.84 after projection, while the
$\theta_{1}$ coordinate is unchanged. The $(\phi,\,\theta_{1})$ coordinates of
the reactant basin $b_{2}$ are also changed from $(-1.68,\,-0.18)$ to
$(-1.25,\,0.0)$. $\theta_{1}$ of the ridge $d_{3}$ is changed from $-0.18$ to
$0.0$, while $\phi$ is unchanged. The fourth prominent peak becomes
undetectable after projection. The persistence diagram (Fig. 8b) is also
significantly different: The third peak is much less prominent with reduced
persistence compared to that in Fig. 6b (see also SI Table S2).

Figure 8: The ($\phi$-$\theta_{1}$)-projection of the 5-d dynamic probability
surface, its topological structure, and distributions of committor values. (a)
The 5-d dynamic probability surface projected onto the ($\phi-\theta_{1}$)
plane. Red and blue dots are locations of probability peaks and ridges,
respectively. (b) The persistent diagram recording the birth and death
probabilities $p(b_{i})$ and $p(d_{i})$ of the peaks in $y$ and$x$,
respectively. (c-e) Distributions of committor values $p_{B}$ for trajectories
from peaks and ridges $b_{3}$, $d_{2}$, and $d_{3}$, respectively. (c)
Transition state conformations are at $b_{3}$.

We then carry out committor test on conformations from locations of these
topological features. The distribution of $p_{B}$ sampled from trajectories
from $b_{3}$ is centered around 0.5 (Fig 8c) and exhibits significant
enrichment of transition state conformations. This is similar to peak $b_{3}$
on the original surface over the 5-d configuration space (Fig 6c), although
the distribution has a broader width.

The committor values of trajectories from the product basin $b_{1}$ are all
1.0, as they all fall back to the product basin. Similarly, trajectories from
$b_{2}$ all fall back to the reactant basin with a $p_{B}$ value of 0.0. The
committor values for the bridges at $d_{2}$ and $d_{3}$ follow one sided
distributions, but only a small amount of conformations have $p_{B}=0.5$ (Fig.
8)d and e.

Overall, these results demonstrate that when projected to the 2d-plane of
($\phi$-$\theta_{1}$), which is formed by the two dominant reaction
coordinates, the dynamic probability surface retain essential dynamic
properties of the transition state surface, and contain rich information such
that the transition state conformations can be recovered.

#### Projecting 5-d dynamic probability surface to the ($\phi-\psi$) plane.

$\phi$ and $\psi$ angles are the standard parameters to describe protein
secondary structures. After projection to the ($\phi-\psi$) plane, there are
three significant probability peaks (Fig. 10a, red and blue dots, and SI Table
S2). The most and the next most prominent peaks $b_{1}$ and $b_{2}$ as shown
in the persistent diagram (Fig. 10b) are the product and reactant basins,
respectively. Similar to projection to the ($\phi$-$\theta_{1}$) plane,
locations of both basins are altered upon projection (SI Table S2). Peak 3
becomes very minor after projection to ($\phi-\psi$) plane, and its location
is also altered from that of the 5-d surface. This can be seen in Fig. 9b,
where the location of peak 3 on the 5-d surface projected onto ($\phi-\psi$)
plane (red point), and the changed location of the new peak after projecting
onto the ($\phi-\psi$) plane (blue dot) are shown.

Figure 9: The 5-d probability surface projected onto two different 2-d planes.
(a) The probability surface projected onto the ($\phi-\theta_{1}$) plane. The
red dot is the location of the third probability peak after reactant and
product peaks, which is where the transition state conformations are located.
(b) The probability surface projected onto the ($\phi-\psi$) plane.
Conformations with the correct transition state value of $\phi$ are at a non-
peak location, which is on a slope below the new peak shown in blue. (c) The
process of projecting the 5-d probability surface onto the 2-d ($\phi-\psi$)
plane for the blue and red dots shown in (b). While the probability at the
location of the red dot is larger on the 5-d probability surface (left bars),
after projection onto the ($\phi-\psi$) plane, the probability at the red dot
is smaller than that at the blue dot (right bar plots). As a consequence, the
projected probability surface on the ($\phi-\psi$) plane does not capture the
actual location of the peak in the 5-d surface.

These observations demonstrate that with projection, locations of topological
features of probability peaks may change, and their prominence as measured by
persistence may also change dramatically. Fig. 9c explains why the projection
to ($\phi$-$\psi$) results in the peak on the 5-d surface (red dot in Fig. 9a)
changing location (Fig. 9b, peak location shown as a blue dot, red dot no
longer at the peak). First, as shown in Fig. 9c (beginning bar plot), the
probability at the location of the red dot is larger than that at the blue dot
on the 5-d surface. The probability $p(\phi,\psi)$ at each point on the
($\phi$-$\psi$) plane (Fig. 9b) is the sum of all points on the 5-d surface
with the same $\phi$ and $\psi$ but with different values in any of the other
three coordinates ($\theta_{1},\alpha$ and $\beta$). Thus, the probability of
each point on the ($\phi-\psi$) plane of Fig. 9b is a 3-d hyper-surface of
$p(\theta_{1},\alpha,\beta)$. When we sum up the 3-d hyper-surface along one
direction (e.g., $\alpha$), we obtain a 2-d probability surface (e.g.
$p(\theta_{1},\beta))$. The 2-d surfaces for the red and blue dots are shown
in the middle panel of Fig. 9c. Reducing the dimension further, we sum up the
2-d probability surface over $\beta$. The resulting distributions are 1-d
probability distributions along the $\theta_{1}$ direction, shown in Fig. 9c
for the blue and red dots. From these 1-d distributions, we can see that when
details in $\theta_{1}$ direction are retained, the probability at the red dot
(actual peak in 5-d) is still higher than the probability at the blue dot. We
next sum up the probability along the $\theta_{1}$ direction (Fig. 9c, final
bar plots). The summed values are the probability values over the 2-d
($\phi-\psi$) plane shown in Fig. 9b. As shown by the final bar plot in Fig.
9c, the total probability mass for the red dot is now less than that for the
blue dot, even though the red dot has a higher probability on the original 5-d
surface. Hence, the location of peak 3 changes when projecting onto
($\phi-\psi$) plane.

We then carried out committor tests (Fig. 10c). None of the topological
features are where the transition state ensemble are located: All trajectories
starting from $b_{2}$, $b_{3}$, and $d_{3}$ go to the reactant basin
($p_{B}=0)$, and all trajectories starting at $b_{1}$ go to the product basin
($p_{B}=1.0$). The committor values for ridge $b_{2}$ follow a one-sided
distribution around $p_{B}=0$ (Fig. 10c). Trajectories starting there mostly
fall back to the reactant basin. Overall, none of the topological features
after projection to the ($\phi$-$\psi$) plane retain the dynamic properties of
the transiton state conformations as the original 5-d surface.

Figure 10: The ($\phi$-$\psi$)-projection of the 5-d dynamic probability
surface, its topological structure, and distribution of committor values. (a)
The 5-d dynamic probability surface projected onto the ($\phi$-$\psi$) plane.
Red and blue dots are locations of probability peaks and ridges, respectively.
(b) The persistent diagram recording the birth and death probabilities
$p(b_{i})$ and $p(d_{i})$ of the peaks in $y$ and $x$, respectively. (c)
Distribution of committor values $p_{B}$ for peak $b_{2}$.

###  3.4 Dimension reduction by PCA destroys dynamic properties inherent in
surface topology

Principal Component Analysis (PCA) is a widely used technique for dimension
reduction. It has found broad applications in molecular simulations 61, 62,
63, 64, 65. However, whether such reduction retains the essential dynamics of
the activation process and whether the surface topology on the PCA space can
uncover the transition state conformations are not known. Here we assess the
dynamic properties of probability surfaces after PCA dimension reduction.

#### Projection of $p(\phi,\psi,\theta_{1},\alpha,\beta,)$ onto
$(PC_{1},\,PC_{2})$ by dPCA.

We first applied dihedral principal component analysis (dPCA) to the 5-d
probability surface 64, 65. dPCA is widely used for dimension reduction on
periodic dimensions. It first maps each periodic dimension of circular angle
to two new dimensions using the $\sin$ and $\cos$ functions. Regular PCA is
then applied for dimension reduction. We use the dPCA procedure and obtain the
first two principal components from the variance matrix of the $1.5\times
10^{11}$ conformations. Collectively, they account for $80.4\%$ of the
variance.

Figure 11: The dynamic probability surface on PCA space, its topological
structure, and committor values on principal components space. (a) The
projection of the 5-d dynamic probability surface
$p(\phi,\,\psi,\,\theta_{1},\,\alpha,\,\beta)$ to the plane of
($PC_{1}$-$PC_{2}$). (b) The persistent diagram exhibits six probability
peaks. (c-e) Distribution of committor values $p_{B}$ for trajectories from
bridge $d_{2}$, peak $b_{6}$, and bridge $d_{6}$, respectively.

The dynamic probability surface after projection to the ($PC1$-$PC2$) plane is
dramatically more complex (Fig. 11a and SI Table S3) than the surfaces when
projected to either the ($\phi$-$\theta_{1}$) plane (Fig. 8) or the
($\phi$-$\psi$) plane (Fig. 10). The persistent diagram (Fig. 11b) is also
very different from that of the original 5-d probability surface (Fig. 6b).

The committor tests show that all committor values for conformations at
$b_{1}$ are 1.0 with trajectories going to the product basin. All committor
values for conformations at $b_{2-5}$ and $d_{3-5}$ are 0.0: trajectories from
there all go to the reactant basin. Committor values at ridge $d_{2}$ and peak
$b_{6}$ follow a one-sided distribution at $p_{B}=0$ (Fig. 11c). Committor
values at ridge $d_{6}$ has higher values both at $p_{B}=0.0$ and $p_{B}=1.0$.
Conformations at this location are a mixture of those close to the reactant
basin and those close to the product basins. There are few conformations from
the transition state ensemble.

Overall, our results show that the dPCA procedure for dimension reduction
removes dynamics relevant information from the topological features of the
probability surface. No conformations in the transition state of this active
process are captured by the topological features of the PCA surface.

#### Projection from 39-d angular space to 5-d dPCA subspace.

We also applied dPCA to the original full-dimensional dynamic probability
surface. After removal of the 21 bond lengths, we apply dPCA to reduce the
remaining 39-dimensional configuration space of angles to 5 principal
components. The contour plot of the 5-d PCA surface projected onto the first
two principal components is shown in Fig. 12a. Persistent homology analysis
identifies 16 peaks (Fig. 12b), each located in a 5-d cube. The location and
probability of the first 6 most dominant peaks and ridges connecting them are
listed in SI Table S3. This 5-d persistent diagram is very different from that
shown in Fig. 6b, exhibiting a significantly more complex surface topology.

Figure 12: The dynamic probability surface on 5-d principal components space
reduced by dPCA from the 39-d configuration space, its topological structure,
and committor values. (a) the 5-d dynamic probability surface
$p(PC_{1},\cdots,PC_{5})$ shown on to the ($PC_{1},PC_{2}$) plane. (b) The
persistent diagram exhibits 16 peaks. (c-d) Distribution of committor $p_{B}$
values for trajectories started at ridges $d_{9}$ and $d_{10}$, respectively.

The committor tests show that all committor values for conformations located
at the peaks and the ridges are either 0 or 1.0. The two exceptions are ridges
$d_{9}$ and $d_{10}$ (Fig. 12c-d). None of the topological features in the
dPCA reduced 5-d dynamic probability surface retain the essential dynamics of
the transition state ensembles. Overall, our results demonstrate that
dimension reduction by dPCA destroy dynamic properties inherent in the
original surface topology of the dynamic probability surface.

#### Direct PCA projections.

Results using direct PCA to project
$p(\phi,\,\psi,\,\theta_{1},\,\alpha,\,\beta)$ onto ($PC1,\,PC2$), as well as
projection from the 39-d space to the 5-d direct PCA subspace are similar and
none retain the dynamic properties in their topological features (see SI for
details).

##  4 Discussion

In this study, we have introduced a novel approach for characterizing the
exact topological features of dynamic probability surfaces. Instead of
examining critical points and Morse indexes, ours is based on homology groups
of a series of superlevel sets of the probability surface. With quantification
of the scales of these topological features by persistent homology, we are
able to uncover the relationship between the topology of the dynamic
probability surface and the dynamics of the activation process of the alanine-
dipeptide isomerization reaction.

This approach allows us to define the topological properties of the high-
dimensional dynamic probability surface that is associated with the transition
state conformations. The probability surface over the transition state region
is the most prominent peak after the reactant and product basins. Instead of a
Morse index of $1$ as conventionally thought 5, 66, transition state ensemble
is on the top of a dynamic probability peak and goes downhill in all
directions. As seen when projected to the ($\phi,\,\theta_{1}$)-plane (Fig
9a), it appears as a small peak rather than a saddle point commonly associated
with transition state on multi-dimensional free energy surface. This is
because the system undergoes certain amount of correlated wandering motions at
the barrier top, before it goes down towards the product basin. Our finding is
against the conventional wisdom that the $C_{7eq}\rightarrow C_{7ax}$
transition is a ballistic process, as it is a small peptide and the transition
occurs in vacuum.

The dynamic probability surface was constructed from naturally occurring
reactive trajectories connecting the reactant and product basins. These
trajectories are unbiased and faithfully reflect how the $C_{7eq}\rightarrow
C_{7ax}$ transition occurs. They contain all the relevant information about
the dynamic process of the activation. Unlike the free energy surface commonly
used in examining the mechanism of an activated process, this probability
surface contains additional information that reflect the non-equilibrium
nature of the transition dynamics.

A common practice in the studies of protein conformational dynamics is to
extract mechanistic insights from the geometry of two-dimensional free energy
surface of a double-well along certain collective variables, which are often
chosen based on heuristics or by intuition. For example, one would associate a
saddle region with the transition states. Our results illustrate the caveats
of such procedures. All three probability surfaces shown in Figs. 6, 8, and 10
exhibit the canonical double-well feature. In Fig. 8, the 5-d surface on the
($\phi,\,\theta_{2}$)-plane has both the product and reactant peaks and the
third peak laying out along the $\phi$-direction alone, indicating that
$\theta_{2}$ is not a reaction coordinate. This is indeed verified by the
committor test: configurations corresponding to the peak in the saddle region
all share the correct value of $\phi$ but samples randomly along the
$\theta_{1}$ direction, as illustrated by the flat committor distribution in
Fig. 7c. Detialed examination shows that the 5-d cube containing the
transition state conformations (red dot in Fig. 13) also contain conformations
with other $\theta_{1}$ values, which are not at the transition state: some go
to the product basin (e.g., green dot, Fig. 13b), and others to the reactant
basin (e.g., blue dot, Fig. 13b).

Figure 13: The horizontal band of the probability surface of Fig. 7a centered
at $\theta_{2}=0$ ($-0.08<\theta_{2}<0.08$) expanded to show distribution in
$\theta_{1}$. (a) This band containing all peaks in Fig. 7a is expanded in the
$\theta_{1}$ direction. The ($\phi$-$\theta_{2}$)-square containing peak
$b_{3}$ in Fig. 7a is expanded in $\theta_{1}$ and is shown as a vertical
strip between the two dashed lines. This strip contains a mixture of
conformations. The red dot shows the location of the transition state
conformations. (b) The distributions of committor values for conformations at
the blue, red, and green dots in (a), respectively. Trajectories from the blue
and green dots fall back to the reactant and product basins, respectively.

In contrast, projection of the 5-d surface on ($\phi,\,\theta_{1}$)-plane (Fig
6) has the two basins along the $\phi$-direction alone, but the transition
region aligned along both $\phi$ and $\theta_{1}$. This second feature is
consistent with the importance of $\theta_{1}$ in determining the barrier
crossing dynamics 1, 2, 3, 4, 17, 21, 15.

Interestingly, the ($\phi,\,\psi$)-surface (Fig 10) has basins and the saddle
region arranged along both directions. Conventional wisdom would have led to
the conclusion that $\psi$ is important in defining both reactant and product
basins as well as the barrier crossing process. However, the double-well
structure exhibited on the ($\phi$-$\psi$) plane is profoundly misleading.
None of the topological features on the surface over the ($\phi$-$\psi$) plane
retain the dynamic properties of the original 5-d surface, as the committor
test showed that configurations corresponding to the peak at the transition
region completely fall into the reactant basin. This demonstrates that the
correlation between $\phi$ and $\psi$, due to minor roles of $\psi$ to the
transition process 1, 2, distorted the probability distribution along $\phi$,
such that the ridge/saddle extended into the reactant peak, leading to
incorrect $\phi$ value that marks the location of the peak in the transition
region. In contrast, $\theta_{2}$ did not impact the distribution of $\phi$,
so the peak in the transition region on the ($\phi$-$\theta_{2}$) plane still
bears the correct $\phi$ value for the transition states.

In general, the dynamic properties of topological features of the probability
surface are very sensitive to the subspace of projection. This is illustrated
by the different locations of the transition state conformations, which are at
a peak location on the ($\phi$-$\theta_{1}$) plane (Fig 9a, red dot), but are
at a slope location below a new peak when projected to the ($\phi$-$\psi$)
plane (Fig 9b, blue dot). While the probability at the correct
($\phi$-$\theta_{1}$) square for the transition state conformations (Fig 9c,
left, red) is at a peak, a different location in the ($\phi$-$\psi$) plane
(Fig. 9b, blue dot) has higher probability, as in this projection the
probability mass distributed along the dimension of $\theta_{1}$ is integrated
over all $\theta_{1}$ values (Fig. 9c, higher bar, right), resulting in the
concentration of probability peak at this new location.

Our results show that without the inclusion of the correct reaction
coordinates, which can be identified by the energy flow theory 29, 17,
probability surface of the same dimension no longer correctly characterizes
the dynamic properties of the active process. While the exact topological
features of the surface will be captured, they no longer correspond to the
location of the transition state ensemble. Without $\theta_{1}$, the 5-d
dynamic probability surface over ($\phi$, $\theta_{2}$, $\psi$, $\alpha$,
$\beta$) fail to capture the dynamics of this active process.

Together, our results showed that intuition-based projection (such as
$\phi$-$\psi$) or other arbitrary projection cannot be relied upon for
understanding the dynamic properties of activated processes. without rigorous
examination such as the committor test, directly assigning mechanistic
significance to features of free energy surface is prone to mistakes,
misinterpretations, and misunderstanding.

Finally, our results show that there are dramatic changes in the topological
properties of the probability surface after dimension reduction, when
techniques such as dPCA are applied. While the simple probability surface on
the properly constructed 2-d ($\phi-\theta_{1}$) plane contains rich dynamic
information and is sufficient to uncover the transition state conformations,
the topological features on PCA-reduced surfaces can become more complex and
no longer reflect essential dynamics and cannot be used to identify the
transition state conformations.

The approach of homology group and the technique of analyzing the persistent
homology of the filtration of the superlevel sets of high-dimensional
probability surfaces introduced here are general. It also works in principle
for systems with more than two reaction coordinates, as long as they are part
of the coordinate system, and the cubic complexes can be computed. We envision
they can be applied to investigate topology of high-dimensional probability
surfaces encountered in other physical problems of activated process.

Additional information about the birth and death places and probabilities of
significant peaks for each situation is available in supporting information.
In addition, supporting information provides results from the dimension
reduction by direct PCA without using dihedral PCA.

We thank Drs. Herbert Edelsbrunner and Hubert Wagner for discussion and for
generous help in extending the cubic complex algorithm. This work is supported
by grants NIH R35 GM127084, NIH R01 GM086536, and NSF CHE-1665104.

## Conflict of Interest Statement

There are no conflict of interests.

## References

  * Chandler 1978 Chandler, D. Statistical mechanics of isomerization dynamics in liquids and the transition state approximation. _J. Chem. Phys._ 1978, _68_ , 2959–2970
  * Kramers 1940 Kramers, H. Brownian motion in a field of force and the diffusion model of chemical reactions. _Physica_ 1940, _7_ , 284 – 304
  * Pechukas 1976 Pechukas, P. _Dynamics of molecular collisions: Part B_ ; Springer US: Boston, MA, 1976; pp 269–322
  * Wigner 1938 Wigner, E. The transition state method. _Trans. Faraday Soc._ 1938, _34_ , 29–41
  * Hänggi et al. 1990 Hänggi, P.; Talkner, P.; Borkovec, M. Reaction-rate theory: fifty years after Kramers. _Rev. Mod. Phys._ 1990, _62_ , 251–341
  * Berne et al. 1988 Berne, B. J.; Borkovec, M.; Straub, J. E. Classical and modern methods in reaction rate theory. _J. Phys. Chem._ 1988, _92_ , 3711–3725
  * Pollak and Talkner 2005 Pollak, E.; Talkner, P. Reaction rate theory: what it was, where is it today, and where is it going? _Chaos_ 2005, _15_ , 026116
  * Bolhuis et al. 2002 Bolhuis, P. G.; Chandler, D.; Dellago, C.; Geissler, P. L. Transition path sampling: throwing ropes over rough mountain passes, in the dark. _Annu. Rev. Phys. Chem._ 2002, _53_ , 291–318, PMID: 11972010
  * Du et al. 1998 Du, R.; Pande, V. S.; Grosberg, A. Y.; Tanaka, T.; Shakhnovich, E. S. On the transition coordinate for protein folding. _J. Chem. Phys._ 1998, _108_ , 334–350
  * Li and Ma 2014 Li, W.; Ma, A. Recent developments in methods for identifying reaction coordinates. _Mol. Simul._ 2014, _40_ , 784–793
  * Schramm and Schwartz 2018 Schramm, V. L.; Schwartz, S. D. Promoting vibrations and the function of enzymes. emerging theoretical and experimental convergence. _Biochem._ 2018, _57_ , 3299–3308, PMID: 29608286
  * Schwartz and Schramm 2009 Schwartz, S. D.; Schramm, V. L. Enzymatic transition states and dynamic motion in barrier crossing. _Nat. Chem. Biol._ 2009, _5_ , 551–558
  * Ryter 1987 Ryter, D. On the eigenfunctions of the Fokker-Planck operator and of its adjoint. _Physica A_ 1987, _142_ , 103 – 121
  * Onsager 1938 Onsager, L. Initial recombination of ions. _Phys. Rev._ 1938, _54_ , 554–557
  * Bolhuis et al. 2000 Bolhuis, P. G.; Dellago, C.; Chandler, D. Reaction coordinates of biomolecular isomerization. _Proc. Natl. Acad. Sci._ 2000, _97_ , 5877–5882
  * Li and Ma 2020 Li, H.; Ma, A. Kinetic energy flows in activated dynamics of biomolecules. _J. Chem. Phys._ 2020, _153_ , 094109
  * Li and Ma 2016 Li, W.; Ma, A. Reaction mechanism and reaction coordinates from the viewpoint of energy flow. _J. Chem. Phys._ 2016, _144_ , 114103
  * Best and Hummer 2005 Best, R. B.; Hummer, G. Reaction coordinates and rates from transition paths. _Proc. Natl. Acad. Sci._ 2005, _102_ , 6732–6737
  * Antoniou and Schwartz 2011 Antoniou, D.; Schwartz, S. D. Toward identification of the reaction coordinate directly from the transition state ensemble using the kernel PCA method. _J. Phys. Chem. B_ 2011, _115_ , 2465–2469, PMID: 21332236
  * Hu et al. 2008 Hu, J.; Ma, A.; Dinner, A. R. A two-step nucleotide-flipping mechanism enables kinetic discrimination of DNA lesions by AGT. _Proc. Natl. Acad. Sci._ 2008, _105_ , 4615–4620
  * Ma and Dinner 2005 Ma, A.; Dinner, A. R. Automatic method for identifying reaction coordinates in complex systems. _J. Phys. Chem. B_ 2005, _109_ , 6769–6779, PMID: 16851762
  * Peters and Trout 2006 Peters, B.; Trout, B. L. Obtaining reaction coordinates by likelihood maximization. _J. Chem. Phys._ 2006, _125_ , 054108
  * Antoniou and Schwartz 2009 Antoniou, D.; Schwartz, S. D. The stochastic separatrix and the reaction coordinate for complex systems. _J. Chem. Phys._ 2009, _130_ , 151103
  * Jung et al. 2019 Jung, H.; Covino, R.; Hummer, G. Artificial intelligence assists discovery of reaction coordinates and mechanisms from molecular dynamics simulations. 2019
  * Sidky et al. 2020 Sidky, H.; Chen, W.; Ferguson, A. L. Machine learning for collective variable discovery and enhanced sampling in biomolecular simulation. _Mol. Phys._ 2020, _118_ , e1737742
  * Bonati et al. 2019 Bonati, L.; Zhang, Y.-Y.; Parrinello, M. Neural networks-based variationally enhanced sampling. _Proc. Nat. Acad. Sci._ 2019, _116_ , 17641–17647
  * Wang et al. 2019 Wang, Y.; Ribeiro, J.; Marcelo, L.; Pratyush, T. Past–future information bottleneck for sampling molecular reaction coordinate simultaneously with thermodynamics and kinetics. _Nat. Commun._ 2019, _10_ , 3573
  * Wang et al. 2020 Wang, Y.; Marcelo, J. L.; Tiwary, P. Machine learning approaches for analyzing and enhancing molecular dynamics simulations. _Curr. Opin. Struct. Biol._ 2020, _61_ , 139 – 145, Theory and Simulation - Macromolecular Assemblies
  * Li and Ma 2016 Li, W.; Ma, A. A benchmark for reaction coordinates in the transition path ensemble. _J. Chem. Phys._ 2016, _144_ , 134104
  * Matsumoto et al. 2002 Matsumoto, Y.; Hudson, K.; Society, A. M.; Saito, M. _An introduction to morse theory_ ; Europe and Central Asia Poverty Reduction and Economic Manag; American Mathematical Society, 2002
  * Liang et al. 1998 Liang, J.; Woodward, C.; Edelsbrunner, H. Anatomy of protein pockets and cavities: Measurement of binding site geometry and implications for ligand design. _Protein Sci._ 1998, _7_ , 1884–1897
  * Liang et al. 1998 Liang, J.; Edelsbrunner, H.; Fu, P.; Sudhakar, P. V.; Subramaniam, S. Analytical shape computation of macromolecules: I. molecular area and volume through alpha shape. _Proteins_ 1998, _33_ , 1–17
  * Edelsbrunner et al. 1998 Edelsbrunner, H.; Facello, M.; Liang, J. On the definition and the construction of pockets in macromolecules. _Discret. Appl. Math._ 1998, _88_ , 83 – 102, Computational Molecular Biology DAM - CMB Series
  * Liang et al. 1998 Liang, J.; Edelsbrunner, H.; Fu, P.; Sudhakar, P. V.; Subramaniam, S. Analytical shape computation of macromolecules: II. Inaccessible cavities in proteins. _Proteins_ 1998, _33_ , 18–29
  * Binkowski et al. 2003 Binkowski, T. A.; Adamian, L.; Liang, J. Inferring functional relationships of proteins from local sequence and spatial surface patterns. _J. Mol. Biol._ 2003, _332_ , 505–526
  * Binkowski et al. 2005 Binkowski, T. A.; Joachimiak, A.; Liang, J. Protein surface analysis for function annotation in high-throughput structural genomics pipeline. _Protein Sci._ 2005, _14_ , 2972–2981
  * Tseng and Liang 2006 Tseng, Y. Y.; Liang, J. Estimation of amino acid residue substitution rates at local spatial regions and application in protein function inference: a Bayesian Monte Carlo approach. _Mol. Biol. Evol._ 2006, _23_ , 421–436
  * Tseng et al. 2009 Tseng, Y. Y.; Dundas, J.; Liang, J. Predicting protein function and binding profile via matching of local evolutionary and geometric surface patterns. _J. Mol. Biol._ 2009, _387_ , 451–464
  * Perez-Rathke et al. 2018 Perez-Rathke, A.; Fahie, M. A.; Chisholm, C.; Liang, J.; Chen, M. Mechanism of OmpG pH-dependent gating from loop ensemble and single channel studies. _J. Am. Chem. Soc._ 2018, _140_ , 1105–1115
  * Tian et al. 2018 Tian, W.; Chen, C.; Lei, X.; Zhao, J.; Liang, J. CASTp 3.0: computed atlas of surface topography of proteins. _Nucleic Acids Res._ 2018, _46_ , W363–W367
  * Xia and Wei 2014 Xia, K.; Wei, G.-W. Persistent homology analysis of protein structure, flexibility, and folding. _International journal for numerical methods in biomedical engineering_ 2014, _30_ , 814–844
  * Xia et al. 2015 Xia, K.; Feng, X.; Tong, Y.; Wei, G. W. Persistent homology for the quantitative prediction of fullerene stability. _Journal of computational chemistry_ 2015, _36_ , 408–422
  * Xia and Wei 2015 Xia, K.; Wei, G.-W. Multidimensional persistence in biomolecular data. _Journal of computational chemistry_ 2015, _36_ , 1502–1520
  * Caiani et al. 1997 Caiani, L.; Casetti, L.; Clementi, C.; Pettini, M. Geometry of dynamics, Lyapunov exponents, and phase transitions. _Phys. Rev. Lett._ 1997, _79_ , 4361
  * Angelani et al. 2005 Angelani, L.; Ruocco, G.; Zamponi, F. Relationship between phase transitions and topological changes in one-dimensional models. _Phys. Rev. E_ 2005, _72_ , 016122
  * Kastner 2008 Kastner, M. Phase transitions and configuration space topology. _Rev, Mod. Phys._ 2008, _80_ , 167
  * Wales 2018 Wales, D. J. Exploring energy landscapes. _Annu. Rev. Phys. Chem._ 2018, _69_ , 401–425
  * Cimasoni and Delabays 2019 Cimasoni, D.; Delabays, R. The topological hypothesis for discrete spin models. _J. Stat. Mech. Theory Exp._ 2019, _2019_ , 033216
  * Kastner and Mehta 2011 Kastner, M.; Mehta, D. Phase transitions detached from stationary points of the energy landscape. _Phys. Rev. lett._ 2011, _107_ , 160602
  * Edelsbrunner and Harer 2010 Edelsbrunner, H.; Harer, J. L. _Computational topology: an introduction_ ; American Mathematical Society: Providence, RI, 2010
  * Edelsbrunner et al. 2002 Edelsbrunner, H.; Letscher, D.; Zomorodian, A. Topological persistence and simplification. _Discrete Comput. Geom._ 2002,
  * Carlsson 2009 Carlsson, G. Topology and data. _Bull. Am. Math. Soc._ 2009, _46_ , 255–308
  * Munkres 2018 Munkres, J. _Elements Of algebraic topology_ ; CRC Press, 2018
  * Hatcher 2005 Hatcher, A. _Algebraic topology_ ; Tsinghua University Press Co., Ltd.; Tsinghua University Press, 2005
  * Cohen-Steiner et al. 2005 Cohen-Steiner, D.; Edelsbrunner, H.; Harer, J. Stability of Persistence Diagrams. Proceedings of the Twenty-First Annual Symposium on Computational Geometry. New York, NY, USA, 2005; p 263–271
  * Wagner et al. 2012 Wagner, H.; Chen, C.; Vuçini, E. _Topological methods in data analysis and visualization II_ ; Springer, 2012; pp 91–106
  * Kaczynski et al. 2006 Kaczynski, T.; Mischaikow, K.; Mrozek, M. _Computational homology_ ; Springer Science & Business Media, 2006; Vol. 157
  * Tian et al. 2021 Tian, W.; Manuchehrfar, F.; Wagner, H.; Edelsbrunner, H.; Liang, J. Persistent homology and moment of probability landscapes of stochastic reaction networks and their changes. _Prepint_ 2021,
  * Bolhuis et al. 2000 Bolhuis, P. G.; Dellago, C.; Chandler, D. Reaction coordinates of biomolecular isomerization. _Proc. Natl. Acad. Sci._ 2000, _97_ , 5877–5882
  * Hess Berk and Erik 2008 Hess Berk, v. d. S. D., Kutzner Carsten; Erik, L. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. _J. Chem. Theory Comput._ 2008, _4_ , 435–447, PMID: 11972010
  * Levy et al. 1984 Levy, R.; Srinivasan, A.; Olson, W.; McCammon, J. Quasi-harmonic method for studying very low frequency modes in proteins. _Biopolymers_ 1984, _23_ , 1099–1112
  * García 1992 García, A. E. Large-amplitude nonlinear motions in proteins. _Phys. Rev. Lett._ 1992, _68_ , 2696
  * Riccardi et al. 2009 Riccardi, L.; Nguyen, P. H.; Stock, G. Free-energy landscape of RNA hairpins constructed via dihedral angle principal component analysis. _J. Phys. Chem. B_ 2009, _113_ , 16660–16668, PMID: 20028141
  * Mu et al. 2005 Mu, Y.; Nguyen, P. H.; Stock, G. Energy landscape of a small peptide revealed by dihedral angle principal component analysis. _Proteins_ 2005, _58_ , 45–52
  * Altis et al. 2007 Altis, A.; Nguyen, P. H.; Hegger, R.; Stock, G. Dihedral angle principal component analysis of molecular dynamics simulations. _J. Chem. Phys._ 2007, _126_ , 244111
  * Peters et al. 2004 Peters, B.; Heyden, A.; Bell, A. T.; Chakraborty, A. A growing string method for determining transition states: Comparison to the nudged elastic band and string methods. _J. Chem. Phys._ 2004, _120_ , 7877–7886

##  5 Supporting Information

Table S1: Locations of birth and death of the probability peaks and the ridges connecting them for dynamic probability surface projected to ($\phi,\psi,\theta_{1},\alpha,\beta$) and ($\phi,\psi,\alpha,\beta,\theta_{2}$) .  $5$-dimensional subspace of ($\phi,\psi,\theta_{1},\alpha,\beta$) | $5$-dimensional subspace of ($\phi,\psi,\alpha,\beta,\theta_{2}$)
---|---
Label | ($\phi,\psi,\theta_{1},\alpha,\beta$) Coordinate | Probability | Label | ($\phi,\psi,\alpha,\beta,\theta_{2}$) Coordinate | Probability
$b_{1}$ | ($1.25,-0.84,0.01,2.24,2.00$) | $1.95\times 10^{-3}$ | $b_{1}$ | ( $1.25,\,-0.83\,,\,2.23\,,\,2.00\,,\,0.00\,$) | $1.95\times 10^{-3}$
$b_{2}$ | ($\,-1.68,\,0.42\,,\,-0.18\,,\,2.18\,,\,1.94\,$) | $1.00\times 10^{-3}$ | $b_{2}$ | ($\,-1.25\,,\,0.00\,,\,2.18,\,2.00\,,\,0.00\,$) | $1.19\times 10^{-3}$
$d_{2}$ | ($\,-0.42\,,\,-0.42\,,\,0.19\,,\,2.34\,,\,2.15\,$) | $5.08\times 10^{-4}$ | $d_{2}$ | ($\,-0.41\,,\,-0.41\,,\,2.34\,,\,2.15\,,\,0.00\,$) | $4.19\times 10^{-4}$
$b_{3}$ | ($\,0.00\,,\,-0.42\,,\,0.00\,,\,2.24\,,\,2.10\,$) | $9.98\times 10^{-4}$ | $b_{3}$ | ($\,0\,,\,-0.42\,,\,2.34\,,\,2.15\,,\,0.00\,$) | $8.46\times 10^{-4}$
$d_{3}$ | ($\,0.42\,,\,-0.42\,,\,-0.18\,,\,2.34\,,\,2.15\,$) | $7.31\times 10^{-4}$ | $d_{3}$ | ($\,0.42\,,\,0.42\,,\,2.34\,,\,2.15\,,\,0.00\,$) | $6.70\times 10^{-4}$
$b_{4}$ | ($\,-2.52\,,\,-2.94\,,\,0.00\,,\,2.18\,,\,1.94\,$) | $1.16\times 10^{-4}$ | $b_{4}$ | ($\,-2.52\,,\,-2.94\,,\,2.18\,,\,1.94\,,\,0.00\,$) | $1.02\times 10^{-4}$
$d_{4}$ | ($\,-2.52\,,\,-1.26\,,\,-0.18\,,\,2.18\,,\,1.34\,$) | $6.73\times 10^{-6}$ | $d_{4}$ | ($\,-2.51\,,\,-1.68\,,\,2.18\,,\,1.94\,,\,0.00\,$) | $5.82\times 10^{-6}$
Table S2: Locations of birth and death of the probability peaks and the ridges connecting them for the projected 5-d probability landscape to ($\phi-\theta_{1}$) plane and ($\phi-\psi$) plane.  $p(\phi,\psi,\theta_{1},\alpha,\beta)$ projected onto ($\phi-\theta_{1}$) plane |  $p(\phi,\psi,\theta_{1},\alpha,\beta)$ projected onto ($\phi-\psi$) plane
---|---
Label | ($\phi,\theta_{1}$) Coordinate | Probability | Label | ($\phi,\psi$) Coordinate | Probability
$b_{1}$ | ($0.84,0.00$) | $6.17\times 10^{-2}$ | $b_{1}$ | ( $1.25,\,-0.83\,$) | $6.80\times 10^{-2}$
$b_{2}$ | ($\,-1.25,\,0.00\,$) | $3.41\times 10^{-2}$ | $b_{2}$ | ($\,-1.25\,,\,0.00\,$) | $4.02\times 10^{-2}$
$d_{2}$ | ($\,-0.42\,0.19\,$) | $1.06\times 10^{-2}$ | $d_{2}$ | ($\,-0.41\,,\,-0.41\,$) | $8.90\times 10^{-3}$
$b_{3}$ | ($\,0.00\,,\,0.00\,$) | $1.70\times 10^{-2}$ | $b_{3}$ | ($\,-2.52\,,\,-2.93\,$) | $3.71\times 10^{-3}$
$d_{3}$ | ($\,0.42\,,\,0.00\,$) | $1.37\times 10^{-4}$ | $d_{3}$ | ($\,-2.52\,,\,-1.68\,$) | $1.34\times 10^{-4}$
Table S3: Locations of the probability peaks and ridges connecting them for a)
the 5-d probability landscape projected to 2 principal components by dPCA. The
two most important eigenvectors $\nu_{1}$ and $\nu_{2}$ of the covariance
matrix of the conformations are
$\nu_{1}=[-0.251,\,-0.771,\,-0.112,\,0.573,\,-0.003,\,-0.014,\,0.025,\,0.019,\,0.025,0.011]^{T}$
and
$\nu_{2}=[0.339,\,-0.208,\,0.910,\,0.052,-0.003,\,0.053,\,-0.027,\,-0.022,\,-0.063,\,-0.030]^{T}$
for
($\cos(\phi),\,\sin(\phi),\,\cos(\psi),\,\sin(\psi),\,\cos(\theta_{1}),\,\sin(\theta_{1}),\,\cos(\alpha),\,\sin(\alpha),\,\cos(\beta),\,\sin(\beta)$),
with eigenvalues of $r_{1}=1.056$ and $r_{2}=0.283$, respectively.
b) The 39-d probability landscape projected to the 5 principal components by dPCA. $d_{9}$ shown in text is at $(PC_{1},\cdots,PC_{5})=(-0.005,\,-0.886,\,0.958,\,-0.382,\,-0.862)$, and the persistence is $b_{9}-d_{9}=1.38\times 10^{-4}-1.99\times 10^{-5}=1.18\times 10^{-4}$. $d_{10}$ is at $(-0.005,\,0.882,\,0.958,\,-0.382,\,-0.862)$ and the persistence is $b_{10}-d_{10}=1.21\times 10^{-4}-1.80\times 10^{-5}=1.03\times 10^{-4}$.  $p(\phi,\psi,\theta_{1},\alpha,\beta)$ projected onto ($PC_{1},PC_{2}$) plane |  $p({\mbox{\boldmath$x$}}\in\mathbb{R}^{39})$ projected onto ($PC_{i}:\,i\in\\{1,\cdots,\,5\\}$) plane
---|---
Label | ($PC_{1},PC_{2}$) Coordinate | Probability | Label | ($PC_{1},\cdots,PC_{5}$) Coordinate | Probability
$b_{1}$ | ($-1.3,\,0.5$) | $6.80\times 10^{-2}$ | $b_{1}$ | ( $-1.35,0.88,0.95,-0.38,-0.86$) | $3.16\times 10^{-4}$
$b_{2}$ | ($0.7,\,1.1\,$) | $5.69\times 10^{-2}$ | $b_{2}$ | ($\,-1.35,-0.88,0.95,-0.38,-0.86$) | $2.54\times 10^{-4}$
$d_{2}$ | ($\,-0.1\,1.3\,$) | $7.14\times 10^{-3}$ | $d_{2}$ | ($\,-1.35,-0.44,0.95,-0.38,-0.86\,$) | $2.41\times 10^{-5}$
$b_{3}$ | ($\,1.1\,,\,0.7\,$) | $3.05\times 10^{-2}$ | $b_{3}$ | ($\,-1.35,-0.88,-0.38,0.94,-0.86\,$) | $2.37\times 10^{-4}$
$d_{3}$ | ($\,1.1\,,\,0.9\,$) | $2.40\times 10^{-2}$ | $d_{3}$ | ($\,-1.35,-0.88,0.95,0.37,-0.68\,$) | $5.99\times 10^{-5}$
$b_{4}$ | ($\,1.3\,,\,0.3\,$) | $2.36\times 10^{-2}$ | $b_{4}$ | ($\,-1.35,0.88,-0.38,0.94,-0.86\,$) | $2.29\times 10^{-4}$
$d_{4}$ | ($\,1.3\,,\,0.5\,$) | $1.29\times 10^{-2}$ | $d_{4}$ | ($\,-1.35,0.88,1.14,0.56,-0.51\,$) | $5.30\times 10^{-5}$
$b_{5}$ | ($\,1.3\,,\,-0.1\,$) | $1.58\times 10^{-2}$ | $b_{5}$ | ($\,-1.35,-0.88,-0.95,0.37,0.86\,$) | $2.03\times 10^{-4}$
$d_{5}$ | ($\,1.3\,,\,0.1\,$) | $9.40\times 10^{-3}$ | $d_{5}$ | ($\,-1.35,-0.88,0.00,-1.33,-0.51\,$) | $7.03\times 10^{-5}$
$b_{6}$ | ($\,-0.3\,,\,1.3\,$) | $1.45\times 10^{-2}$ | $b_{6}$ | ($\,-1.35,0.88,0.38,-0.95,0.86\,$) | $1.86\times 10^{-4}$
$d_{6}$ | ($\,-0.5\,,\,1.1\,$) | $1.12\times 10^{-2}$ | $d_{6}$ | ($\,-1.35,0.88,0.19,-1.14,-0.68\,$) | $4.66\times 10^{-5}$

### Dimension reduction by PCA

#### Projection of $p(\phi,\psi,\theta_{1},\alpha,\beta)$ onto ($PC1,PC2$) by
PCA.

Here we applied direct PCA without using dihedral PCA to the 5-d probability
surface. We first shift each direction to its mean value and then perform PCA
to obtain the first two principal components from the covariance matrix of the
$1.5\times 10^{11}$ conformations.

Figure S1: The dynamic probability surface on PCA space, its topological
structure, and committor values on principal components space. (a) The
projection of the 5-d dynamic probability surface
$p(\phi,\,\psi,\,\theta_{1},\,\alpha,\,\beta)$ to the plane of
($PC_{1}$-$PC_{2}$). (b) The persistent diagram exhibits six probability
peaks. (c-d) Distribution of committor values $p_{B}$ for trajectories from
bridge $d_{2}$, peak $d_{6}$, respectively.

The dynamic probability surface after projection to the ($PC1,PC2$) plane is
shown in Fig. S1a. This probability surface is dramatically different than
that of the 5-d surface when projected onto ($\phi-\theta_{1}$) plane or
($\phi-\psi$) plane. The persistence diagram also shows different peaks with
different birth and death places (Fig. S1b). Committor tests show that all
committor values for the conformations at $b_{2-6}$, and $d_{3-5}$ are 0.0
with trajectories going to the reactant basin. All committor values for
conformations at $b_{1}$ are 1.0: Trajectories from the location of $b_{1}$
all go to the product basin. Committor values for ridges $d_{2}$ and $d_{6}$
follow a one-sided distribution. Trajectories starting from these locations
mostly fall back to the reactant basin. Overall, non of the topological
feature after PCA retain dynamic properties of the transition state
conformations as the oroginal 5-d surface.

#### Projecting from 39-d space to 5-d PCA subspace.

We applied direct PCA to the original full dimensional dynamic probability
surface, after removal of 21 bond lengths, to reduce the remaining
39-dimensional subspace to 5 principal components. The contour plot of the 5-d
PCA surface projected onto the first two principal components is shown in Fig.
S2a. Persistence homology analysis identifies 9 peaks (Fig. S2b), each located
in a 5-d cube. These 5-d persistent diagram is very different from the
persistence diagram shown for the original 5-d subspace.

Figure S2: The dynamic probability surface on 5-d principal components space
reduced by dPCA from the 39-d configuration space, its topological structure,
and committor values. (a) the 5-d dynamic probability surface
$p(PC_{1},\cdots,PC_{5})$ shown on to the ($PC_{1},PC_{2}$) plane. (b) The
persistent diagram exhibits 9 peaks. (c-e) Distribution of committor $p_{B}$
values for trajectories started at ridges $d_{6}$, $d_{7}$, and $d_{8}$,
respectively.

The committor tests show that all committor values for conformations located
at the peaks $b_{6-8}$ are 0.0. That is all the trajectories starting from
these locations fall back to the reactant basin. In addition, all committor
values for the conformations located at the peaks $b_{1-5},b_{9}$, and the
ridges $d_{2-5},d_{9}$ are 1.0. Trajectories from these locations move forward
to the product basin. Committor values for ridges $d_{6-8}$ follow one-sided
distributions with the peak at $p_{B}=0$ (Fig. S2c-e). Trajectories starting
from these ridges mostly fall back to the reactant basin.

Overall, These results demonstrate that dimension reduction by PCA destroys
dynamic properties inherent on the original surface topology of the dynamic
probability surface.