# Finding dominant reaction pathways via global optimization of action

Juyong Lee  Laboratory of Computational Biology, National Heart, Lung, and
Blood Institute (NHLBI), National Institutes of Health (NIH), Bethesda,
Maryland 20892, U.S.A.  In-Ho Lee  Korea Research Institute of Standards and
Science, Daejeon 34113, Republic of Korea  InSuk Joung  School of
Computational Sciences, Korea Institute of Advanced Study, Seoul 02455, Korea
Jooyoung Lee<EMAIL_ADDRESS>School of Computational Sciences, Korea
Institute of Advanced Study, Seoul 02455, Korea  Bernard R. Brooks
<EMAIL_ADDRESS>Laboratory of Computational Biology, National Heart, Lung, and
Blood Institute (NHLBI), National Institutes of Health (NIH), Bethesda,
Maryland 20892, USA

###### Abstract

We present a new computational approach to sample multiple reaction pathways
with fixed initial and final states through global optimization of the
Onsager-Machlup action using the conformational space annealing method. This
approach successfully samples not only the most dominant pathway but also many
other possible paths without initial guesses on reaction pathways. Pathway
space is efficiently sampled by crossover operations of a set of paths and
preserving the diversity of sampled pathways. The sampling ability of the
approach is assessed by finding pathways for the conformational changes of
alanine dipeptide and hexane. The benchmarks demonstrate that the rank order
and the transition time distribution of multiple pathways identified by the
new approach are in good agreement with those of long molecular dynamics
simulations. We also show that the folding pathway of the mini-protein FSD-1
identified by the new approach is consistent with previous molecular dynamics
simulations and experiments.

###### pacs:

87.15.ap,87.15.hm,82.30.Qt,87.15.hp

Finding multiple reaction pathways between two end states remains a
challenging problem in computational biophysics Elber (2016). For this
purpose, performing a long-time molecular dynamics (MD) simulation is a
commonly used approach. Despite recent progress in the methodologies of MD,
this approach still suffers from a timescale problem. Many biological
reactions such as protein folding and conformational transitions occur in the
microsecond or millisecond ranges, which are hard to be performed with
conventional computers. In addition, MD simulations starting from one end
state are not guaranteed to reach the other end state of interest. Thus
developing an efficient computational method to find multiple possible
reaction pathways connecting two end states is necessary. There are currently
no methods that can efficiently produce the multiple dominant pathways
connecting two well-defined end point states in a complex system. The
objective of this work is to present such a method. Other approaches using a
conformational driving force do not sample alternatives. Methods that are
robust, such as transition path sampling Dellago  _et al._ (1998); Bolhuis
_et al._ (2002), are very expensive to use for complex systems with and
multiple steps and barriers.

Various chain-of-state methods have been suggested based on the assumption
that a dominant transition pathway between two states follows the minimum
energy pathway Czerminski and Elber (1990); Henkelman  _et al._ (2000); E  _et
al._ (2002). The limitations of these methods are that they do not consider
the dynamics of a system and find only the nearest local minimum solution from
a given initial pathway Elber (2016). Alternative methods based on the
principle of least action have been suggested Gillilan and Wilson (1992);
Olender and Elber (1996); Elber  _et al._ (1999); Elber and Shalloway (2000);
Passerone and Parrinello (2001); Passerone  _et al._ (2003); Lee  _et al._
(2003a). Passerone and Parrinello suggested the action-derived molecular
dynamics (ADMD) method based on the combination of classical action and a
penalty term that conserves the total energy of a system Passerone and
Parrinello (2001); Passerone  _et al._ (2003). To enhance the convergence of
ADMD calculations, Lee et al. introduced a kinetic energy penalty term based
on the equipartition theorem Lee  _et al._ (2003a, 2005, 2012a, 2013a).
Although ADMD yields physically relevant pathways, it has two practical
limitations Crehuet and Field (2003); Lee  _et al._ (2003a): i) it strongly
depends on the initial guess of a pathway; and ii) it cannot identify the most
dominant pathway when there are multiple pathways because the classical least
action principle is an extremum principle.

For diffusive processes, the second problem can be avoided by using the
Onsager-Machlup (OM) action $S_{\text{OM}}$ Onsager and Machlup (1953);
Machlup and Onsager (1953); Olender and Elber (1996); Eastman  _et al._
(2001); Zuckerman and Woolf (2000); Fujisaki  _et al._ (2010, 2013). Onsager
and Machlup showed that the relative probability to observe a pathway with an
OM action of $S$ is proportional to $e^{-S/k_{\text{B}}T}$ where
$k_{\text{B}}$ is the Boltzmann constant and $T$ is a temperature. Thus the
most dominant pathway corresponds to the one that minimizes $S_{\text{OM}}$
and the same result can be obtained by solving the Fokker-Planck equation
Faccioli  _et al._ (2006); Sega  _et al._ (2007); Beccara  _et al._ (2012).
This property recasts the problem of finding dominant pathways into a global
optimization problem. However, finding the global minimum of $S_{\text{OM}}$
is a numerically challenging task because the minimization of $S_{\text{OM}}$
requires the second derivatives of a potential function, which are
computationally expensive, at best, and wholly unavailable for may quantum
mechanical energy surfaces.

