# Optimal control for sampling the transition path process
and estimating rates

Jiaxin Yuan<EMAIL_ADDRESS>Department of Mathematics, University of Maryland,
College Park, MD 20742, USA Amar Shah<EMAIL_ADDRESS>College of
Letters and Science, University of California, Berkeley, CA 94709, USA
Channing Bentz<EMAIL_ADDRESS>College of Arts and Sciences, University
of Nebraska, Lincoln, NE 68508, USA  Maria Cameron<EMAIL_ADDRESS>Department
of Mathematics, University of Maryland, College Park, MD 20742, USA

###### Abstract

Many processes in nature such as conformal changes in biomolecules and
clusters of interacting particles, genetic switches, mechanical or
electromechanical oscillators with added noise, and many others are modeled
using stochastic differential equations with small white noise. The study of
rare transitions between metastable states in such systems is of great
interest and importance. The direct simulation of rare transitions is
difficult due to long waiting times. Transition path theory is a mathematical
framework for the quantitative description of rare events. Its crucial
component is the committor function, the solution to a boundary value problem
for the backward Kolmogorov equation. The key fact exploited in this work is
that the optimal controller constructed from the committor leads to the
generation of transition trajectories exclusively. We prove this fact for a
broad class of stochastic differential equations. Moreover, we demonstrate
that the committor computed for a dimensionally reduced system and then lifted
to the original phase space still allows us to construct an effective
controller and estimate the transition rate with reasonable accuracy.
Furthermore, we propose an all-the-way-through scheme for computing the
committor via neural networks, sampling the transition trajectories, and
estimating the transition rate without meshing the space. We apply the
proposed methodology to four test problems: the overdamped Langevin dynamics
with Mueller’s potential and the rugged Mueller potential in 10D, the noisy
bistable Duffing oscillator, and Lennard-Jones-7 in 2D.

Keywords: transition path theory, optimal stochastic control, transition path
process, committor, transition rate, Langevin dynamics, overdamped Langevin
dynamics, collective variables, sampling, neural network, Lennard-Jones-7,
Duffing oscillator.

##  1 Introduction

The study of rare events in stochastic systems is crucial for understanding
natural phenomena such as conformal changes in biomolecules and clusters of
interacting particles, protein folding, noise-driven transitions in nonlinear
oscillator systems, genetic switches, and many others. Often rare events in
such systems are associated with transitions between metastable states
separated by high energetic barriers. Direct simulations of rare events are
difficult due to long waiting times. Deterministic approaches based on solving
partial differential equations are hampered by the high dimensionality of
phase space or other numerical issues. In this work, we propose an approach
for sampling transition trajectories between metastable states based on
optimal control and use them to calculate transition rates. This approach is
inspired by a remarkable fact demonstrated in a recent work by Zhang, Sahai,
and Marzouk [1] that a highly effective controller for a broad class of
stochastic systems can be obtained using a rough approximation to the solution
of an appropriate partial differential equation. The theoretical foundation of
this approach is motivated by work [2] of Lu and Nolen and recent work [3] of
Gao, Li, Li, and Liu detailing the case of overdamped Langevin dynamics.

###  1.1 An overview

Numerical methods for the study of rare events can be divided into two large
classes: deterministic and stochastic.

Stochastic methods include direct simulation suitable only for the case where
the noise is relatively large [4] as well as various enhanced sampling
algorithms. These include methods aimed at sampling rare transitions e.g.
transition path sampling [5, 6, 7], milestoning [8], weighted ensemble [9],
and adaptive splitting methods [10, 11], and methods for the exploration of
the configurational space e.g. metadynamics [12, 13] and machine learning-
assisted techniques – see [14, 15] and references therein.

The class of deterministic methods can be subdivided into several categories.
First, there are methods aiming at finding maximum likelihood transition paths
via a gradient descent in the path space [16, 17, 18, 19, 20, 21] and via
using control theory and a shooting approach [22, 23]. While they are
computationally cheap and suitable for high dimensions, these methods produce
a single transition path along which the transition flux is focused as the
noise amplitude tends to zero, and are only able to give an asymptotic scaling
for the exponential factor of the transition rate in the limit of noise
coefficient approaching zero. Second, methods for computing the
_quasipotential_ on a mesh are useful for visualization of the effective
potential for two- or three-dimensional systems with nongradient drifts – see
[24, 25, 26, 27]. Third, there are methods using the framework of transition
path theory [28, 29] where the key component is the numerical solution of the
_committor problem_ , a boundary-value problem for the stationary backward
Kolmogorov equation. Novel techniques developed for accomplishing this task in
dimensions higher, and even in some cases, much higher than three, include
those based on training neural networks [30, 31, 32], crafting diffusion maps
[33, 34, 35], or representing the solution by tensor trains [36].

###  1.2 Optimal control

The first two categories of deterministic methods, i.e. methods for finding
maximum likelihood transition paths in the vanishing noise limit and methods
for computing the quasipotential, are connected via the _deterministic optimal
control_. For example, consider a system evolving according to an SDE of the
form

$dX_{t}=b(X_{t})dt+\sigma(X_{t})\sqrt{\epsilon}dW_{t},$ (1)

where the drift field $b$ is smooth and has a finite number of attractors
lying within a ball of a finite radius around the origin, the matrix function
$\sigma$ is smooth, and $\epsilon$ is a small parameter. Let $\mathcal{A}$ be
an attractor of the corresponding ODE $\dot{x}=b(x)$. The escape problem from
the attractor $\mathcal{A}$ can be viewed as an optimal control problem where
_an optimal realization of the Brownian motion $W_{t}$ driving the process out
of the basin of $\mathcal{A}$ is sought_. Therefore, the controlled ODE with a
controller $u$ is

$\dot{x}=b(x)+\sigma(x)u,\quad x(0)\in\mathcal{A}.$ (2)

The cost functional is derived in the large deviations theory [37, 38]. If
$\sigma$ is nonsingular everywhere, the cost functional defined for all
absolutely continuous paths $\phi(\cdot)$ is given by

$S_{T}(\phi)=\frac{1}{2}\int_{0}^{T}\|u(\phi(t))\|^{2}dt\equiv\frac{1}{2}\int_{0}^{T}\|\sigma^{-1}(\phi)[\dot{\phi}-b(\phi)]\|^{2}dt.$
(3)

The last expression is exactly the Freidlin-Wentzell action functional for SDE
(1) [38]. The optimal controller $u$ is given by

$u=\sigma^{\top}\nabla U,\quad{\rm where}\quad
U(x)=\inf_{\phi,T}\left\\{S_{t}(\phi)~{}|~{}\phi(0)\in\mathcal{A},\phi(T)=x\right\\}$
(4)

is the _quasipotential_. The infimum in (4) taken over all paths and all final
times is always achieved at $T=\infty$ since $\phi(0)\in\mathcal{A}$ [38].
Plugging this optimal controller to ODE (2) results in the ODE governing the
optimal escape path from the attractor $\mathcal{A}$ [39, 40]

$\dot{x}=b(x)+\sigma(x)\sigma(x)^{\top}\nabla U(x).$ (5)

In practice, if the quasipotential is found, one can find the optimal escape
path from the basin of $\mathcal{A}$ by integrating ODE (5) backward in time
starting at the point at the boundary of the basin of $\mathcal{A}$ where the
quasipotential $U$ is minimal. Therefore, _the quasipotential determines the
optimal controller for finding the most probable escape path in the zero noise
limit_ , and the escape path is governed by ODE (5).

A similar connection via _stochastic optimal control_ exists between
transition path sampling and methods relying on solving the committor problem.
These two approaches address the case where the noise is small but finite. The
study of stochastic optimal control problems in the context of molecular
dynamics applications was started by Hartmann and collaborators in 2012 [41].
In contrast to seeking an optimal realization of the Brownian motion as in the
deterministic optimal control problem outlined above, the setting of the
stochastic optimal control problem leaves the noise term unchanged. Instead,
it aims at _finding a minimal modification to the drift term of the governing
SDE that would make all trajectories accomplish the desired transition_. For
example, suppose a system is governed by the overdamped Langevin dynamics

$dX_{t}=-\nabla V(X_{t})dt+\sqrt{2\beta^{-1}}dW_{t},$ (6)

where $V(x)$ is a smooth and coercive potential with a finite number of
isolated minima and $\beta^{-1}$ is a small parameter often interpreted in
chemical physics applications as the temperature. Let $A$ and $B$ be open
disjoint sets surrounding local minima $x_{A}$ and $x_{B}$ of $V(x)$, and let
$q(x)$ be the committor function, i.e., the probability that the process
governed by (6) and starting at $x$ will first reach $B$ rather than $A$. The
_committor function determines the optimal controller for stochastic
dynamics_. Precisely, the dynamics of transition paths from $A$ to $B$ are
governed by [2]

$dX_{t}=\left[-\nabla V(X_{t})+{2\beta^{-1}}\nabla\log
q(x)\right]dt+\sqrt{2\beta^{-1}}dW_{t},\quad X_{0}\in\partial A.$ (7)

The function $2\beta^{-1}\nabla\log q(x)$ is the optimal controller with
respect to the cost functional [42, 3]

$C(u)=\mathbb{E}_{P}\left[\frac{1}{2}\int_{0}^{\tau_{AB}}\|u(X_{t})\|^{2}dt+g(X_{\tau_{AB}})\right],\quad{\rm
where}\quad g(x)=\begin{cases}+\infty,&x\in\partial A\\\ 0,&x\in\partial
B\end{cases}$ (8)

is the exit cost, $\tau_{AB}=\inf\\{t>0~{}|~{}X_{t}\in\bar{A}\cup\bar{B}\\}$
is the stopping time, and $P$ is probability measure on the path space of SDE
(6).

###  1.3 Applications

In this work, we are especially interested in applications coming from
chemical physics and mechanical engineering.

In chemical physics, molecular motion is often modeled by Langevin dynamics in
$\mathbb{R}^{6N}$ where $N$ is the number of atoms. To alleviate the problem
of high dimensionality and make results more interpretable, physically
motivated collective variables are often introduced. Collective variables
(CVs) are functions of atomic coordinates effectively capturing the main
dynamical modes. Common choices of CVs are dihedral angles along the backbone
of a studied biomolecule, interatomic distances between particular key atoms,
etc. The dynamics in collective variables $Z_{t}\in\mathbb{R}^{d}$, $d\ll 6N$,
are often modeled by the _overdamped Langevin dynamics in collective
variables_

${dZ_{t}=\left[-M(Z_{t})\nabla F(Z_{t})+\beta^{-1}\nabla\cdot
M(Z_{t})\right]dt+\sqrt{2\beta^{-1}}M^{1/2}(Z_{t})dW_{t}}$ (9)

where $M(Z_{t})$ and $F(Z_{t})$ are the diffusion tensor and the free energy
found by standard techniques using molecular dynamics (MD) data [19]. The
computation of $M$ and $F$ is detailed in Appendix A of [34]. We note that the
dynamics of SDE (9) do not necessarily accurately represent the dynamics of
the CVs evaluated along the trajectory governed by the original SDE even if
the original SDE is merely the overdamped Langevin dynamics (6) and with only
one collective variable [43]. It is shown in [43] that _the level sets of the
collective variable should be normal to the manifold_ along which the dynamics
of the original system are focused for an accurate estimate of residence time
near metastable states using the reduced system, i.e., SDE (9) with $X_{t}$
being a scalar function.

Among mechanical engineering models, we are interested in nonlinear
oscillators with small added noise:

$\begin{cases}dX_{t}&=m^{-1}P_{t}dt\\\ dP_{t}&=\left[-\gamma P_{t}-\nabla
V(X_{t})\right]dt+\sqrt{2\gamma m\epsilon}dW_{t},\end{cases}$ (10)

where $\epsilon$ is the parameter regulating the noise amplitude, $\gamma$ is
the friction coefficient, and $V$ is the potential energy function.

Note that SDEs (9) and (10) are of the form of SDE (1) with a matrix function
$\sigma(\cdot)$ being ${d\times r}$ where $d$ is the dimension of the phase
space and $r\leq d$ is the rank of $\sigma$, i.e., the matrix function
$\sigma$ has linearly independent columns. In this case, $W_{t}$ is the
standard $r$-dimensional Brownian motion.

###  1.4 Goals and summary of main results

The goal of this work is two-fold. The first objective is to establish the
connection between the committor and the optimal control for a broad class of
SDEs. The second objective is to develop a methodology based on optimal
control and transition path theory for sampling transition paths from a
metastable region $A$ to a metastable region $B$ and finding the transition
rate from $A$ to $B$.

We develop a workflow that allows one to generate the transition trajectories
and compute the transition rate without meshing the ambient space. Our Python
codes implementing this workflow are posted on GitHub [44, 45].

1\. Theoretical result: the solution to the optimal control problem. We have
proven a theorem (Theorem 3.1) that established the relationship between the
committor and the optimal control for SDE (1) with $\sigma(\cdot)$ being
${d\times r}$, ${\sf rank}(\sigma)=r$. The optimally controlled dynamics are
found to be of the form

$dX_{t}=\left[b(X_{t})+\sigma\sigma^{\top}v^{*}(X_{t})\right]dt+\sigma(X_{t})dW_{t},$
(11)

where the optimal control $\sigma^{\top}v^{*}$ satisfies

$\sigma^{\top}v^{*}=\sigma^{\top}\nabla\log q^{+}$ (12)

where $q^{+}$ is the forward committor (17). This is a generalization of
Theorem 3.3 in [3] and is related to the results in [2] and [1].

2\. Compute the committor. In MD applications, it is challenging to compute
the committor $q(x)$ accurately for the dynamics of interest due to high
dimensionality. Therefore, we compute it for the reduced dynamics in CVs (9)
and then lift it to the original phase space assuming that the dynamics in it
is overdamped Langevin (6). It is also difficult to obtain an accurate
solution to the committor problem for SDE (10) due to the degeneracy of the
elliptic PDE. Therefore, in both cases, we expect to have an approximate
solution to the committor problem.

Neural network-based solvers for the committor problem have several
advantages. First, they find a globally defined smooth solution function whose
gradient needed for the controller $v^{*}$ is readily accessible via automatic
differentiation. Second, they do not require artificial boundary conditions on
the outer boundary of the computational domain unlike finite difference and
finite element methods. Finally, they do not require meshing the space which
makes them more amenable to promotion to higher dimensions. Our neural
network-based solver for the committor problem for SDE (9) is similar to the
one by Li, Lin, and Ren [31] that exploits the variational formulation and
sets up a solution model that automatically satisfies the boundary conditions.
The committor problem for SDE (10) does not admit a variational formulation.
Therefore, we use the PINNs framework by Raissi, Perdikaris, and Karniadakis
[46].

3\. Sample transition trajectories. The transition path theory framework
allows one to compute the transition rate once the committor is available.
However, if the committor is inexact e.g., due to suboptimal or insufficient
set of CVs, the transition rate determined in this way is likely to be highly
inaccurate [43, 35]. On the other hand, even a rough approximation to the
solution of the backward Kolmogorov equation yields a very good controller.
Therefore, we use the found committor to construct a controller according to
(11)–(12) and sample transition trajectories.

4\. Estimate the transition rate. The transition rate $\nu_{AB}$ is defined as
the average number of transitions from $A$ to $B$ observed per unit time. We
propose to estimate it as

$\nu_{AB}=\frac{\rho_{AB}}{\mathbb{E}[\tau_{AB}]},$ (13)

where $\rho_{AB}$ is the probability of a trajectory being reactive i.e. on
its way from $A$ to $B$, and $\mathbb{E}[\tau_{AB}]$ is the _expected
crossover time_ from $A$ to $B$ found by simulating the controlled process
(11)–(12). The probability $\rho_{AB}$ is estimated using the computed
committor.

5\. Validation. We apply the proposed methodology to four test problems: the
overdamped Langevin equation with Mueller’s potential in 2D and with the
rugged Mueller potential in 10D as in [31], a single bistable Duffing
oscillator as in [1], and Lennard-Jones cluster of 7 particles in 2D (LJ7) as
in [47, 34]. We assess the accuracy and demonstrate the efficacy of the
proposed methodology. In all test cases, the estimates of the transition rate
by formula (13) are consistent with those found by brute force even if the
committor is not very accurate as in the case of LJ7. On the other hand, the
transition rate estimated directly using the computed committor is notably
less accurate. An explanation for this phenomenon is offered.

The paper is organized as follows. Section 2 provides the necessary background
on transition path theory and the transition path process. Section 3 and 4
contains our theoretical results. Section 5 describes the numerical methods
used in this work. Section 6 presents the application of the proposed
methodology to three benchmark test problems. Section 7 summarizes the results
and gives perspectives for future work. The proof of the main theorem (Theorem
3.1) as well as a number of technical aspects are elaborated in appendices.

##  2 Background

In this section, we will provide the necessary theoretical background on
transition path theory (TPT) and the transition path process.

###  2.1 Transition Path Theory

Transition Path Theory (TPT) is a celebrated mathematical framework for the
quantitative description of rare transitions in stochastic systems [28, 29].
Suppose a system is evolving according to SDE

$dX_{t}=b(X_{t})dt+\sigma(X_{t})dW_{t},\quad
X_{t}\in\Omega\subseteq\mathbb{R}^{d},\quad x\in\Omega\subset\mathbb{R}^{d}.$
(14)

Throughout this work, we will adopt the following assumptions about SDE (14).

######  Assumption 1.

The domain $\Omega$ is either $\mathbb{R}^{d}$ or a manifold without boundary
with metric being locally Euclidean.

For example, $\Omega$ can be a $d$-dimensional “flat” torus, $\mathbb{T}^{d}$,
or a direct product $\mathbb{T}^{k}\times\mathbb{R}^{d-k}$, $1\leq k\leq d-1$.

######  Assumption 2.

The drift field $b:\Omega\rightarrow\mathbb{R}^{d}$ is a smooth vector
function. The corresponding ODE $\dot{x}=b(x)$ has a finite number of
attractors and all its trajectories approach one of the attractors as
$t\rightarrow\infty$.

######  Assumption 3.

The matrix function $\sigma:\Omega\rightarrow\mathbb{R}^{d\times r}$, $r:={\sf
rank}(\sigma)$, is smooth. The entries of $\sigma(x)$ are bounded, and ${\sf
rank}(\sigma)$ is the same for all $x\in\Omega$. The singular values of
$\sigma$ are bounded from above and from below in $\Omega$.

######  Assumption 4.

There exists a unique invariant density $\mu(x)$, and the system is ergodic.

The infinitesimal generator $\mathcal{L}$ of the process governed by SDE (14)
is defined as

$\mathcal{L}f(x)=b\cdot\nabla f+\frac{1}{2}{\sf
tr}\left(\sigma\sigma^{T}\nabla\nabla f\right).$ (15)