In this work, we propose an efficient computational method that finds not only
the most dominant pathway but also multiple suboptimal pathways without second
derivative calculations. For global optimization of $S_{\text{OM}}$, we used
an efficient global optimization method called conformational space annealing
(CSA) based on a combination of genetic algorithm, simulated annealing, and
Monte Carlo minimization Lee  _et al._ (1997, 1999). The CSA method has been
demonstrated to be extremely efficient in solving various global optimization
problems including finding low energy conformations of Lennard-Jones clusters
Lee  _et al._ (2003b), protein structure prediction Lee  _et al._ (1999,
2003b); Steinbach (2004); Joo  _et al._ (2007); Lee  _et al._ (2008); Joo  _et
al._ (2009); Lee  _et al._ (2011), multiple sequence alignment Joo  _et al._
(2008), and community detection in networks Lee  _et al._ (2012b); Lee and Lee
(2013); Lee  _et al._ (2013b). We extend the CSA approach to examining
pathways, preserving all features that make it robust and efficient, by
applying it to sets of entire pathways represented as a chain-of-states. From
benchmark simulations using alanine dipeptide, we observed that our method
finds multiple transition pathways, which are consistent with long-time
Langevin dynamics (LD) simulations. In addition, the rank order statistics and
transition time distributions of the multiple transition pathways are in good
agreement with those of the LD results. We will call the method Action-CSA.

Here, we briefly review the theoretical background behind Action-CSA. If a
system with $N$ atoms with a potential energy $V$ follows the overdamped
Langevin dynamics,

$\gamma\mathbf{m}\mathbf{\dot{x}}=-\frac{\partial
V}{\partial\mathbf{x}}+\mathbf{R},$ (1)

where $\mathbf{m}$ is a diagonal mass matrix, $\mathbf{x}$ is a $3N$
dimensional coordinate vector, $\gamma$ is collision frequency, and
$\mathbf{R}$ is a Gaussian random force, the relative probability of finding a
final state $\mathbf{x}_{\textup{f}}$ at a time $t$ from an initial state
$\mathbf{x}_{\textup{i}}$ via diffusive trajectories $\mathbf{x}(t)$ is
determined by using the path integral approach and OM action
$S_{\text{OM}}[\mathbf{x}(t)]$ Onsager and Machlup (1953); Machlup and Onsager
(1953):

$P(\mathbf{x}_{\textup{f}}|\mathbf{x}_{\textup{i}};t)=\int_{0}^{t}\mathcal{D}\mathbf{x}(t)\;\text{exp}\left(-\frac{S_{\text{OM}}[\mathbf{x}(t)]}{k_{\text{B}}T}\right),$
(2)

where $\mathcal{D}\mathbf{x}(t)$ indicates that the integration runs over all
possible pathways $\mathbf{x}(t)$. This relationship suggests that if the
$S_{\text{OM}}$ values of all physically accessible pathways are obtained, one
can determine the relative populations of multiple pathways. Thus,
$S_{\text{OM}}$ is a proper target objective function of global optimization.
The generalized OM action of a pathway $\mathbf{x}(t)$ is defined Onsager and
Machlup (1953); Machlup and Onsager (1953); Hunt and Ross (1981); Adib (2008):

$S_{\text{OM}}[\mathbf{x}(t)]=\frac{\Delta
V}{2}+\frac{1}{4\gamma}\int_{0}^{t}d\tau\left\\{\left[\gamma\mathbf{m}\mathbf{\dot{x}}(\tau)\right]^{2}+\norm{\nabla
V[\mathbf{x}(\tau)]}^{2}-2k_{\text{B}}T\nabla^{2}V[\mathbf{x}(\tau)]\right\\},$
(3)

where $\Delta V=V(\mathbf{x}_{\textup{f}})-V(\mathbf{x}_{\textup{i}})$. The
last term is related to trajectory entropy connected with fluctuations Adib
(2008); Haas  _et al._ (2014a, b) and was not presented in the orignal work by
Onsager and Machlup because the harmonic potential was considered Onsager and
Machlup (1953); Machlup and Onsager (1953). Note that the minimization of
$S_{\text{OM}}$ using analytic local minimization algorithms requires third
derivatives. This makes the direct global optimization of $S_{\text{OM}}$ hard
to be applied to detect transition pathways of biomolecules with all-atom
force fields due to the complexity of implementation and high computational
cost. For numerical calculations based on a chain-of-state representation, the
OM action should be discretized. The method uses the second-order
discretization of the symmetric OM formula, which uses only gradients for
$S_{\text{OM}}$ calculations Miller and Predescu (2007):

$\begin{split}S_{\text{OM}}[\mathbf{x}(t)]=\frac{\Delta
V}{2}+\sum_{i=0}^{P-1}\frac{\Delta
t}{4\gamma}\left\\{\left[\frac{\gamma\mathbf{m}(\mathbf{x}_{i+1}-\mathbf{x}_{i})}{\Delta
t}\right]^{2}+\frac{\norm{\nabla V(\mathbf{x}_{i})}^{2}+\norm{\nabla
V(\mathbf{x}_{i+1})}^{2}}{2}\right.\\\
\left.-\frac{\gamma\mathbf{m}(\mathbf{x}_{i+1}-\mathbf{x}_{i})}{\Delta
t}[\nabla V(\mathbf{x}_{i+1})-\nabla V(\mathbf{x}_{i})]\right\\},\end{split}$
(4)

where $P+1$ is the number of replicas, $\Delta t$ is a time step between
successive replicas, and $t=P\Delta t$ is the total transition time. This
formula is superior to the direct implementation of Eq. (3) since it requires
only the first derivatives of $V$ to evaluate $S_{\text{OM}}$.

Here, we describe the application of CSA to optimize $S_{\text{OM}}$. In
general, a pathway is represented as a chain of $P-1$ replicas with $N$ atoms
for each replica leading to $3N(P-1)$ total degrees of freedom. Each replica
is represented by a sequence of $3N-6$ internal dihedral angles and 6 net
translational/rotational degrees of freedom. An Action-CSA calculation starts
with a set of random pathways on a pathway space. Subsequently, the actions of
the random pathways should be locally optimized. We call this set of pathways
a bank, and update conformations in the bank during the Action-CSA
calculation. As stated previously, direct minimization of $S_{\text{OM}}$
using analytic gradients is computationally challenging.

For a computationally feasible local action optimization, we optimized a
pathway using a modified action from ADMD instead of using $S_{\text{OM}}$.
The discretized classical action is defined:

$S_{\text{classical}}[\mathbf{x}(t)]=\sum_{i=0}^{P-1}L_{i}(\mathbf{x}_{i})\Delta
t=\sum_{i=0}^{P-1}\left[\frac{\mathbf{m}(\mathbf{x}_{i}-\mathbf{x}_{i+1})^{2}}{2\Delta
t^{2}}-V(\mathbf{x}_{i})\right]\Delta t.$ (5)

Physically accessible pathways correspond to the stationary points of
$S_{\text{classical}}$. Finding such pathways is a computationally difficult
task because $S_{\text{classical}}$ is not bounded; $S_{\text{classical}}$ can
be minimized or maximized, and the stationary points of $S_{\text{classical}}$
can be minima, maxima or saddle points. Another practical problem is that the
total energies of pathways satisfying the stationary condition $\delta
S_{\text{classical}}=0$ may not be conserved Passerone and Parrinello (2001).
To find pathways that satisfy the principle of least action and conserve total
energies, a modified action with a penalty term restraining total energy was
suggested Passerone and Parrinello (2001):

$\begin{split}\Theta(\mathbf{x}_{i};E)&=\mu_{A}S_{\text{classical}}+\mu_{E}\sum_{i=0}^{P-1}(E_{i}-E)^{2}\\\
&=\mu_{A}\sum_{i=0}^{P-1}\left[\frac{\mathbf{m}(\mathbf{x}_{i}-\mathbf{x}_{i+1})^{2}}{2\Delta
t^{2}}-V(\mathbf{x}_{i})\right]\Delta
t+\mu_{E}\sum_{i=0}^{P-1}\left\\{\left[\frac{\mathbf{m}(\mathbf{x}_{i}-\mathbf{x}_{i+1})^{2}}{2\Delta
t^{2}}+V(\mathbf{x}_{i})\right]-E\right\\}^{2},\end{split}$ (6)

where $E$ is a targeted total energy of a system, $\mu_{A}$ and $\mu_{E}$ are
the weighting parameters of the classical action, and the restraint term for
energy conservation. The minimization of $\Theta[\mathbf{x}(t);E]$ requires
only the first derivatives of $V$.

We call the set of locally optimized initial random pathways using
$\Theta[\mathbf{x}(t);E]$ the _first bank_. The first bank remains the same
throughout the optimization and is used as the reservoir of partially
optimized pathways to enhance the diversity of pathway search. A copy of the
first bank is generated and called a _bank_. The conformations in the bank are
updated during a calculation while the size of the bank is kept constant. By
using the pathways included in the first bank and the bank, new trial pathways
are generated by performing crossover and random perturbation operations. For
a crossover operation, two pathways, a seed pathway from the bank and a random
pathway either from the bank or the first bank, are selected and random parts
of two selected configurations are swapped. For a random perturbation, a
certain number of degrees of freedom of a seed pathway, up to 5% of total
degrees of freedom, were randomly changed. The generated trial pathways are
locally optimized using $\Theta[\mathbf{x}(t);E]$ to remove any possible
artifacts generated by the crossover and the random perturbation operations.
However, after local minimizations, the bank was updated by comparing the
$S_{\text{OM}}$ values of the existing pathways and the new ones instead of
$\Theta[\mathbf{x}(t);E]$.

A key feature of CSA is a sophisticated bank-update procedure that prevents a
search being trapped in local minima during the optimization and keeps the
diversity of the bank. For a newly obtained configuration $\alpha$, the
distances between $\alpha$ and the existing ones in the bank are calculated.
If the distance $D$ between $\alpha$ and its closest neighbor is less than a
cutoff distance $D_{\text{cut}}$, only the better configuration in terms of
the objective function is selected. If $D>D_{\text{cut}}$, $\alpha$ is
considered a new configuration and it replaces the worst configuration in the
bank. At initial stages of a calculation, $D_{\text{cut}}$ is kept large for
wider sampling. As the calculation proceeds, it gradually decreases for a
refined search near the global minimum. The bank keeps updating until no
pathway with a lower $S_{\text{OM}}$ is found. In this work, a distance
between two pathways was measured by the Fréchet distance Alt and Godau
(1995). More details on a general CSA procedure are described elsewhere Lee
_et al._ (1997, 1999, 2003b); Joo  _et al._ (2008, 2009); Lee  _et al._
(2011).