The invariant density is the solution to $\mathcal{L}^{\ast}\mu=0$,
$\int_{\Omega}\mu dx=1$, where $\mathcal{L}^{\ast}$ is the adjoint to the
generator (15).

Suppose we want to study transitions between two disjoint open sets $A$ and
$B$ in $\Omega$. For example, if we are interested in transitions between
neighborhoods of two distinct attractors of the ODE $\dot{x}=b(x)$, we choose
$A$ and $B$ to be these neighborhoods. Let
$\\{X_{t}~{}|~{}-\infty<t<\infty\\}$ be an infinitely long trajectory of (14).
Due to ergodicity, it will visit $A$ and $B$ infinitely many times. TPT
studies statistics of pieces of such a trajectory that start at $\partial A$
and next hit $\partial B$ without returning to $\bar{A}$111The bar above a set
denotes its closure. in between.  Such pieces are called the _reactive
trajectories_ , and $A$ and $B$ are called the _reactant_ and _product sets_
respectively. Key concepts of TPT are _forward and backward committors_
$q^{+}(x)$ and $q^{-}(x)$. The forward committor $q^{+}(x)$ is defined as the
probability that the process starting at $x$ will first hit $\bar{B}$ rather
than $\bar{A}$. The backward committor is the probability that the process
arriving at $x$ has hit $\bar{A}$ last rather than $\bar{B}$. Specifically,
let $\tau^{+}_{D}$ ($\tau^{-}_{D}$) be the first (last) hitting time of region
$\bar{D}$ for trajectory starting (arriving) at $x$, and $\tau^{-}_{D}$ be the
last hitting time of region $\bar{D}$ for trajectory arriving at $x$.
Equivalently, $\tau^{-}_{D}$ is the first hitting time of $\bar{D}$ for the
time-reversed process of $\hat{X}_{t}=X_{-t}$:

$\begin{cases}\tau^{+}_{D}(x)=\inf\\{t\geq 0:X_{t}\in\bar{D},x(0)=x\\}\\\
\tau^{-}_{D}(x)=\inf\\{t\geq
0:\hat{X}_{t}\in\bar{D},\hat{x}(0)=x\\}.\end{cases}$ (16)

Given two disjoint regions $A$ and $B$, the forward committor function
$q^{+}(x)$ and backward committor function $q^{-}(x)$ are defined as follows

$\begin{cases}q^{+}:\Omega\rightarrow[0,1],q^{+}(x)=\mathbb{P}\\{\tau^{+}_{B}(x)<\tau^{+}_{A}(x)\\}\\\
q^{-}:\Omega\rightarrow[0,1],q^{-}(x)=\mathbb{P}\\{\tau_{A}^{-}(x)<\tau_{B}^{-}(x)\\}.\end{cases}$
(17)

The region $\Omega$ with removed sets $\bar{A}$ and $\bar{B}$ will be denoted
by $\Omega_{AB}$:

$\Omega_{AB}:=\Omega\backslash(\bar{A}\cup\bar{B}).$ (18)

The forward and backward committors $q^{+}$ and $q^{-}$ are the solutions to
the following boundary value problems (BVPs) [28]

$\begin{cases}\mathcal{L}q^{+}=0&x\in\Omega_{AB}\\\ q^{+}(x)=0&x\in\partial
A\\\ q^{+}(x)=1&x\in\partial B,\\\ \frac{\partial
q^{+}}{\partial\hat{n}}=0,&x\in\partial\Omega,\end{cases}\hfill\quad\begin{cases}\mathcal{L}^{\dagger}q^{-}=0&x\in\Omega_{AB}\\\
q^{-}(x)=1&x\in\partial A\\\ q^{-}(x)=0&x\in\partial B,\\\ \frac{\partial
q^{-}}{\partial\hat{n}}=0,&x\in\partial\Omega.\end{cases}$ (19)

We will refer to the BVP for the forward committor as the _committor problem_.
Here, $\mathcal{L}$ is the infinitesimal generator (15).
$\mathcal{L}^{\dagger}$ is the infinitesimal generator of the corresponding
time-reversed process [29]

$\mathcal{L}^{\dagger}f(x)=-b\cdot\nabla f+\frac{1}{\mu}{\sf
div}(\sigma\sigma^{T}\mu)+\frac{1}{2}{\sf
tr}\left(\sigma\sigma^{T}\nabla\nabla f\right),$ (20)

where ${\sf div}(\sigma\sigma^{T}\mu)$ is the divergence of the matrix
$\sigma\sigma^{T}\mu$, i.e., a vector with components $\left[{\sf
div}(\sigma\sigma^{T}\mu)\right]_{i}=\sum_{j}\partial_{x_{j}}(\sigma\sigma^{T}\mu)_{ji}$,
$1\leq i\leq d$. The homogeneous Neumann boundary conditions on
$\partial\Omega$ in (19) are relevant if $\Omega$ has a reflecting boundary;
$\hat{n}$ is the external unit normal to $\partial\Omega$. Note that for both
SDEs of our interest, (9) and (10), the invariant density $\mu$ is known. For
overdamped Langevin dynamics in CVs (9),

$\mu(x)=Z_{F}e^{-\beta F(x)},~{}~{}Z_{F}=\int_{\Omega}e^{-\beta
F(x)}dx\quad\text{for SDE \eqref{sde4}}$ (21)

and for Langevin dynamics (10),

$\displaystyle\mu(x,p)$ $\displaystyle=Z_{H}e^{-H(x,p)/\epsilon},\quad
H(x,p)=\frac{p^{2}}{2m}+V(x),$ (22) $\displaystyle Z_{H}$
$\displaystyle=\int_{\Omega}e^{-H(x,p)/\epsilon}dxdp\quad\text{for SDE
\eqref{sde5}}.$

The generators for SDEs (9) and (10), respectively, are

$\mathcal{L}f=\beta^{-1}e^{\beta F}\nabla\cdot\left(e^{-\beta F}M\nabla
f\right)$ (23)

and

$\mathcal{L}f=\frac{p}{m}\cdot\nabla_{x}f-\nabla_{x}V\cdot\nabla_{p}f-\gamma
p\cdot\nabla_{p}f+\gamma m\epsilon\Delta_{p}f,$ (24)

where $\nabla_{x}$ and $\nabla_{p}$ denote gradients with respect to the
coordinates and momenta respectively, and $\Delta_{p}$ is the Laplacian with
respect to momenta.

If the governing SDE is time-reversible, e.g. the overdamped Langevin dynamics
(6) and the overdamped Langevin dynamics in collective variables (9), the
backward committor is readily found from the forward committor:
$q^{-}(x)=1-q^{+}(x)$. Langevin dynamics (10) is not time-reversible, however,
there is a nice relationship between the forward and backward committors [29]:
$q^{-}(x,p)=1-q^{+}(x,-p)$.

The time-reversible SDE (9) admits a variational formulation for the committor
problem [19] that motivates the construction of the loss function for neural
network-based committor solvers [30, 31]:

$q^{+}(x)=\operatorname*{arg\,min}_{f\in\mathcal{Q}}\int_{\Omega_{AB}}\nabla
f(x)^{\top}M(x)\nabla f(x)\mu(x)dx,$ (25)

where $\mu$ is the invariant density given by (21) and $\mathcal{Q}$ is the
set of all continuously differentiable functions $f$ satisfying the boundary
conditions $f(x)=0$, for $x\in\partial A$, and $f(x)=1$, for $x\in\partial B$.

The probability density of reactive trajectories is given by

$\mu_{AB}=\mu q^{+}q^{-}.$ (26)

The integral of $\mu_{AB}$ over $\Omega_{AB}$,

$\rho_{AB}=\int_{\Omega_{AB}}\mu q^{+}q^{-}dx$ (27)

is the probability that a stochastic trajectory at a randomly picked time is
reactive, i.e., is on its way from $A$ to $B$.

The transition rate from $A$ to $B$, $\nu_{AB}$, is defined as

$\nu_{AB}=\lim_{T\rightarrow\infty}\frac{N_{AB}}{T},$ (28)

where $N_{AB}$ is the number of transitions from $A$ to $B$ observed during
the time interval $[0,T]$. The escape rate from $A$ is defined as

$k_{AB}=\lim_{T\rightarrow\infty}\frac{N_{AB}}{T_{A}}\equiv\frac{\nu_{AB}}{\rho_{A}},$
(29)

where $T_{A}$ is the total time within $[0,T]$ during which the system last
hit $A$ and $\rho_{A}$ is the probability that an infinitely long trajectory
at any randomly picked moment of time has last hit $\bar{A}$. The probability
$\rho_{A}$ is equal to

$\rho_{A}=\int_{\Omega}\mu(x)q^{-}(x)dx.$ (30)

The reactive current is a vector field such that its flux through any surface
separating $A$ and $B$ is the transition rate $\nu_{AB}$. For SDEs (9) and
(10), it is given, respectively, by

${J}_{AB}=\beta^{-1}Z^{-1}e^{-\beta F}M\nabla q^{+}$ (31)

and

${J}_{AB}=Z^{-1}_{H}e^{-H/\epsilon}q^{+}q^{-}\begin{pmatrix}p\\\
-\nabla_{x}V\end{pmatrix}+\epsilon\gamma
Z^{-1}_{H}e^{-H/\epsilon}\begin{pmatrix}0\\\ q^{-}\frac{\partial
q^{+}}{\partial p}-q^{+}\frac{\partial q^{-}}{\partial p}\end{pmatrix}.$ (32)

A natural choice of a surface separating $A$ and $B$ is an isocommittor
surface $\Sigma_{\alpha}:=\\{x\in\Omega~{}|~{}q^{+}(x)=\alpha\\}$ for any
$\alpha\in[0,1]$. The transition rate can be expressed as (Proposition 6 in
[28]):

$\nu_{AB}=\frac{1}{2}\int_{\Omega_{AB}}\mu[\nabla
q^{+}]^{\top}\sigma\sigma^{\top}\nabla
q^{+}dx=\frac{1}{2}\int_{\Omega_{AB}}\mu[\nabla
q^{-}]^{\top}\sigma\sigma^{\top}\nabla q^{-}dx.$ (33)

In particular, for SDEs (9) and (10), the transition rates, respectively, are

$\nu_{AB}=\beta^{-1}Z_{F}^{-1}\int_{\Omega_{AB}}(\nabla q^{+})^{\top}M\nabla
q^{+}e^{-\beta F}dx$ (34)

and

$\nu_{AB}=\epsilon\gamma
Z^{-1}_{H}\int_{\Omega_{AB}}\left[\nabla_{p}q^{+}\right]^{\top}m\left[\nabla_{p}q^{+}\right]e^{-H/\epsilon}dxdp,$
(35)

where $m$ is the diagonal mass matrix.

###  2.2 Transition path process

TPT received an interesting development in the paper “Reactive trajectories
and the transition path process” by Lu and Nolen (2015) [2]. It is proven in
[2] that the dynamics of the reactive trajectories for SDE (14), with matrix
$\sigma\sigma^{\top}$ satisfying the strong ellipticity condition, are
governed by the SDE

$d{X}_{t}=\left[b({X}_{t})+\sigma({X}_{t})\sigma^{\top}({X}_{t})\nabla\log
q^{+}({X}_{t})\right]dt+\sigma({X}_{t})dW_{t},$ (36)

and this expression was suggested by the Doob $h$-transform. The process
governed by (36) is called _the transition path process_. Though we do not
assume strong ellipticity condition for $\sigma\sigma^{\top}$, most results in
[2] can be extended to the processes considered in this paper i.e., evolving
according to SDE (14) with Assumptions 1–4.

The equilibrium unnormalized density of reactive trajectories on $\partial A$
is

$\eta_{A}=\mu(x)\mathcal{L}q^{+}(x)=-\mu(x)\hat{n}(x)^{\top}\sigma(x)\sigma^{\top}(x)\nabla
q^{+}(x),\quad x\in\partial A,$ (37)

where $\hat{n}(x)$ is the unit normal pointing inside $A$. Note that while
$\mathcal{L}q^{+}=0$ in
$\Omega_{AB}\equiv\Omega\backslash(\bar{A}\cup\bar{B})$, it does not need to
be zero on $\partial A$. Equation (37) should be understood as the rate at
which transition paths exit $A$. The equilibrium unnormalized density of
reactive trajectories on $\partial B$ is

$\eta_{B}=-\mu(x)\hat{n}(x)^{\top}\sigma(x)\sigma^{\top}(x)\nabla
q^{-}(x),\quad x\in\partial B,$ (38)

where $\hat{n}(x)$ is the unit normal pointing inside $B$.

The expected crossover time $\mathbb{E}[\tau_{AB}]$ of the reactive
trajectories is defined as

$\mathbb{E}[\tau_{AB}]=\lim_{N_{AB}\rightarrow\infty}\frac{1}{N_{AB}}\sum_{k=0}^{N_{AB}-1}\left(\tau_{B,k}^{+}-\tau_{A,k}^{-}\right),$
(39)

where $N_{AB}$ is the number of transitions from $A$ to $B$ and
$\tau_{B,k}^{+}$ and $\tau_{A,k}^{-}$ are $k$th entrance time to $B$ and $k$th
exit time from $A$, respectively, registered for a long trajectory of SDE (14)
initiated at $x\in A$. The _expected crossover time_ is the ratio of the
probability $\rho_{AB}$ (27) and the transition rate [2] (also see A)

$\mathbb{E}[\tau_{AB}]=\frac{\rho_{AB}}{\nu_{AB}}.$ (40)

This relationship is very important for this work as it will be used for
restoring the transition rate using sampled transition trajectories (see
Section 4).

The transition path process governed by SDE (36) can be thought of as a
process in which the trajectories are killed as they reach $\partial B$ and
reintroduced at $\partial A$ after a waiting time at the rate $\nu_{AB}$.
Moreover, the trajectories entering $\Omega_{AB}$ through the boundary of $A$
are distributed according to (37).

The probability density of reactive trajectories $\mu_{AB}$ is an invariant
measure of the transition path process [28, 2]. Indeed, the backward and
forward Kolmogorov operators for (36) are, respectively,

$\displaystyle\mathcal{L}_{AB}f$
$\displaystyle=\left[b+\sigma\sigma^{\top}\nabla\log q^{+}\right]\cdot\nabla
f+\frac{1}{2}{\sf tr}\left(\sigma\sigma^{\top}\nabla\nabla f\right)\quad{\rm
and}$ (41) $\displaystyle\mathcal{L}^{*}_{AB}f$
$\displaystyle=-\nabla\cdot\left[(b+\sigma\sigma^{\top}\nabla\log
q^{+})f-\frac{1}{2}{\sf div}\left(\sigma\sigma^{\top}f\right)\right].$ (42)

One can check that

$\displaystyle\mathcal{L}^{*}_{AB}\mu_{AB}$
$\displaystyle=-\nabla\cdot\left[(b+\sigma\sigma^{\top}\nabla\log q^{+})\mu
q^{+}q^{-}-\frac{1}{2}{\sf div}\left(\sigma\sigma^{\top}\mu
q^{+}q^{-}\right)\right]$
$\displaystyle=-\nabla\cdot\left[\left(b\mu-\frac{1}{2}{\sf
div}\left(\sigma\sigma^{\top}\mu\right)\right)q^{+}q^{-}+\frac{1}{2}\mu\sigma\sigma^{\top}\left(q^{-}\nabla
q^{+}-q^{+}\nabla q^{-}\right)\right]$ (43) $\displaystyle=-\nabla\cdot
J_{AB}=0.$ (44)

The symbol $J_{AB}$ is the reactive current [29]. A calculation showing that
$-\nabla\cdot J_{AB}=0$ in $\Omega_{AB}$ is detailed in B. In order to make
the transition path process an equilibrium process we assume that the
transition trajectories that have reached $\partial B$ are transported back to
$\partial A$ as shown in Fig. 1. A similar construction was used for Markov
jump processes in [48].

Figure 1: An illustration for the transition path process. Transition
trajectories appear at the boundary of the region $A$, travel to the region
$B$ without returning to $A$, and disappear at the boundary of $B$. To
complete the dynamics of transition trajectories to an equilibrium process,
the trajectories absorbed at $\partial B$ must be transported back to
$\partial A$. The pipe in the figure symbolizes this reverse transport.

The probability to find the system at $s_{AB}$ at a randomly picked moment of
time is $1-\rho_{AB}$. As a result, the invariant measure of the transition
path process becomes

$\mu_{AB}=\begin{cases}\mu q^{+}q^{-},&x\in\Omega_{AB}\\\
1-\rho_{AB},&x=s_{AB}\end{cases}.$ (45)

It follows from (43) and (19) that the reactive current $J_{AB}$ on the
boundaries $\partial A$ and $\partial B$ of $s_{AB}$ is [2]

$J_{AB}=\begin{cases}\frac{1}{2}\mu\sigma\sigma^{\top}\nabla
q^{+},&x\in\partial A\\\ -\frac{1}{2}\mu\sigma\sigma^{\top}\nabla
q^{-},&x\in\partial B\end{cases}.$ (46)

Since $\nabla\cdot J_{AB}=0$ in $\Omega_{AB}$ we have

$\int_{\partial A}J_{AB}\cdot\hat{n}ds+\int_{\partial
B}J_{AB}\cdot\hat{n}ds=0$

where $\hat{n}$ is the outer unit normal and $ds$ is the surface element.
Equation (46) is consistent with (37) and (38).

##  3 Optimally controlled dynamics

In this section, we show that SDE (36) governing the transition path process
can be obtained as the solution to a stochastic optimal control problem.
Precisely, we generalize Theorem 3.3 proven in [3] for the case of the
overdamped Langevin dynamics.

###  3.1 The general Ito SDE

We consider the dynamics governed by the general Ito SDE (14) where the drift
field $b$ and the diffusion matrix $\sigma$ are as described at the beginning
of Section 2. The controlled dynamics are set to be of the form

$dY_{t}=\left[b(Y_{t})+\sigma(Y_{t})\sigma^{\top}(Y_{t})v(Y_{t})\right]dt+\sigma(Y_{t})dW_{t}$
(47)

where $v:[0,\infty)\rightarrow\Omega\subset\mathbb{R}^{d}$ is to be chosen in
an optimal manner. The form of the modification to the drift,
$\sigma\sigma^{\top}v$, is borrowed from [1] and is motivated as follows. We
are seeking a drift term that will make all trajectories of the resulting
process obey the same statistics as the transition trajectories of the
original process (14). This means that all trajectories of the controlled
process can be observed in the original process at particular noise
realizations. Hence, the span of the modification to the drift at every point
should match the span of the noise term which is the column space of the
matrix-valued function $\sigma(\cdot)$. That’s why there is the factor
$\sigma$ on the left of $\sigma\sigma^{\top}v$. The factor $\sigma^{\top}$ is
used for convenience so that the control $v$ is of the same dimension as the
process $X_{t}$.

We have chosen the cost functional to be

$C_{x}[v(\cdot)]=\mathbb{E}_{P}\left[\frac{1}{2}\int_{0}^{\tau_{AB}}\|\sigma^{\top}(Y_{s})v(Y_{s})\|^{2}ds+g(Y_{\tau})~{}|~{}X_{0}=x\right],$
(48)

where $P$ is the probability measure on the path space of SDE (14),

$\tau_{AB}=\inf\\{t>0~{}|~{}Y_{t}\in(\bar{A}\cup\bar{B})\\}$ (49)

is the stopping time, and $g$ is the exit cost defined by

$g(x)=\begin{cases}+\infty,&x\in\overline{A}\\\
0,&x\in\overline{B}\end{cases}.$ (50)

Cost functional (48) gives finite cost only if the trajectory leaves
$\Omega_{AB}$ via the boundary of $B$. The function $\sigma^{\top}v$ is the
standard form in which the optimal control is sought [49]. The optimal control
problem is to find the function $v(\cdot)$ that minimizes the cost functional
(48). Its solution is given in the following theorem.

######  Theorem 3.1.

Let $X_{t}$ be a process governed by SDE (14) satisfying Assumptions 1–4 and
let $Y_{t}$ be the corresponding controlled process governed by SDE (47). In
addition, we assume that $\Omega$ is compact with a reflecting boundary and
$\partial A$ and $\partial B$ are smooth. The infimum of the cost functional
(48) is given by

$c^{*}(x):=\inf_{v\in\mathcal{V}}C_{x}[v(\cdot)]=-\log q^{+}(x),$ (51)

where $\mathcal{V}$ is the set of admissible controls

$\mathcal{V}:=\overline{\left\\{v\in
C^{1}(\Omega_{AB})~{}|~{}\mathbb{E}_{P}\left[\exp\left(\int_{0}^{\tau{AB}}\frac{1}{2}\|\sigma^{T}(Y_{s})v(Y_{s})\|^{2}ds\right)\right]<\infty\right\\}},$
(52)

$P$ is the measure on the path space of (14), and $q^{+}$ is the forward
committor for SDE (14). The corresponding optimal control $v^{*}$ satisfies

$\sigma^{\top}v^{*}=\sigma^{\top}\nabla\log q^{+}.$ (53)

The proof of Theorem 3.1 combines ideas from Gao et al. ([3], the proof of
Theorem 3.3) and L. C. Evans’s notes on the control theory [50]. It is found
in C.

######  Remark 1.

Theorem 3.1 shows that the optimally controlled dynamics are governed by SDE
(36) for the transition path process. Therefore, all facts about the
transition path process stated in Section (2.2) are valid for the optimally
controlled dynamics.

######  Remark 2.

The requirement that the domain $\bar{\Omega}$ is compact is often implemented
in numerical simulations: the computational domain is always bounded.
Deterministic techniques always use meshes or point clouds of finite size. In
stochastic simulations, particles are often put in a box. Therefore, the
assumption that $\Omega$ is compact is not practically restrictive.

######  Remark 3.

Theorem 3.1 suggests a new approach for finding the committor via minimization
of the cost functional (51). We leave the investigation into the viability of
this approach for the future.

###  3.2 Overdamped Langevin equation in collective variables

Theorem 3.1 has an immediate application to a practical scenario: overdamped
Langevin equation in collective variables (9). The diffusion matrix $M(x)$ in
(9) is symmetric positive definite everywhere in $\Omega$. Applying Theorem
3.1 to SDE (9) results in the following corollary.

######  Corollary 3.1.1.

For the overdamped Langevin equation in collective variables (9) the optimal
controller is

$v^{*}(x)=\nabla\log q^{+}(x),$ (54)

and the controlled process is

$\displaystyle d{X}_{t}$ $\displaystyle=[-M({X}_{t})\\{\nabla
F({X}_{t})-2\beta^{-1}\nabla\log q^{+}\\}+\beta^{-1}\nabla\cdot M({X}_{t})]dt$
$\displaystyle+\sqrt{2\beta^{-1}}M({X}_{t})^{\frac{1}{2}}dW_{t}.$ (55)

###  3.3 Full Langevin equations

The Langevin dynamics (10) can be written as follows

$d\begin{bmatrix}X_{t}\\\ P_{t}\end{bmatrix}=\begin{bmatrix}m^{-1}P_{t}\\\
-(\nabla U(X_{t})+\gamma P_{t})\end{bmatrix}dt+\sqrt{2\gamma\epsilon
m}\begin{pmatrix}0\\\ I\end{pmatrix}dW_{t}$ (56)

where the diffusion matrix is

$\sigma(X_{t})=\sqrt{2\gamma\epsilon m}\begin{pmatrix}0\\\
I\end{pmatrix}\in\mathbb{R}^{2d\times d}.$ (57)

The application of Theorem 3.1 to SDE (56) gives the following controlled
process.

######  Corollary 3.1.2.

For the Langevin dynamics (56), the optimal controller $v^{*}$ can be chosen
to be

$v^{*}(x,p)=\nabla_{p}\log q^{+}(x,p)$ (58)

and the corresponding controlled process is

$\begin{cases}d{X}_{t}=m^{-1}{P}_{t}dt\\\ d{P}_{t}=\left(-\nabla
U({X}_{t})-\gamma{P}_{t}+2\gamma\epsilon m\nabla_{p}\log
q\right)dt+\sqrt{2\gamma\epsilon m}dW_{t}.\end{cases}$ (59)

The form of the optimal controller $v^{*}$ follows from (53) and (57). Note
that (53) and (57) do not define $v^{*}$ uniquely as $v^{*}$ can have
arbitrary first $d$ components corresponding to the coordinate subspace.
However, the control in SDE (59) is of the form $\sigma\sigma^{\top}v^{*}$
where

$\sigma\sigma^{\top}=2\gamma\beta^{-1}m\left[\begin{array}[]{cc}0&0\\\
0&1\end{array}\right].$

Hence the components of $v^{*}$ in the coordinate subspace are eliminated in
SDE (59).

##  4 Estimation of the transition rate

The problem of estimating the transition rate between metastable states has
been one of the central problems addressed by chemists, physicists, and
mathematicians working on quantifying rare events and remains a subject of
active research [51, 52]. Practical methods for finding transition rates can
be roughly divided into two categories, splitting and reweighting.

Splitting methods, e.g. transition interface sampling [53] and forward flux
sampling [54], stratify $\Omega_{AB}$ using level sets of a reaction
coordinate. The level sets are denoted by $\lambda_{i}$, $i=0,\ldots,n$, where
$\lambda_{0}=\partial A$ and $\lambda_{n}=\partial B$. Next, the transition
probabilities $\mathbb{P}(\lambda_{i+1}|\lambda_{i})$ to reach $\lambda_{i+1}$
starting from $\lambda_{i}$ before returning to $A$ are estimated. Then the
escape rate from $A$ to $B$ is calculated according to the formula

$k_{AB}=\frac{\nu_{A,\lambda_{1}}}{\rho_{A}}\prod_{i=1}^{n-1}\mathbb{P}(\lambda_{i+1}|\lambda_{i}),$
(60)

where $\nu_{A,\lambda_{1}}$ is the average number of transitions from $A$ to
$\lambda_{1}$ per unit time and $\rho_{A}$ is the probability that the system
has last visited $A$ rather than $B$. It was shown in Ref. [47] that these
methods can suffer from an unfortunate choice of the reaction coordinate.

Reweighting methods use enhanced sampling techniques and restore the
statistics of unbiased transition paths using an appropriate reweighting
scheme. These include weighted ensemble [9] and the Girsanov reweighting [55].
The Girsanov reweighting was used in [1] in combination with optimal control
with a fixed stopping time.

The settings in which we need to determine the transition rate are different
than those in the works mentioned above. We plan to compute the committor for
the system under consideration or for its reduced model. This means that we
can compute the transition rate using (33) provided that the invariant density
$\mu$ is known.

However, the rate computed in this manner may be inaccurate due to

  * •

an inaccurate estimate of the normalization constant for the invariant density
(see Section 6.1) and/or

  * •

a suboptimal choice set of collective variables when model reduction is used.

The issue with the normalization constant for the invariant density can be
eliminated if the escape rate from the set $A$ is computed instead – see
equation (63) below. The problem with the choice of collective variables can
be hard to overcome if the system is complicated. If the underlying dynamics
are time-reversible, i.e. given by (9), and $\xi(x)$ is the set of collective
variables then the transition rate $\tilde{\nu}_{AB}$ estimated in the space
of collective variables $\xi$ is always exaggerated and relates to the
original transition rate $\nu_{AB}$ via (Proposition 6 in Zhang, Hartmann,
Schuette (2016)[56])

$\nu_{AB}\leq\tilde{\nu}_{AB}=\nu_{AB}+\frac{1}{\beta}\int_{\Omega_{AB}}\nabla[q(x)-\tilde{q}(\xi(x))]^{\top}M(x)\nabla[q(x)-\tilde{q}(\xi(x))]\mu(x)dx,$
(61)

where $q$ and $\tilde{q}$ are the committors computed for the original and
reduced systems respectively. Equation (61) shows that if the set of the
collective variables $\xi$ were perfect, i.e., if $q(x)=\tilde{q}(\xi(x))$,
then $\nu_{AB}=\tilde{\nu}_{AB}$. In particular, this means that the lowering
dimensionality per se does not lead to an error in the transition rate.
Otherwise, there will be a model reduction error in the transition rate.

We propose the following scheme utilizing the fact that the controlled SDE
(47) with the optimal controller (51) exactly matches SDE (36) that governs
the transition trajectories of SDE (14). In particular, this means that the
expected crossover time $\mathbb{E}[\tau_{AB}]$ for trajectories of SDE (47)
with (51) is the same as that for the transition trajectories of (14).
Therefore, first one needs to generate a set of transition trajectories using
the optimally controlled dynamics (47) with (51). This allows us to compute
the expected crossover time $\mathbb{E}[\tau_{AB}]$ and find the transition
rate $\nu_{AB}$:

$\nu_{AB}=\frac{\rho_{AB}}{\mathbb{E}[\tau_{AB}]}.$ (62)

The expected escape time from $A$ can be readily found as well using (29) and
(62):

$\mathbb{E}[\tau_{A}]=\mathbb{E}[\tau_{AB}]\frac{\rho_{A}}{\rho_{AB}}.$ (63)

The probabilities $\rho_{A}$ and $\rho_{AB}$ will be estimated using the
computed committors and formulas (27) and (30). The error in their estimates
due to model reduction via imperfect collective variables is less impaired
than the error in the rate. The reason is that the formula for the rate uses
the gradient of the committor while the formulas for $\rho_{A}$ and
$\rho_{AB}$ involve only the committors. The error in the gradient of the
committor is amplified due to the differentiation. This effect can be
illustrated using an example from [43] found in D. Moreover, the estimate of
the expected crossover time $\mathbb{E}[\tau_{AB}]$ from the controlled
dynamics (47) remains reasonably accurate even if the estimate of the
committor is rough. This issue is explored in E.

##  5 Numerical solution to the committor problem

This section offers descriptions of numerical methods that we have used for
finding the committors for the test problems reported in Section 6. As
mentioned in Section 1.4, neural network-based (NN-based) solvers have several
advantages. First, they yield a globally defined smooth solution whose
derivative is readily available due to automatic differentiation. Second, they
do not require artificial boundary conditions on the outer boundary of
$\Omega_{AB}$. Finally, they do not require meshing the space which makes them
more amenable for promotion to higher dimensions. The finite element method
(FEM) is used for validation of the committors computed using NN-based
methods. Its implementation for the time-reversible dynamics (9) and for the
Langevin dynamics is detailed in F.

###  5.1 NN-solver based on the variational form of the committor problem

Two NN-based solvers for the committor problem (19) for the overdamped
Langevin dynamics (6) based on the variational formulation (25) were proposed
by Khoo, Lu, and Ying (2018) [30] and Li, Lin, and Ren (2019) [31]. These
solvers both use the loss function motivated by the minimizing property of the
committor (25) and require only the first derivatives of the committor. They
use different solution models for the committor and enforce the boundary
conditions at $\partial A$ and $\partial B$ in different ways. In [30], the
solution model involves Green’s function for Laplace’s equation and the
boundary conditions are implemented via penalty terms in the loss function. In
[31], the solution model is designed similarly to the first NN-based PDE
solver (Lagaris et al. (1998) [57]) so that it automatically satisfies the
boundary conditions. In this work, we chose to use the NN-based solver by Li,
Lin, and Ren [31] as it is simpler and its extension to the committor problem
for the overdamped Langevin dynamics in collective variables is
straightforward.

Since the overdamped Langevin dynamics (6) and the overdamped Langevin
dynamics in collective variables (9) are time-reversible, the forward and
backward committor are related via $q^{-}=1-q^{+}$. Therefore, for brevity, we
use the notation $q$ for the forward committor whenever the underlying
dynamics is time-reversible.

Following [31], we use the following solution model to the committor problem
(19)

$q(x;\theta)=(1-\chi_{A}(x))[(1-\chi_{B}(x))\mathcal{N}(x;\theta)+\chi_{B}(x)],\quad
x\in\Omega_{AB},$ (64)

where $\mathcal{N}(x;\theta)$ is the output of a fully connected neural
network (NN) and $\chi_{A}(x)$ and $\chi_{B}(x)$ are smooth approximations to
the indicator functions of $\partial A$ and $\partial B$. In this work, we use
fully connected neural networks with $L$ hidden layers with ${\sf tanh}$
activation functions and the outer layer with the sigmoid function ${\sf
sigmoid}(x)=(1+e^{-x})^{-1}$. For example, for $L=1$ and $L=2$, the neural
networks are

$\begin{cases}\mathcal{N}(x;\theta)={\sf sigmoid}\left[W_{2}{\sf
tanh}(W_{1}x+b_{1})+b_{2}\right],&L=1,\\\ \mathcal{N}(x;\theta)={\sf
sigmoid}\left(W_{3}{\sf tanh}\left[W_{2}{\sf
tanh}(W_{1}x+b_{1})+b_{2}\right]+b_{3}\right),&L=2.\\\ \end{cases}$ (65)

The argument $\theta$ comprises all entries of the matrices $W_{j}$ and the
shift vectors $b_{j}$.

The loss function is derived from the minimizing property of the committor
(25):

$q(x)=\arg\min\int_{\Omega_{AB}}\nabla f(x)^{\top}M(x)\nabla f(x)e^{-\beta
F(x)}dx,$ (66)

where the minimum is taken among all functions
$f\in\mathcal{H}^{1}(\Omega_{AB})$ (the Sobolev space) such that $f(\partial
A)=0$ and $f(\partial B)=1$. The integral in (66) is the expectation of
$\nabla f(x)^{\top}M(x)\nabla f(x)$ where $x$ is a random variable with
invariant density proportional to $e^{-\beta F(x)}$ supported in
$\bar{\Omega}_{AB}$:

$\int_{\Omega_{AB}}\nabla f(x)^{\top}M(x)\nabla f(x)e^{-\beta
F(x)}dx=\mathbb{E}_{\begin{subarray}{c}x\in\Omega_{AB}\\\ x\sim e^{-\beta
F}\end{subarray}}\left[\nabla f^{\top}M\nabla f\right].$ (67)

If we want $x$ to be distributed according to a density $\rho$, we rewrite
this expectation as

$\mathbb{E}_{\begin{subarray}{c}x\in\Omega_{AB}\\\ x\sim e^{-\beta
F}\end{subarray}}\left[\nabla f^{\top}M\nabla
f\right]=\mathbb{E}_{\begin{subarray}{c}x\in\Omega_{AB}\\\
x\sim\rho\end{subarray}}\left[\nabla f^{\top}M\nabla f\frac{e^{-\beta
F}}{\rho}\right].$ (68)

The last expectation can be approximated as a sample mean. Hence, if the
training points $x_{k}\in\Omega_{AB}$, $1\leq k\leq K$, are sampled from a
density $\rho$, then the loss function is defined by

${\sf Loss}(\theta)=\frac{1}{K}\sum_{k=1}^{K}\left[\nabla
q(x_{k};\theta)^{\top}M(x_{k})\nabla q(x_{k};\theta)\frac{e^{-\beta
F(x_{k})}}{\rho(x_{k})}\right].$ (69)

We found that it is advantageous to create the set of training points in two
stages. First, a large point cloud is generated using the metadynamics
algorithm [12] (also see [31]). Then this point cloud is rarefied into a
delta-net [58, 59], i.e., a spatially quasiuniform set of points obtained as
follows. Let $\\{x_{i}\\}_{i=1}^{N}$ be the generated point cloud and $\delta$
be the desired distance between the nearest neighbors in the training set.
Initially, we assign labels 0 to all points. We take $x_{1}$, compute the
distance from $x_{1}$ to all other points in the point cloud, assign label
keep to $x_{1}$ and labels discard to all points at distance less than
$\delta$. Then we find the point with the smallest index that has label 0,
compute the distance from it to all points with label 0, and change its label
to keep and labels of all points at distance less than $\delta$ from it to
discard. And so we continue until there are no more points with labels 0. All
points labeled as keep will form the training set. The resulting set of
training points is spatially quasiuniform. Therefore we set $\rho(x_{k})=1$ in
the loss function (69).

The loss function is minimized using the stochastic Adam optimizer [60]. Our
codes are written in Python and use the PyTorch library for the implementation
of neural networks and automatic differentiation [44].

###  5.2 PINNs for Langevin dynamics

The Langevin SDE (10) is not time-reversible and its generator (24) is not
self-adjoint. As a result, the variational formulation of the committor
problem is not available. In this case, we opt to use the NN-based solver
called the _physics-informed neural networks_ (PINNs) proposed by Raissi,
Perdikaris, and Karniadakis (2019) [46] (also see [61]). In this solver, the
loss function is the sum of the mean squared discrepancy between the left- and
right-hand side of the PDE and the mean squared error at the Dirichlet
boundary of the domain:

$\displaystyle{\sf Loss}(\theta)$
$\displaystyle=\frac{1}{K}\sum_{(x_{k},p_{k})\in\Omega_{AB}}\left|\mathcal{L}\mathcal{N}(x_{k},p_{k};\theta)\right|^{2}$
(70) $\displaystyle+\frac{1}{K_{\partial A}}\sum_{(x_{k},p_{k})\in{\partial
A}}|\mathcal{N}(x_{k},p_{k})|^{2}+\frac{1}{K_{\partial
B}}\sum_{(x_{k},p_{k})\in{\partial B}}|\mathcal{N}(x_{k},p_{k})-1|^{2}.$

Here, $K$, $K_{\partial A}$, and $K_{\partial B}$ are the numbers of training
points in $\Omega_{AB}$, $\partial A$, and $\partial B$ respectively,
$\mathcal{L}$ is the generator given by (24), and
$\mathcal{N}(x_{k},p_{k};\theta)$ is a neural network defined similar to (65).
As before, the loss function is minimized using the stochastic Adam optimizer.