To verify that Action-CSA successfully finds multiple transition pathways and
allows one to determine the rank order of the pathways based on their
optimized $S_{\text{OM}}$ values, we applied our method to investigate the
conformational transition of alanine dipeptide from $C7_{\text{eq}}$ to
$C7_{\text{ax}}$ in the vacuum. Here, we used the polar hydrogen
representation in the PARAM19 force field Neria  _et al._ (1996) and
dielectric constant was set to 1.0 Loncharich  _et al._ (1992). We performed
Action-CSA simulations with different transition times, $t$ in Eqs. (4) and
(6), ranging from 0.2 ps to 2.0 ps with an interval of 0.1 ps. The numbers of
replicas were adjusted with $t$ to keep the time step between replicas $\Delta
t$ = 5 fs. All simulations were performed at temperature $T$ = 350 K with a
collision frequency $\gamma$ = 1.0 ps-1. The reference total energy $E$ in Eq.
(5) was obtained by adding the initial potential energy
$V(\mathbf{x}_{i})=-43.3$ kcal/mol and a kinetic energy of 12.5 kcal/mol
estimated by $3Nk_{\text{B}}T/2$ with the number of atoms $N$ = 12. The
weighting parameters $\mu_{A}$ and $\mu_{E}$ in Eq. (5) were set to $-1.0$ and
$1.0$, respectively. For comparison purposes, we performed 500 $\mu$s LD
simulations of alanine dipeptide under the same condition and counted the
number of the $C7_{\text{eq}}$ $\to$ $C7_{\text{ax}}$ transitions.

Now we will show that the Action-CSA identifies not only the most dominant
pathway but also multiple possible pathways. We identified 8 different
pathways for the $C7_{\text{eq}}$ $\to$ $C7_{\text{ax}}$ transition by
clustering all pathways sampled from the Action-CSA simulations (Figure 1A).
From the $S_{\text{OM}}$ values obtained with different transition times
(Figure 1B), it is clear that the pathway that passes barrier B has the lowest
$S_{\text{OM}}$ values along all transition times tested, which indicates that
it is the most dominant pathway regardless of transition time. This is
consistent with the LD simulation results. To compare the Action-CSA result
with LD, we performed 5000 independent 100 ns LD simulations amounting to 500
$\mu$s trajectories (Table 1). From the simulations, we observed 1350
transitions from $C7_{\text{eq}}$ to $C7_{\text{ax}}$ and categorized them by
finding the nearest neighbor from the 8 pathways obtained with CSA. Among
them, the pathway that crosses barrier B was identified as the most dominant
one with all transition times This demonstrates that Action-CSA correctly
identified the minimum OM action pathway and that it matches the most dominant
pathway observed in LD simulations.

Figure 1: Upper panel: eight different pathways for the $C7_{\text{eq}}$ $\to$ $C7_{\text{ax}}$ transition and the potential energy surface for the $\Phi$ and $\Psi$ angles with the PARAM19 force field (in units of kcal/mol). Potential energy barriers are labeled in order of their heights (from A to F). Lower panel: the $S_{\text{OM}}$ values of 6 pathways for the $C7_{\text{eq}}$ $\to$ $C7_{\text{ax}}$ transitions of alanine dipeptide along different transition times.  Table 1: The frequences of pathways observed from 500 $\mu$s Langevin dynamics simulations. Path ID | Frequency
---|---
Path1 | 1183
Path2 | 116
Path3 | 25
Path4 | 7
Path5 | 4
Path6 | 4
Path7 | 10
Path8 | 1

In addition, it is also identified that the Action-CSA simulations can provide
information on the transition times of various pathways. Until $t<0.8$ ps, the
pathway that crosses barrier C (Path2) has the second lowest $S_{\text{OM}}$
and the lowest $S_{\text{OM}}$ value was observed at 0.4 ps. These are
consistent with the LD results in which all 118 transitions that crossed
barrier C occurred within 1.1 ps and their most probable transition time was
0.7 ps (the inset of Fig 1B). However, when $t>0.8$ ps, Path3, which passes
the fully extended conformation region ($\Phi$, $\Psi$) = (-180∘, 180∘) and
barrier A and B becomes the pathway with the second lowest $S_{\text{OM}}$.
From the LD simulations, when $t>0.9$ ps, 25 pathways similar to Path3 were
identified, which makes them the second dominant pathway. These results
demonstrate that the profile of $S_{\text{OM}}$ values is consistent with the
distributions of transition times obtained from the LD simulations. Note that
the most probable transitions times observed from the LD simulations are
longer than the minimum action transition times obtained from the CSA
simulations. This is because high-frequency motions due to thermal
fluctuations are filtered out in the minimum action pathways Olender and Elber
(1996); Elber  _et al._ (1999); Elber (2016). This means that the dwell time
is well filtered out in the simulation, where a physically sufficient sampling
time is assumed.

The second example is finding possible pathways for the conformational change
of hexane from the all-gauche(-) (g-g-g-) to the all-gauche(+) state (g+g+g+).
We assessed the sampling ability of Action-CSA by investigating how many
pathways are found. If we assume that dihedral angles do not cross a high
barrier around the _cis_ state, all possible transition pathways between the
two all-gauche states can be enumerated (Table LABEL:tab:hexane_pathway_list).
For this reaction, there exist 44 possible pathways in total. If the
symmetries of dihedral angles and atomic order are considered, all 44 pathways
can be categorized into 14 unique pathway types. We repeated the Action-CSA
calculation of the reaction 40 times by using 200 trial pathways consisting of
100 replicas and a transition time of 3 ps.