##  6 Test problems

In this section, we demonstrate the effectiveness of the proposed methodology
on the following four test problems:

  1. 1.

the overdamped Langevin dynamics with Mueller’s potential, a common 2D test
problem in chemical physics,

  2. 2.

the overdamped Langevin dynamics with the rugged Mueller potential in 10D with
settings as in Ref. [31],

  3. 3.

the bistable Duffing oscillators with added white noise, and

  4. 4.

the overdamped Langevin dynamics in collective variables for the Lennard-
Jones-7 system in 2D.

The codes for these examples are available on GitHub:

  * •

neural network-based committor solvers and sampling algorithms for the
transition path process [44];

  * •

Finite element committor solvers [45].

The dynamics in test problems 1, 2, and 4 are time-reversible. Therefore, as
mentioned in Section 5.1, the forward and backward committors are related via
$q^{-}=1-q^{+}$. Hence, it suffices to compute only the forward committor in
these problems. For brevity, we will denote the forward committor simply by
$q$ in Sections 6.1, 6.2, and 6.4 containing test problems 1, 2, and 4
respectively.

###  6.1 Mueller’s potential

We first consider an overdamped Langevin equation with Mueller’s potential
(see Fig. 2)

$V(x_{1},x_{2})=\sum\limits_{i=1}^{4}D_{i}\exp\left\\{a_{i}(x_{1}-X_{i})^{2}+b_{i}(x_{1}-X_{i})(x_{2}-Y_{i})+c_{i}(x_{2}-Y_{i})^{2}\right\\}$
(71)

where

$\displaystyle[a_{1},a_{2},a_{3},a_{4}]$ $\displaystyle=[-1,-1,-6.5,0.7]$
$\displaystyle[b_{1},b_{2},b_{3},b_{4}]$ $\displaystyle=[0,0,11,0.6]$
$\displaystyle[c_{1},c_{2},c_{3},c_{4}]$ $\displaystyle=[-10,-10,-6.5,0.7]$
$\displaystyle[D_{1},D_{2},D_{3},D_{4}]$ $\displaystyle=[-200,-100,-170,15]$
$\displaystyle[X_{1},X_{2},X_{3},X_{4}]$ $\displaystyle=[1,0,-0.5,-1]$
$\displaystyle[Y_{1},Y_{2},Y_{3},Y_{4}]$ $\displaystyle=[0,0.5,1.5,1]$

####  6.1.1 Computing the committor

The two deepest minima of Mueller’s potential are located near
$a=(-0.558,1.441)$ and $b=(0.623,0.028)$. Following [31], the sets $A$ and $B$
are chosen to be the balls centered at $a$ and $b$ respectively with radius
$r=0.1$, and the smoothed indicator function functions of $\partial A$ and
$\partial B$ are defined as

$\displaystyle\chi_{A}(x)$
$\displaystyle=\frac{1}{2}-\frac{1}{2}\tanh\left[1000(\|x-a\|^{2}-(r+0.02)^{2})\right],$
$\displaystyle\chi_{B}(x)$
$\displaystyle=\frac{1}{2}-\frac{1}{2}\tanh\left[1000(\|x-b\|^{2}-(r+0.02)^{2})\right].$

The temperature is set to be $\beta^{-1}=10$ as in [31]. At this temperature,
the transitions between $A$ and $B$ are rare.

We compute the committor using FEM (see F.1) and the NN-based approach
employing the variational formulation of the committor problem (_variational
NN_ – see Section 5.1). For FEM, the domain $\Omega$ is defined as

$\Omega=\\{x\in\mathbb{R}^{2}~{}|~{}V(x)\leq 250\\}$ (72)

and triangulated using the DistMesh algorithm [62]. We also discretized
$\Omega$ using mesh2d [63] and found that the difference between the committor
computed on these two meshes was about $10^{-4}$ in the max norm. The
committor computed using FEM is displayed in Fig. 2(a).

(a) FEM

(b) Variational NN

Figure 2: The committor for the overdamped Langevin dynamics with Mueller’s
potential is computed by (a) FEM and (b) NN at $\beta^{-1}=10$. The level sets
of Mueller’s potential are superimposed. The computational domain is larger
than the shown box in both cases.

At $\beta=0.1$, sampling from the invariant Gibbs density leaves the
transition region severely underresolved. Therefore, a set of training points
for the variational NN-based solver was generated using a standard enhanced
sampling algorithm called _metadynamics_ [12] with settings used in [31].
Metadynamics is implemented as follows. The overdamped Langevin dynamics (1)
is simulated with the time step $\Delta t=10^{-5}$, and Gaussian functions of
the form

$g_{j}(x)=w\exp\left(-\frac{(x_{1}-x_{1}(t_{j}))^{2}}{2\sigma_{1}^{2}}-\frac{(x_{2}-x_{2}(t_{j}))^{2}}{2\sigma_{2}^{2}}\right)$

with height $w=5$ and $\sigma_{1}=\sigma_{2}=0.05$ are added to the potential
at times $t_{j}=500j\Delta t$, $1\leq j\leq N_{g}=2000$. Then, the overdamped
Langevin dynamics in the modified potential $V+\sum_{j}g_{j}$ is simulated
with the same time step $\Delta t=10^{-5}$, and the initial set of points is
recorded. Finally, the obtained set of points is converted into a spatially
quasi-uniform set, a _delta-net_ with $\delta=0.015$, as described in Section
5.1. The resulting training set contains a total of $N_{\sf train}=15466$
points.

The solution model is given by (64) with a neural network (65) with $L=2$
hidden layers and $N=40$ neurons in each layer. The neural network was trained
for 1000 epochs at learning rate $\eta=10^{-4}$. The resulting solution is
shown in Fig. 2(b).

To assess numerical errors in the computed committor, we use error measures
weighted by the probability density of transition trajectories: the weighted
mean absolute error (wMAE) and the weighted root mean squared error (wRMSE):

wMAE $\displaystyle=\sum_{i=1}^{N_{\sf test}}w(x^{i})\left|q_{\sf
nn}(x^{i})-q_{{\sf fem}}(x^{i})\right|$ (73) wRMSE
$\displaystyle=\sqrt{\sum_{i=1}^{N_{\sf test}}w(x^{i})\left|q_{\sf
nn}(x^{i})-q_{\text{fem}}(x^{i})\right|^{2}}$ (74)

where $q_{\sf fem}$ and $q_{\sf nn}$ are the forward committors computed by
FEM and the NN-based solver respectively, $x^{i}$, $1\leq i\leq N_{\sf test}$
are the test points, and the weights $w(x^{i})$ are defined so that they are
proportional to $\mu_{AB}$ and their sum is one:

$w(x^{i})=\frac{q_{\sf fem}(x^{i})(1-q_{\sf
fem}(x^{i}))\mu(x^{i})}{\sum_{j=1}^{N_{\sf test}}q_{\sf fem}(x^{j})(1-q_{\sf
fem}(x^{j}))\mu(x^{j})}.$ (75)

The subset of the nodes of the FEM mesh lying within the box
$[-1.5,1]\times[-0.5,2]$ was used as the test point set.

Table 1 shows the wMAE and wRMSE for the variational NN solver with $L=2$
hidden layers and $W=40$ neurons per layer.

Temperature | NN structure | wMAE | wRMSE
---|---|---|---
$\beta^{-1}=10$ |  $L=2$, $W=40$ | 2.6e-3 | 4.1e-3
Table 1: Errors wMAE and wRMSE in the forward committor for the overdamped
Langevin dynamics with Mueller’s potential at temperature $\beta^{-1}=10$
computed using the variational NN solver with a spatially quasiuniform set of
approximately $16000$ training points. The parameter $L$ indicates the number
of hidden layers (see (65)), and $W$ indicates the number of neurons per
hidden layer.

####  6.1.2 Estimation of the transition rate using the controlled process

The transition rate $\nu_{AB}$ is computed by equation (62). The expected
crossover time $\mathbb{E}[\tau_{AB}]$ is calculated by averaging crossover
times of 250 sampled transition trajectories governed by the controlled
process

$dX_{t}=\left(-\nabla V(X_{t})+2\beta^{-1}\frac{\nabla q_{\sf
nn}(X_{t})}{q_{\sf nn}(X_{t})}\right)dt+\sqrt{2\beta^{-1}}dW_{t}.$ (76)

The initial points of these trajectories are sampled according to (37) as
follows. First, $N_{\partial A}$ points $x_{j}$ equispaced on a circle of
radius $r+\delta r$ where $r$ is the radius of $A$ and $\delta r$ is a small
positive number. Then weights are assigned $w_{j}$ to these points according
to

$w_{j}=\frac{e^{-\beta V(x_{j})}|\hat{n}(x_{j})\cdot\nabla q_{\sf
nn}(x_{j})|}{\sum_{k=1}^{N_{\partial A}}e^{-\beta
V(x_{k})}|\hat{n}(x_{k})\cdot\nabla q_{\sf nn}(x_{k})|}$ (77)

where $\hat{n}(x_{j})$ is the outer unit normal to $B_{r+\delta r}(a)$ at
$x_{j}$. Then the points $x_{j}$ are sampled according to their probability
weights $w_{j}$ visualized in Fig. 3. We used $r=0.1$, $\delta r=10^{-3}$, and
$N_{\partial A}=1000$.

Figure 3: The probability distribution of the starting points of transition
trajectories on the boundary of $A$ for the overdamped Langevin dynamics in
Mueller’s potential is shown by color code. The gray curves are level sets of
Mueller’s potential.

Three samples of trajectories of (76) are displayed in Fig. 4(a). Three
trajectories of

$dX_{t}=-\nabla V(X_{t})dt+\sqrt{2\beta^{-1}}dW_{t}$ (78)

with the same realizations of the Brownian motion are shown in Fig. 4(b) for
comparison.

(a) Controlled SDE (76)

(b) Uncontrolled SDE (78)

Figure 4: Comparison of three controlled (a) and uncontrolled (b) processes
with the same initial positions and noise realizations at $\beta^{-1}=10$.
Each trajectories consist of 5000 time steps, with time step $\Delta
t=10^{-5}$. The sets $A$ and $B$ are marked by blue and orange circles of
radius $r=0.1$ respectively.

Given the FEM committor $q_{\sf fem}$ and the triangulated domain
$\Omega_{AB}$, $\rho_{AB}$ is computed directly from (27). Given $q_{\sf nn}$
and the set of training points $x_{j}$ quasiuniformly distributed in
$\Omega_{AB}$, $\rho_{AB}$ is obtained by means of Monte Carlo integration:

$\rho_{AB}=\frac{1}{N_{\sf train}}\sum_{j=1}^{N_{\sf
train}}Z_{V}^{-1}e^{-\beta V(x_{j})}q_{\sf nn}(x_{j})(1-q_{\sf nn}(x_{j})).$
(79)

The normalization constant $Z_{V}$ is computed as

$Z_{V}=\frac{1}{N_{\sf train}}\sum_{j=1}^{N_{\sf train}}e^{-\beta V(x_{j})}.$
(80)

We first generated a total of $10^{5}$ points by running metadynamics [12] and
then subsampling it into a spatially quasi-uniform delta-net with
$\delta=0.005$. The resulting set of $56322$ points is shown in Fig. 5.

Figure 5: The set of points for Mueller’s potential used for finding the
normalization constant for the invariant density and the probability
$\rho_{AB}$.

We also calculate the transition rate $\nu_{AB}$ using (33) adapted for the
overdamped Langevin dynamics (78)

$\nu_{AB}=\beta^{-1}Z_{V}^{-1}\int_{\Omega_{AB}}\|\nabla q(x)\|^{2}e^{-\beta
V(x)}dx$ (81)

with the committors $q_{\sf fem}$ and $q_{\sf nn}$. The results are presented
in Table 2. The computation of the 95% confidence interval is detailed in G.

The discrepancy between the estimates for $\rho_{AB}$ obtained using the
committors $q_{\sf nn}$ and $q_{\sf fem}$ and equations (79) and (27) is about
3%. The difference between the transition rates $\nu_{AB}$ obtained using
$q_{\sf nn}$ with Monte Carlo integration and $q_{\sf fem}$ with formula (28)
is less then 2%. The expected crossover time $\mathbb{E}[\tau_{AB}]$ computed
using the controlled process (76) and $q_{\sf nn}$ is used for computing the
transition rate using $\rho_{AB}$ obtained using (40) with $q_{\sf nn}$ and
Monte Carlo integration. The result differs from the other rate values by less
than 4%. Thus, we conclude that

  * •

the values for the transition rate obtained in three different ways are all
consistent; all of them fall into the 95% confidence interval of the rate
value computed using (40);

  * •

the method for estimating $\rho_{AB}$ and $\nu_{AB}$ using the neural network
solver and Monte Carlo integration can be promoted to higher dimensions.

Finally, we remark that the estimate for $\nu_{AB}$ is sensitive to the
normalization constant $Z_{V}$ in (80). It is important to estimate it
accurately. For example, Mueller’s potential has rather large region in which
the potential energy $V$ is relatively low – see Fig. 5. Sampling points for
determining $Z_{V}$ from a smaller region (a lower sublevel set of $V$) leads
to a notable discrepancy in the transition rate estimate. At the same time,
the estimate for $\rho_{AB}$ is much less sensitive to the accuracy of
$Z_{V}$. The reason in the difference of sensitivity is that $\nu_{AB}$ has
the gradient of the committor in its integral, while $\rho_{AB}$ has the
committor itself.

| Simulations, optimal control | TPT, NN | TPT, FEM
---|---|---|---
$\rho_{AB}$ | NA | 2.43e-4 | 2.36e-4
$\mathbb{E}[\tau_{AB}]$ |  5.05e-2 $\pm$ 0.42e-2 | NA | NA
$\nu_{AB}$ | 4.80e-3, [4.43e-3, 5.23e-3] | 4.99e-3 | 4.93e-3
Table 2: Estimates for $\rho_{AB}$, $\mathbb{E}[\tau_{AB}]$ and $\nu_{AB}$ for
the overdamped Langevin dynamics in Mueller’s potential at $\beta^{-1}=10$
computed using three schemes. The value for $\nu_{AB}=4.80e-3$ is obtained
with $\rho_{AB}=2.43e-4$. The 95% confidence interval is obtained using the
error estimate for $\mathbb{E}[\tau_{AB}]$.

###  6.2 The rugged Mueller potential in 10D

The test problem with Mueller’s potential can be upgraded by making it
10-dimensional and perturbing its energy landscape with an oscillatory
function:

$V(x)=V_{0}(x_{1},x_{2})+\gamma\sin(2k\pi x_{1})\sin(2k\pi
x_{2})+\frac{1}{\sigma^{2}}\sum_{i=3}^{10}x_{i}^{2}.$ (82)

Here, $V_{0}(x_{1},x_{2})$ is Mueller’s potential (71) and
$\gamma=9,k=5,\sigma=0.05$ as in [31]. Following [31], the set $A$ and $B$ are
chosen to be cylinders centered at $a=(-0.558,1.441)$ and $b=(0.623,0.028)$
with radius $r=0.1$. The exact solution to the committor problem for with $V$
given by (82) and such sets $A$ and $B$ is independent of
$x_{3},\ldots,x_{10}$. This allows us to use the FEM solver in 2D to test the
solution computed using the variational NN-based solver in 10D.

####  6.2.1 Computing the committor

We compute the committor using the same procedure as detailed in Section
6.1.1. For the FEM solver, the computational domain is
$\Omega=\\{x\in\mathbb{R}^{2}~{}|~{}V(x)\leq 250\\}$. For the variational NN-
based solver, a training set of $N_{\sf train}=64882$ points is generated by
sampling $2.1\times 10^{5}$ points in 10D and rarefying them into a delta-net
with $\delta=0.005$. The neural network in the solution model (64) has $L=3$
hidden layers and $N=10$ neurons in each layer. The committors computed by the
FEM and variational NN-based solvers are shown in Fig. 6(a) and (b)
respectively.

(a) FEM

(b) Variational NN

Figure 6: The committor for the overdamped Langevin dynamics with the rugged
Mueller potential (82) is computed by (a) FEM in 2D and (b) by variational NN
in 10D. The temperature is $\beta^{-1}=10$. The level sets of the potential
are superimposed.

We compute wMAE (73) and wRMSE (74) to assess numerical errors. A set of
$N_{\sf test}=44938$ test points in 10D is generated using metadynamics and
delta-net postprocessing. The variational NN solution is evaluated at these
test points and projected onto the 2D space $(x_{1},x_{2})$ to compare with
the FEM solution. The resulting wMAE and wRMSE are reported in Table 3.

Temperature | NN structure | wMAE | wRMSE
---|---|---|---
$\beta^{-1}=10$ |  $L=3$, $W=10$ | 1.08e-2 | 3.13e-2
Table 3: Errors wMAE and wRMSE in the forward committor for the overdamped
Langevin dynamics with rugged Mueller’s potential in 10D at temperature
$\beta^{-1}=10$ computed using the variational NN solver with a spatially
quasiuniform set of $64882$ training points. The parameter $L$ indicates the
number of hidden layers, and $W$ indicates the number of neurons per hidden
layer.

####  6.2.2 Estimation of the transition rate using the controlled process

The transition rate $\nu_{AB}$ is found using equation (62).

The expected crossover time $\mathbb{E}[\tau_{AB}]$ is estimated by averaging
crossover times of 250 transition trajectories sampled using the controlled
process (76) in 10D. The initial points of the trajectories are sampled
according to (37) as in the previous test problem. First, $N_{\partial A}$
points are equispaced on the circle of radius $r+\delta r$ centered at $a$
lying in the subspace $(x_{1},x_{2})$. The weights of these points assigned
according to (77) are shown in Fig. 7. Then these points are sampled based on
their probability weight. Three samples of trajectories of the controlled and
uncontrolled processes in 10D with the same initial state and the same
realizations of the Brownian motion projected onto the
$(x_{1},x_{2})$-subspace are visualized in Fig. 8 for comparison.

Figure 7: The probability distribution of the starting points of transition
trajectories on the boundary of $A$ for the overdamped Langevin dynamics in
the rugged Mueller potential (82) is shown by color code. The gray curves are
level sets of the potential.

(a) Controlled SDE (76)

(b) Uncontrolled SDE (78)

Figure 8: Comparison of three controlled (a) and uncontrolled (b) trajectories
for the rugged Mueller potential in 10D with the same initial positions and
the same noise realizations at $\beta^{-1}=10$ projected onto the
$(x_{1},x_{2})$-subspace. Each trajectory consists of 5000 time steps, with
time step $\Delta t=10^{-5}$. The sets $A$ and $B$ are marked by blue and
orange circles of radius $r=0.1$ respectively.

The probability $\rho_{AB}$ is computed using (27) for both committors $q_{\sf
fem}$ and $q_{\sf nn}$. Monte Carlo integration is used with $q_{\sf nn}$ over
the set of $N_{\sf test}=44938$ test points as described in Section 6.1.2.

The transition rate $\nu_{AB}$ is also estimated using (81) and the committors
$q_{\sf fem}$ and $q_{\sf nn}$ for comparison. The results are presented in
Table 4.

The discrepancy between the estimates for $\rho_{AB}$ obtained using the
committors $q_{\sf nn}$ in 10D and $q_{\sf fem}$ in 2D and equations (79) and
(27) is about 4.5%. The difference between the transition rates $\nu_{AB}$
obtained using $q_{\sf nn}$ with Monte Carlo integration and $q_{\sf fem}$
with formula (28) is approximately 13%. On the other hand, the transition rate
computed via (62) uses the expected crossover time and $\rho_{AB}$ acquired
from $q_{\sf nn}$, resulting in a transition rate that differs from the FEM
rate by 6.5%. Hence we conclude that our proposed scheme for estimating
$\rho_{AB}$ and $\nu_{AB}$ with a neural network solver yields a reasonable
accuracy in higher dimensions.

| Simulations, optimal control | TPT, NN | TPT, FEM
---|---|---|---
$\rho_{AB}$ | NA | 2.56e-4 | 2.45e-4
$\mathbb{E}[\tau_{AB}]$ |  5.94e-2 $\pm$ 0.46e-2 | NA | NA
$\nu_{AB}$ | 4.31e-3, [4.0e-3,4.67e-3] | 5.23e-3 | 4.61e-3
Table 4: Estimates for $\rho_{AB}$, $\mathbb{E}[\tau_{AB}]$ and $\nu_{AB}$ for
the overdamped Langevin dynamics in the rugged Mueller potential (82) at
$\beta^{-1}=10$ in $\mathbb{R}^{10}$ computed using three schemes. The value
for $\nu_{AB}=4.31e-3$ is obtained with $\rho_{AB}=2.56e-4$. The 95%
confidence interval is obtained using the error estimate for
$\mathbb{E}[\tau_{AB}]$.

###  6.3 Duffing Oscillator in 1D

Now we test the proposed methodology on the bistable Duffing oscillator with
mass $m=1$, friction coefficient $\gamma=0.5$, and the potential energy
function $V(x)=0.25(x^{2}-1)^{2}$. The dynamics are governed by the Langevin
SDE

$d\begin{bmatrix}X_{t}\\\ P_{t}\end{bmatrix}=\begin{bmatrix}P_{t}\\\
-X_{t}(X_{t}^{2}-1)-\frac{1}{2}P_{t}\end{bmatrix}dt+\sqrt{\epsilon}\begin{bmatrix}0\\\
dW_{t}\end{bmatrix}.$ (83)

The system has two stable equilibria $a=(-1,0)$ and $b=(1,0)$ and an unstable
equilibrium at the origin. The full energy of the system is
$H(x,p)=\tfrac{1}{2}p^{2}+V(x)$. The invariant probability density is
$\mu(x,p)=Z_{H}^{-1}\exp(-H(x,p)/\epsilon)$. The sets $A$ and $B$ are chosen
to be ellipses with radii $r_{x}=0.3$ and $r_{y}=0.4$ centered at $a$ and $b$
respectively. Two values of the noise coefficient are used: $\epsilon=0.1$ and
$\epsilon=0.05$.

####  6.3.1 Computation of forward and backward committor functions

Since the dynamics (83) are not time-reversible, we implement the Physics-
Informed Neural Network (PINN) approach detailed in Section 5.2 to solve the
committor problem. A uniform grid with a total of 16000 points in the
rectangle $[-2.5,2.5]\times[-2,2]$ is taken as training data. The architecture
of the neural network is as in equation (65) with a single hidden layer,
$L=1$, and $W=40$ neurons in it. The Adam optimizer is used with the learning
rate $10^{-3}$ for 500 epochs. We also compute the committors using FEM as
described in F.2. The computed forward and backward committors, $q^{+}(x,p)$
and $q^{-}(x,p)$, for $\epsilon=0.1$ and $\epsilon=0.05$ are displayed in
Figs. 9 and 10 respectively. The theoretical relationship between them is
$q^{-}(x,p)=1-q^{+}(x,-p)$. However, we still computed $q^{-}$ using FEM
because the FEM mesh is not symmetric.

(a) FEM (solid) & PINN (dashed), $\epsilon=0.1$

(b) FEM (solid) & PINN (dashed), $\epsilon=0.05$

Figure 9: Comparison of contours of forward committor computed using FEM and
PINNs, marked in solid and dashed lines respectively for Duffing Oscillator at
$\epsilon=0.1$ and $\epsilon=0.05$.

(a) FEM (solid) & PINN (dashed), $\epsilon=0.1$

(b) FEM (solid) & PINN (dashed), $\epsilon=0.05$

Figure 10: Comparison of contours of backward committor computed using FEM and
PINNs, marked in solid and dashed lines respectively for Duffing Oscillator at
$\epsilon=0.1$ and $\epsilon=0.05$.

We call the discrepancies between the FEM and PINN solutions computed by (73)
the weighted mean absolute difference (wMAD) and the weighted root mean square
difference (wRMSD) with weights at the training points $(x_{i},p_{i})$ given
by

$w(x_{i},p_{i})=\frac{{q}^{+}(x_{i},p_{i}){q}^{-}(x_{i},p_{i})e^{-H(x_{i},p_{i})/\epsilon}}{\sum_{j=1}^{N_{\sf
train}}{q}^{+}(x_{j},p_{j}){q}^{-}(x_{j},p_{j})e^{-H(x_{j},p_{j})/\epsilon}},\quad
i=1,\ldots,N_{\sf train}.$ (84)

Table 5 shows the wMAD and wRMSD for forward and backward committors computed
using PINN and FEM for $\epsilon=0.1$ and $\epsilon=0.05$.

$\epsilon$ | NN structure |  $q^{+}$, wMAD |  $q^{+}$, wRMSD |  $q^{-}$, wMAD |  $q^{-}$, wRMSD
---|---|---|---|---|---
$\epsilon=0.1$ |  $L=1$, $W=40$ | 1.6e-2 | 2.0e-2 | 1.8e-2 | 2.2e-2
$\epsilon=0.05$ |  $L=1$, $W=40$ | 1.3e-2 | 2.0e-2 | 1.3e-2 | 2.0e-2
Table 5: Comparison of the numerical solutions for the forward and backward
committors computed using PINNs and FEM. Metrics wMAD and wRMSD are used. The
parameter $L$ indicates the number of hidden layers (see (65)), and $W$
indicates the number of neurons per hidden layer.

####  6.3.2 Estimation of the transition rate using the controlled process

The optimally controlled process for the Langevin dynamics is governed by

$\begin{cases}dX_{t}=P_{t}dt\\\
dP_{t}=-\left[X_{t}(X_{t}^{2}-1)+\tfrac{1}{2}P_{t}-\epsilon\frac{\partial}{\partial
p}\log(q_{\sf nn}^{+})\right]dt+\sqrt{\epsilon}dW.\end{cases}$ (85)

Fig. 11 shows three sampled trajectories with and without the influence of the
optimal control starting at the same initial position and the same
realizations of the Brownian motion for $\epsilon=0.1$ (Fig. 11(a,b)) and
$\epsilon=0.05$ (Fig. 11(c,d)). The trajectories governed by the original
Langevin dynamics stay near region $A$. In contrast, the trajectories governed
by the optimally controlled dynamics (85) leave $A$ and reach region $B$.

(a) $\epsilon=0.1$, w/o o/c

(b) $\epsilon=0.1$, o/c

(c) $\epsilon=0.05$, w/o o/c

(d) $\epsilon=0.05$, o/c

Figure 11: Comparison of three uncontrolled and controlled processes with the
same initial random position for Duffing oscillator at $\epsilon=0.1$ and
$\epsilon=0.05$. (a) and (c) shows three trajectories under the original
dynamics at $\epsilon=0.1$ and $\epsilon=0.05$ respectively. (b) and (d) shows
three trajectories under the original dynamics for the two cases respectively.
5000 time steps are sampled with $\Delta t=10^{-3}$. Region $A$ and $B$ are
ellipses colored in blue and orange.

Next, we find the transition rate $\nu_{AB}$ at $\epsilon=0.1$ and
$\epsilon=0.05$ in four ways. The results are summarized in Tables 6 and 7.

  1. 1.

_Simulations with optimal control._ The expected crossover time
$\mathbb{E}[\tau_{AB}]$ is averaged over 250 trajectories governed by (85).
The distribution of the starting points of these trajectories defined by (37)
is obtained as described in Section 6.1.2. It is displayed in Fig. 12. The
probability $\rho_{AB}$ is found by (27) similarly to how it is done in the
test problem with Mueller’s potential. The PINN committors and the PINN
training points are used for Monte Carlo integration. The normalization
constant for the invariant density is also found using Monte Carlo
integration. Then formula (62) is used to find $\nu_{AB}$.

  2. 2.

_Simulations without optimal control._ Direct simulations of the uncontrolled
Langevin dynamics (83) are used to find $\nu_{AB}$. Ten simulations of
$10^{7}$ time steps each with timestep $\Delta t=5\times 10^{-3}$ were
performed.

  3. 3.

_TPT, NN._ The rate $\nu_{AB}$ was found by (35) using the gradient of the
committor computed using PINNs.

  4. 4.

_TPT, FEM._ Likewise, except for the FEM committor was used.

Figure 12: The probability distribution of the starting points of transition
trajectories on the boundary of $A$ for the Duffing oscillator is shown by
color code. The gray curves are level sets of Hamiltonian.  Duffing oscillator
$\epsilon=1/10$
---
| Simul., o/c | Simul., w/o o/c | TPT, NN | TPT, FEM
$\rho_{AB}$ | NA |  4.31e-2 $\pm$ 0.12e-2 | 3.97e-2 | 4.04e-2
$\mathbb{E}[\tau_{AB}]$ | 6.88 $\pm$ 0.34 | 7.32 $\pm$ 0.14 | NA | NA
$\nu_{AB}$ | [5.50e-3,6.07e-3] | [5.76e-3,6.01e-3] | 4.53e-3 | 5.74e-3
Table 6: Comparison of the estimates for $\rho_{AB}$, $\mathbb{E}[\tau_{AB}]$
and $\nu_{AB}$ for Duffing oscillator at $\epsilon=0.1$.  Duffing oscillator
$\epsilon=1/20$
---
| Simul., o/c | Simul., w/o o/c | TPT, NN | TPT, FEM
$\rho_{AB}$ | NA |  4.5e-3 $\pm$ 0.4e-3 | 4.23e-3 | 4.07e-3
$\mathbb{E}[\tau_{AB}]$ | 7.34 $\pm$ 0.33 | 7.48 $\pm$ 0.49 | NA | NA
$\nu_{AB}$ | [5.53e-4,6.06e-4] | [5.49e-4,6.51e-4] | 4.72e-4 | 5.49e-4
Table 7: Comparison of the estimates for $\rho_{AB}$, $\mathbb{E}[\tau_{AB}]$
and $\nu_{AB}$ for Duffing oscillator at $\epsilon=0.05$.

The results in tables 6 and 7 show that the 95% confidence intervals for the
transition rate $\nu_{AB}$ at $\epsilon=0.1$ and $\epsilon=0.05$ estimated by
means of simulations with and without optimal control largely overlap. The 95%
confidence intervals expected crossover times $\mathbb{E}[\tau_{AB}]$ largely
overlap for $\epsilon=0.05$ and slightly overlap for $\epsilon=0.1$. The
estimates for probability $\rho_{AB}$ that a trajectory at a random time is
reactive obtained by direct simulations of uncontrolled dynamics and TPT&NN
and TPT&FEM are all consistent for $\epsilon=0.05$ and both TPT-based
estimates for $\rho_{AB}$ are smaller than those by direct simulations for
$\epsilon=0.1$. At both values of $\epsilon$, the TPT estimates for
$\rho_{AB}$ obtained using the FEM and PINN committors are consistent, while
there is a notable discrepancy between the estimates for $\nu_{AB}$ by TPT&NN
and TPT&FEM. This discrepancy must be caused by the fact that $\nu_{AB}$ uses
the gradient of the committors while $\rho_{AB}$ involves the committors
themselves as shown in D. At both values of $\epsilon$, the TPT&FEM estimate
for $\nu_{AB}$ falls into the 95% confidence intervals obtained using
simulations, controlled or uncontrolled, while the TPT&NN seems to
underestimate the transition rate.

###  6.4 Lennard-Jones-7 in 2D

Finally, we apply the proposed methodology to estimate the transition rate
between the trapezoidal and the hexagonal configurations of the Lennard-
Jones-7 cluster (LJ7) in a plane. This is a popular test problem in chemical
physics [64, 65, 66, 67]. In this example, we will compute the committor for
the reduced model 2D and use it to construct an approximation to the optimal
controller in the original 14D model. The expected crossover time for the
original 14D model and the estimate for $\rho_{AB}$ for the 2D model will be
used to determine the transition rate. The result will be compared with those
obtained via brute force simulations of the original uncontrolled dynamics in
14D.

We consider seven two-dimensional particles interacting according to the
Lennard-Jones pair potential

$V_{\sf
pair}(r)=4a\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]$
(86)

where $\sigma>0$ and $a>0$ are parameters controlling the range and strength
of interparticle interaction respectively. We set $\sigma=1$ and $a=1$. The
potential energy of the system

$V_{\sf LJ}(x)=\sum_{\begin{subarray}{c}i,j=1\\\ i<j\end{subarray}}^{7}V_{\sf
pair}(\|x_{i}-x_{j}\|)$ (87)

has four geometrically distinct local minima denoted by
$C_{0},C_{1},C_{2},C_{3}$ shown in Fig. 13.

Figure 13: The free energy of LJ7 in 2D in CVs second and third central moment
of coordination number $\mu_{2},\mu_{3}$ and the four geometrically distinct
local minima.

We assume that the original system is evolving according to the overdamped
Langevin dynamics (6). This choice is dictated by our wish to construct a
controller using the committor for the reduced model in collective variables.
Since the collective variables are functions only of $x$, the committor for
the reduced model lifted to the original phase space will depend only on $x$
and not on the momenta $p$. Therefore, it cannot give an approximation to the
optimal control for the Langevin dynamics – see equation (59). We set
$\beta=5$ as in [34].

####  6.4.1 The reduced model

Following [68, 69, 34], we pick the 2nd and 3rd central moments of the
coordination numbers as the collective variables (CVs) for LJ7. The coordinate
number of particle $i$ is a smooth function approximating the number of
nearest neighbors of $i$:

$c_{i}(x)=\sum_{i\neq
j}\frac{1-(\frac{r_{ij}}{1.5\sigma})^{8}}{1-(\frac{r_{ij}}{1.5\sigma})^{16}},\quad
r_{ij}=\|x_{i}-x_{j}\|.$ (88)

The $k$-th central moment of $c_{i}(x)$ is defined by

$\mu_{k}(x)=\frac{1}{7}\sum_{i=1}^{7}(c_{i}(x)-\bar{c}(x))^{k},\quad{\rm
where}\quad\bar{c}(x)=\frac{1}{7}\sum_{j=1}^{7}c_{j}(x).$ (89)

The reduced model is governed by the overdamped Langevin dynamics in
collective variables (9)

$\displaystyle d\left[\begin{array}[]{c}\mu_{2}\\\ \mu_{3}\end{array}\right]$
$\displaystyle=\left[-\nabla
M(\mu_{2},\mu_{3})F(\mu_{2},\mu_{3})+\beta^{-1}\cdot\nabla
M(\mu_{2},\mu_{3})\right]dt$ (92)
$\displaystyle+\sqrt{2\beta^{-1}}M^{1/2}(\mu_{2},\mu_{3})\left[\begin{array}[]{c}dw_{1}\\\
dw_{2}\end{array}\right].$ (95)

The corresponding generator is given by (23). Fig. 13 displays the free
energy111We thank Luke Evans for sharing with us the free energy and the
diffusion matrix in CVs $\mu_{2},$ and $\mu_{3}$. $F(\mu_{2},\mu_{3})$. The
diffusion matrix $M(\mu_{2},\mu_{3})$ varies significantly throughout the
accessible free energy region (see Fig. 8 in [34]). The computation of
$F(\mu_{2},\mu_{3})$ and $M(\mu_{2},\mu_{3})$ is detailed in Appendix A in
[34].

Regions $A_{\sf CV}$ and $B_{\sf CV}$ are chosen around minima $C_{3}$ (the
trapezoid) and $C_{0}$ (the hexagon) respectively. We use the subscript to CV
to indicate that these regions are defined in the set of collective variables.
Region $A_{\sf CV}$ is a circle centered at $(0.5526,-0.0935)$ of radius
$r=0.1$ and while region $B_{\sf CV}$ is a tilted ellipse defined by the
equation

$\frac{\left((x-c_{x})\cos\theta+(y-c_{y})\sin\theta\right)^{2}}{r_{x}^{2}}+\frac{\left((x-c_{x})\sin\theta+(y-c_{y})\cos\theta\right)^{2}}{r_{y}^{2}}=1$
(96)

where $(c_{x},c_{y})=(0.7184,1.1607),r_{x}=0.15,r_{y}=0.03$ and
$\theta=5\pi/12$.

####  6.4.2 Computation of the committor for the reduced model

The committor for the reduced model is computed in two ways: using FEM and the
variational NN-based solver described in Section 5.1. For the variational NN,
the neural network (65) with $L=2$ hidden layers and $W=10$ neurons per layer
has been used to minimize the loss (69). The training points were $10^{4}$
trajectory data projected to the space $(\mu_{2},\mu_{3})$ and assumed to be
distributed according to the invariant density e.g., $\rho\sim e^{\beta F}$
where $F$ is the free energy in $(\mu_{2},\mu_{3})$ – see Section 4.2.2 in
Ref. [34] for more details. The resulting loss to be minimized hence becomes

${\sf Loss}(\theta)=\frac{1}{K}\sum_{i=k}^{K}\left[\nabla
q(x_{k};\theta)^{\top}M(x_{k})\nabla q(x_{k};\theta)\right].$ (97)

The results are displayed in Fig. 14. The wMAE and wRMSE are given in Table 8.

(a)  FEM

(b) Variational NN

Figure 14: The forward committor for LJ7 in 2D computed in the space of collective variables at $\beta=5$ using (a) FEM and (b) the variational NN-based solver. The gray curves in (b) are level sets of the free energy. Temperature | NN structure | wMAE | wRMSE
---|---|---|---
$\beta=5$ |  $L=2$, $W=10$ | 1.5e-2 | 2.4e-2
Table 8: The errors wMAE and wRMSE of the committor function for LJ7 computed
using the variational NN-based solver for the two-dimensional reduced model.
The committor computed using FEM is taken as the ground truth.