In all simulations, the 6 lowest action pathways, C+C+, C-C-, T+C+, T+C-,
C+M+, and C+M-, were found consistently. The other higher action pathways
except for the highest action pathway, M+XM-, were found in at least 29 out of
40 CSA simulations. Only M+XM- was found in 9 simulations. On average, a
single CSA simulation sampled 12 out of 14 unique path types and 26 out of 44
possible pathways. These results show that Action-CSA assuredly samples a
number of lowest action, most dominant, pathways. The majority of the
remaining pathways with higher actions can also be found with a tendency that
relatively lower action pathways are more likely to be found. The sampling
ability of Action-CSA can be improved by using a larger bank size. The
potential energy landscape of the C+C+ pathway corresponding to the least
$S_{\text{OM}}$ shows that hexane crosses 6 energy barriers (Fig. 2).

Figure 2: The changes of potential energy and Onsager-Machlup action along the
lowest action pathway between the all-gauche(-) to the all-gauche(+)
conformations of hexane in the vacuum, the C+C+ pathway.

The third example is finding the folding pathway of FSD-1, a 28-residue mini-
protein that has been widely investigated as a model system for studying
protein folding Jang  _et al._ (2002); Lei and Duan (2004); Lei  _et al._
(2006); Wu and Shea (2010); Lei  _et al._ (2009); Lee  _et al._
(2012a).Folding pathways of FSD-1 from the fully extended conformation to the
native structure were represented by using 100 replicas, a total folding time
of 10 ps, and a temperature of 300 K. The protein was represented by the
PARAM19 force field Neria  _et al._ (1996) and solvation effects were
considered using the FACTS implicit solvent model Haberthür and Caflisch
(2008).

The characteristics of the identified lowest action folding pathway are
consistent with experiments where the N-terminal $\beta$-hairpin is more
flexible than the C-terminal $\alpha$-helix Feng  _et al._ (2009). A
comparison of the RMSD values indicates that the $\alpha$-helix folds first.
Afterward, the folding of $\beta$-hairpin and the formation of hydrophobic
core occur concurrently (the upper panel of Fig. 3). The potential energy
landscape of FSD-1 folding shows that the potential energy decreases quickly
after the 80 step suggesting that the step may be the transition state of
folding (the lower panel of Fig. 3). The conformation at the 80 step shows
that the helix is almost formed while the C-terminal region is not folded and
the hydrophobic core is partially exposed, which is similar to transition
states observed in previous conventional MD simulations Jang  _et al._ (2002);
Lei and Duan (2004); Lei  _et al._ (2006); Wu and Shea (2010).

Figure 3: Upper panel: The RMSD values of the entire FSD-1 (blue), the
C-terminal $\alpha$-helix (green, residue 14-28), and the N-terminal
$\beta$-hairpin (red, residue 1-13) from the native structure along the
folding pathway. Lower panel: The evolutions of potential energy (blue) and
the Onsager-Machlup action (red) of FSD-1 along the folding pathway.

In conclusion, we demonstrated that efficient global optimization of Onsager-
Machlup action reveals multiple reaction pathways successfully. In this work,
we introduced a new computational method that samples not only the most
dominant pathway but also other possible pathways by optimizing Onsager-
Machlup action using the CSA method. The advantages of our method over
existing pathway sampling methods are the fact that it samples multiple
pathways regardless of the quality of initial guesses on pathways; it requires
only the calculation of first derivatives; and its sampling ability is not
limited by high energy barriers separating pathways, which is a major limiting
factor of previous MD-based pathway sampling methods in exploring pathway
space Eastman  _et al._ (2001); Faccioli  _et al._ (2006); Mathews and Case
(2006); Sega  _et al._ (2007); Fujisaki  _et al._ (2010); Beccara  _et al._
(2012); Fujisaki  _et al._ (2013). In addition, we identified that the profile
of minimum Onsager-Machlup actions found with different transition time
parameters provide kinetic information on pathways. In terms of
implementation, Action-CSA calculation is massively parallel because the local
minimization of each trial pathway is independent of each other. Thus, pathway
samplings for larger systems are readily possible with help of a large cluster
system.