####  6.4.3 Estimation of the transition rate using the reduced model and the
controlled process

We set up the controlled process in the original 14-dimensional coordinate
space as

$dX_{t}=\left(-\nabla V_{\sf LJ}(X_{t})+2\beta^{-1}\nabla_{x}\ln q_{\sf
nn}(\mu_{2}(X_{t}),\mu_{3}(X_{t}))\right)dt+\sqrt{2\beta^{-1}}dW_{t},$ (98)

where $V_{\sf LJ}$ is defined by (87) and $q_{\sf nn}$ is the committor
computed for the reduced model using the variational NN-based solver. The
Metropolis-Adjusted Langevin Algorithm (MALA) [70] with the time step $\Delta
t=5\times 10^{-5}$ has been used for time integration to prevent very large
moves of the system that can occur due to extremely strong repulsive forces.

The sets $A$ and $B$ in the original coordinate space $\mathbb{R}^{14}$ are
defined by lifting the sets $A_{\sf CV}$ and $B_{\sf CV}$:

$\displaystyle A$
$\displaystyle:=\\{x\in\mathbb{R}^{14}~{}|~{}(\mu_{2}(x),\mu_{3}(x))\in A_{\sf
CV}\\},$ $\displaystyle B$
$\displaystyle:=\\{x\in\mathbb{R}^{14}~{}|~{}(\mu_{2}(x),\mu_{3}(x))\in B_{\sf
CV}\\}.$

Three trajectories of the controlled process (98) and three trajectories of
the uncontrolled overdamped Langevin dynamics (6) in the 14D with the same
three realizations of the Brownian motion projected to the space
$(\mu_{2},\mu_{3})$ are displayed in Fig. 16 (a) and (b) respectively.

The expected crossover time $\mathbb{E}[\tau_{AB}]$ is averaged over 254
trajectories of SDE (98). The starting points for these trajectories near the
boundary of $A$ are sampled as follows. First, the probability weights of the
points on the boundary of $A_{\sf CV}$ are computed as described in Section
6.1.2 (see Fig. 15). Then, points on $\partial A_{\sf CV}$ are sampled
according to these weights and lifted to the original coordinate space by
running biased simulations as described in Appendix A of [34] (see equations
(A.3) and (A.4) there).

Figure 15: The probability distribution of the starting points of transition
trajectories on the boundary of $A$ for LJ7 is shown by color code. The gray
curves are level sets of the free energy.

(a) $\beta=5$, o/c

(b) $\beta=5$, w/o o/c

Figure 16: Three trajectories of (a) the controlled process (98) and (b) the
uncontrolled process projected into the space of collective variables, i.e.,
the overdamped Langevin dynamics (6) in the 14D projected into the space of
collective variables with the same starting points and the same realizations
of Brownian motion.

We compute the transition rate $\nu_{AB}$ in four ways analogous to those for
the Duffing oscillator (Section 6.3.2).

  1. 1.

_Simulations, optimal control, 14D._ The rate is calculated using (62). The
expected crossover time is estimated as described above using the controlled
process (98) in 14D. The probability $\rho_{AB}$ is obtained for the reduced
2D model using (79)–(80) with the committor computed by the variational NN-
based solver and a uniform grid of points rather than the actual training
points for the neural network.

  2. 2.

_Simulations without optimal control, 14D._ Ten runs of direct simulations of
the uncontrolled overdamped Langevin dynamics (6) in 14D of $10^{8}$ timesteps
with $\Delta t=5\cdot 10^{-5}$ were executed. The numbers of transitions from
$A$ to $B$ that occurred in these runs were 28, 138, 160, 146, 93, 63, 158,
165, 171, and 160.

  3. 3.

_TPT, NN._ The rate $\nu_{AB}$ was found by (28) for the reduced 2D model
using the gradient of the committor computed using the variational NN-based
solver.

  4. 4.

_TPT, FEM._ Likewise, except for the FEM committor was used.

The results are displayed in Table 9. The following observations can be made.

  1. 1.

For the transition rates obtained using controlled and uncontrolled
simulations in the 14D, the relative error is about 14%. There is a large
overlap of the 95% confidence intervals.

  2. 2.

There is a discrepancy in the expected crossover time $\mathbb{E}[\tau_{AB}]$
computed using controlled and uncontrolled simulations in the 14D: the
expected crossover time for the controlled process exceeds that for the
uncontrolled process by approximately 35%. This discrepancy is primarily
caused by the fact that the controller obtained by lifting the committor
computed for the reduced model is not optimal. It still drives the
trajectories away from the boundary of $A$ but, unlike the optimal controller,
somewhat affects the statistics for the transition trajectories.

  3. 3.

The rates obtained using the reduced 2D model are highly exaggerated (by the
factor of approximately four) as one can expect given the Zhang-Hartmann-
Schuette rate formula (61). Indeed, the collective variables $\mu_{2}$ and
$\mu_{3}$ are chosen due to their ability to separate the four geometrically
distinct local minima of LJ7 while there is no indication that they are
supposed to represent the dynamics accurately.

  4. 4.

The values of $\rho_{AB}$ and $\nu_{AB}$ computed for the reduced 2D model
using the FEM and the variational NNs are in good agreement with each other.

| Simul., o/c | Simul., w/o o/c | TPT, NN | TPT, FEM
---|---|---|---|---
Dimension | 14D | 14D | 2D | 2D
$\rho_{AB}$ | NA |  0.080 $\pm$ 0.024 | 0.106 | 0.108
$\mathbb{E}[\tau_{AB}]$ |  4.88 $\pm$ 0.48 |  3.16 $\pm$ 0.21 | NA | NA
$\nu_{AB}$ |  0.022, [0.020, 0.024] |  0.025, [0.019, 0.033] | 0.097 | 0.086
Table 9: Comparison of the estimated $\rho_{AB}$, $\mathbb{E}(\tau_{AB})$ and
$\nu_{AB}$ for LJ7 at $\beta=5$.

##  7 Conclusion

In this work, we have proposed a methodology for sampling transition
trajectories and estimating transition rates in systems governed by SDEs using
optimal control and perhaps model reduction.

Our main theoretical contribution is the proof of Theorem 3.1 establishing the
optimality of the control obtained from the committor via the Doob
$h$-transform for a broad class of processes including the Langevin dynamics
and the overdamped Langevin dynamics in collective variables.

We have elaborated on a number of practical aspects related to the use of
neural network-based solvers, finite element methods, and sampling reactive
trajectories.

We have conducted in-depth case studies of three benchmark systems. In
particular, we have demonstrated that the optimal control and the estimate for
the probability of a trajectory to be reactive at a random moment of time
obtained for the reduced model result is a reasonably good estimate of the
transition rate even if the collective variables do not represent the dynamics
accurately.

Our codes are published on GitHub [44, 45].

Further improvement of the proposed methodology can be done in the following
two directions. First, the design of collective variables is important for an
accurate representation of the dynamics. Autoencoders (see e.g. [15] and
references therein) with an appropriate choice of the loss function seem to be
a promising tool. Second, the neural network-based techniques for solving the
committor problem are promotable to higher dimensions. In this work, we
intentionally calculated all required quantities for the use of the transition
path theory without meshing the space. We did not attempt, though, to use
these techniques in higher dimensions. We are leaving these research topics
for future work.

##  8 Acknowledgements

We thank Dr. Luke Evans for providing us with the free energy and diffusion
matrix data for the Lennard-Jones-7 test problem. We also thank the UMD REU
students Luke Triplett, Dmitry Pinchuk, Prisca Calkins, and William Clark for
their investigation into neural network-based committor solvers. This work was
partially supported by AFOSR MURI grant FA9550-20-1-0397 and by NSF REU grant
DMS-2149913.

##  Appendix A  Proof of equation (40):
$\nu_{AB}={\rho_{AB}}/{\mathbb{E}[\tau_{AB}]}$

###### Proof.

Let $X_{t}$, $0\leq t\leq T$, be a long trajectory. We decompose the interval
$[0,T]$ into two subsets $[0,T]=\mathcal{I}_{A}\cup\mathcal{I}_{B}$ where

$\begin{cases}\mathcal{I}_{A}([0,T]):=\\{t\in[0,T]~{}|~{}\tau_{A}^{-}(X_{t})<\tau_{B}^{-}(X_{t})\\},\\\
\mathcal{I}_{B}([0,T]):=\\{t\in[0,T]~{}|~{}\tau_{B}^{-}(X_{t})<\tau_{A}^{-}(X_{t})\\}.\end{cases}$
(A-1)

In words, $\mathcal{I}_{A}([0,T])$ is the set of all moments of time $t$ in
the interval $[0,T]$ such that the trajectory at time $t$, $X_{t}$, last
visited $\bar{A}$ rather than $\bar{B}$. The set $\mathcal{I}_{B}([0,T])$ is
described likewise. Let $T_{A}$ and $T_{B}$ be the total lengths of
$\mathcal{I}_{A}([0,T])$ and $\mathcal{I}_{B}([0,T])$ respectively. The
probabilities $\rho_{A}$ and $\rho_{B}$ that the trajectory at a randomly
picked moment of time last visited $A$ or $B$ are, respectively,

$\rho_{A}=\lim_{T\rightarrow\infty}\frac{T_{A}}{T},\quad\rho_{B}=\lim_{T\rightarrow\infty}\frac{T_{B}}{T}.$
(A-2)

The set $\mathcal{I}_{A}([0,T])$ is further decomposed into two subsets
$\mathcal{I}_{A}([0,T])=\mathcal{I}_{AA}([0,T])\cup\mathcal{I}_{AB}([0,T])$ of
total lengths $T_{AA}$ and $T_{AB}$ respectively where

$\begin{cases}\mathcal{I}_{AA}([0,T]):=\\{t\in[0,T]~{}|~{}\tau_{A}^{-}(X_{t})<\tau_{B}^{-}(X_{t})~{}\&~{}\tau_{A}^{+}(X_{t})<\tau_{B}^{+}(X_{t})\\},\\\
\mathcal{I}_{AB}([0,T]):=\\{t\in[0,T]~{}|~{}\tau_{A}^{-}(X_{t})<\tau_{B}^{-}(X_{t})~{}\&~{}\tau_{B}^{+}(X_{t})<\tau_{A}^{+}(X_{t})\\}.\end{cases}$
(A-3)

I.e., $\mathcal{I}_{AA}([0,T])$ is the set of moments of times $t\in[0,T]$
such that the trajectory at time $t$, $X_{t}$, last visited $\bar{A}$ rather
than $\bar{B}$ and going to hit next $\bar{A}$ rather than $\bar{B}$, while
$\mathcal{I}_{AB}([0,T])$ is the subset of moments of time $t\in[0,T]$ such
that the trajectory is reactive. Respectively, the probability $\rho_{AA}$
that a trajectory at a randomly picked time $t$ last hit $\bar{A}$ rather than
$\bar{B}$ and is not reactive and the probability $\rho_{AB}$ that a
trajectory at a randomly picked time $t$ is reactive are given by

$\rho_{AA}=\lim_{T\rightarrow\infty}\frac{T_{AA}}{T},\quad\rho_{AB}=\lim_{T\rightarrow\infty}\frac{T_{AB}}{T}.$
(A-4)

Now we recall the definitions of the transition rate $\nu_{AB}$ and the
expected crossover time $\mathbb{E}[\tau_{AB}]$:

$\nu_{AB}=\lim_{T\rightarrow\infty}\frac{N_{AB}}{T},\quad\mathbb{E}[\tau_{AB}]=\lim_{T\rightarrow\infty}\frac{T_{AB}}{N_{AB}}.$
(A-5)

Hence the expected crossover time $\mathbb{E}[\tau_{AB}]$ can be written as

$\nu_{AB}=\lim_{T\rightarrow\infty}\frac{N_{AB}}{T}=\lim_{T\rightarrow\infty}\frac{N_{AB}}{T_{AB}}\frac{T_{AB}}{T}=\frac{\rho_{AB}}{\mathbb{E}[\tau_{AB}]}.$
(A-6)

∎

##  Appendix B Proof that $-\nabla\cdot J_{AB}=0$ in $\Omega_{AB}$

Let us show that the divergence of the reactive current, or, equivalently, the
stationary current of the transition path process, vanishes in $\Omega_{AB}$.
We will need a formula for the divergence of a matrix-vector product. It can
be checked directly that for any $A\in\mathbb{R}^{d\times d}$ and any
$y\in\mathbb{R^{d}}$,

$\nabla\cdot(Ay)={\sf div}A\cdot y+{\sf tr}\left(A\nabla y\right),\quad{\rm
where}\quad\nabla y=\left[\begin{array}[]{c}\nabla y_{1}^{\top}\\\ \vdots\\\
\nabla y_{d}^{\top}\end{array}\right].$ (B-7)

Using (B-7) we calculate:

$\displaystyle\nabla\cdot J_{AB}$
$\displaystyle=-\nabla\cdot\left[\left(b\mu-\frac{1}{2}{\sf
div}\left(\sigma\sigma^{\top}\mu\right)\right)q^{+}q^{-}+\frac{1}{2}\mu\sigma\sigma^{\top}\left(q^{-}\nabla
q^{+}-q^{+}\nabla q^{-}\right)\right]$
$\displaystyle=q^{+}q^{-}\nabla\cdot\left(b\mu-\frac{1}{2}{\sf
div}\left(\sigma\sigma^{\top}\mu\right)\right)$
$\displaystyle+\left(b\mu-\frac{1}{2}{\sf
div}\left(\sigma\sigma^{\top}\mu\right)\right)\cdot\left(q^{+}\nabla
q^{-}+q^{-}\nabla q^{+}\right)$ $\displaystyle+\frac{1}{2}{\sf
div}\left(\sigma\sigma^{\top}\mu\right)\cdot\left(q^{+}\nabla
q^{-}-q^{-}\nabla q^{+}\right)+\frac{1}{2}{\sf
tr}\left(\sigma\sigma^{\top}\nabla\nabla q^{+}\right)\mu q^{-}+$
$\displaystyle+\frac{1}{2}{\sf tr}\left(\sigma\sigma^{\top}\mu\nabla
q^{-}[\nabla q^{+}]^{\top}\right)-\frac{1}{2}{\sf
tr}\left(\sigma\sigma^{\top}\mu\nabla q^{+}[\nabla
q^{-}]^{\top}\right)-\frac{1}{2}{\sf tr}\left(\sigma\sigma^{\top}\nabla\nabla
q^{-}\right)\mu q^{+}$ $\displaystyle=q^{+}q^{-}\nabla\cdot J+\mu
q^{-}\mathcal{L}q^{+}-\mu q^{+}\mathcal{L}^{\dagger}q^{-}=0.$

In the last expression, $J$ is the stationary current for the invariant
density $\mu$, and hence $\mathcal{L}^{*}\mu=-\nabla\cdot J=0$ in $\Omega$.
The last two terms are zero as by (19).

##  Appendix C Proof of Theorem 3.1

###### Proof.

This proof combines ideas from Gao et al. ([3], the proof of Theorem 3.3) and
L. C. Evans’s notes on the control theory [50]

Step 1. Regularization. We first consider a regularized optimal control
problem in which the exit cost (50) is replaced with a finite exit cost

$g_{N}(x)=\begin{cases}N,&x\in\overline{A}\\\
0,&x\in\overline{B}\end{cases},\quad{\rm where}~{}N~{}\text{is a large
number}.$ (C-8)

Let $c_{N}^{\ast}(x)$ be the infimum of the cost functional
$C_{x}[N;v(\cdot)]$ with the regularized exit cost (C-8) among all admissible
controls. Note that the admissible set $\mathcal{V}$ is not empty because the
time $\tau_{AB}<\infty$ almost surely since the system is ergodic and the
domain $\Omega$ is compact. Furthermore, $c_{N}^{\ast}(x)<\infty$ as $v\equiv
0$ is an admissible control and the corresponding cost functional is
$C_{x}[N;v(\cdot)]=N(1-q^{+}(x))<\infty$. Indeed, if $v\equiv 0$, then the
process hits $\partial B$ first with probability $q^{+}(x)$ and scores zero
and hits $\partial A$ first with probability $1-q^{+}(x)$ and scores $N$.

We also define a regularized forward committor as the solution to the
following boundary-value problem

$\begin{cases}\mathcal{L}q_{N}^{+}=0,&x\in\Omega_{AB},\\\
q_{N}^{+}=e^{-N},&x\in\partial A\\\ q_{N}^{+}=1,&x\in\partial B\\\
\frac{\partial q_{N}^{+}}{\partial\hat{n}}=0,&x\in\partial\Omega\end{cases},$
(C-9)

where $\mathcal{L}$ is the generator for (14). It is easy to check that

$q_{N}^{+}=q^{+}+(1-q^{+})e^{-N}.$ (C-10)

Step 2. Show that $c_{N}^{*}\geq-\log q_{N}^{+}$. The regularized forward
committor can be written as

$q_{N}^{+}(x)=\mathbb{E}_{P}\left[e^{-g_{N}(X_{\tau_{AB}})}~{}|~{}X_{0}=x\right]\equiv\mathbb{E}_{P,x}\left[e^{-g_{N}(X_{\tau_{AB}})}\right].$
(C-11)

Indeed, the process $X_{t}$ governed by (14) with $X_{0}=x$ reaches $\partial
A$ at time the stopping time $\tau_{AB}$ with probability $(1-q^{+})$ and
scores $e^{-N}$, and reaches $\partial B$ at $\tau_{AB}$ and scores $1$. This
results in the expectation given by the right-hand side of (C-10) which is
equal to $q_{N}^{+}(x)$.