Table 2: List of 14 unique pathway types and 44 non redundant pathways for conformational change of hexane from g-g-g- to g+g+g+. Unique path type |  Non redundant path
---|---
C+C+ |  g-g-g- $\rightarrow$ tg-g- $\rightarrow$ tg-t $\rightarrow$ ttt $\rightarrow$ tg+t $\rightarrow$ tg+g+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-g-t $\rightarrow$ tg-t $\rightarrow$ ttt
$\rightarrow$ tg+t $\rightarrow$ g+g+t $\rightarrow$ g+g+g+
C+C- |  g-g-g- $\rightarrow$ g-g-t $\rightarrow$ tg-t $\rightarrow$ ttt $\rightarrow$ tg+t $\rightarrow$ tg+g+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ tg-g- $\rightarrow$ tg-t $\rightarrow$ ttt
$\rightarrow$ tg+t $\rightarrow$ g+g+t $\rightarrow$ g+g+g+
T+C+ |  g-g-g- $\rightarrow$ tg-g- $\rightarrow$ tg-t $\rightarrow$ ttt $\rightarrow$ ttg+ $\rightarrow$ tg+g+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-g-t $\rightarrow$ tg-t $\rightarrow$ ttt
$\rightarrow$ g+tt $\rightarrow$ g+g+t $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-g-t $\rightarrow$ g-tt $\rightarrow$ ttt
$\rightarrow$ tg+t $\rightarrow$ g+g+t $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ tg-g- $\rightarrow$ ttg- $\rightarrow$ ttt
$\rightarrow$ tg+t $\rightarrow$ tg+g+ $\rightarrow$ g+g+g+
T+C- |  g-g-g- $\rightarrow$ g-g-t $\rightarrow$ tg-t $\rightarrow$ ttt $\rightarrow$ ttg+ $\rightarrow$ tg+g+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ tg-g- $\rightarrow$ tg-t $\rightarrow$ ttt
$\rightarrow$ g+tt $\rightarrow$ g+g+t $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-g-t $\rightarrow$ g-tt $\rightarrow$ ttt
$\rightarrow$ tg+t $\rightarrow$ tg+g+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ tg-g- $\rightarrow$ ttg- $\rightarrow$ ttt
$\rightarrow$ tg+t $\rightarrow$ g+g+t $\rightarrow$ g+g+g+
C+M+ |  g-g-g- $\rightarrow$ tg-g- $\rightarrow$ tg-t $\rightarrow$ ttt $\rightarrow$ ttg+ $\rightarrow$ g+tg+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-g-t $\rightarrow$ tg-t $\rightarrow$ ttt
$\rightarrow$ g+tt $\rightarrow$ g+tg+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-tg- $\rightarrow$ ttg- $\rightarrow$ ttt
$\rightarrow$ tg+t $\rightarrow$ g+g+t $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-tg- $\rightarrow$ g-tt $\rightarrow$ ttt
$\rightarrow$ tg+t $\rightarrow$ tg+g+ $\rightarrow$ g+g+g+
C+M- |  g-g-g- $\rightarrow$ g-g-t $\rightarrow$ tg-t $\rightarrow$ ttt $\rightarrow$ ttg+ $\rightarrow$ g+tg+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ tg-g- $\rightarrow$ tg-t $\rightarrow$ ttt
$\rightarrow$ g+tt $\rightarrow$ g+tg+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-tg- $\rightarrow$ g-tt $\rightarrow$ ttt
$\rightarrow$ tg+t $\rightarrow$ g+g+t $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-tg- $\rightarrow$ ttg- $\rightarrow$ ttt
$\rightarrow$ tg+t $\rightarrow$ tg+g+ $\rightarrow$ g+g+g+
T+T+ |  g-g-g- $\rightarrow$ tg-g- $\rightarrow$ ttg- $\rightarrow$ ttt $\rightarrow$ ttg+ $\rightarrow$ tg+g+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-g-t $\rightarrow$ g-tt $\rightarrow$ ttt
$\rightarrow$ g+tt $\rightarrow$ g+g+t $\rightarrow$ g+g+g+
T+T- |  g-g-g- $\rightarrow$ tg-g- $\rightarrow$ ttg- $\rightarrow$ ttt $\rightarrow$ g+tt $\rightarrow$ g+g+t $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-g-t $\rightarrow$ g-tt $\rightarrow$ ttt
$\rightarrow$ ttg+ $\rightarrow$ tg+g+ $\rightarrow$ g+g+g+
T+M+ |  g-g-g- $\rightarrow$ g-g-t $\rightarrow$ g-tt $\rightarrow$ ttt $\rightarrow$ g+tt $\rightarrow$ g+tg+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ tg-g- $\rightarrow$ ttg- $\rightarrow$ ttt
$\rightarrow$ ttg+ $\rightarrow$ g+tg+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-tg- $\rightarrow$ g-tt $\rightarrow$ ttt
$\rightarrow$ g+tt $\rightarrow$ g+g+t $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-tg- $\rightarrow$ ttg- $\rightarrow$ ttt
$\rightarrow$ ttg+ $\rightarrow$ tg+g+ $\rightarrow$ g+g+g+
T+M- |  g-g-g- $\rightarrow$ g-g-t $\rightarrow$ g-tt $\rightarrow$ ttt $\rightarrow$ ttg+ $\rightarrow$ g+tg+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ tg-g- $\rightarrow$ ttg- $\rightarrow$ ttt
$\rightarrow$ g+tt $\rightarrow$ g+tg+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-tg- $\rightarrow$ g-tt $\rightarrow$ ttt
$\rightarrow$ ttg+ $\rightarrow$ tg+g+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-tg- $\rightarrow$ ttg- $\rightarrow$ ttt
$\rightarrow$ g+tt $\rightarrow$ g+g+t $\rightarrow$ g+g+g+
M+M- |  g-g-g- $\rightarrow$ g-tg- $\rightarrow$ g-tt $\rightarrow$ ttt $\rightarrow$ ttg+ $\rightarrow$ g+tg+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-tg- $\rightarrow$ ttg- $\rightarrow$ ttt
$\rightarrow$ g+tt $\rightarrow$ g+tg+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-tg- $\rightarrow$ g-tt $\rightarrow$ ttt
$\rightarrow$ g+tt $\rightarrow$ g+tg+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-tg- $\rightarrow$ ttg- $\rightarrow$ ttt
$\rightarrow$ ttg+ $\rightarrow$ g+tg+ $\rightarrow$ g+g+g+
T+XT- |  g-g-g- $\rightarrow$ g-g-t $\rightarrow$ g-tt $\rightarrow$ g-tg+ $\rightarrow$ ttg+ $\rightarrow$ tg+g+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ tg-g- $\rightarrow$ ttg- $\rightarrow$ g+tg-
$\rightarrow$ g+tt $\rightarrow$ g+g+t $\rightarrow$ g+g+g+
TXM |  g-g-g- $\rightarrow$ g-tg- $\rightarrow$ g-tt $\rightarrow$ g-tg+ $\rightarrow$ ttg+ $\rightarrow$ tg+g+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ tg-g- $\rightarrow$ ttg- $\rightarrow$ g+tg-
$\rightarrow$ g+tt $\rightarrow$ g+tg+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-g-t $\rightarrow$ g-tt $\rightarrow$ g-tg+
$\rightarrow$ ttg+ $\rightarrow$ g+tg+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-tg- $\rightarrow$ ttg- $\rightarrow$ g+tg-
$\rightarrow$ g+tt $\rightarrow$ g+g+t $\rightarrow$ g+g+g+
M+XM- |  g-g-g- $\rightarrow$ g-tg- $\rightarrow$ g-tt $\rightarrow$ g-tg+ $\rightarrow$ ttg+ $\rightarrow$ g+tg+ $\rightarrow$ g+g+g+
|  g-g-g- $\rightarrow$ g-tg- $\rightarrow$ ttg- $\rightarrow$ g+tg-
$\rightarrow$ g+tt $\rightarrow$ g+tg+ $\rightarrow$ g+g+g+