Let $Y_{t}$ be the controlled process governed by (47) with a control $\sigma
v$, $v\in C^{1}(\Omega_{AB}$, satisfying Novikov’s condition

$\mathbb{E}_{P}\left[\exp\left(\int_{0}^{\tau{AB}}\frac{1}{2}\|\sigma^{T}(Y_{s})v(Y_{s})\|^{2}ds\right)\right]<\infty$
(C-12)

and $P_{v}$ be the probability measure on the path space of this process.
According to the Girsanov theorem (Theorem 8.6.5, p.158 in [71]),

$\mathbb{E}_{P,x}\left[e^{-g_{N}(X_{\tau_{AB}})}\right]=\mathbb{E}_{P_{v},x}\left[e^{-g_{N}(Y_{\tau_{AB}})}\right]=\mathbb{E}_{P,x}\left[e^{-g_{N}(Y_{\tau_{AB}})}\frac{dP_{v}}{dP}\right],$
(C-13)

where the Radon-Nikodym derivative

$\frac{dP_{v}}{dP}=\exp\left\\{-\int_{0}^{\tau_{AB}}\sigma^{\top}v(Y_{s})\cdot
dW_{s}-\frac{1}{2}\int_{0}^{\tau_{AB}}\|\sigma^{\top}v(Y_{s})\|^{2}ds\right\\}>0~{}~{}{\rm
P-a.s.}$ (C-14)

Therefore,

$\mathbb{E}_{P,x}\left[e^{-g_{N}(X_{\tau_{AB}})}\right]=\mathbb{E}_{P,x}\left[e^{-g_{N}(Y_{\tau_{AB}})-\int_{0}^{\tau_{AB}}\sigma^{\top}v(Y_{s})\cdot
dW_{s}-\frac{1}{2}\int_{0}^{\tau_{AB}}\|\sigma^{\top}v(Y_{s})\|^{2}ds}\right].$
(C-15)

By Jensen’s inequality, for any smooth convex function $\phi$ and a random
variable $Z$ we have $\phi(\mathbb{E}[Z])\leq\mathbb{E}[\phi(Z)]$. Applying it
to the right-hand side of (C-15) we get

$e^{-\mathbb{E}_{P,x}\left[g_{N}\left(Y_{\tau_{AB}}\right)+\int_{0}^{\tau_{AB}}\sigma^{\top}v(Y_{s})dW_{s}+\frac{1}{2}\int_{0}^{\tau_{AB}}\|\sigma^{\top}v(Y_{s})\|^{2}ds\right]}\leq\mathbb{E}_{P,x}\left[e^{-g_{N}\left(X_{\tau_{AB}}\right)}\right].$
(C-16)

Since the expectation of the Ito stochastic integral is zero, i.e.,

$\mathbb{E}_{P,x}\left[\int_{0}^{\tau_{AB}}\sigma^{\top}v(Y_{s})\cdot
dW_{s}\right]=0,$

and the cost functional $C_{x}[N;v]$ defined in (48) is exactly

$C_{x}[N;v]\equiv\mathbb{E}_{P,x}\left[g_{N}\left(Y_{\tau_{AB}}\right)+\frac{1}{2}\int_{0}^{\tau_{AB}}\|\sigma^{\top}v(Y_{s})\|^{2}ds\right],$

and recalling (C-11) we get

$e^{-C_{x}[N;v]}\leq\mathbb{E}_{P,x}\left[e^{-g_{N}\left(X_{\tau_{AB}}\right)}\right]\equiv
q_{N}^{+}(x).$ (C-17)

Taking logarithms of the left- and right-hand side of (C-17) and multiplying
the result by $-1$ we obtain the following lower bound for the cost
functional: for any control $\sigma^{\top}v$ satisfying Novikov’s condition
(C-12),

$C_{x}[N;v]\geq-\log q_{N}^{+}(x).$ (C-18)

Since the admissible set $\mathcal{V}$ is closed, the bound (C-18) holds for
any admissible $v$. This means that for any $v\in\mathcal{V}$,

$c^{*}_{N}(x)\geq-\log q_{N}^{+}(x).$ (C-19)

Step 3. Derive the Hamilton-Jacobi-Bellman equation for the minimal cost
$c^{*}_{N}$. This upper bound will be derived via the Hamilton-Jacobi-Bellman
equation. Let $v(x)\in C^{1}(\Omega_{AB})$ be such that Novikov’s condition
(C-12) holds and let $h$ be a small positive number. Then for the process
$Y_{t}$ governed by the controlled SDE (47) with the control $\sigma^{\top}v$
we have the following upper bound:

$c_{N}^{\ast}(x)\leq\mathbb{E}_{P,x}\left[\frac{1}{2}\int_{0}^{h\land\tau_{AB}}\|\sigma^{T}(Y_{s})v(Y_{s})\|^{2}ds+c_{N}^{\ast}(Y_{h\land\tau_{AB}})\right].$
(C-20)

The equality is reached if $v$ is an optimal controller. We observe that if
$v\equiv 0$ than $C_{x}[N;0]=N(1-q^{+}(x))<N$ for $x\in\Omega_{AB}$.
Therefore, $c_{N}^{\ast}(x)\leq C_{x}[N;0]<N<\infty$. Therefore, we subtract
$c_{N}^{\ast}(x)$ from both sides of the inequality (C-20) and get

$\displaystyle 0$
$\displaystyle\leq\mathbb{E}_{P,x}\left[\frac{1}{2}\int_{0}^{h\land\tau_{AB}}\|\sigma^{T}(Y_{s})v(Y_{s})\|^{2}ds+c_{N}^{\ast}(Y_{h\land\tau_{AB}})-c_{N}^{\ast}(x)\right]$
(C-21)
$\displaystyle=\mathbb{E}_{P,x}\left[\frac{1}{2}\int_{0}^{h\land\tau_{AB}}\|\sigma^{T}(Y_{s})v(Y_{s})\|^{2}ds\right]+\mathbb{E}_{P,x}\left[c_{N}^{\ast}(Y_{h\land\tau_{AB}})\right]-c_{N}^{\ast}(x).$
(C-22)

Dividing by $h$ and letting $h\rightarrow 0$ we obtain:

$0\leq\frac{1}{2}\|\sigma^{\top}v(x)\|^{2}+\lim_{h\rightarrow
0}\frac{\mathbb{E}_{P,x}\left[c_{N}^{\ast}(Y_{h})\right]-c_{N}^{\ast}(x)}{h}.$
(C-23)

Here we took into account that the optimal control is continuously
differentiable in $\Omega_{AB}$ and $b$ and $\sigma$ are smooth. Therefore the
drift and the diffusion in (47) are finite and hence the probability that
$\tau_{AB}<h$ tends to zero as $h\rightarrow 0$. Furthermore, we note that

$\lim_{h\rightarrow
0}\frac{\mathbb{E}_{P,x}\left[c_{N}^{\ast}(Y_{h})\right]-c_{N}^{\ast}(x)}{h}=\lim_{h\rightarrow
0}\frac{\mathbb{E}_{P_{v},x}\left[c_{N}^{\ast}(Y_{h})\right]-c_{N}^{\ast}(x)}{h}=\mathcal{L}_{v}c^{*}_{N}(x),$
(C-24)

where $\mathcal{L}_{v}$ is the generator of the controlled process (47). The
equality (C-24) follows from the fact that

$\mathbb{E}_{P_{v},x}\left[c_{N}^{\ast}(Y_{h})\right]=\mathbb{E}_{P,x}\left[\frac{dP_{v}}{dP}c_{N}^{\ast}(Y_{h})\right],$

where

$\frac{dP_{v}}{dP}=\exp\left\\{-\int_{0}^{h}\sigma^{\top}v(Y_{s})\cdot
dW_{s}-\frac{1}{2}\int_{0}^{h}\|\sigma^{\top}v(Y_{s})\|^{2}ds\right\\}\rightarrow
1~{}~{}{\rm a.s.}~{}~{}{\rm as}~{}~{}h\rightarrow 0.$

Therefore, (C-23) is equivalent to

$\frac{1}{2}\|\sigma^{\top}v\|^{2}+\left[b+\sigma\sigma^{\top}v\right]\cdot\nabla
c^{*}_{N}+\frac{1}{2}{\sf tr}\left(\sigma\sigma^{\top}\nabla\nabla
c_{N}^{*}\right)\geq 0.$ (C-25)

Furthermore, the equality is reached if and only if the control $v$ is
optimal, i.e.,

$\inf_{v\in\mathcal{V}}\left[\frac{1}{2}\|\sigma^{\top}v\|^{2}+\sigma\sigma^{\top}v\cdot\nabla
c^{*}_{N}\right]+b\cdot\nabla c^{*}_{N}+\frac{1}{2}{\sf
tr}\left(\sigma\sigma^{\top}\nabla\nabla c_{N}^{*}\right)=0.$ (C-26)

The function in the square brackets in (C-26) is convex quadratic in $v$. To
minimize it, we take its gradient and set it to zero:

$\nabla_{v}\left[\sigma\sigma^{\top}v\cdot\nabla
c_{N}^{*}+\frac{1}{2}v^{\top}\sigma\sigma^{\top}v\right]=\sigma\sigma^{\top}\nabla
c_{N}^{*}+\sigma\sigma^{\top}v=0.$ (C-27)

Therefore, a minimizer $v_{N}^{*}$ must satisfy $\sigma\sigma^{\top}(\nabla
c^{*}_{N}+v_{N}^{*})=0$. Since columns of $\sigma$ are linearly independent,
this condition is equivalent to

$\sigma^{\top}(\nabla c_{N}^{*}+v_{N}^{*})=0\quad{\rm
or}\quad\sigma^{\top}v_{N}^{*}=-\sigma^{\top}\nabla c_{N}^{*}.$ (C-28)

Plugging this into (C-26) we obtain the following equation for the minimal
cost $c^{*}_{N}$:

$\displaystyle 0$ $\displaystyle=\frac{1}{2}{\sf
tr}\left(\sigma\sigma^{\top}\nabla\nabla c_{N}^{*}\right)+b\cdot\nabla
c_{N}^{*}-\sigma\sigma^{\top}\nabla c_{N}^{*}\cdot\nabla
c_{N}^{*}+\frac{1}{2}(\nabla c_{N}^{*})^{\top}\sigma\sigma^{\top}\nabla
c_{N}^{*}$ $\displaystyle=\frac{1}{2}{\sf
tr}\left(\sigma\sigma^{\top}\nabla\nabla c_{N}^{*}\right)+b\cdot\nabla
c_{N}^{*}-\frac{1}{2}(\nabla c_{N}^{*})^{\top}\sigma\sigma^{\top}\nabla
c_{N}^{*}.$ (C-29)

Step 4. Show that $c^{*}_{N}=-\log q_{N}^{+}(x)$ is the solution to the HJB
equation. Plugging $c_{N}^{*}=-\log q_{N}^{+}$ into (C-29) we get

$\displaystyle-$ $\displaystyle\frac{1}{2q_{N}^{+}}{\sf
tr}\left(\sigma\sigma^{\top}\nabla\nabla
q_{N}^{+}\right)+\frac{1}{2(q_{N}^{+})^{2}}\left\|\sigma^{\top}\nabla
q_{N}^{+}\right\|^{2}-b\cdot\nabla\log
q_{N}^{+}-\frac{1}{2(q_{N}^{+})^{2}}\left\|\sigma^{\top}\nabla
q_{N}^{+}\right\|^{2}=$ $\displaystyle-$
$\displaystyle\frac{1}{2q_{N}^{+}}{\sf
tr}\left(\sigma\sigma^{\top}\nabla\nabla q_{N}^{+}\right)-b\cdot\nabla\log
q_{N}^{+}=$ $\displaystyle-$
$\displaystyle\frac{1}{q_{N}^{+}}\mathcal{L}q_{N}^{+}=0.$

The last equality follows from the fact that $\mathcal{L}q_{N}^{+}=0$ in
$\Omega_{AB}$. The boundary conditions for $c_{N}^{*}=-\log q_{N}^{+}$ are
readily checked: $-\log q_{N}^{+}=N$ on $\partial A$, $-\log q_{N}^{+}=0$ on
$\partial B$, and

$\frac{\partial}{\partial\hat{n}}\left(-\log
q^{+}_{N}\right)=-\frac{1}{q^{+}_{N}}\frac{\partial
q_{N}^{+}}{\partial\hat{n}}=0,\quad x\in\partial\Omega.$

The optimal control associated with $c^{*}_{N}=-\log q_{N}^{+}(x)$ given by
(C-28) is

$\sigma^{\top}v_{N}^{*}=-\sigma^{\top}\nabla c_{N}^{*}=\sigma^{\top}\nabla\log
q^{+}_{N}.$ (C-30)

Step 5. Show that the control $\sigma^{\top}v^{*}_{N}=\sigma^{\top}\log
q^{+}_{N}$ is admissible. Equation (C-10) implies that

$\sigma^{\top}v^{*}_{N}=-\sigma^{\top}\nabla c^{*}_{N}=\sigma^{\top}\nabla\log
q_{N}^{+}=\sigma^{\top}\left[\frac{\nabla
q^{+}\left(1-e^{-N}\right)}{q^{+}\left(1-e^{-N}\right)+e^{-N}}\right].$ (C-31)

Hence

$\left\|\sigma^{\top}v^{*}_{N}\right\|\leq
e^{N}\max_{x\in\Omega_{AB}}\left\|\nabla q^{+}_{N}(x)\right\|<\infty.$ (C-32)

The stopping time $\tau_{AB}<\infty$ a.s. as the system is ergodic and the
domain is compact. Therefore

$\frac{1}{2}\int_{0}^{\tau_{AB}}\|\sigma^{\top}v^{*}_{N}(Y_{s})\|^{2}ds<\infty~{}~{}{\rm
a.s.}$

and hence

$\mathbb{E}_{P,x}\left[e^{\frac{1}{2}\int_{0}^{\tau_{AB}}\|\sigma^{\top}v^{*}_{N}(Y_{s})\|^{2}ds}\right]<\infty,$
(C-33)

i.e. $v_{N}^{*}$ is admissible.

Step 6. Take the limit $N\rightarrow\infty$. Letting $N\rightarrow\infty$ in
(C-19) and taking into account the explicit expression (C-10) for
$q^{+}_{N}(x)$ we conclude that

$c^{*}(x)\geq-\log q^{+}(x).$ (C-34)

On the other hand, as we have shown in Step 3, (C-19) is actually an equality,
and the corresponding optimal control satisfies

$\sigma^{\top}v^{*}_{N}=-\sigma^{\top}\nabla c^{*}_{N}=\sigma^{\top}\nabla\log
q_{N}^{+}=\sigma^{\top}\left[\frac{\nabla
q^{+}\left(1-e^{-N}\right)}{q^{+}\left(1-e^{-N}\right)+e^{-N}}\right].$ (C-35)

Taking limit $N\rightarrow\infty$ we obtain

$\sigma^{\top}v^{*}=\sigma^{\top}\frac{\nabla
q^{+}}{q^{+}}=\sigma^{\top}\nabla\log q^{+}.$ (C-36)

Since the admissible set $\mathcal{V}$ is closed, $v^{*}=\nabla\log
q^{+}\in\mathcal{V}$. One can readily check that the corresponding solution
the Hamilton-Jacobi-Bellman equation (C-26) with the boundary conditions
$c^{*}=+\infty$ on $\partial A$, $c^{*}=0$ on $\partial B$, and
$\tfrac{\partial c^{*}}{\partial\hat{n}}=0$ on $\partial\Omega$ is

$c^{*}(x)=-\log q^{+}(x).$ (C-37)

This completes the proof of Theorem (3.1). ∎

##  Appendix D Errors due to model reduction: an example

We will illustrate the error due to model reduction in the transition rate
$\nu_{AB}$ as well as in the probabilities $\rho_{A}$ and $\rho_{AB}$ on the
example used in [43]. A system is evolving according to the overdamped
Langevin dynamics (6) with the potential given by

$V(x,y)=(x^{2}-1)^{2}+\epsilon^{-1}(y+x^{2}-1)^{2},$ (D-38)

where $\epsilon$ is a small parameter. The second term in (D-38) effectively
restricts the dynamics to a small neighborhood of the parabola $y=1-x^{2}$. It
is shown in [43] that $x$ is a suboptimal choice of a collective variable
because the gradient of $x$ with respect to $(x,y)$ is not orthogonal to the
normal vector to the manifold near which the dynamics live.

Let us consider the signed arclength parameter along the parabola $y=1-x^{2}$

$s(x)=\int_{0}^{x}\sqrt{1+4z^{2}}dz=\frac{1}{2}x\sqrt{1+4x^{2}}+\frac{1}{4}\log(2x+\sqrt{1+4x^{2}})$
(D-39)

as a collective variable. The function $s(x)$ is monotone and hence
invertible. In the limit $\epsilon\rightarrow 0$, the dynamics are one-
dimensional and governed by

$ds=-\frac{dV_{0}(x(s))}{ds}dt+\sqrt{2\beta^{-1}}dw,$ (D-40)

where $V_{0}:=(x^{2}-1)^{2}$. We set $\beta=3$ choose the sets $A$ and $B$ as
in [43]:

$A=\\{x<a\\},\quad B=\\{x>b\\},\quad a=-0.5,\quad b=0.5.$ (D-41)

We calculate the committor $q(s)$ and $\tilde{q}(x)$ using the exact formula
for the one-dimensional case:

$q(s)=\frac{\int_{s(a)}^{s}e^{\beta
V_{0}(x(s^{\prime}))}ds^{\prime}}{\int_{s(a)}^{s(b)}e^{\beta
V_{0}(x(s^{\prime}))}ds^{\prime}},\quad\tilde{q}(x)=\frac{\int_{a}^{x}e^{\beta
V_{0}(x^{\prime})}dx^{\prime}}{\int_{a}^{b}e^{\beta
V_{0}(x^{\prime})}dx^{\prime}}.$ (D-42)

Using the inverse of $s(x)$, we obtain $\tilde{q}(x(s))$. The plots of $q(s)$
and $\tilde{q}(x(s))$ and their derivatives in $s$ are displayed in Fig.
17(left). It is evident that the difference between their derivatives is
notably larger.

Figure 17: An illustration to D.

Next, we use $q(s)$ and $\tilde{q}(x)$ to calculate the transition rate
$\nu_{AB}$ from $A$ to $B$ via (33) and the probability $\rho_{AB}$ via (27).
The notation with tilde will indicate the results obtained using $x$ as a
collective variable. We get:

$\displaystyle\nu_{AB}$ $\displaystyle=9.46\cdot 10^{-3},$
$\displaystyle\rho_{AB}=9.47\cdot 10^{-3},$ $\displaystyle\tilde{\nu}_{AB}$
$\displaystyle=2.19\cdot 10^{-2},$ $\displaystyle\tilde{\rho}_{AB}=7.98\cdot
10^{-3}.$

The transition rate estimated using $x$ as a collective variable exceeds to
true rate by the factor of approximately 2.3, while the error in the estimate
of the probability to be reactive is about 16%.

##  Appendix E Robustness of the crossover time: an example.

In this appendix, we will examine how the quality of the committor estimate
affects the crossover time in the controlled dynamics on the example of the
bistable Duffing oscillator (83) with $\epsilon=0.05$. The controlled dynamics
of the Duffing oscillator are governed by SDE (85) where $q^{+}_{\sf nn}$ is
the estimate to the forward committor computed using PINN.

The PINN committor solver computes the committor via an optimization process
called _training the neural network_. The number of training steps is measured
in epochs222One epoch comprises the number of iterations necessary to use all
training data $N_{\sf train}$ one time. For example, if the optimization
algorithm is deterministic and all training data are used for computing the
gradient of the objective function in each iteration, then one epoch is equal
to one iteration. If an optimizer is stochastic and a subset, a _batch_ , of
$N_{\sf batch}$ training data is used at each iteration to evaluate the
direction of the step, then one epoch consists of ${\sf round}(N_{\sf
train}/N_{\sf batch})$ iterations. .  Stopping the training process too early
results in a rough approximation to the forward committor. The results
reported in Section 6.3 are obtained as a result of training the neural
network for 500 epochs. Hence, we form a sequence of approximations to the
forward committor by evaluating the solution model after 100, 125, 150, 300,
and 500 epochs of training visualized in the figures in Table 10 using dashed
contour plots. The solution becomes progressively closer to the final solution
$q_{\sf nn}^{+}$ evaluated at 500 epochs, which, in turn, in close to the FEM
solution depicted using solid contour plots. The discrepancies MAD and RMSD
between the FEM forward committor and the approximations to it progressively
shrink.

For each of these “undertrained” solutions, we evaluate the expected crossover
time $\mathbb{E}[\tau_{AB}]$ by averaging the crossover times of 250
transition trajectories governed by the controlled SDE (85) with the
corresponding “undertrained” forward committor. The results are shown in the
last column of Table 10.

|  wMAD |  wRMSD |  Visual comparison of $q^{+}$ of PINNs and FEM |  $\mathbb{E}[\tau_{AB}]$
---|---|---|---|---
Case 0 (epoch 500) |  1.3e-2 |  2.0e-2 |  |  7.34 $\pm$ 0.33
Case 1 (epoch 300) |  4.9e-2 |  6.1e-2 |  |  7.88 $\pm$ 0.36
Case 2 (epoch 150) |  6.0e-2 |  7.7e-2 |  |  8.48 $\pm$ 0.84
Case 3 (epoch 125) |  7.9e-2 |  10.2e-2 |  |  12.18 $\pm$ 2.49
Case 4 (epoch 100) |  13.6e-2 |  16.7e-2 |  FEM(solid), PINN(dashed) |  48.86 $\pm$ 5.71
Table 10: Estimates of wMAD, wRMSD and $\mathbb{E}[\tau_{AB}]$ computed for
five different approximations to the forward committor (dashed contours) of
the bistable Duffing oscillator (83) with $\epsilon=0.05$ obtained at various
stages of training the neural network-based solution model with the PINN loss
function. The FEM forward committor (solid contours) and the corresponding
$\mathbb{E}[\tau_{AB}]=\mathbf{7.48\pm 0.49}$ are treated as the ground truth.

Table 10 shows that as the accuracy of the committor decreases, the expected
crossover time $\mathbb{E}[\tau_{AB}]$ increases. Nonetheless, the relative
increment in $\mathbb{E}[\tau_{AB}]$ is notably smaller compared to the
discrepancy in the committors as evident from Figure 18.

Figure 18: Test case: the bistable Duffing oscillator with $\epsilon=0.05$.
The growth of the relative error in the expected crossover time
$\mathbb{E}[\tau_{AB}]$ versus the relative growth of the weighted mean
absolute difference (wMAD) of the “undertrained” approximations to the forward
committors and the FEM committor. The approximate committors $q_{\sf nn}^{+}$
are obtained by evaluating the solution model after $100$, $125$, $150$,
$300$, and $500$ epochs of training.

In summary, this investigation suggests that the estimate of the expected
crossover time obtained by means of sampling controlled trajectories remains
reasonably accurate even if the approximation to the forward committor is
rough.

##  Appendix F Finite element method for the committor problem

###  F.1 Time-reversible dynamics

If the governing SDE is time-reversible as it is in the case of the overdamped
Langevin dynamics (6) or the overdamped Langevin dynamics in collective
variables (9), the committor problem (19) with the generator (23) is self-
adjoint. In this case, we proceed in the standard way detailed in [72]. First,
we decompose the committor $q$ into $q=q_{1}+q_{0}$ where is a prescribed
function such that $q_{1}=1$ on $\partial B$ and $q_{1}=0$ outside a small
neighborhood of $\partial B$ and $q_{0}$ needs to be found. The boundary value
problem for $q_{0}$ is

$\displaystyle\mathcal{L}q_{0}$ $\displaystyle=-\mathcal{L}q_{1},\quad
x\in\Omega_{AB},$ (F-43) $\displaystyle q_{0}$ $\displaystyle=0,\quad
x\in\partial A\cup\partial B,$ (F-44) $\displaystyle\frac{\partial
q_{0}}{\partial\hat{n}}$ $\displaystyle=0,\quad x\in\partial\Omega.$ (F-45)

Second, we multiply multiply (F-43) by a test function $w\in
H^{1}_{0}(\bar{\Omega})$ where the subscript $0$ means that $w=0$ on $\partial
A\cup\partial B$, integrate over $\Omega_{AB}$ and apply the generalized
divergence theorem to both parts. The result is the following integral
equation for $q_{0}$ that must hold for all $w\in C^{1}_{0}(\bar{\Omega})$:

$\int_{\Omega_{AB}}e^{-\beta F(x)}w(x)^{\top}M(x)\nabla
q_{0}=-\int_{\Omega_{AB}}e^{-\beta F(x)}w(x)^{\top}M(x)\nabla q_{1}.$ (F-46)

Next, we triangulate $\Omega_{AB}$ and denote the associated finite element
space by $S_{h}$ with the standard piecewise-linear basis
$\\{v_{j}(x)\\}_{i\in\mathcal{I}}$ where $\mathcal{I}$ is the set of vertices
of the triangles [72]. The subset of vertices that do not belong to $\partial
A\cup\partial B$ is denoted by $\mathcal{I}_{\rm free}$. We choose $q_{1}\in
S_{h}$ so that $q_{1}=1$ only at the nodes lying on $\partial B$ and $q_{1}=0$
at all other nodes. We seek the finite element solution for $q_{0}$ of the
form

$q_{0}=\sum_{k\in\mathcal{I}_{\rm free}}v_{k}(x)(q_{0})_{k}$ (F-47)

where the vector $\\{(q_{0})_{k}\\}_{k\in\mathcal{I}_{free}}$ is the solution
to the linear system

$\sum_{k\in\mathcal{I}_{\rm free}}{\sf A}_{jk}(q_{0})_{k}=-\sum_{k\in\partial
B}{\sf A}_{jk}$ (F-48)

with the matrix elements ${\sf A}_{jk}$ given by

${\sf A}_{jk}=\int_{\Omega_{AB}}e^{-\beta F(x)}\nabla
v_{j}(x)^{\top}M(x)\nabla v_{k}(x)dx.$ (F-49)

The integral in (F-49) is the sum of the integrals over all triangles. In each
triangle, the gradients of the basis functions are constant, and $F(x)$ and
$M(x)$ are approximated by their values at the center of mass. Finally, the
finite element solution $q_{\sf fem}$ is found at the sum $q_{\sf
fem}=q_{0}+q_{1}$.

###  F.2 The Langevin dynamics

For the Langevin dynamics (10), the committor problem (19) with the generator
(24) is hypoelliptic, and the application of the finite element method (FEM)
requires care. We design a FEM solver for this case motivated by the article
by Morton on FEM for non-self-adjoint problems [73] that suggests to make the
problem as close to self-adjoint as possible. Since in the case of Langevin
dynamics FEM is practical only if the space $(x,p)$ is two-dimensional, $x$
and $p$ will be one-dimensional in the presentation below.

As in F.1, we start by decompositing $q^{+}$ into $q^{+}=q^{+}_{0}+q^{+}_{1}$.
The boundary value problem for $q^{+}_{0}$ is of the form (F-43)–(F-45) except
for $x$ is replaced with $(x,p)$. Then we multiply the PDE for $q^{+}_{0}$ by
$\epsilon^{-1}\exp\left(-\tfrac{p^{2}}{2\epsilon}\right)$ and get:

$\frac{e^{-\frac{p^{2}}{2\epsilon}}}{\epsilon}\mathcal{L}q^{+}_{0}=\frac{e^{-\frac{p^{2}}{2\epsilon}}}{\epsilon}\left[\frac{p}{m}\frac{dq_{0}}{dx}-V^{\prime}(x)\frac{dq^{+}_{0}}{dp}\right]+\gamma\frac{d}{dp}\left(e^{-\frac{p^{2}}{2\epsilon}}\frac{dq^{+}_{0}}{dp}\right)=-\frac{e^{-\frac{p^{2}}{2\epsilon}}}{\epsilon}\mathcal{L}q^{+}_{1}.$
(F-50)

We denote $(x,p)$ by $z$ and the gradient with respect to $z$ by $\nabla$ and
rewrite (F-50) in a matrix form:

$\gamma\nabla\cdot\left(e^{-\frac{p^{2}}{2\epsilon}}\left[\begin{array}[]{cc}0&0\\\
0&m\end{array}\right]\nabla
q^{+}_{0}\right)+\frac{e^{-\frac{p^{2}}{2\epsilon}}}{\epsilon}\left[\begin{array}[]{c}p/m\\\
V^{\prime}(x)\end{array}\right]\cdot\nabla
q^{+}_{0}=-\frac{e^{-\frac{p^{2}}{2\epsilon}}}{\epsilon}\mathcal{L}q^{+}_{1}.$
(F-51)

Then we follow the steps in F.1. We multiply (F-51) by a test function
$w(z)\in H_{0}^{1}(\bar{\Omega}_{AB})$, integrate over $\Omega_{AB}$ and apply
the generalized divergence theorem. This results in the following integral
equation for $q^{+}_{0}$ that must hold for all $w(z)\in
H_{0}^{1}(\bar{\Omega}_{AB})$:

$\displaystyle\gamma\int_{\Omega_{AB}}e^{-\frac{p^{2}}{2\epsilon}}\nabla
w^{\top}\left[\begin{array}[]{cc}0&0\\\ 0&m\end{array}\right]\nabla
q^{+}_{0}dx-\frac{e^{-\frac{p^{2}}{2\epsilon}}}{\epsilon}\int_{\Omega_{AB}}w\left[\begin{array}[]{c}p/m\\\
V^{\prime}(x)\end{array}\right]\cdot\nabla q^{+}_{0}$ (F-56) $\displaystyle=$
$\displaystyle-\gamma\int_{\Omega_{AB}}e^{-\frac{p^{2}}{2\epsilon}}\nabla
w^{\top}\left[\begin{array}[]{cc}0&0\\\ 0&m\end{array}\right]\nabla
q^{+}_{1}dx+\frac{e^{-\frac{p^{2}}{2\epsilon}}}{\epsilon}\int_{\Omega_{AB}}w\left[\begin{array}[]{c}p/m\\\
V^{\prime}(x)\end{array}\right]\cdot\nabla q^{+}_{1}dx.$ (F-61)

Then we triangulate $\Omega_{AB}$, introduce the standard FEM basis
$\\{v_{i}\\}_{i\in\mathcal{I}}$ in $S_{h}$, and represent $q^{+}_{0}(z)$ as a
linear combination of the basis functions associated with the nodes not in
$\partial A\cup\partial B$, and obtain the following linear system for the
coefficients $(q^{+}_{0})_{k}$:

$\sum_{k\in\mathcal{I}_{\rm free}}{\sf
A}_{jk}(q^{+}_{0})_{k}-\sum_{k\in\mathcal{I}_{\rm free}}{\sf
B}_{jk}(q^{+}_{0})_{k}=-\sum_{k\in\partial B}{\sf A}_{jk}+\sum_{k\in\partial
B}{\sf B}_{jk}.$ (F-62)

The matrix elements in (F-62) are given by

${\sf A}_{jk}=\gamma\int_{\Omega_{AB}}e^{-\frac{p^{2}}{2\epsilon}}\nabla
v_{j}^{\top}\left[\begin{array}[]{cc}0&0\\\ 0&m\end{array}\right]\nabla
v_{k}dx,\quad{\sf
B}_{jk}=\frac{e^{-\frac{p^{2}}{2\epsilon}}}{\epsilon}\int_{\Omega_{AB}}v_{j}\left[\begin{array}[]{c}p/m\\\
V^{\prime}(x)\end{array}\right]\cdot\nabla v_{k}dx.$ (F-63)

Computing the integrals in (F-63) over each triangle, all nonlinear functions
are approximated by their values at the centers of mass of the triangle.
Finally, $q^{+}_{\sf fem}=q^{+}_{0}+q^{+}_{1}$. The backward committor is
readily found by $q^{-}_{\sf fem}(x,p)=1-q^{+}_{\sf fem}(x,-p)$.

##  Appendix G Confidence interval

To compute confidence intervals for simulated transition time from $A$ to $B$,
we first compute mean $\bar{\tau}_{AB}$ and standard error of the mean
$\text{se}_{\tau}$ using scipy.stats.sem, which is the sample standard
deviation divided by square root of the sample size:
$\text{se}_{\tau}=\frac{\sigma}{\sqrt{n}}$. The confidence interval then is
obtained using t distribution:

$\text{confidence
interval}=[\bar{\tau}_{AB}-\text{se}_{\tau}t(\alpha),\bar{\tau}_{AB}+\text{se}_{\tau}t(\alpha)]$
(G-64)

where $t(\alpha)$ satisfies $\int_{t(\alpha)}^{\infty}f(x)dx=\alpha$, and can
be found using scipy.stats.t.ppf. For all examples, $95\%$ confidence
intervals are used.

## References

  * [1] B. J. Zhang, T. Sahai, Y. M. Marzouk, A koopman framework for rare event simulation in stochastic differential equations, Journal of Computational Physics 456 (2022) 111025.  doi:https://doi.org/10.1016/j.jcp.2022.111025.
  * [2] J. Lu, J. Nolen, Reactive trajectories and the transition path process, Probability Theory and Related Fileds (2015).
URL https://link.springer.com/article/10.1007/s00440-014-0547-y

  * [3] Y. Gao, T. Li, X. Li, J.-G. Liu, Transition path theory for langevin dynamics on manifolds: Optimal control and data-driven solver, Multiscale Modeling & Simulation 21 (1) (2023) 1–33.  doi:10.1137/21M1437883.
  * [4] T. Breunung, B. Balachandran, Noise color influence on escape times in nonlinear oscillators - experimental and numerical results, Theoretical and Applied Mechanics Letters (2022) 100420doi:https://doi.org/10.1016/j.taml.2022.100420.
URL https://www.sciencedirect.com/science/article/pii/S2095034922001003

  * [5] C. Dellago, P. Bolhuis, F. Csajka, D. Chandler, Transition path sampling and the calculation of rate constants, J. Chem. Phys. 108 (1997) 1964–1977.  doi:10.1063/1.475562.
  * [6] C. Dellago, P. G. Bolhuis, D. Chandler, On the calculation of reaction rate constants in the transition path ensemble, The Journal of Chemical Physics 110 (14) (10 1998).  doi:https://doi.org/10.1063/1.478569.
  * [7] P. G. Bolhuis, D. Chandler, C. Dellago, P. L. Geissler, Transition path sampling: Throwing ropes over rough mountain passes, in the dark, Annual Review of Physical Chemistry 53 (1) (2002) 291–318.  doi:10.1146/annurev.physchem.53.082301.113146.
  * [8] R. Elber, Milestoning: An efficient approach for atomically detailed simulations of kinetics in biophysics, Annual Review of Biophysics 49 (1) (2020) 69–85.
  * [9] D. Aristoff, J. Copperman, G. Simpson, R. J. Webber, D. M. Zuckerman, Weighted ensemble: Recent mathematical developments, The Journal of Chemical Physics 158 (2023) 014108.
  * [10] F. Cérou, A. Guyader, Adaptive multilevel splitting for rare event analysis, Stochastic Analysis and Applications 25 (2) (2007) 417–443.  doi:10.1080/07362990601139628.
  * [11] F. Cérou, A. Guyader, M. Rousset, Adaptive multilevel splitting: Historical perspective and recent results, Chaos: An Interdisciplinary Journal of Nonlinear Science 29 (4) (2019) 043108.  doi:10.1063/1.5082247.
  * [12] A. Laio, M. Parrinello, Escaping free-energy minima, Proceedings of the National Academy of Sciences 99 (20) (2002) 12562–12566.  doi:10.1073/pnas.202427399.
  * [13] A. Barducci, G. Bussi, M. Parrinello, Well-tempered metadynamics: a smoothly converging and tunable free-energy method, Physical review letters 100 (2) (2008) 020603.
  * [14] Y. Wang, J. Marcelo, L. Ribeiro, P. Tiwary, Machine learning approaches for analyzing and enhancing molecular dynamics simulations, Current Opinion in Structural Biology 61 (2020) 139–145.
  * [15] Z. Belkacemi, P. Gkeka, T. Lelièvre, G. Stoltz, Chasing collective variables using autoencoders and biased trajectories, Journal of chemical theory and computation 18 (2022) 59–78.
  * [16] H. Jonsson, G. Mills, J. K. W., Nudged Elastic Band Method for Finding Minimum Energy Paths of Transitions, World Scientific, Singapore, 1998, Ch. 16, pp. 385–404.
  * [17] W. E, W. Ren, E. Vanden-Eijnden, String method for the study of rare events, Phys. Rev. B 66 (2002) 052301.  doi:10.1103/PhysRevB.66.052301.
URL https://link.aps.org/doi/10.1103/PhysRevB.66.052301

  * [18] W. E, W. Ren, E. Vanden-Eijnden, Minimum action method for the study of rare events, Communications on Pure and Applied Mathematics 57 (2004).
URL https://api.semanticscholar.org/CorpusID:123082141

  * [19] L. Maragliano, A. Fischer, E. Vanden-Eijnden, G. Ciccotti, String method in collective variables: Minimum free energy paths and isocommittor surfaces, The Journal of chemical physics 125 (2) (2006) 024106.
  * [20] W. E, W. Ren, E. Vanden-Eijnden, Simplified and improved string method for computing the minimum energy paths in barrier-crossing events., The Journal of chemical physics 126 16 (2007) 164103.
  * [21] E. Vanden-Eijnden, M. Heymann, The geometric minimum action method for computing minimum energy paths, The Journal of Chemical Physics 128 (6) (2008) 061103.  arXiv:https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.2833040/13625043/061103\\_1\\_online.pdf, doi:10.1063/1.2833040.
URL https://doi.org/10.1063/1.2833040

  * [22] S. Beri, R. Mannella, D. G. Luchinsky, A. N. Silchenko, P. V. E. McClintock, Solution of the boundary value problem for optimal escape in continuous stochastic systems and maps, Phys. Rev. E 72 (2005) 036131.  doi:10.1103/PhysRevE.72.036131.
URL https://link.aps.org/doi/10.1103/PhysRevE.72.036131

  * [23] L. Cilenti, M. Cameron, B. Balachandran, Most probable escape paths in periodically driven nonlinear oscillators, Chaos 32 (2022) 083140.
  * [24] S. Yang, S. F. Potter, M. K. Cameron, Computing the quasipotential for nongradient sdes in 3d, Journal of Computational Physics 379 (2019) 325–350.
  * [25] M. Cameron, S. Yang, Computing the quasipotential for highly dissipative and chaotic sdes an application to stochastic lorenz-63, Communications in Applied Mathematics and Computational Science 14 (2) (2019) 207–246.