###### Acknowledgements.

The authors wish to acknowledge helpful discussions with Attila Szabo. The
research was supported by the Intramural Research Program of the NIH, NHLBI.
Computational resources and services used in this work were provided by the
LoBoS cluster of the National Institutes of Health.

## References

  * Elber (2016) R. Elber, The Journal of Chemical Physics 144, 060901 (2016).
  * Dellago  _et al._ (1998) C. Dellago, P. G. Bolhuis, F. S. Csajka, and D. Chandler, Journal of Chemical Physics 108, 1964 (1998).
  * Bolhuis  _et al._ (2002) P. G. Bolhuis, D. Chandler, C. Dellago, and P. L. Geissler, Annual Review of Physical Chemistry 53, 291 (2002).
  * Czerminski and Elber (1990) R. Czerminski and R. Elber, International Journal of Quantum Chemistry 186, 167 (1990).
  * Henkelman  _et al._ (2000) G. Henkelman, B. P. Uberuaga, and H. Jónsson, Journal of Chemical Physics 113, 9901 (2000).
  * E  _et al._ (2002) W. E, W. Ren, and E. Vanden-Eijnden, Physical Review B 66, 5 (2002).
  * Gillilan and Wilson (1992) R. E. Gillilan and K. R. Wilson, Journal of Chemical Physics 97, 1757 (1992).
  * Olender and Elber (1996) R. Olender and R. Elber, The Journal of Chemical Physics 105, 9299 (1996).
  * Elber  _et al._ (1999) R. Elber, J. Meller, and R. Olender, The Journal of Physical Chemistry B 103, 899 (1999).
  * Elber and Shalloway (2000) R. Elber and D. Shalloway, The Journal of Chemical Physics 112, 5539 (2000).
  * Passerone and Parrinello (2001) D. Passerone and M. Parrinello, Physical Review Letters 87, 108302 (2001).
  * Passerone  _et al._ (2003) D. Passerone, M. Ceccarelli, and M. Parrinello, Journal of Chemical Physics 118, 2025 (2003).
  * Lee  _et al._ (2003a) I.-H. Lee, J. Lee, and S. Lee, Physical Review B 68, 064303 (2003a).
  * Lee  _et al._ (2005) I.-H. Lee, S.-Y. Kim, and J. Lee, Chemical Physics Letters 412, 307 (2005).
  * Lee  _et al._ (2012a) I.-H. Lee, S.-Y. Kim, and J. Lee, The Journal of Physical Chemistry. B 116, 6916 (2012a).
  * Lee  _et al._ (2013a) I.-H. Lee, S.-Y. Kim, and J. Lee, International Journal of Molecular Sciences 14, 16058 (2013a).
  * Crehuet and Field (2003) R. Crehuet and M. J. Field, Physical Review Letters 90, 089801; author reply 089802 (2003), arXiv:0104204v1 [arXiv:cond-mat] .
  * Onsager and Machlup (1953) L. Onsager and S. Machlup, Physical Review 91, 1505 (1953).
  * Machlup and Onsager (1953) S. Machlup and L. Onsager, Physical Review 91, 1512 (1953).
  * Eastman  _et al._ (2001) P. Eastman, N. Grønbech-Jensen, and S. Doniach, Journal of Chemical Physics 114, 3823 (2001).
  * Zuckerman and Woolf (2000) D. M. Zuckerman and T. B. Woolf, Physical Review. E, Statistical, nonlinear, and soft matter physics 63, 016702 (2000), arXiv:0005073v2 [arXiv:physics] .
  * Fujisaki  _et al._ (2010) H. Fujisaki, M. Shiga, and A. Kidera, Journal of Chemical Physics 132, 134101 (2010), arXiv:1003.0971 .
  * Fujisaki  _et al._ (2013) H. Fujisaki, M. Shiga, K. Moritsugu, and A. Kidera, The Journal of Chemical Physics 139, 054117 (2013), arXiv:1304.7342 .
  * Faccioli  _et al._ (2006) P. Faccioli, M. Sega, F. Pederiva, and H. Orland, Physical Review Letters 97, 108101 (2006).
  * Sega  _et al._ (2007) M. Sega, P. Faccioli, F. Pederiva, G. Garberoglio, and H. Orland, Physical Review Letters 99, 118102 (2007).
  * Beccara  _et al._ (2012) S. Beccara, T. Skrbic, R. Covino, and P. Faccioli, Proceedings of the National Academy of Sciences U.S.A. 109, 2330 (2012), arXiv:1111.3518v1 .
  * Lee  _et al._ (1997) J. Lee, H. Scheraga, and S. Rackovsky, Journal of Computational Chemistry 18, 1222 (1997).
  * Lee  _et al._ (1999) J. Lee, A. Liwo, and H. Scheraga, Proceedings of the National Academy of Sciences U.S.A. 96, 2025 (1999).
  * Lee  _et al._ (2003b) J. Lee, I.-H. Lee, and J. Lee, Physical Review Letters 91, 080201 (2003b), arXiv:0307690 [cond-mat] .
  * Steinbach (2004) P. J. Steinbach, Proteins: Structure, Function and Genetics 57, 665 (2004).
  * Joo  _et al._ (2007) K. Joo, J. Lee, S. Lee, J.-H. Seo, S. J. Lee, and J. Lee, Proteins 69, 83 (2007).
  * Lee  _et al._ (2008) J. Lee, K. Joo, Kim, Seung-Yeon, and J. Lee, Journal of Computational Chemistry 29, 2479 (2008).
  * Joo  _et al._ (2009) K. Joo, J. Lee, J.-H. Seo, K. Lee, B.-G. Kim, and J. Lee, Proteins: Structure, Function, and Bioinformatics 75, 1010 (2009).
  * Lee  _et al._ (2011) J. Lee, J. Lee, T. N. Sasaki, M. Sasai, C. Seok, and J. Lee, Proteins: Structure, Function, and Bioinformatics 79, 2403 (2011).
  * Joo  _et al._ (2008) K. Joo, J. Lee, I. Kim, S. J. Lee, and J. Lee, Biophysical Journal 95, 4813 (2008).
  * Lee  _et al._ (2012b) J. Lee, S. P. Gross, and J. Lee, Physical Review E 85, 056702 (2012b).
  * Lee and Lee (2013) J. Lee and J. Lee, PLOS ONE 8, e60372 (2013).
  * Lee  _et al._ (2013b) J. Lee, S. P. Gross, and J. Lee, Scientific Reports 3, 2197 (2013b).
  * Hunt and Ross (1981) K. L. C. Hunt and J. Ross, The Journal of Chemical Physics 75, 976 (1981).
  * Adib (2008) A. B. Adib, The Journal of Physical Chemistry. B 112, 5910 (2008), arXiv:0712.1255 .
  * Haas  _et al._ (2014a) K. R. Haas, H. Yang, and J. W. Chu, Journal of Physical Chemistry Letters 5, 999 (2014a).
  * Haas  _et al._ (2014b) K. R. Haas, H. Yang, and J. W. Chu, Journal of Physical Chemistry B 188, 8099 (2014b).
  * Miller and Predescu (2007) T. F. Miller and C. Predescu, Journal of Chemical Physics 126, 144102 (2007).
  * Alt and Godau (1995) H. Alt and M. Godau, International Journal of Computational Geometry & Applications 05, 75 (1995).
  * Neria  _et al._ (1996) E. Neria, S. Fischer, and M. Karplus, The Journal of Chemical Physics 105, 1902 (1996).
  * Loncharich  _et al._ (1992) R. J. Loncharich, B. R. Brooks, and R. W. Pastor, Biopolymers 32, 523 (1992).
  * Jang  _et al._ (2002) S. Jang, S. Shin, and Y. Pak, Journal of the American Chemical Society 124, 4976 (2002).
  * Lei and Duan (2004) H. Lei and Y. Duan, Journal of Chemical Physics 121, 12104 (2004).
  * Lei  _et al._ (2006) H. Lei, S. G. Dastidar, and Y. Duan, Journal of Physical Chemicstry B 110, 22001 (2006).
  * Wu and Shea (2010) C. Wu and J. E. Shea, PLoS Computational Biology 6, e1000998 (2010).
  * Lei  _et al._ (2009) H. Lei, Z. X. Wang, C. Wu, and Y. Duan, Journal of Chemical Physics 131, 165105 (2009).
  * Haberthür and Caflisch (2008) U. Haberthür and A. Caflisch, Journal of Computational Chemistry 29, 701 (2008).
  * Feng  _et al._ (2009) J. A. Feng, J. Kao, and G. R. Marshall, Biophysical Journal 97, 2803 (2009).
  * Mathews and Case (2006) D. H. Mathews and D. A. Case, Journal of Molecular Biology 357, 1683 (2006).