A primer on noise-induced transitions in applied dynamical systemsE. Forgoston
and R.O. Moore

# A primer on noise-induced transitions in applied dynamical systems

Eric Forgoston Department of Mathematical Sciences, Montclair State
University, Montclair, NJ 07043, USA<EMAIL_ADDRESS>Richard
O. Moore Department of Mathematical Sciences, New Jersey Institute of
Technology, Newark, NJ 07102, USA<EMAIL_ADDRESS>
###### Abstract

Noise plays a fundamental role in a wide variety of physical and biological
dynamical systems. It can arise from an external forcing or due to random
dynamics internal to the system. It is well established that even weak noise
can result in large behavioral changes such as transitions between or escapes
from quasi-stable states. These transitions can correspond to critical events
such as failures or extinctions that make them essential phenomena to
understand and quantify, despite the fact that their occurrence is rare. This
article will provide an overview of the theory underlying the dynamics of rare
events for stochastic models along with some example applications.

###### keywords:

Dynamical systems, escape, extinction, mean exit time, large deviation theory,
mode-locked lasers, noise-induced transitions, optimal path, rare events,
stochasticity, epidemiology, bit error ratios.

##  1 Introduction

Recent years have seen a dramatic rise in the use of stochastic systems to
model a wide variety of important physical and biological phenomena. Studies
of sub-cellular processes and tissue dynamics [85], large-scale population
dynamics [74], genetic switching [7], magnetic devices [57], optical
devices[51, 81, 75, 20], Josephson junctions [37], fluid mechanics[34, 33,
48], weather and climate[87], and geosciences[82] have included numerous
investigations into how noise affects physical and biological phenomena at a
wide variety of scales. One often sees rare transition events in these systems
that are induced by noise which may be internal or external to the system.
These noise-induced rare events may be associated with a desirable outcome,
such as the extinction of an infectious disease outbreak [6, 11] or
eradication of a pest [71], or an undesirable outcome, such as the sudden
clustering of cancerous cells [55], the outbreak of an infectious disease
[72], or a bit error in a digital communication system [29]. Another important
class of systems are those that exhibit rare events due to intermittently
activated finite-time instabilities [64, 30]. These types of rare events
include freak water waves in the ocean [17] and ship rolling and capsizing
[65]. The need to understand these phenomena has fueled a demand for methods
that quantify the impact of noise on complex systems.

In these stochastic systems, noise can affect the system in a variety of ways.
Assessing the full impact of noise is rarely possible due to the “curse of
dimensionality”. Analysis and computations must often concentrate on the most
important noise-induced events, which include spontaneous switching between
coexisting stable states, critical failure of a system, and nucleation of
coherent structures. One important feature of interest when studying noise-
induced transitions is the optimal transition pathway of escape from a
metastable state either to another metastable state or to a stable (absorbing)
state. We define the optimal path as the path that is most likely to occur
among all possible paths, recognizing that this path may not be unique or even
exist (i.e., belong to the set of admissible paths). For systems out of
equilibrium such that detailed balance is not satisfied, knowledge of the
optimal path is required for the computation of the mean switching time
between states or the mean time to exit.

This article provides a tutorial overview of noise-induced rare events. By
providing theoretical and numerical background along with several example
applications of recently developed methods to compute the most likely noise
pathways to critical events, we hope to enable the reader to quickly apply the
methods to original applications of interest.

###  1.1 Noise sources

Stochasticity manifests itself as either external or internal noise. External
noise comes from a source outside the system being considered (e.g. population
growth under the influence of climatic effects, or a random signal fed into a
transmission line), and often is modeled by replacing an external parameter
with a random process. Internal noise is inherent to the system itself and is
caused by the random interactions of discrete particles (e.g. individuals in a
population, or molecules undergoing a chemical reaction) [44, 39, 86, 84].
Both types of noise can lead to a rare switching event between metastable
states or a rare escape event from a metastable state.

Figure 1(a) shows a snapshot in time of a particle moving in a quartic
potential well under the influence of external noise. Figure 1(b) shows a
single realization of the particle position in time. One sees multiple rare
stochastic switches from one well to another.

Figure 1: (a)Particle moving in a potential well under the influence of
external noise. (b)Stochastic realization of the particle position in time.

Figure 2 shows a single realization of the population size in time for an
Allee effect [2] problem in population biology. In this example, the noise is
internal and arises from the random interactions of individuals in a
population. One sees the population fluctuating about the carrying capacity
for a long period of time until the rare extinction event occurs.

Figure 2: Stochastic realization of the population size in time where the
local dynamics of the population exhibit the Allee effect. The green, solid
line denotes the carrying capacity, while the red, dashed line denotes the
Allee threshold.

There are many possible escape/switching paths, but there is a path along
which switching or escape is most likely to occur. We call this most likely
path of escape or switching the optimal path. It is of great importance in a
variety of applied problems to determine this optimal path since knowledge of
the path then enables the determination of the mean time to escape from a
metastable state or to switch from one metastable state to another metastable
state.

Mathematically, the effect of external noise is often described using a
Langevin equation or the associated Fokker-Planck equation (though the
dynamics of external noise may sometimes be described by a master equation
[78]). Feynman famously pointed out that each noise realization corresponds to
a particular trajectory of the system, and therefore the probability density
of realizations of trajectories is determined by the probability density of
noise realizations [31]. This idea can be used to formulate a variational
problem to find the optimal path that ultimately reduces to considering
trajectories of an auxiliary Hamiltonian dynamical system. One can solve for
the Hamiltonian dynamics, either analytically or numerically, for the most
probable (i.e., optimal) path of escape or switching [36, 25].

The effect of internal noise due to the random interactions of individuals
within the system is described mathematically using a master equation. The
master equation is a large, or even infinite, set of differential equations
that cannot in general be solved analytically. However, the optimal
fluctuation path is found to solve an eikonal equation obtained via Wentzel-
Kramers-Brillouin-Jeffreys (WKBJ) analysis of the master equation in the limit
of weak internal noise [60, 38, 25, 28, 54, 6, 32, 80, 71, 11, 72]. As in the
case of external noise, the optimal path is found by transforming the original
stochastic problem into a new deterministic system described by a Hamiltonian
$\mathcal{H}(x,\lambda)$. The dimension of the Hamiltonian is twice the
dimension of the original system due to the conjugate momentum variables
$\lambda$. An optimal path starting at a metastable state is identified with
$\mathcal{H}(x,\lambda)=0$; therefore, the method amounts to finding a zero-
energy trajectory of an effective mechanical system, and at least one of the
solutions to the zero-energy Hamiltonian is the optimal path. There may be
other escape, switching, or extinction paths, but the optimal path is the path
that maximizes the probability of escape, switching, or extinction [10].

###  1.2 Infinite-dimensional equations and large deviation theory

Even though it is often straightforward to extend the framework discussed in
this article to infinite-dimensional systems with states defined on Hilbert
spaces [18], we generally restrict our discussion to finite-dimensional
examples. The exception to this is the formal introduction of stochastic PDEs
used to model mode-locked lasers and communication systems introduced in
Sections 3.1 and 3.2. The simulations performed in these sections,
particularly those using importance sampling applied to (95), effectively
limited the bandwidth of noise that was otherwise taken to be delta-correlated
in both time and space. To the authors’ knowledge no proof of well-posedness
exists for this equation, and indeed when $c_{2}=0$ in (95), there is evidence
to suggest that the stochastic PDE is ill-posed [14, 15].

Finally, we note that large deviation results for stochastic PDEs are becoming
more common in the literature, particularly in the context of parabolic
equations with damping at high wavenumbers [46, 45] Large deviation principles
for the nonlinear Schrödinger equation with additive and multiplicative noise
have been shown when the noise is colored in space [19, 40, 41].

##  2 Transitions, exits, and extinctions

###  2.1 Mean first exit time due to external noise

The effect of noise on a physical system can be described in a variety of
ways, all of which aim to capture how a single deterministic solution has been
transformed by noise into a distribution of solutions and what the
characteristics are of that distribution. This can take the form of a boundary
value problem with random data or random coefficients, as one might find in
the forward (scattering) problem associated with imaging of the earth or of
living tissue, or it can take the form of an initial value problem with random
data, as one might find in numerical weather prediction. We will focus here on
initial value problems with time-dependent random forcing, specifically on
stochastic dynamical systems described by a well-defined state evolving in
time, where the evolution is probabilistic due to the presence of noise. This
section provides some definitions to set the context of what follows. Although
we will consider a one-dimensional (1D) framework, one can naturally extend
the results to the multi-dimensional case.

We begin by considering the stochastic differential equation (SDE) described
by

(1) $\dot{x}=f(x)+\eta(t),$

where $\eta(t)$ expresses an additive randomly distributed noise term
influencing the state $x(t)$. Whether resulting from truly random (e.g.,
quantum) phenomena or simply an assortment of physical phenomena not included
in the model for reasons of scale or uncertainty, $\eta(t)$ can often be
equated with the formal time derivative of a Brownian motion, resulting in
paths $x(t)$ that are continuous and Markovian (i.e., “memoryless”). Since
these paths, referred to as diffusion processes, are generally not
differentiable, it is helpful to consider the statistics of finite increments
$\Delta\eta=\eta(t+\Delta t)-\eta(t)$. These increments are drawn from a
Gaussian distribution with mean and variance given by

(2)
${\mathbb{E}}[\Delta\eta]=0\quad\mbox{and}\quad{\mathbb{E}}[(\Delta\eta)^{2}]=D\Delta
t,$

where ${\mathbb{E}}[\cdot]$ denotes the expectation or expected value of a
quantity with respect to its distribution. A natural way to think of (1) is
therefore as the limit of the Euler-Maruyama method [73], which is simply the
expression of the comments above in the form of a numerical integration method
to obtain $\Delta x=x(t+\Delta t)-x(t)$:

(3) $\Delta x=f(x(t))\Delta t+\Delta\eta.$

Then (1) is understood as the limiting process as $\Delta t\rightarrow 0$,
such that

(4)
${\mathbb{E}}[\eta(t)]=0\quad\mbox{and}\quad{\mathbb{E}}[\eta(t)\eta(t^{\prime})]=D\delta(t-t^{\prime}),$

where $D$ is the noise intensity. Details of the Euler-Maruyama method can be
found in Appendix A.

Some of the SDEs considered in the following sections have state-dependent
functions multiplying the random term, i.e.,

(5) $\dot{x}=f(x)+\sigma(x)\eta(t).$

This expression is subject to interpretation depending on how one determines
the limit described above. Two standard interpretations of this stochastic
integral predominate, including the prepoint or Itô interpretation,

(6) $\int_{t}^{t+\Delta
t}\sigma(x(s))\eta(s)\,ds\approx\sigma(x(t))\Delta\eta,$

and the midpoint or Stratonovich interpretation,

(7) $\int_{t}^{t+\Delta
t}\sigma(x(s))\eta(s)\,ds\approx\frac{1}{2}\left[\sigma(x(t))+\sigma(x(t+\Delta
t))\right]\Delta\eta.$

The examples to follow will use only the Itô interpretation for simplicity,
although for many physical models the Stratonovich interpretation is more
appropriate [86].

Reflecting the fact that (1) is best understood as a distribution of paths
determined by possible realizations of the random driving term $\eta(t)$, its
solution can be alternatively expressed as a transition density $p(x,t|y,s)$,
such that the probability of finding $x(t)$ in a measurable set $\Omega$
conditioned on $x(s)=y$ with $s\leq t$ is given by

(8) $P(x(t)\in\Omega|x(s)=y)=\int_{\Omega}p(x,t|y,s)\,dx.$

The time evolution of this transition density satisfies the forward Fokker-
Planck equation,

(9) $\frac{\partial
p}{\partial{t}}=-\frac{\partial}{\partial{x}}(fp)+\frac{1}{2}D\frac{\partial^{2}p}{\partial{x}^{2}}$

with initial condition $p(x,s|y,s)=\delta(x-y)$. Conditional expectations can
then be computed against this density using

(10) $\overline{g}(y,t)={\mathbb{E}}[g(x(t))|y,s]=\int g(x)p(x,t|y,s)\,dx,$

or the quantity $\overline{g}(y,t)$ can be evolved directly using the backward
Fokker-Planck equation,

(11)
$\frac{\partial\overline{g}}{\partial{t}}=f(y)\frac{\partial\overline{g}}{\partial{y}}+\frac{1}{2}D\frac{\partial^{2}\overline{g}}{\partial{y}^{2}}$

with initial condition $\overline{g}(y,s)=g(y)$.

An important application of the Fokker-Planck equation is in computing exit
times of state $x$ from a domain $\Omega$. If we define $g$ to be the
indicator function associated with set $\Omega$,
$g(x)={\mathbbm{1}}_{x\in\Omega}$, then $\overline{g}(y,t)$ is the probability
of being in $\Omega$ at time $t$ having started at $x(0)=y$ (we set $s=0$
without loss of generality for time-homogenous processes). If we focus only on
the first exit from $\Omega$, we note that this probability satisfies (11)
with absorbing boundary condition $\overline{g}|_{\partial\Omega}=0$ [39, 76].
Now we note that over a time interval $(t,t+\Delta t)$, the probability of
exit from a trajectory starting at $y$ is given by

(12) $\overline{g}(y,t)-\overline{g}(y,t+\Delta
t)\approx-\frac{\partial\overline{g}}{\partial{t}}\Delta t.$

In other words, $-\partial\overline{g}/\partial t$ is the exit time density.
Under the assumption that all trajectories exit with probability one, the
density can be used to compute the mean first exit time $\tau(y)$ for any
$y\in\Omega$:

(13)
$\tau(y)=\int_{0}^{\infty}t(-\frac{\partial\overline{g}}{\partial{t}})\,dt=\int_{0}^{\infty}\overline{g}(y,t)\,dt,$

where we have integrated by parts and assumed that $t\overline{g}\rightarrow
0$ as $t\rightarrow\infty$. Integrating (11) thus yields a boundary value
problem referred to as Dynkin’s Equation for the mean first exit time [39],

(14)
$f(y)\frac{\partial\tau}{\partial{y}}+\frac{1}{2}D\frac{\partial^{2}\tau}{\partial{y}^{2}}=-1,\qquad\tau\left|\right._{\delta\Omega}=0.$

###  2.2 Large deviation theory, action and optimal paths

Section 2.1 poses two broad alternatives for computing expectations of
functions of a diffusion process. The first is to solve the Fokker-Planck
equation (9) using, for example, eigenfunction expansions (with an inner
product defined with respect to a weighting that renders the operator self-
adjoint). This approach is viable in one dimension as represented in (9), but
quickly becomes numerically intractable if the state is high- or infinite-
dimensional. This is partly due to the fact that, in the limit of small noise
where $D\rightarrow 0$, the differential operator in (9) is singular. This
feature implies the existence of boundary layers, rendering the computation
numerically stiff. Analytically, the existence of a small parameter suggests
the use of asymptotic methods such as singular perturbation or WKBJ theory to
approximate the expectation [12, 63].

The second broad approach is to simulate the random walk in (1) multiple times
and to compute the expectations empirically using, for example, a Monte Carlo
method. The difficulty in using this approach for rare events in high-
dimensional systems stems primarily from the fact that the desired expectation
is dominated by a small set of events of interest. This problem is exacerbated
in the small-noise limit due to an exponential drop-off in likelihood as the
random walker moves away from this set. In this scenario of high dimension and
small noise, it is thus critically important to bias the Monte Carlo
simulations such that a substantial number of simulated trajectories coincide
with those that sample from the set of interest with high likelihood. This
search for a good bias function is made all the more important by the inherent
inefficiency of the Monte Carlo method, which produces estimators with
relative error that scales as $1/\sqrt{N}$ where $N$ is the sample size of
simulations. The importance of this issue has led to the development of many
techniques appropriate for rare-event sampling in a variety of physical
contexts [47, 88, 90, 1].

Regardless of whether one uses asymptotic methods to approximate an
expectation or uses importance sampling to concentrate Monte Carlo simulations
on states of interest, both approaches lead to the theory of large deviations
described in detail in the text by Freidlin and Wentzell [36]. We will now
present a description of this theory, focusing on providing a practical
introduction to exit problems on finite and infinite time horizons.

We return to (1),

(15) $\dot{x}=f(x)+\eta(t),$

where the noise $\eta(t)$ is as described in (4). Each noise realization
$\eta(t)$ produces an associated trajectory $x(t)$. Therefore, the probability
of generating a given path $x(t)$ is proportional to that of generating the
corresponding noise realization $\eta(t)$ [31]. The formal probability density
for this Gaussian noise process satisfies

(16)
$P[\eta(t)]\propto\exp{\left(-\frac{1}{2D}\int_{t_{i}}^{t_{f}}\eta^{2}\,dt\right)}.$

Substitution of (15) into (16) gives [31, 24]

(17)
$P[x(t)]\propto\exp{\left(-\frac{1}{2D}\int_{t_{i}}^{t_{f}}(\dot{x}-f(x))^{2}\,dt\right)}.$

Since we are concerned with maximizing this probability density over a set of
eligible paths satisfying, for example, boundary conditions $x(t_{i})=x_{i}$
and $x(t_{f})=x_{f}$, we express (17) using the conventions of variational
calculus, where

(18)
$P[x(t)]\propto\exp{\left(-\frac{1}{2D}\int[\dot{x}-f(x)]^{2}\,dt\right)}\propto\exp{\left(-\mathcal{S}[x(t)]/D\right)},$

with action ${\mathcal{S}}$ and Langrangian ${\mathcal{L}}$ [35, 24, 36, 62]
defined by

(19)
$\mathcal{S}[x(t)]=\int\limits_{t_{i}}^{t_{f}}\mathcal{L}(\dot{x},x;t)\,dt,\qquad\mathcal{L}(\dot{x},x;t)=\frac{1}{2}[\dot{x}-f(x)]^{2}.$

Next, we calculate the probability of a transition between two states. In the
limit as $D\rightarrow 0$, this probability will be increasingly dominated by
the path that minimizes ${\mathcal{S}}[x(t)]$ so that $P[x(t)]$ is a maximum.
This optimal fluctuational path $x_{\rm opt}$ is found as a solution to the
variational problem $\delta\mathcal{S}[x]=0$, where

$\displaystyle\delta\mathcal{S}[x]$ $\displaystyle=\mathcal{S}[x+\delta
x]-\mathcal{S}[x]=\int\limits_{t_{i}}^{t_{f}}\mathcal{L}(\dot{x}+\delta\dot{x},x+\delta
x;t)\,dt-\int\limits_{t_{i}}^{t_{f}}\mathcal{L}(\dot{x},x;t)\,dt$ (20)
$\displaystyle=\int\limits_{t_{i}}^{t_{f}}\left[\frac{\partial\mathcal{L}}{\partial
x}-\frac{d}{dt}\left(\frac{\partial\mathcal{L}}{\partial\dot{x}}\right)\right]\delta
x\,dt=0.$

Since the variation $\delta x$ is arbitrary one is left with the Euler-
Lagrange equation

(21) $\frac{\partial\mathcal{L}}{\partial
x}-\frac{d}{dt}(\frac{\partial\mathcal{L}}{\partial\dot{x}})=f(x)f^{\prime}(x)-\ddot{x}=0,$

which is then solved for the optimal path $x_{\rm opt}$. It is worth noting
that the optimal path has been found by transforming the original stochastic
problem into a new deterministic system described by the Euler-Lagrange
equation (21). The dimensions of the Euler-Lagrange equation are twice the
dimensions of the original system.

Following the standard progression of classical mechanics, it is sometimes of
use to move from the Lagrangian ${\mathcal{L}}$ expressed in terms of the
state $x$ and associated velocity $\dot{x}$ to the Hamiltonian ${\mathcal{H}}$
formulated as a function of $x$ and the conjugate momentum $\lambda$ defined
as

(22) $\lambda=\frac{\partial\mathcal{L}}{\partial\dot{x}}=\dot{x}-f(x).$

Note that we are denoting the conjugate momentum variable as $\lambda$ rather
than the usual $p$ to avoid conflict with our use of $p$ for probability
density and to emphasize the connection between conjugate momentum and the
costate (i.e., Lagrange multiplier) from the optimal control formulation
below. One way to move from the Lagrangian ${\mathcal{L}}$ to the Hamiltonian
${\mathcal{H}}$ is through the Legendre transformation

(23) $\mathcal{L}=\lambda\dot{x}-\mathcal{H}(x,\lambda),$

where $\dot{x}$ is found as a function of $x$ and $\lambda$ from
$\lambda=\frac{\partial\mathcal{L}}{\partial\dot{x}}$ using the inverse
function theorem. Starting with (1), the Legendre transformation yields

(24)
$\dot{x}=f(x)+\lambda,\quad\mbox{and}\quad{\mathcal{H}}=\frac{\lambda^{2}}{2}+\lambda
f(x),$

illustrating the connection between the classical momentum $\lambda$ and the
optimal fluctuation $\eta$.

One can then write the action $\mathcal{S}$ as

(25)
$\mathcal{S}[x(t)]=\int\limits_{t_{i}}^{t_{f}}\mathcal{L}(\dot{x},x;t)\,dt=\int\limits_{t_{i}}^{t_{f}}\left(\lambda\dot{x}-\mathcal{H}(x,\lambda)\right)\,dt.$

The optimal path is found as a solution to the variational problem
$\delta\mathcal{S}[x,\lambda]=0$, where variations in both $x$ and $\lambda$
must be considered. This leads to

$\displaystyle\delta\mathcal{S}[x,\lambda]=\mathcal{S}[x+\delta
x,\lambda+\delta\lambda]-\mathcal{S}[x,\lambda]$
$\displaystyle=\int\limits_{t_{i}}^{t_{f}}\left[(\lambda+\delta\lambda)(\dot{x}+\delta\dot{x})-\mathcal{H}(x+\delta
x,\lambda+\delta\lambda)\right]\,dt-\int\limits_{t_{i}}^{t_{f}}\left[\lambda\dot{x}-\mathcal{H}(x,\lambda)\right]\,dt$
(26)
$\displaystyle=\int\limits_{t_{i}}^{t_{f}}\left[\left(\dot{x}-\frac{\partial\mathcal{H}}{\partial\lambda}\right)\delta\lambda-\left(\dot{\lambda}+\frac{\partial\mathcal{H}}{\partial
x}\right)\delta x\right]\,dt=0.$

As with the Lagrangian formulation, both variations are arbitrary. Therefore
one is left with Hamilton’s equations

(27)
$\dot{x}=\frac{\partial\mathcal{H}}{\partial\lambda},\qquad\mbox{and}\qquad\dot{\lambda}=-\frac{\partial\mathcal{H}}{\partial
x}.$

To avoid potential inverse function theorem issues, one could find the
Hamiltonian from the outset. Thus,

(28)
$P[\eta(t)]\propto\exp\left(-\frac{1}{2D}\int_{t_{i}}^{t_{f}}\eta^{2}\,dt\right)=\exp(-\mathcal{R}_{\eta}/D)\quad\mbox{and}\quad\mathcal{R}_{\eta}[\eta(t)]=\frac{1}{2}\int_{t_{i}}^{t_{f}}\eta^{2}\,dt,$

where we regard $R_{\eta}$ as an objective function. To find the optimal
trajectory, one must minimize the functional

(29)
$\mathcal{R}[x,\eta,\lambda]=\mathcal{R}_{\eta}[\eta(t)]+\int_{t_{i}}^{t_{f}}\lambda[\dot{x}-f(x)-\eta(t)]\,dt,$

where $\lambda$ is the Lagrange multiplier enforcing the dynamic constraint
(1). Therefore

(30) $P[x(t)]\propto\exp(-R/D),$

where the value $R$ is obtained from

(31) $R={\rm{min}}\mathcal{R}[x,\eta,\lambda].$

We extremize $\mathcal{R}[x,\eta,\lambda]$ with respect to $\eta$, $\lambda$,
and $x$ as follows:

(32)
$\frac{\delta\mathcal{R}}{\delta\eta}=\mathcal{R}[x,\eta+\delta\eta,\lambda]-\mathcal{R}[x,\eta,\lambda]=\int_{t_{i}}^{t_{f}}(\eta-\lambda)\delta\eta\,dt=0,$

so that

(33) $\lambda=\eta;$

(34)
$\frac{\delta\mathcal{R}}{\delta\lambda}=\mathcal{R}[x,\eta,\lambda+\delta\lambda]-\mathcal{R}[x,\eta,\lambda]=\int_{t_{i}}^{t_{f}}(\dot{x}-f(x)-\eta(t))\delta\lambda\,dt=0,$

so that

(35) $\dot{x}=f(x)+\eta(t);\qquad\mbox{and}$

(36) $\frac{\delta\mathcal{R}}{\delta x}=\mathcal{R}[x+\delta
x,\eta,\lambda]-\mathcal{R}[x,\eta,\lambda]=\int_{t_{i}}^{t_{f}}(-\dot{\lambda}-\lambda\frac{\partial
f}{\partial x})\delta x\,dt=0,$

so that

(37) $\dot{\lambda}=-\lambda\frac{\partial f}{\partial x}.$

Taken together we have

(38) $\dot{x}=f(x)+\lambda,\quad\dot{\lambda}=-\lambda\frac{\partial
f}{\partial
x},\quad\mbox{and}\quad\mathcal{H}(x,\lambda)=\frac{\lambda^{2}}{2}+\lambda
f(x).$

By comparison with the Lagrangian formulation, it is clear that the Lagrange
multiplier plays the role of the conjugate momentum.

As mentioned in Sec. 1.1, in both formulations the optimal path is found by
transforming the original stochastic problem into a new deterministic system
described by a Hamiltonian with $\mathcal{H}(x,\lambda)=0$ when $x_{i}$ is a
metastable state and $t_{i}=-\infty$. The dimensions of the Hamiltonian are
twice the dimensions of the original system due to the conjugate momenta
$\lambda$. The method amounts to finding a zero-energy trajectory of an
effective mechanical system, and at least one of the solutions to the zero-
energy Hamiltonian is in the set of optimal paths.

###  2.3 Escape in finite and infinite time

In Sec. 2.2 we formulated a criterion satisfied by optimal (i.e., most likely)
diffusion paths produced by (1). These paths must also satisfy time
constraints and boundary conditions reflecting, for example, a transition
between quasi-stable states either in finite time or over an arbitrarily long
time. The minimal action and associated paths strongly depend on the time
horizon relevant to the computations. Computing mean first exit times from
$\Omega$ generated by (5) requires minimizing the following action over all
exit times $t_{\mathrm{f}}$:

(39)
${\mathcal{S}}_{t_{\mathrm{f}}}=\frac{1}{2}\int_{t_{\mathrm{i}}}^{t_{\mathrm{f}}}(\dot{x}-f(x))^{\dagger}a^{-1}(x)(\dot{x}-f(x))\,dt,$

where we generalize the presentation in Sec. 2.1 to higher dimension so that
$x\in{\mathbb{R}}^{n}$ and $a(x)=\sigma(x)\sigma^{\dagger}(x)$. Performing
this minimization over paths from points $x_{\mathrm{i}}$ to $x_{\mathrm{f}}$
yields the function

(40)
$Q(x_{\mathrm{i}},x_{\mathrm{f}})=\inf_{t_{\mathrm{f}}}{\mathcal{S}}_{t_{\mathrm{f}}}[x(t)|x(t_{\mathrm{i}})=x_{\mathrm{i}},x(t_{\mathrm{f}})=x_{\mathrm{f}}],$

known as the quasi-potential for its similarity in form to a potential defined
between any two points. Indeed, in the case of overdamped Langevin dynamics
[58, 77, 39, 86] when (5) is a noise-driven gradient flow so that

(41) $f(x)=-\left(\frac{\partial U}{\partial{x}}\right)^{\dagger},$

the quasi-potential is intimately related to the potential $U$, and under
additional assumptions
$Q(x_{\mathrm{i}},x_{\mathrm{f}})=2(U(x_{\mathrm{f}})-U(x_{\mathrm{i}}))$.
Regardless of whether $f$ is a gradient, the quasi-potential determines the
scaling law of exit rates $\tau(x):=\inf\\{t:X_{t}\notin\Omega|X_{i}=x\\}$ as
$D\rightarrow 0$, with

${\mathbb{E}}[\tau(x)]\sim\exp\left(\inf_{y\notin\Omega}Q(x,y)/D\right).$

SDEs with sufficient complexity require a numerical approach to compute the
quasi-potential [13], starting either from the optimal control formulation in
(29) or from the Hamiltonian framework presented in (27). To include the
matrix $a(x)$ in (39), assumed positive-definite to reflect a non-degenerate
noise process in (5), the Hamiltonian is changed slightly to

(42)
${\mathcal{H}}(x,\lambda)=\frac{1}{2}\lambda^{\dagger}a\lambda+\lambda^{\dagger}f.$

Finite-time paths can then be solved using a shooting method to solve the two-
point boundary value problem expressed by

(43) $\displaystyle\dot{x}=f(x)+a\lambda\quad\mbox{and}$ (44)
$\displaystyle\dot{\lambda}=-\left(\frac{\partial
f}{\partial{x}}\right)^{\dagger}\lambda-\frac{1}{2}\left[\frac{\partial}{\partial{x}}(\lambda^{\dagger}a\lambda)\right]^{\dagger},$

where we have used the convention of representing vectors as columns and
derivatives with respect to those vectors as rows. Given a typical exit
problem with fixed initial condition $x_{\mathrm{i}}$ and exit criterion
forcing $x(t_{\mathrm{f}})$ to lie on the boundary $\partial\Omega$ of set
$\Omega$, for instance, a straightforward shooting method would start with an
initial guess for $\lambda$ (e.g., $\lambda(t)\equiv 0$), integrate (43),
adjust $\lambda(t_{\mathrm{f}})$ and integrate (44) backwards, then iterate
this process until the terminal condition is satisfied [16, 28].

Different alternatives exist for finding optimal paths in the case of gradient
systems, including the nudged elastic band method [52] and the string method
[27]. A powerful method for more general dynamical systems is the minimum
action method (MAM), which can be adapted via a rescaling in time to
accommodate infinite-time optimal paths, e.g., where either $t_{i}=-\infty$ in
(39) or (40) or the minimization over $t_{f}$ in those expressions yields
$t_{f}=\infty$. It sets up the minimization as the gradient descent of the
action ${\mathcal{S}}_{t_{\mathrm{f}}}$. This can be accomplished in a number
of ways, including by differentiating out the dependence of Hamilton’s
equations on the quasi-momentum $\lambda$ and advancing a PDE whose steady
state is the optimal path sought. The adaptation of this method to infinite-
time problems, referred to as the geometric minimum action method (GMAM) [49],
reparameterizes the optimal path to depend on arclength rather than time,
removing the obvious difficulty of enforcing boundary conditions at infinity.
The iterative action minimization method (IAMM) is a similar method that
employs a direct, fully explicit iterative scheme [61]. Basic details of both
the GMAM and the IAMM can be found in Appendix B.

####  2.3.1 Escape from potential well

As an example of an infinite-time calculation, consider a particle in a
potential well under the influence of external Gaussian noise. The noise
causes the particle to escape from the well in which it resides, an event that
grows increasingly rare as the strength of the noise decreases. To determine
the mean time to escape (MTE) in this limit, we consider again the Langevin
equation

(45) $\dot{x}=-\frac{\partial U}{\partial x}+\eta(t),$

where $x$ is the particle’s position, $U$ determines the potential well, and
$\eta$ is additive Gaussian noise as described by (4). The associated Fokker-
Planck equation is

(46) $\frac{\partial p(x,t)}{\partial t}=\frac{\partial}{\partial
x}\left[\frac{\partial U(x)}{\partial
x}p(x,t)\right]+\frac{1}{2}D\frac{\partial^{2}p(x,t)}{\partial x^{2}},$

where $p$ is the probability density. The first term on the right-hand side is
the drift or transport term, and the second term on the right-hand side is the
diffusion or fluctuation term. It is possible to rewrite (46) as

(47) $\frac{\partial p}{\partial t}=-\frac{\partial J}{\partial x},\qquad
J=-\frac{\partial U}{\partial x}p-\frac{1}{2}D\frac{\partial p}{\partial x}.$

Therefore, $\frac{\partial p}{\partial t}+\frac{\partial J}{\partial x}=0$ is
a continuity equation for the probability density $p$, and $J$ is interpreted
as a probability current. For a stationary process, $J=\,$constant.

Figure 3 shows a schematic of a potential well $U(x)$. We assume the well
height is much larger than the noise intensity so that $\Delta U/D\gg 1$. We
will determine the escape rate/time for particles sitting in a deep well near
$x=x_{\rm min}$.

Figure 3: Schematic of a potential well $U(x)$. The local minimum of the well
is located at $x=x_{\rm min}$, while the top of the barrier is located at
$x=x_{\rm max}$.

The probability current given by (47) can be rewritten as

(48) $J=-\frac{1}{2}D\exp{(-2U(x)/D)}\frac{\partial}{\partial
x}\left[\exp{(2U(x)/D)}p(x,t)\right].$

Integration from $x_{\rm min}$ to $A$ gives, under the assumption of quasi-
stationarity (i.e., $J$ is constant),

(49) $J=\frac{D\exp{(2U(x_{\rm min})/D)}p(x_{\rm min},t)}{2\int\limits_{x_{\rm
min}}^{A}\exp{(2U(x)/D)}\,dx},$

where $x=A$ is the location of an arbitrary point far to the right of the
location of the top of the well barrier (see Fig. 3). For a high barrier the
quasi-stationary distribution of the Fokker-Planck equation satisfies

(50) $0\approx\frac{\partial}{\partial x}\left[\frac{\partial U(x)}{\partial
x}p(x,t)\right]+\frac{1}{2}D\frac{\partial^{2}p(x,t)}{\partial x^{2}},$

where the time dependence of $p(x,t)$ arises from a slow leakage of
probability from the well. The quasi-stationary distribution function near
$x_{\rm min}$ is then

(51) $p(x,t)=p(x_{\rm min},t)\exp{(-2[U(x)-U(x_{\rm min})]/D)}.$

The probability $P(t)$ to find the particle in the well is

(52) $P=\int\limits_{x_{1}}^{x_{2}}p(x,t)\,dx=p(x_{\rm min},t)\exp{(2U(x_{\rm
min})/D)}\int\limits_{x_{1}}^{x_{2}}\exp{(-2U(x)/D)}\,dx,$

where $x=x_{1}$ and $x=x_{2}$ are the locations of arbitrary points in the
potential well found respectively to the left and the right of $x=x_{\rm min}$
(see Fig. 3). Since probability $P$ times the escape rate $r$ is the
probability current $J$, the inverse escape rate (escape time) is

(53)
$\tau=\frac{1}{r}=\frac{P}{J}=\frac{2}{D}\int\limits_{x_{1}}^{x_{2}}\exp{(-2U(x)/D)}\,dx\int\limits_{x_{\rm
min}}^{A}\exp{(2U(x)/D)}\,dx.$

The first integral is largest near $x_{\rm min}$. We expand

(54) $U(x)\approx U(x_{\rm min})+\frac{1}{2}U^{\prime\prime}(x_{\rm
min})(x-x_{\rm min})^{2}.$

The second integral is largest near $x_{\rm max}$. Again, we expand

(55) $U(x)\approx U(x_{\rm max})-\frac{1}{2}|U^{\prime\prime}(x_{\rm
max})|(x-x_{\rm max})^{2}.$

By extending the integral boundaries to $\pm\infty$ in both directions, one
obtains the expression for the MTE as

(56) $\tau=\frac{2\pi}{\sqrt{U^{\prime\prime}(x_{\rm
min})|U^{\prime\prime}(x_{\rm max})|}}\exp{(2[U(x_{\rm max})-U(x_{\rm
min})]/D)}=K\exp{(2\Delta U/D)}.$

It is important to note in (56) that the escape time increases exponentially
with increasing barrier height and decreasing noise intensity, and the
prefactor $K$ depends on the curvatures at $x_{\rm min}$ and $x_{\rm max}$.

A specific example is shown in Fig. 4(a)-(b). Figure 4(a) shows the potential
well given by $U(x)=-x^{3}+\frac{3}{4}x$ with associated barrier height
$\Delta U=\frac{1}{2}$. We numerically integrate (45) using the Euler-Maruyama
method (Appendix A), and the time needed for a particle to escape from the
well is determined. By performing the computation for 10,000 particles, the
mean time to escape is computed for a variety of noise intensities. By taking
the logarithm of (56), one has $\ln{\tau}=\ln{K}+\frac{2\Delta U}{D}$. By
plotting $\ln{\tau}$ versus $2/D$, the analytical prediction is that the data
should lie along a line with slope $m=\Delta U$. Figure 4(b) shows the
numerical results for the potential shown in Fig. 4(a) along with a line of
best fit through the data. One can see excellent agreement between the
numerical and analytical results.

Figure 4: (a) Graph of the potential well given by $U(x)=-x^{3}+\frac{3}{4}x$
with associated well height $\Delta U=\frac{1}{2}$. (b) Logarithm of the
numerically computed mean time to extinction (black circles) as a function of
$2/D$ with a line of best fit (red line) passing through the data. The
equation of the best fit line is $y=0.5x+0.59$, whose slope agrees perfectly
with the analytical prediction of $\Delta U=\frac{1}{2}$.

Returning to the SDE description, the optimal fluctuation equations given by
(27) are easily shown in the gradient case under consideration to give

(57) $\ddot{x}=\frac{\partial
U}{\partial{x}}\frac{\partial^{2}U}{\partial{x}^{2}},$

which integrates once to give

(58) $\dot{x}^{2}=\left(\frac{\partial U}{\partial{x}}\right)^{2}+c.$

In the case of infinite-time exit, the constant $c=0$, leading to an
asymptotic dynamic behavior where the conjugate momentum $\lambda$ is either
zero or acts against the gradient, i.e.,

(59) $\dot{x}=-\frac{\partial U}{\partial{x}}+\lambda=\mp\frac{\partial
U}{\partial{x}},$

depending on whether or not the uncontrolled dynamics (i.e., the deterministic
dynamics obtained by setting $\eta\equiv 0$ in (45)) leads $x$ in the
direction of an exit. If $x(t_{\mathrm{i}})<-1/2$, the uncontrolled dynamics
takes infinite time to descend into the stable fixed point, after which the
controlled dynamics follows the trajectory

(60) $x(t)=\frac{1}{2}\tanh(\frac{3}{2}t),$

taking infinitely long both to rise out of the stable fixed point and to
approach the saddle, thereby effecting an exit. The minimum action for any
initial condition to the left of the stable fixed point is
${\mathcal{S}}=2\Delta U=1$.

To see how this changes when a finite time horizon is imposed, note that in
this case (57) no longer implies (59).

Figure 5: Optimal paths computed from the stable fixed point to the saddle for
final times $t_{\mathrm{f}}=1/2$ (blue), $t_{\mathrm{f}}=1$ (red),
$t_{\mathrm{f}}=2$ (orange), $t_{\mathrm{f}}=4$ (purple), $t_{\mathrm{f}}=8$
(green) and $t_{\mathrm{f}}=\infty$ (black). (a) Top left: Optimal paths
rescaled to unit exit time. (b) Top right: Optimal paths with time shifted so
that passages through zero coincide. (c) Bottom left: Phase portrait
($\dot{x}$ vs. $x$). (d) Bottom right: Minimum action ${\mathcal{S}}$ vs.
$t_{\mathrm{f}}$.

Figures 5(a)-(d) show the optimal paths and their phase portraits associated
with $t_{\mathrm{f}}\in\\{1/2,1,2,4,8\\}$ and compares them to the infinite-
time optimal path. Whereas for short time $t_{\mathrm{f}}=1/2$ the optimal
path is near a straight line path (more generally, a geodesic with respect to
the metric imposed by $a(x)$), the optimal path quickly approaches the
infinite-time path as $t_{\mathrm{f}}$ increases. Once these paths are
indistinguishable, i.e., by $t_{\mathrm{f}}=4$, increases in $t_{\mathrm{f}}$
simply add to the time spent vanishingly close to either fixed point. By this
value of $t_{\mathrm{f}}$, the action has also converged to its quasi-
potential value of ${\mathcal{S}}=Q(-1/2,0)=1$. The phase portrait depicted in
Fig. 5(c) demonstrates again that for short times $t_{\mathrm{f}}$ the optimal
path tends to a constant velocity that generates the straight-line path
depicted in Fig. 5(b); however, for large times the path through phase space
collapses onto the anti-gradient flow described above.

Finally, we note that the reference to “infinite-time optimal paths” implies
passage from, to, or through a hyperbolic fixed point. The optimization
suggested by (39) performed on paths connecting two points that are not
hyperbolic fixed points of the dynamical system and that do not pass through
one, e.g., $x(t_{\mathrm{i}})=-1/4$ and $x(t_{\mathrm{f}})=1/4$, yields an
action determined by the quasi-potential $Q=2\Delta U=11/16$ and an optimal
path time of

(61) $\Delta t=\frac{4}{3}\tanh^{-1}\frac{1}{2}=0.73.$

###  2.4 Master equation formalism for internal noise

Just as we considered the mean exit time due to external noise in Sec. (2.1),
we would like to determine the mean time to extinction (MTE) and the optimal
path for internal noise. Assuming again that we can model stochastic
population dynamics using a Markov process, we describe the evolution of the
probability density using the master equation

(62)
$\frac{\partial\rho(X,t)}{\partial{t}}=\sum\limits_{r}\left[W_{r}(X-r)\rho(X-r,t)-W_{r}(X)\rho(X,t)\right],$

where $\rho(X,t)$ is the probability of finding $X$ individuals at time $t$,
and $W_{r}(X)$ is the transition rate from $X$ to $X+r$, where $r$ is a
positive or negative integer increment. The first term of (62) is the gain to
state $X$ from state $X-r$. The second term of (62) is the loss from state $X$
to other states.

The master equation is a large, or even infinite, set of differential
equations. A diffusion approximation leading to a Fokker-Planck equation as in
(9) is often used but is very (exponentially) inaccurate [42, 21, 54, 5].
Assuming that the typical size of the population $X\sim K$ as determined by
the maximum of the quasi-stationary probability distribution, we can use a
WKBJ approximation of the master equation [60, 38, 25, 28, 54]. If the mean-
field equation associated with the original stochastic problem is
$n$-dimensional, then the WKBJ method transforms the problem into one of
classical mechanics in $2n$-dimensional phase space where the doubling of
dimension is due to the conjugate momenta. It is important to note that these
deterministic dynamics are not the system size $K\to\infty$ mean-field
dynamics of the original stochastic model, but rather provide auxiliary
dynamics that describe processes not captured by the mean-field dynamics.

For large populations, the time to extinction is long. The mean time to
extinction is determined by the probability flux (extinction rate) into the
extinct state, and is determined by the tail of the quasi-stationary
probability distribution where $\frac{\partial\rho}{\partial t}\approx 0$.
Therefore

(63) $0=\sum\limits_{r}\left[W_{r}(X-r)\rho(X-r)-W_{r}(X)\rho(X)\right].$

Let $X$ be scaled by $K$, the typical population size in the metastable state.
Using $x=X/K$, the transition rate $W_{r}(X)=W_{r}(Kx)$ can be represented as
the following expansion in $K$,

(64) $W_{r}(Kx)=Kw_{r}(x)+u_{r}(x)+\mathcal{O}(1/K),$

where $x$ and the scaled transition rates $w_{r}$ and $u_{r}$ are
$\mathcal{O}(1)$. Additionally we can write,

(65) $\rho(X)\equiv\rho(Kx)=\pi(x).$

Then the scaled master equation is

(66)
$\sum\limits_{r}\left[w_{r}(x-\frac{r}{K})\pi(x-\frac{r}{K})-w_{r}(x)\pi(x)\right]=0.$

For $K\gg 1$ we approximate the scaled master equation using the WKB
approximation. To apply the WKB approximation, we assume that

(67) $\pi(x)=A(x)\exp{[(-K\mathcal{S}(x))(1+\mathcal{O}(1/K))]}.$

We substitute the WKB ansatz along with
$\mathcal{S}(x-\frac{r}{K})=\mathcal{S}(x)-\frac{r}{K}\mathcal{S}^{\prime}(x)+\frac{r^{2}}{2K^{2}}\mathcal{S}^{\prime\prime}(x)+\mathcal{O}(1/K^{3})$
and $A(x-\frac{r}{K})=A(x)-\frac{r}{K}A^{\prime}(x)+\mathcal{O}(1/K^{2})$ into
the master equation, and at leading order one obtains a Hamilton-Jacobi
equation with Hamiltonian

(68)
$\mathcal{H}(x,\lambda)=\sum\limits_{r}w_{r}(x)[\exp{(r\lambda)-1}]=0,\qquad\lambda=\frac{\partial\mathcal{S}}{\partial
x}.$

Hamilton’s equations are

(69a) $\displaystyle\dot{x}=$
$\displaystyle\frac{\partial\mathcal{H}}{\partial\lambda}=\sum\limits_{r}rw_{r}(x)\exp{(r\lambda)},$
(69b) $\displaystyle\dot{\lambda}=$
$\displaystyle-\frac{\partial\mathcal{H}}{\partial
x}=-\sum\limits_{r}[\exp{(r\lambda)}-1]\frac{\partial w_{r}(x)}{\partial x}.$

In summary, the WKB solution amounts to finding a zero-energy trajectory of an
effective mechanical system, and at least one of the solutions to the zero-
energy Hamiltonian is the optimal path. The $x$ dynamics along the $\lambda=0$
deterministic line are described by

(70)
$\dot{x}=\frac{\partial\mathcal{H}(x,\lambda)}{\partial\lambda}\rvert_{\lambda=0}$

which is the rescaled mean-field rate equation associated with the original
deterministic problem. For a single step process, this simplifies to
$\dot{x}=w_{1}(x)-w_{-1}(x)$.

In Sec. 3.3, we will discuss the application of this master equation formalism
to two examples involving extinction in population biology. For the simple
models discussed in Sec. 3.3, the deterministic steady states are nodes. It is
easy to show that the WKB method transforms these steady state nodes in the
original 1D deterministic system into steady state saddle points in the 2D set
of Hamilton’s equations. This allows for noise-induced escape from a
metastable state and provides a path to extinction that did not exist in the
original deterministic model.

###  2.5 Mean Time to Extinction

In a single step process, such as the two examples considered in Sec. 3.3, the
optimal path $\lambda_{opt}(x)$ will always have the general form

(71) $\lambda_{opt}(x)=-\ln{\left(w_{+1}(x)/w_{-1}(x)\right)}.$

Using the definition of the conjugate momentum $\lambda=d{\cal S}/dx$, the
action ${\cal S}_{opt}$ along the optimal path $\lambda_{opt}(x)$ is given by

(72) ${\cal S}_{opt}=\int^{x}\lambda_{opt}(x)dx,$

where the limits of integration are determined by the appropriate equilibrium
points of the specific problem. Therefore, the mean time to extinction (MTE)
can be approximated by

(73) $\tau=B\exp({K{\cal S}_{opt}}),$

where $B$ is a prefactor that depends on the system parameters and on the
population size. An accurate approximation of the MTE depends on obtaining
$B$. The specific form of the prefactor differs for the two types of
extinction problems considered in Sec. 3.3, and are provided in that section.

###  2.6 Sampling for exits

In the simple gradient flow discussed in Sec. 2.3.1 it is possible to obtain
or approximate not only the exponential scaling law in $D$ of the exit
probability over finite times, but also the normalization constant that
provides the pre-exponential factor $K$ in (56). In flows that are more
complex, including non-potential flows that do not satisfy detailed balance
[86] and those of higher (or infinite) dimension, this is no longer possible
analytically. Nor is it practical computationally, due to the formation of
boundary layers with the singular perturbation introduced as $D\rightarrow 0$
in Dynkin’s equation (14) or the time-dependent Fokker-Planck equation (9).

An alternate method to obtain the normalization constant uses sampling through
Monte Carlo simulations. This too is subject to the “curse of dimensionality”
in cases of high dimension such as discretizations of stochastic PDEs;
however, variance reduction techniques can be used to improve the numerical
efficiency of obtaining an estimate. We will focus here on importance
sampling, a variance reduction method that uses information from the
controlled dynamics of the deterministic system to inform a bias applied to
the distribution from which the Monte Carlo simulations are drawn.

Suppose we wish to calculate the probability that, after time
$t_{\mathrm{f}}$, a stochastic process $x(t)$ evolving according to (5)
terminates outside of set $\Omega$, i.e., with
$x(t_{\mathrm{f}})\notin\Omega$, where we assume for concreteness that the
deterministic evolution $\dot{x}=f(x,t)$ yields $x(t)\in\Omega$ for
$x(t_{\mathrm{i}})=x_{\mathrm{i}}$ and $t\leq t_{\mathrm{f}}$. In other words,
we wish to calculate the probability that the stochastic forcing term
$\sigma(x,t)\eta$ drives the state some distance away from its deterministic
trajectory. The standard Monte Carlo estimator for this probability $P$ is
given by

(74)
$\hat{P}_{\mathrm{MC}}=\frac{1}{N}\sum_{k=1}^{N}I(X_{t_{\mathrm{f}}}^{(k)})$

where $\\{X_{t_{\mathrm{f}}}^{(k)},k=1,\dots,N\\}$ are $N$ independent
simulations of (1) and $I(x)$ is an indicator function evaluating to 0 if
$x\in\Omega$ and 1 otherwise. Thus, $I(X_{t_{\mathrm{f}}})$ forms a Bernoulli
random variable, implying that the estimator $\hat{P}_{\mathrm{MC}}$ has
variance given by

(75) $\mathbb{V}\left[\hat{P}_{\mathrm{MC}}\right]=\dfrac{P(1-P)}{N},$

with a relative error (coefficient of variation) of,

(76)
$\mathbb{C}_{\mathrm{var}}\left[\hat{P}_{\mathrm{MC}}\right]=\frac{\sqrt{\mathbb{V}\left[\hat{P}_{\mathrm{MC}}\right]}}{\mathbb{E}\left[\hat{P}_{\mathrm{MC}}\right]}=\dfrac{\sqrt{1-P}}{\sqrt{NP}}.$

Ensuring that (74) is an accurate estimator for the true probability requires
that $\mathbb{C}_{\mathrm{var}}\left[\hat{P}_{\mathrm{MC}}\right]\ll 1$.
However, if the set $\Omega$ and noise strength $D$ are such that $P\ll 1$
(i.e., $X_{t_{\mathrm{f}}}\notin\Omega$ is a rare event), this requirement is
approximately expressed by $N\gg 1/P$, which implies that a very large (and
often unattainable) number of Monte Carlo runs is required to obtain an
accurate estimate.

The idea behind importance sampling for diffusions is to replace the original
dynamics given by (5) with

(77) $\dot{x}=f(x,t)+\sigma(x,t)(b+\eta),$

so that the noise term $\eta$ is biased in such a way that the events of
interest, $x({t_{\mathrm{f}}})\notin\Omega$, are much more likely. The
particular form of bias suggested by (77) is mean-biasing; other forms of
biasing the noise are possible but less common. Mean biasing effectively
replaces increments $\Delta\eta$ drawn from the original distribution (e.g.,
Gaussian) $p(\Delta\eta)$ by shifted increments

(78) $\Delta\tilde{\eta}=\Delta\eta+b\Delta t,$

which is equivalent to drawing $\Delta\tilde{\eta}$ from biased distribution
$\tilde{p}(\Delta\tilde{\eta})=p(\Delta\tilde{\eta}-b\Delta t)$. Despite the
biased dynamics described by (77), an unbiased estimator can be recovered from
Monte Carlo simulations of this equation through an appropriate weighting
referred to as the likelihood ratio

(79)
$l(\Delta\tilde{\eta})=\frac{{\tilde{p}}(\Delta\tilde{\eta})}{p(\Delta\tilde{\eta})}$

with the likelihood ratio for the $k^{\mathrm{th}}$ simulated path given by
the product of these increments. The unbiased estimator is then

(80)
$\hat{P}_{\mathrm{IS}}=\frac{1}{N}\sum_{k=1}^{N}I({\tilde{X}}_{t_{\mathrm{f}}}^{(k)})l^{(k)},\quad\mathrm{with}\quad
l^{k}=\prod_{j}l(\Delta\tilde{\eta}^{(k)}_{j})$

and with variance

(81)
$\mathbb{{\tilde{V}}}\left[\hat{P}_{\mathrm{IS}}\right]=\dfrac{\mathbb{{\tilde{E}}}\left[I({\tilde{X}}_{t_{\mathrm{f}}})l^{2}\right]-P^{2}}{N}\approx\dfrac{\mathbb{E}\left[I({\tilde{X}}_{t_{\mathrm{f}}})l\right]}{N},$

where ${\mathbb{\tilde{E}}}$ and ${\mathbb{\tilde{V}}}$ denote expectation and
variance with respect to the biased noise process.

A “good” importance sampling density is one that produces a small variance,
i.e., one that satisfies $\tilde{p}\gg p$ for the events of interest,
$x({t_{\mathrm{f}}})\notin\Omega$. In the present example, the biased density
is determined by the control term in (77). In general, determining the control
is a difficult task and one that depends sensitively on the deterministic
dynamics and noise covariance structure. To motivate a particularly effective
choice of $b$, we appeal again to the optimal control problem described above
arising in the theory of large deviations. In particular, note that (77)
essentially just adds noise to (43) with control $b=\sigma^{\dagger}\lambda$.
Thus, building a control that satisfies (44) provides paths that form the most
likely routes to the rare events of interest in (1).

##  3 Applications

###  3.1 Phase slips in mode-locked lasers

Figure 6:  Linewidth analysis of a mode-locked laser, demonstrating the
dominant role played by timing and phase dynamics in determining linewidth.
Reproduced from Ref. [9].

Mode-locked lasers provide an important example of systems subject to failure
caused by large noise-induced phase excursions. These lasers produce pulses of
light whose electric field envelope’s phase is locked to the phase of the
underlying carrier wave. In frequency space, this corresponds to a highly
regular comb of laser lines with broad bandwidth and narrow linewidth [51].
Important applications of the temporal and spectral properties of mode-locked
lasers depicted in Fig. (7) include laser manipulation of ultrafast chemical
reactions and extremely precise measurement of time intervals [89].

A typical analysis of the effect of noise on laser operation is restricted to
the induced linewidth obtained by linearizing the model equations for the
laser about their desired operational state to produce Ornstein-Uhlenbeck
dynamics; see for example [9],

$\dot{x}=-Ax+\eta,\quad x=(\Delta g,\Delta
w,\Delta\omega,\Delta\tau,\Delta\theta)^{\dagger},$

where $\Delta g$, $\Delta w$, $\Delta\omega$, $\Delta\tau$, and $\Delta\theta$
represent fluctuations in saturated gain, pulse energy, central frequency,
central pulse timing, and phase, respectively. The conclusion suggested by
Fig. 6 is that the timing and phase dynamics provide the dominant mechanisms
by which noise affects linewidth. As explained below, however, this analysis
neglects another important source of uncertainty.

The generation of a frequency comb with narrow individual lines over a broad
spectral region requires precise control over the carrier-envelope phase
offset, depicted as $\Delta\phi$ in Fig. 7. This control is implemented
through interferometric feedback that typically only measures $\Delta\phi$
modulo full rotations of $2\pi$. If a sequence of noise events successfully
pushes $\Delta\phi$ to $\pm\pi$ (assuming without loss of generality that
$\Delta\phi=0$ for the noise-free operational state) then the feedback drives
$\Delta\phi$ to the neighboring equilibrium value of $\pm 2\pi$. This
situation is referred to as a phase slip and is an additional source of
uncertainty in the laser’s output. A noise-driven damped pendulum provides a
good conceptual framework for these two sources of uncertainty, where the
linewidth is provided by the stationary measure obtained in the small-
amplitude limit of a damped linear oscillator, while phase slips correspond to
the occasional full rotations induced in the pendulum by very unlikely
sequences of noise increments.

Figure 7:  Sketch of electric field output from a mode-locked laser, in time
(top) and in frequency space (bottom). A phase slip occurs when $\Delta\phi$
undergoes a noise-driven rotation of $2\pi$. Reproduced from Ref. [51].

To account for these rare events, consider the heuristic mode-locked laser
model with active modulation introduced in Ref. [22]:

(82) $i\frac{\partial
u}{\partial{t}}+\frac{1}{2}\frac{\partial^{2}u}{\partial{\xi}^{2}}+|u|^{2}u=-c_{0}\cos(\omega\xi)u-ic_{1}u+ic_{2}\frac{\partial^{2}u}{\partial{\xi}^{2}}+id_{1}|u|^{2}u-id_{2}|u|^{4}u+i\eta(\xi,t),$

where $u$ is the normalized complex electric field envelope, $\xi$ is the
averaged transverse variable (e.g., time, for temporal pulses as seen on an
oscilloscope), and $t$ is the longitudinal variable (e.g., a continuous
approximation of the cavity round-trip index). Constants $c_{1}$, $c_{2}$,
$d_{1}$, and $d_{2}$ respectively represent coefficients of linear cavity
loss, spectral filtering, nonlinear gain, and gain saturation. Constants
$c_{0}$ and $\omega$ are the amplitude and inverse width of active modulation
intended to stabilize the pulse with its center position at the origin, and
$\eta(\xi,t)$ is space-time white noise, with

(83)
${\mathbb{E}}[\eta^{\dagger}(\xi,t)\eta(\xi^{\prime},t^{\prime})]=D\delta(\xi-\xi^{\prime})\delta(t-t^{\prime}).$

The stabilization term in (82) provides a restoring force on the pulse
parameters that mitigates the impact of disturbances, including noise. This
effect is characterized below through a linearized analysis of the reduced
system of pulse parameters. The stabilization has a secondary effect, however,
of producing additional equilibria that become dynamically accessible upon the
introduction of noise. Transitions to these equilibria represent phase slips
in the context of mode-locked lasers, and their rate of occurrence is an
important quantity. We demonstrate below how this rate is computed
approximately using the quasi-potential, and how the optimal paths provide
insight into how phase slips arise.

When all coefficients are trivial, so that
$c_{0}=c_{1}=c_{2}=d_{1}=d_{2}=D=0$, (82) is the integrable nonlinear
Schrödinger equation, which supports a four-parameter family of soliton
solutions with functional form

(84) $u_{s}(\xi,t)=A(t){\rm sech}\hskip
0.72229pt[A(t)(\xi-\Xi(t))]e^{i\Theta(t)+i\xi\Omega(t)},$

where $\dot{A}=\dot{\Omega}=0$, $\dot{\Xi}=\Omega$, and
$\dot{\Theta}=\frac{1}{2}(A^{2}-\Omega^{2})$. In the perturbative limit where
the coefficients are small, the dynamics of (82) can be approximated well for
short times by a diffusion process in the four parameters induced by the noise
term $\eta$ [50]. Related approximations made outside of the perturbative
limit are not rigorously justifiable yet often agree well with numerics [56].
We employ an approach based on averaging the Lagrangian density

(85) $\quad
L[u_{s},u_{s\xi},u_{st}]:=\Im(u_{st}^{*}u_{s})-\frac{1}{2}|u_{s\xi}|^{2}+\frac{1}{2}|u_{s}|^{4}+c_{0}\cos(\omega\xi)|u_{s}|^{2}.$

of the variational (i.e., non-dissipative) terms in (82) over the ansatz
expressed by (83). The dissipative terms in (82) are included through their
projections against the tangent space of (83) [3, 8], i.e.,

$\displaystyle\nabla_{y}{{\mathcal{L}}}-$
$\displaystyle\frac{\partial}{\partial{t}}\left(\nabla_{\dot{y}}{{\mathcal{L}}}\right)=$
(86) $\displaystyle 2\mathop{\rm Re}\int
i[-c_{1}u_{s}+c_{2}u_{s\xi\xi}+d_{1}|u_{s}|^{2}u_{s}-d_{2}|u_{s}|^{4}u_{s}+\eta(\xi,t)]\nabla_{y}{u_{s}^{*}}\,d\xi$

where ${\mathcal{L}}=\int L\,d\xi$ and $y=(A,\Omega,\Xi,\Theta)^{T}$. Note
that this Lagrangian ${\mathcal{L}}$ and its density $L$ are associated with
the variational model reduction technique used in this application, and are
not related to the Lagrangian associated with the action minimization
described in Sec. 2.2.

Figure 8:  Left axis: Expected phase slip rate according to large deviation
theory. Right axis: Square root of diagonal elements of covariance $\Sigma$,
representing the root-mean-squared linewidth based on a conventional
linearized analysis. Minimum escape probability is achieved at a different
value of $\omega$ than minimum position jitter $\Sigma_{\Xi}^{1/2}$.

Keeping only the first three components of $y$, i.e., $x_{j}=y_{j}$ for
$j=1,\dots,3$, we have

(87) $\dot{x}=f(x)+\sigma(x)\eta(t)$

with drift

(88)
$f(x)=\begin{pmatrix}-2c_{1}x_{1}+(\frac{4}{3}d_{1}-\frac{2}{3}c_{2})x_{1}^{3}-\frac{16}{15}d_{2}x_{1}^{5}-2c_{2}x_{1}x_{2}^{2}\\\
-\frac{4}{3}c_{2}x_{1}^{2}x_{2}-\pi c_{0}\omega^{2}\mathop{\rm
csch}\nolimits(\pi\omega/2x_{1})\sin(\omega x_{3})/2x_{1}^{3}\\\
x_{2}\end{pmatrix}$

and diffusivity

(89) $\sigma(x)=\begin{pmatrix}\sqrt{x_{1}}&0&0\\\ 0&\sqrt{x_{1}/3}&0\\\
-x_{3}/\sqrt{x_{1}}&0&\sqrt{\pi^{2}/12x_{1}^{3}+x_{3}^{2}/x_{1}}\end{pmatrix}.$

The phase dynamics for $y_{4}=\Theta$ are slaved to the other independent
variables, with

$\displaystyle\dot{\Theta}$ $\displaystyle=-\frac{\pi\omega
c_{0}}{A^{3}}\cos(\omega\Xi)\mathop{\rm
csch}\nolimits\left(\frac{\pi\omega}{2A}\right)\left(1+\frac{\pi\omega}{2A^{2}}\coth\left(\frac{\pi\omega}{2A}\right)\right)$
(90)
$\displaystyle-\Xi+\frac{1}{2}(A^{2}-\Omega^{2})+\left(\frac{\sqrt{12+\pi^{2}}}{6\sqrt{A}}\right)\eta_{4}.$

The noise term $\eta(t)\in{\mathbb{R}}^{3}$ in (87) satisfies

(91) ${\mathbb{E}}[\eta(t)\eta(t^{\prime})^{T}]=DI\delta(t-t^{\prime})$

and is obtained from $\eta(\xi,t)$ in (82) using the same averaged Lagrangian
method, and essentially represents that portion of the stochastic driving term
that directly influences the parameters in $x$. Matrix $I$ is the $3\times 3$
identity.

Figure 9: (a) Minimum action connecting stable fixed point at origin ($n=0$)
to saddle ($n=1$), plotted vs. parameter $\omega$ from (82). (b) Laser pulse
(black) and active feedback for $\omega=0$ (blue), $\omega=1.9$ (red), and
$\omega=10$ (green), corresponding to the squares in (a).

In the deterministic ($D=0$) limit, the nontrivial equilibria are
$x^{(n)}=(A_{0},\;0,\;n\pi/\omega)$ for $n\in{\mathbb{Z}}$, with

$A_{0\pm}^{2}=\frac{5}{16d_{2}}\left(2d_{1}-c_{2}\pm\sqrt{(2d_{1}-c_{2})^{2}-\frac{96}{5}c_{1}d_{2}}\right).$

Of these, the equilibria with $A_{0+}$ are stable nodes or spirals for $n$
even and saddles for $n$ odd. The equilibria with $A_{0-}$ are always saddles.
To repeat the analysis from Ref. [9], linearization about $u_{0}$ gives steady
state covariance

(92) $\Sigma=\frac{D}{2\pi}\int
ds\left[(isI-M)^{-1}\sigma(x^{0})\right]\left[(isI-M)^{-1}\sigma(x^{0})\right]^{{\dagger}},$

where

(93) $M=\begin{pmatrix}8c_{1}-\frac{4}{3}(2d_{1}-c_{2})A_{0}^{2}&0&0\\\
0&-\frac{4}{3}c_{2}A_{0}^{2}&-\frac{(-1)^{n}\pi
c_{0}\omega^{3}}{2A_{0}^{3}}\mathop{\rm
csch}\nolimits(\frac{\pi\omega}{2A_{0}})\\\ 0&1&0\end{pmatrix}.$

As discussed above, however, the feedback mechanism requires that we consider
two sources of uncertainty in the laser output, the first associated with
laser linewidth and the second associated with the expected rate of phase
slips. A complete analysis of the phase slip rate requires a consideration of
the phase dynamics expressed by (87) and (90), including contributions from
the dynamics of all other parameters. For simplicity of presentation, we
consider here transitions between the equilibria described above, noting that
the drift term in (90) changes by $2\pi/\omega$ between neighboring stable
equilibria. In appropriate parameter regimes, this is the dominant
contributing factor to phase slips.

As described in Sec. 2.3, the expected time for transition between neighboring
equilibria is $\exp(Q(x^{0},x^{2})/D)$, where $Q$ is the quasi-potential.
Figure 8 plots these distinct types of uncertainty against one of the active
feedback parameters, $\omega$, where the reciprocal of the transition time has
been expressed as a rate. It is clear from the figure that one cannot in
general expect to be able to minimize both linewidth and phase slip rate
simultaneously; practical considerations inform the balance that must be
struck between these sources of uncertainty [66].

Although the action $Q$ provided by the GMAM computation is the quantity that
determines the mean exit rate estimate, it is also instructive to examine the
paths associated with minimum action for various choices of physical
constants.

Figure 10:  Minimum action paths from stable fixed point ($n=0$) to saddle
($n=1$) for (a) $\omega=0$, (b) $\omega=1.9$, and (c) $\omega=10$,
corresponding to squares in Fig. 9(a).

Figure 11: (a) Minimum action path connecting stable fixed point at origin
($n=0$) to saddle ($n=1$), plotted vs. parameter $c_{0}$ from (82). Minimum
action paths for squares on (a) are plotted in (b), (c), and (d).

Figure 9 graphs the minimum action as a function of the inverse width of the
active feedback trap, represented by $\omega$ in (82). For small values of
$\omega$, the trap is so wide as to be essentially constant over the support
of the pulse. This leaves (82) effectively invariant with respect to
translations of the pulse offering no restoring force to noise that moves the
pulse’s position from the origin. The action for a change in position along
this neutrally stable manifold is therefore zero. This is reflected in Fig.
10(a), where the most likely transition path involves simple translation in
$\xi$. Conversely, for large values of $\omega$, as seen in the green curve,
the trap’s wavelength is much shorter than the pulse’s width, so that its
ability to stop the pulse from moving is limited. Neighboring equilibria are
also very close together in this limit. The action therefore approaches zero
and the most likely transition path between equilibria is again a linear
translation, as shown in Fig. 10(c). Between these values of $\omega$,
however, is a range of values where variations of the trap are commensurate in
scale with the support of the soliton, such that the trap offers maximal
resistance to large excursions and generates a large value of action. Figure
10(b) shows that this increased action is associated with an exit path that
exploits the internal dynamics of the SDE much more significantly. This
picture is qualitatively similar to the impact of varying $\omega$ on the
standard deviation, or jitter, in pulse position, but the minimum is achieved
at a different value of $\omega$, as illustrated in Fig. 8.

Figure 11 illustrates a different dependence of action on parameter value
$c_{0}$, the amplitude of the active feedback trap for fixed $\omega=1$. In
this case, the action increases monotonically with $c_{0}$, and the minimum
action paths are seen to become more and more tortuous as the trap becomes
more and more effective at blocking the transition mechanism of simple
translation.

###  3.2 Importance sampling for large soliton walks

In the laser dynamics example in Sec. 3.1, the minimum action calculated for
passage from a stable fixed point to a neighboring saddle is used to
approximate the average rate at which noise is expected to drive the laser
from its desired operating point. The paths associated with these minimum
action computations are informative in their own right since they illustrate
the dynamic paths taken on these most probable excursions. In fact, these
paths provide the information necessary to implement importance sampling in
Monte Carlo simulations designed to compute correct expectations of exit
probabilities in finite time as described in Sec. 2.6. Figure 12 demonstrates
that these paths change significantly as the time horizon for an exit is
varied, from a geodesic path with respect to a norm dictated by the diffusion
tensor to the limiting path corresponding to the quasi-potential calculations
from Sec. 3.1.

Figure 12: Minimum action paths for the SDE given by (87) as the horizon time
$t_{\mathrm{f}}$ is varied.

For exit probabilities on finite times, traditional large deviation theory
provides an asymptotic scaling law,

(94) $\lim_{D\rightarrow 0}D\ln
P\left(x(t_{\mathrm{f}})\notin\Omega\right)=-\inf_{x(t_{\mathrm{f}})\notin\Omega}{\mathcal{S}}_{t_{\mathrm{f}}}[x],$

where the probability $P$ contains an undetermined prefactor. For the reduced
model provided by (87), computing this prefactor is at least possible using
more substantial analysis of the Fokker-Planck equation [63]. Given that the
model reduction from the stochastic PDE (82) is itself imperfect, one would
ideally seek to perform Monte Carlo simulations on the original model to
verify the accuracy of estimates obtained through the reduced model. To
compute exit probabilities in the context of small noise, a naive Monte Carlo
approach is computationally prohibitive, requiring the use of importance-
sampled Monte Carlo simulations based on the paths computed above [67, 68,
83].

Figure 13:  Probability density $p$ associated with the diffusion of position
$\Xi$ obtained from importance-sampled Monte Carlo simulations of the
stochastic PDE (95) (black squares) and from unbiased simulations (brown
circles). The black curve is the Gaussian probability density function
predicted by linearizing the reduced ODE given by (96)-(98).

For example, a simplified version of (82) relevant to optical communications
includes spectral filtering, compensatory gain ($c_{1}<0$), and white noise:

(95) $\displaystyle i\frac{\partial u}{\partial{t}}+$
$\displaystyle\frac{1}{2}\frac{\partial^{2}u}{\partial{\xi}^{2}}+|u|^{2}u=-ic_{1}u+ic_{2}\frac{\partial^{2}u}{\partial{\xi}^{2}}+i\eta(\xi,t).$

Filters have been used as a mechanism for reducing jitter in the position of
solitons by adding a restoring force to the spatial frequency, which is tied
to the soliton’s speed through the dispersion relation [50]. Damping
fluctuations in spatial frequency is therefore an effective way to suppress
position jitter, as illustrated in Fig. 13. The same reduction method used
above produces the SDE given by

(96) $\dot{x}=f(x)+\sigma(x)\eta$

with $x=(A,\,\Omega,\,\Xi)^{T}$, where the drift is

(97)
$f(x)=\begin{pmatrix}-2c_{1}x_{1}-\frac{2}{3}c_{2}x_{1}^{3}-2c_{2}x_{1}x_{2}^{2}\\\
-\frac{4}{3}c_{2}x_{1}^{2}x_{2}\\\ x_{2}\end{pmatrix},$

and diffusivity is

(98) $\sigma(x)=\begin{pmatrix}\sqrt{x_{1}}&0&0\\\ 0&\sqrt{x_{1}/3}&0\\\
-x_{3}/\sqrt{x_{1}}&0&\sqrt{\pi^{2}/12x_{1}^{3}+x_{3}^{2}/x_{1}}\\\
\end{pmatrix},$

and where we have again omitted the dynamics in $\Theta$. The GMAM method
presented in Sec. 2.3 can be applied to compute the optimal path from a pulse
centered at the origin with amplitude, frequency, and position coordinates
$(A,\Omega,\Xi)=(1,0,0)$ to a final position $\Xi(t_{\mathrm{f}})=\Xi_{f}$
with the other two final coordinates left unspecified. By picking different
values of $\Xi_{f}$, one can reconstruct a histogram of probabilities from
which the entire probability density function (pdf) for the final position
$\Xi$ can be computed well down into the low-probability tails. Figure 13
illustrates such a reconstructed pdf, where 50,000 importance-sampled Monte
Carlo simulations of the stochastic PDE (95) have been used to compute pdf
tail probabilities on the order of $10^{-12}$ and below. Also shown is the
Gaussian pdf predicted by a linearization of the ODE reduction given by
(96)-(98). This curve and simulations of the nonlinear ODEs (not shown)
demonstrate that the low-order reduction fails to accurately compute
probabilities in the tails of the pdf, but Fig. 13 clearly shows that it is
effective at generating appropriately biased simulations of the original PDE
model.

###  3.3 Extinction in biological populations

We now consider internal noise and how it may induce the rare event of
extinction of a species or of an infectious disease. Extinction can be good
(disease) or bad (species); either way, it is important to understand why and
how often it can be expected to happen for different model parameters. We
consider two examples: (1) the Susceptible-Infectious-Susceptible (SIS)
epidemic model, and (2) the Allee effect model.

$x_{0}$Extinctstate$x_{1}$Endemicstate(a)

$x_{0}$Extinctstate$x_{2}$Carryingcapacity$x_{1}$Alleethreshold(b)

Figure 14: (a) Topology of the deterministic SIS model. There are two fixed
points: a stable endemic state and an unstable extinct state. (b) Topology of
the deterministic Allee effect model. There are three fixed points: a stable
carrying capacity, a stable extinct state, and an unstable Allee threshold.

The deterministic SIS model has two fixed points. One is an extinct state
where no infectious individuals are present, while the other is an endemic
state where the infection is maintained. The stability of these two fixed
points is determined by the value of the reproductive number $R_{0}$. The
reproductive number can be thought of as the average number of new infectious
individuals that one infectious individual generates over the course of the
infectious period, in an entirely susceptible population. For $R_{0}>1$, the
extinct state is unstable while the endemic state is stable, as shown in Fig.
14(a). Note that since the model is deterministic, a population at the
attracting endemic state can never go extinct.

The deterministic Allee effect model has three fixed points. One can see in
Fig. 14(b) that the extinct state $x_{0}$ is stable, as is the carrying
capacity $x_{2}$. The Allee threshold $x_{1}$ is unstable, so when initial
values lie between $x_{1}$ and $x_{2}$, the deterministic solution will
increase to $x_{2}$, while for initial values less than $x_{1}$ the
deterministic solution decreases to the extinct state. Similar to the
deterministic SIS, a population at the attracting carrying capacity can never
go extinct.

Figure 15:  Stochastic realization of an SIS system showing extinction after
fluctuating about the endemic state for a long time.

To capture extinction events, we must consider a stochastic model with
internal noise that represents the random interactions of individuals in the
population. Figure 15 shows a time series of infectious individuals for a
stochastic SIS model. One can see fluctuations about the endemic state for a
long period of time until eventually the noise-induced rare event (extinction)
occurs. Similarly, Fig. 2 shows a time series of individuals for a stochastic
Allee effect model. One sees the population fluctuating about the carrying
capacity for a long period of time until the rare extinction event occurs.
Details of the Monte Carlo method used in both simulations can be found in
Appendix A.

The two examples considered in this section have different mean-field
topologies as seen in Fig. 14. For the SIS model, the extinct state $x_{0}$ is
a repelling point of the deterministic mean-field equation, while the endemic
state $x_{1}$ is an attracting point. In this extinction scenario, the optimal
path to extinction is a heteroclinic trajectory with non-zero momentum that
connects the equilibrium point $(x,\lambda)=(x_{1},0)$ with a fluctuational
extinct state $(x,\lambda)=(0,\lambda_{f})$. This new fluctuational extinct
state is an equilibrium point of Hamilton’s equations given by (69). Figure 16
shows the optimal path topology in the expanded 2D space.

For the Allee effect model, the extinct state $x_{0}$ is an attracting point
of the deterministic mean-field equation. Additionally, there is an
intermediate repelling point $x_{1}$ that lies between the extinct state and
another attracting state $x_{2}$. In this extinction scenario, the optimal
path to extinction, is composed of two segments. The first segment is a
heteroclinic trajectory with non-zero momentum that connects the equilibrium
point $(x,\lambda)=(x_{2},0)$, with the intermediate equilibrium point
$(x,\lambda)=(x_{1},0)$. The second segment consists of the segment along
$\lambda=0$ from $x_{1}$ to the extinct state $x_{0}$. Figure 17 shows the
optimal path topology in the expanded 2D space.

$(x_{1},0)$$x$$(0,0)$$(0,\lambda_{f})$$\lambda_{opt}$$\lambda$ Figure 16:
Steady states of Hamilton’s equations (69) and zero-energy trajectories of the
Hamiltonian (68) for the SIS model as well as other models with deterministic
mean-field topology of the type shown in 14(a). The extinct state $x_{0}=0$ is
a repelling point of the deterministic mean-field equation, while the endemic
state $x_{1}$ is an attracting point. In this extinction scenario, the optimal
path to extinction $\lambda_{opt}$ is a heteroclinic trajectory with non-zero
momentum that connects the equilibrium point $(x,\lambda)=(x_{1},0)$ with a
fluctuational extinct state $(x,\lambda)=(0,\lambda_{f})$.
$(0,0)$$(x_{1},0)$$x$$(x_{2},0)$$\lambda_{opt}$$\lambda$ Figure 17:  Steady
states of Hamilton’s equations (69) and zero-energy trajectories of the
Hamiltonian (68) for the Allee effect model as well as other models with
deterministic mean-field topology of the type shown in 14(b). The extinct
state $x_{0}=0$ is an attracting point of the deterministic mean-field
equation. Additionally, there is an intermediate repelling point $x_{1}$ that
lies between the extinct state and another attracting state $x_{2}$. In this
extinction scenario, the optimal path to extinction, is composed of two
segments. The first segment is a heteroclinic trajectory with non-zero
momentum that connects the equilibrium point $(x,\lambda)=(x_{2},0)$, with the
intermediate equilibrium point $(x,\lambda)=(x_{1},0)$. The second segment
consists of the segment along $\lambda=0$ from $x_{1}$ to the extinct state
$x_{0}$.

####  3.3.1 SIS epidemic model

The SIS problem has mean-field equations given as

(99a) $\displaystyle\dot{S}$ $\displaystyle=\mu K-\frac{\beta}{K}SI+\gamma
I-\mu S,$ (99b) $\displaystyle\dot{I}$ $\displaystyle=\frac{\beta}{K}SI-\gamma
I-\mu I,$

where $\beta$ is the contact rate, $\gamma$ is the recovery rate, and $\mu$ is
the birth/death rate. Using a constant population assumption $S+I=K$ leads to
[59, 69, 4, 70]

(100) $\dot{I}=\frac{\beta}{K}(K-I)I-\gamma I-\mu I.$

Steady states of (100) are $I=0$, and $I=K\left(1-\frac{1}{R_{0}}\right)$,
where $R_{0}=\frac{\beta}{\mu+\gamma}$ is the reproductive number. The endemic
state is stable for $R_{0}>1$, as seen in Fig. 14(a). Therefore,
deterministically, there is no way for the disease to go extinct as mentioned
previously. However, as shown in Fig. 15, the internal noise can in fact
induce a large fluctuation which brings the population into the extinct state.
Employing the theory described in Secs. 2.4 and 2.5, one can find the optimal
path to extinction which connects the deterministic endemic state to a new
fluctuational extinct state.

We constrain the population size so that $S=K-I$. This allows one to consider
the dynamics of the constrained SIS model in terms of only the infectious
individuals $I$. Rescaling time by $\mu+\gamma$, the mean-field equation of
the constrained SIS model becomes

(101) $\dot{I}=\frac{R_{0}}{K}(K-I)I-I.$

In rescaling time by $\mu+\gamma$, the corresponding stochastic population
model is represented by the transition processes of infection, and removal due
to recovery [21, 26, 79]. The associated rates $W(X;r)$ are given as

Infection: |  $W(I;1)=(R_{0}/K)I(K-I)$
---|---
Recovery: |  $W(I;-1)=I$

and the master equation is then

(102) $\displaystyle\frac{\partial\rho(I,t)}{\partial t}$
$\displaystyle=[(I+1)\rho(I+1,t)-I\rho(I,t)]$ (103)
$\displaystyle+\frac{R_{0}}{K}[(I-1)(K-(I-1))\rho(I-1,t)-I(K-I)\rho(I,t)].$

The scaled transition rates in (64) are given as

(104) $\begin{array}[]{lll}w_{+1}(x)=R_{0}(1-i)i,&{}{}{}{}{}&w_{-1}(x)=i,\\\
&&\\\ u_{+1}(x)=0,&{}{}{}{}{}&u_{-1}(x)=0,\end{array}$

where $i=I/K$ is the fraction of infectious individuals in the population.

Substitution of (104) into (68) leads to the Hamiltonian given as

(105) $\mathcal{H}(I,\lambda)=R_{0}(1-i)i(e^{\lambda}-1)+i(e^{-\lambda}-1).$

Solutions to $\mathcal{H}(I,\lambda)=0$ are

(106) $i=0,\quad\lambda=0,\quad{\rm and}\quad\lambda(i)=-\ln{(R_{0}(1-i))}.$

The third solution is $\lambda_{opt}$ and can also be found using (71). Taking
derivatives of (105) with respect to $x$ and $\lambda$ (see (69)) leads to the
following system of Hamilton’s equations:

(107a) $\displaystyle\dot{i}$
$\displaystyle=\frac{\partial\mathcal{H}}{\partial\lambda}=R_{0}(1-i)ie^{\lambda}-ie^{-\lambda},$
(107b) $\displaystyle\dot{\lambda}$
$\displaystyle=-\frac{\partial\mathcal{H}}{\partial
i}=-R_{0}(1-2i)(e^{\lambda}-1)-(e^{-\lambda}-1).$

This system of Hamilton’s equations has three steady states given by an
extinct state $(i,\lambda)=(0,0)$, an endemic state
$(i,\lambda)=(1-1/R_{0},0)$, and a fluctuational extinct state
$(i,\lambda)=(0,-\ln{(R_{0})})$. These steady states along with the zero
energy curves of the Hamiltonian are shown in Fig. 18.

Figure 18:  Steady states of Hamilton’s equations (107) and zero-energy curves
of the SIS model Hamiltonian (105) for $R_{0}=1.5$. The extinct state is
located at $(0,0)$ (black dot), the endemic state is located at $(1/3,0)$ (red
dot), and the fluctuational extinct state is located at $(0,-\ln{(1.5)})$
(green dot). The optimal path $\lambda_{opt}=-\ln{(1.5(1-i))}$ (blue curve)
runs from the endemic state to the fluctuational extinct state.

The action along the optimal path given by (72) is

(108) $\mathcal{S}_{\rm
opt}=\int\limits_{1-\frac{1}{R_{0}}}^{0}\lambda(i)\,di=\ln{(R_{0})}-1+\frac{1}{R_{0}}.$

The general form of the MTE for single-step problems with a similar topology
to that shown in Fig. 14(a) is

(109) $\tau=\frac{\sqrt{2\pi
R}\exp\left(\int_{x_{0}}^{x_{1}}\left(\frac{u_{+1}(x)}{w_{+1}(x)}-\frac{u_{-1}(x)}{w_{-1}(x)}\right)dx\right)}{(R-1)w_{+1}(x_{1})\sqrt{K\lambda^{\prime}_{opt}(x_{1})}}\exp\left(K\int_{x_{0}}^{x_{1}}\ln\left(\frac{w_{+1}(x)}{w_{-1}(x)}\right)dx\right),$

where $R=w_{+1}^{\prime}(0)/w_{-1}^{\prime}(0)$. The SIS problem therefore has
a MTE given as

(110) $\tau=B\exp{(K\mathcal{S}_{\rm
opt})}=\frac{R_{0}}{(R_{0}-1)^{2}}\sqrt{\frac{2\pi}{N}}\exp{\left[K\left(\ln{(R_{0})}-1+\frac{1}{R_{0}}\right)\right]}.$

The analytical MTE is confirmed using numerical simulations. By numerically
computing thousands of stochastic realizations and the associated extinction
times, one can calculate the MTE. Figure 19 shows a comparison of the
analytical and numerical mean time to extinction as a function of reproductive
number $R_{0}$. There is excellent agreement except at low $R_{0}$ where the
quasi-stationary assumption breaks down.

Figure 19:  Comparison of analytical (red curve) and numerical (black circles)
mean time to extinction versus $R_{0}$ for SIS model.

####  3.3.2 Allee effect model

The Allee effect problem has a mean-field equation given as

(111) $\dot{x}=\frac{\alpha x^{2}}{2}-\mu x-\frac{\sigma x^{3}}{6},$

where $\mu$ is the death rate of low-density populations, $\alpha$ is the
growth rate when the population is large enough, and $\sigma$ is a negative
growth rate for an overcrowded population.

The steady states of (111) are $x_{0}=0$, and
$x_{1,2}=\frac{3\alpha\mp\sqrt{9\alpha^{2}-24\sigma\mu}}{2\sigma}$, where
$x_{0}$ is stable, $x_{1}$ is unstable, and $x_{2}$ is stable as shown in Fig.
14(b) [6, 71]. Therefore, deterministically, there is no way for the disease
to go extinct, as mentioned previously. However, as shown in Fig. 2 the
internal noise can in fact induce a large fluctuation which brings the
population to a vicinity of the repelling fixed point of the deterministic
rate equation. From there, the population travels essentially
deterministically to the extinct state. Employing the theory described in
Secs. 2.4 and 2.5, one can find the optimal path to extinction.

Figure 20:  Steady states of Hamilton’s equations (115) and zero-energy curves
of the Hamiltonian (113) for $\mu=0.2$, $\alpha=1.5$, and $\sigma=3$. The
extinct state is located at $(0,0)$ (black dot), the Allee threshold state is
located at $(0.35,0)$ (green dot), and the carrying capacity state is located
at $(1.15,0)$ (red dot). The optimal path
$\lambda_{opt}=-\ln{[(1.2+3x^{2})/4.5x]}$ (blue curve) runs from the carrying
capacity state to the Allee threshold state, and from there continues
deterministically to the extinct state.

The corresponding stochastic population model is represented by the following
transition processes and associated rates $W(X;r)$.

${\displaystyle X\overset{\mu}{\longrightarrow}\varnothing}$ |  $\mu X$,
---|---
${\displaystyle 2X\overset{\alpha/K}{\longrightarrow}3X}$ |  $\alpha\frac{X(X-1)}{2K}$,
${\displaystyle 3X\overset{\sigma/K^{2}}{\longrightarrow}2X}$ |  $\sigma\frac{X(X-1)(X-2)}{6K^{2}}$.

The first two transitions involving the death rate $\mu$ and the growth rate
$\alpha$ are required to capture the Allee effect. The negative growth rate
$\sigma$ allows for an overcrowded population to decline to the carrying
capacity $K$ [6, 71].

The scaled transition rates in (64) are given as

(112) $\begin{array}[]{lll}w_{+1}(x)=\frac{\alpha
x^{2}}{2},&{}{}{}{}{}&w_{-1}(x)=\mu x+\frac{\sigma x^{3}}{6},\\\ &&\\\
u_{+1}(x)=-\frac{\alpha x}{2},&{}{}{}{}{}&u_{-1}(x)=-\frac{\sigma
x^{2}}{2}.\end{array}$

Substitution of (112) into (68) leads to the Hamiltonian given as

(113) $\mathcal{H}(x,\lambda)=\frac{\alpha x^{2}}{2}(e^{\lambda}-1)+\left(\mu
x+\frac{\sigma x^{3}}{6}\right)(e^{-\lambda}-1).$

Solutions to $\mathcal{H}(x,\lambda)=0$ are

(114) $x=0,\quad\lambda=0,\quad{\rm
and}\quad\lambda(x)=\ln{\left(\frac{6\mu+\sigma x^{2}}{3\alpha x}\right)}.$

The third solution is $\lambda_{opt}$ and can also be found using (71). Taking
derivatives of (113) with respect to $x$ and $\lambda$ (see (69)) lead to the
following system of Hamilton’s equations:

(115) $\displaystyle\dot{x}$
$\displaystyle=\frac{\partial\mathcal{H}}{\partial\lambda}=\frac{\alpha
x^{2}}{2}e^{\lambda}-\left(\mu x+\frac{\sigma x^{3}}{6}\right)e^{-\lambda},$
(116) $\displaystyle\dot{\lambda}$
$\displaystyle=-\frac{\partial\mathcal{H}}{\partial x}=-\alpha
x(e^{\lambda}-1)-\left(\mu+\frac{\sigma x^{2}}{2}\right)(e^{-\lambda}-1).$

This system of Hamilton’s equations has three steady states given as
$(x,\lambda)=(x_{0},0)=(0,0)$, $(x,\lambda)=(x_{1},0)$, and
$(x,\lambda)=(x_{2},0)$ where $x_{0}$, $x_{1}$, and $x_{2}$ are the steady
states of the mean-field equation (111) provided above. These steady states
along with the zero energy curves of the Hamiltonian are shown in Fig. 20.

The general form of the MTE for single-step problems with a similar topology
to that shown in Fig. 14(b) is

(117)
$\tau=\frac{2\pi\exp\left(\int_{x_{1}}^{x_{2}}\left(\frac{u_{+1}(x)}{w_{+1}(x)}-\frac{u_{-1}(x)}{w_{-1}(x)}\right)dx\right)}{w_{1}(x_{2})\sqrt{|\lambda^{\prime}_{opt}(x_{1})|\lambda^{\prime}_{opt}(x_{2})}}\exp\left(K\int_{x_{2}}^{x_{1}}\ln\left(\frac{w_{-1}(x)}{w_{+1}(x)}\right)dx\right).$

It is worth noting that the derivation of (117) involves matching the solution
from $x_{2}$ to $x_{1}$ asymptotically with the deterministic solution from
$x_{1}$ to $x_{0}$. Because this latter solution is associated with
$\lambda=0$, its final form does not involve an integral from $x_{1}$ to
$x_{0}$. Nevertheless, the deterministic contribution is in fact included in
(117). Details regarding the derivation of the prefactor $B$ in (109) and
(117) as well as details regarding the MTE for multi-step problems can be
found in [6].

The analytical MTE for the Allee effect problem is found using (117) and is
confirmed using numerical simulations. By numerically computing thousands of
stochastic realizations and the associated extinction times, one can calculate
the MTE. Figure 21 shows the comparison between the analytical and the
numerical mean time to extinction as a function of $\alpha$ for various
choices of $\sigma$ [71].

Figure 21:  Comparison of analytical and numerical mean time to extinction
versus $\alpha$ for various $\sigma$ for the stochastic Allee effect problem.

##  4 Conclusions

The study of rare events in noise-driven physical and biological models
presents a unique opportunity to combine tools from different arenas of
applied mathematics, including variational calculus, optimal control, singular
perturbation theory, WKBJ theory, asymptotics of integrals, and sampling
techniques, to provide considerable insight into ostensibly intractable
problems. The examples included here are just a sample of the important
questions that can be answered by an intuitive combination of theory and
numerics. The applicability of these tools can only grow with the explosion of
interest in stochastic partial differential equations and geometric growth in
computing power. It is the authors’ hope that this primer facilitates the
entry of new researchers into this exciting and rapidly growing field.

##  Appendix A Numerical Methods for Stochastic Differential Equations

The simplest numerical method to generate a random walk approximating a
diffusion process is the stochastic Euler method, sometimes referred to as the
Euler-Maruyama method. Given the Itô stochastic differential equation
expressed in (5), the stochastic Euler scheme is given by

(118) $X_{n+1}=X_{n}+f(X_{n})\Delta t+\sigma(X_{n})\Delta\eta_{n},$

where the time increment is $\Delta t=t_{n+1}-t_{n}$ and the noise increment
is $\Delta\eta_{n}=\eta_{t_{n+1}}-\eta_{t_{n}}$. To implement the stochastic
Euler scheme, we note that the noise increments are independent Gaussian
random variables with the following first and second moments:

(119)
$\mathbb{E}(\Delta\eta_{n})=0,\qquad\mathbb{E}((\Delta\eta_{n})^{2})=D\Delta
t,$

so that

(120) $\Delta\eta_{n}\sim{\mathcal{N}}(0,D\Delta t).$

In practice, one can solve (118) numerically, using a random number generator
to draw noise values from the distribution given by (120). Alternatively, one
can numerically solve

(121) $X_{n+1}=X_{n}+f(X_{n})\Delta t+\sqrt{D\Delta
t}\;\sigma(X_{n})\;\hat{\eta}_{n},$

where

(122) $\hat{\eta}_{n}\sim{\mathcal{N}}(0,1).$

For the Ornstein-Uhlenbeck process expressed by $f(x)=-x$ with $D\sigma\equiv
1$, this becomes

(123) $X_{n+1}=(1-\Delta t)X_{n}+\sqrt{\Delta
t}\;\eta_{n},\qquad\eta_{n}\sim{\mathcal{N}}(0,1).$

To generate a solution of a stochastic equation, where the noise is internal,
we use the Doob-Gillespie algorithm (also known as the Gillespie algorithm or
Gillespie’s stochastic simulation algorithm (SSA)) [23, 43]. The algorithm is
a type of Monte Carlo method that was originally proposed by Kendall [53] for
simulating birth-death processes and was popularized by Gillespie [43] as a
useful method for simulating chemical reactions based on molecular collisions.
The results of a Gillespie simulation is a stochastic trajectory that
represents an exact sample from the probability function that solves the
master equation. Therefore the method can be used to simulate population
dynamics where molecular collisions are replaced by individual events and
interactions including birth, death, and infection.

Let ${\bf x}=(x_{1},\ldots,x_{n})^{T}$ denote the state variables of a system,
where $x_{i}$ provides the number of individuals in state $x_{i}$ at time $t$.
The first step of the algorithm is to initialize the number of individuals in
the population compartments ${\bf x}_{0}$. For a given state ${\bf x}$ of the
system, one calculates the transition rates (birth rate, death rate, contact
rate, etc.) denoted as $a_{i}({\bf x})$ for $i=1\ldots l$, where $l$ is the
number of transitions. Thus the sum of all transition rates is given by
$a_{0}=\sum\limits_{i=1}^{l}a_{i}({\bf x})$.

Random numbers are generated to determine both the next event to occur as well
as the time at which the next event will occur. One simulates the time $\tau$
until the next transition by drawing from an exponential distribution with
mean $1/a_{0}$. This is equivalent to drawing a random number $r_{1}$
uniformly on $(0,1)$ and computing $\tau=(1/a_{0})\ln{(1/r_{1})}$. During each
random time step exactly one event occurs. The probability of any particular
event taking place is equal to its own transition rate divided by the sum of
all transition rates $a_{i}({\bf x})/a_{0}$. A second random number $r_{2}$ is
drawn uniformly on $(0,1)$, and it is used to determine the transition event
that occurs. If $0<r_{2}<a_{1}({\bf x})/a_{0}$, then the first transition
occurs; if $a_{1}({\bf x})/a_{0}<r_{2}<(a_{1}({\bf x})+a_{2}({\bf x}))/a_{0}$,
then the second transition occurs, and so on. Lastly, both the time step and
the number of individuals in each compartment are updated, and the process is
iterated until the disease goes extinct or until the simulation time has been
exceeded.

##  Appendix B Numerical Methods for Optimal Path Computations

Although the two-point boundary value problem given by (43) and (44) in Sec.
(2.3) naturally suggests a shooting method for numerical solution, an
alternative approach used in the present work is to combine the state and
costate equations into a single second-order differential equation,

$\displaystyle\ddot{x}_{i}$ $\displaystyle=\frac{\partial
f_{i}}{\partial{x_{j}}}\dot{x}_{j}+(\dot{x}_{l}-f_{l})(a^{-1})_{lj}\frac{\partial
a_{ij}}{\partial{x_{k}}}\dot{x}_{k}-a_{ij}\frac{\partial
f_{k}}{\partial{x_{j}}}(a^{-1})_{kl}(\dot{x}_{l}-f_{l})$ (124)
$\displaystyle+\frac{1}{2}a_{ij}(\dot{x}_{k}-f_{k})\frac{\partial(a^{-1})_{kl}}{\partial{x_{j}}}(\dot{x}_{l}-f_{l}),\quad
i=1,\dots,d,$

where $a:=\sigma\sigma^{T}$ and we have used the summation convention over
doubled indices. The supplementary boundary conditions for $\lambda$ are re-
expressed as Robin conditions for $x$ at $t=t_{\mathrm{f}}$. Introduction of
an artificial time $T$ allows a simple initial condition for $x(t,T)$,
generally a linear interpolant, to relax to the optimal path through

$\displaystyle\frac{\partial x_{i}}{\partial{T}}=$
$\displaystyle\frac{\partial^{2}{x}_{i}}{\partial{t}^{2}}-\frac{\partial
f_{i}}{\partial{x_{j}}}\frac{\partial{x}_{j}}{\partial{t}}-\left(\frac{\partial{x}_{l}}{\partial{t}}-f_{l}\right)(a^{-1})_{lj}\frac{\partial
a_{ij}}{\partial{x_{k}}}\frac{\partial{x}_{k}}{\partial{t}}+a_{ij}\frac{\partial
f_{k}}{\partial{x_{j}}}(a^{-1})_{kl}\left(\frac{\partial{x}_{l}}{\partial{t}}-f_{l}\right)$
(125)
$\displaystyle-\frac{1}{2}a_{ij}\left(\frac{\partial{x}_{k}}{\partial{t}}-f_{k}\right)\frac{\partial(a^{-1})_{kl}}{\partial{x_{j}}}\left(\frac{\partial{x}_{l}}{\partial{t}}-f_{l}\right).$

In the present work, this PDE was solved using central differences in space
and a semi-implicit method in time, with the second derivative handled
implicitly and the remaining terms on the right-hand side of (125) handled
explicitly. The PDE was evolved in artificial time $T$ until the maximum of
the residual dropped below threshold.

The relaxation described in 125 occurs over curves parameterized by true time
$t$. This parameterization is not convenient in cases where the true time of
travel over the trajectory is infinite, such as for quasi-potential or mean
exit time computations, or where the speed is highly nonuniform and the action
is therefore concentrated in a particular segment of the path. It is
convenient in such cases to re-parameterize the curve in terms of arclength.
Following Ref. [49], we see that the action

(126) $\displaystyle{\mathcal{S}}_{t_{\mathrm{f}}}$
$\displaystyle=\frac{1}{2}\int_{t_{0}}^{t_{\mathrm{f}}}\|\dot{x}-f(x)\|^{2}\,dt$
(127)
$\displaystyle=\frac{1}{2}\int_{t_{0}}^{t_{\mathrm{f}}}\left((\|\dot{x}\|-\|f\|)^{2}+2\|\dot{x}\|\|f\|-2(\dot{x},f)\right)\,dt$
(128)
$\displaystyle\geq\int_{t_{0}}^{t_{\mathrm{f}}}\left(\|\dot{x}\|\|f\|-(\dot{x},f)\right)\,dt$
(129)
$\displaystyle=2\int_{t_{0}}^{t_{\mathrm{f}}}\|\dot{x}\|\|f\|\sin^{2}\frac{1}{2}\theta(t)\,dt$

where $\theta$ is the angle between $\dot{x}$ and $f(x)$. In the case of
infinite-time paths, time can be rescaled to force

(130) $\|\dot{x}\|=\|f(x)\|,$

so that the above inequality becomes an equality, and

(131)
$V(0,y)=2\inf_{\gamma}\int_{\gamma}\|f(x(\alpha))\|\sin^{2}\frac{1}{2}\theta(\alpha)\,d\alpha,$

i.e., the time-parameterized minimization can be replaced by a minimization
using arclength as a parameter, where in (131) the admissible paths are
absolutely continuous. This approach to the numerical computation of Wentzell-
Freidlin action minimizers is referred to as the geometric minimum-action
method (GMAM), and requires a modification of 125 to simultaneously compute
the time parameterization. The reader is referred to Ref. [49] for more
details.

A similar method that employs a direct, fully explicit iterative scheme is the
iterative action minimization method (IAMM) [61]. The IAMM is useful in the
general situation where a path connecting steady states $C_{a}$ and $C_{b}$
starts at $C_{a}$ at $t=-\infty$ and ends at $C_{b}$ at $t=+\infty$. A time
parameter $t$ exists such that $-\infty<t<\infty$. For this method, we require
a numerical approximation of the time needed to leave the region of $C_{a}$
and arrive in the region of $C_{b}$. Therefore, we define a time
$T_{\epsilon}$ such that $-\infty<-T_{\epsilon}<t<T_{\epsilon}<\infty$.
Additionally, $C(-T_{\epsilon})\approx C_{a}$ and $C(T_{\epsilon})\approx
C_{b}$. In other words, the solution stays very near the equilibrium $C_{a}$
for $-\infty<t\leq-T_{\epsilon}$, has a transition region from
$-T_{\epsilon}<t<T_{\epsilon}$, and then stays near $C_{b}$ for
$T_{\epsilon}<t<+\infty$. The interval $[-T_{\epsilon},T_{\epsilon}]$ is
discretized into $n$ segments using a uniform step size $h=(2T_{\epsilon})/n$
or a suitable non-uniform step size $h_{k}$. The corresponding time series is
$t_{k+1}=t_{k}+h_{k}$.

The derivative of the function value $\mathbf{q}_{k}$ is approximated using
central finite differences by the operator $\delta_{h}$ given as

(132)
$\frac{d}{dt}\mathbf{q}_{k}\approx\delta_{h}\mathbf{q}_{k}\equiv\frac{h^{2}_{k-1}\mathbf{q}_{k+1}+(h^{2}_{k}-h^{2}_{k-1})\mathbf{q}_{k}-h^{2}_{k}\mathbf{q}_{k-1}}{h_{k-1}h^{2}_{k}+h_{k}h^{2}_{k-1}},\quad
k=0,\ldots,n.$

Clearly, if a uniform step size is chosen then Eq. (132) simplifies to the
familiar form given as

(133)
$\frac{d}{dt}\mathbf{q}_{k}\approx\delta_{h}\mathbf{q}_{k}\equiv\frac{\mathbf{q}_{k+1}-\mathbf{q}_{k-1}}{2h},\quad
k=0,\ldots,n.$

Thus, one can develop the system of nonlinear algebraic equations

(134) $\delta_{h}\mathbf{x}_{k}-\frac{\partial
H(\mathbf{x}_{k},\mathbf{p}_{k})}{\partial\mathbf{p}}=0,\quad\delta_{h}\mathbf{p}_{k}+\frac{\partial
H(\mathbf{x}_{k},\mathbf{p}_{k})}{\partial\mathbf{x}}=0,\quad k=0,\ldots,n,$

which is solved using a general Newton’s method. We let

(135)
$\mathbf{q}_{j}(\mathbf{x,p})=\\{\mathbf{x}_{1,j}...\mathbf{x}_{n,j},\mathbf{p}_{1,j}...\mathbf{p}_{n,j}\\}^{T}$

be an extended vector of $2nN$ components that contains the $j^{\rm th}$
Newton iterate, where $N$ is the number of populations. When $j=0$,
$\mathbf{q}_{0}(\mathbf{x,p})$ provides the initial “guess” as to the location
of the path that connects $C_{a}$ and $C_{b}$. Given the $j^{\rm th}$ Newton
iterate $\mathbf{q}_{j}$, the new $\mathbf{q}_{j+1}$ iterate is found by
solving the linear system

(136)
$\mathbf{q}_{j+1}=\mathbf{q}_{j}-\frac{\mathbf{F}\left(\mathbf{q}_{j}\right)}{\mathbf{J}\left(\mathbf{q}_{j}\right)},$

where $\mathbf{F}$ is the function defined by Eq. (134) acting on
$\mathbf{q}_{j}$, and $\mathbf{J}$ is the Jacobian. Equation (136) is solved
using LU decomposition with partial pivoting.

## Acknowledgments

We warmly thank Lora Billings for her insight and Michael Khasin for reading
and commenting on an early version of the manuscript. We also gratefully
acknowledge support from the National Science Foundation awards CMMI-1233397
(EF), DMS-1418956 (EF), and DMS-1109278 (ROM).

## References

  * [1] D. A. Adams, L. M. Sander, and R. M. Ziff, The barrier method: A technique for calculating very long transition times, The Journal of Chemical Physics, 133 (2010), p. 124103.
  * [2] W. C. Allee, Animal Aggregations, a Study in General Sociology, Univ. Chicago Press, 1931.
  * [3] D. Anderson, M. Lisak, and A. Berntson, A variational approach to nonlinear evolution equations in optics, Pramana, 57 (2001), pp. 917–936.
  * [4] H. Andersson and B. Djehiche, A threshold limit theorem for the stochastic logistic epidemic, J. App.. Prob., 35 (1998), pp. 662–670.
  * [5] M. Assaf and B. Meerson, Spectral theory of metastability and extinction in a branching-annihilation reaction, Phys. Rev. E, 75 (2007), p. 031122.
  * [6] , Extinction of metastable stochastic populations, Phys. Rev. E, 81 (2010), p. 021116.
  * [7] M. Assaf, E. Roberts, and Z. Luthey-Schulten, Determining the stability of genetic switches: Explicitly accounting for mRNA noise, Phys. Rev. Lett., 106 (2011), p. 248102.
  * [8] B. G. Bale and J. N. Kutz, Variational method for mode-locked lasers, Journal of the Optical Society of America B, 25 (2008), p. 1193.
  * [9] C. Bao, A. C. Funk, C. Yang, and S. T. Cundiff, Pulse dynamics in a mode-locked fiber laser and its quantum limited comb frequency uncertainty, Optics Letters, 39 (2014), p. 3266.
  * [10] M. Bauver, E. Forgoston, and L. Billings, Computing the optimal path in stochastic dynamical systems, Chaos, 26 (2016), p. 083101.
  * [11] L. Billings and E. Forgoston, Seasonal forcing in stochastic epidemiology models, Ricerche di Matematica, accepted (2017).
  * [12] B. Z. Bobrovsky and Z. Schuss, A Singular Perturbation Method for the Computation of the Mean First Passage Time in a Nonlinear Filter, SIAM Journal on Applied Mathematics, 42 (1982), pp. 174–187.
  * [13] A. E. Bryson and Y.-C. Ho, Applied optimal control: optimization, estimation, and control, Taylor & Francis, Levittown, PA, 2001.
  * [14] D. S. Cargill and R. O. Moore, Incorporating Radiation in Noise-Induced Phase Evolution of Optical Solitons, SIAM Journal on Applied Dynamical Systems, 15 (2016), pp. 1025–1061.
  * [15] D. S. Cargill, R. O. Moore, and C. J. McKinstrie, Noise bandwidth dependence of soliton phase in simulations of stochastic nonlinear Schrödinger equations, Opt. Lett., 36 (2011), pp. 1659–1661.
  * [16] A. Chernykh and M. Stepanov, Large negative velocity gradients in burgers turbulence, Physical Review E, 64 (2001), p. 026306.
  * [17] W. Cousins and T. P. Sapsis, Reduced-order precursors of rare events in unidirectional nonlinear water waves, Journal of Fluid Mechanics, 790 (2016), pp. 368–388.
  * [18] G. Da Prato and J. Zabczyk, Stochastic equations in infinite dimensions, no. 45 in Encyclopedia of mathematics and its applications, Cambridge Univ. Press, Cambridge, digitally printed version (with corrections) ed., 2008.  OCLC: 612518037.
  * [19] A. de Bouard and A. Debussche, The Stochastic Nonlinear Schrodinger Equation in H1, Stochastic Analysis & Applications, 21 (2003), p. 97.
  * [20] E. Desurvire, Erbium-doped fiber amplifiers: principles and applications, Wiley-Interscience, Hoboken, N.J, 2002.
  * [21] C. R. Doering, K. V. Sargsyan, and L. M. Sander, Extinction times for birth-death processes: Exact results, continuum asymptotics, and the failure of the Fokker-Planck approximation, Multiscale Model. Simul., 3 (2005), pp. 283–299.
  * [22] G. M. Donovan and W. L. Kath, Rare event simulation of the performance of an actively mode-locked fiber laser model, in Photonic Applications Systems Technologies Conference, 2007.
  * [23] J. L. Doob, Topics in the theory of markoff chains, Transactions of the American Mathematical Society, 52 (1942), pp. 37–64.
  * [24] M. I. Dykman and M. A. Krivoglaz, Theory of fluctuational transitions between stable states of a nonlinear oscillator, Sov. Phys. JETP, 50 (1979), pp. 30–37.
  * [25] M. I. Dykman, E. Mori, J. Ross, and P. M. Hunt, Large fluctuations and optimal paths in chemical-kinetics, J. Chem. Phys., 100 (1994), pp. 5735–5750.
  * [26] M. I. Dykman, I. B. Schwartz, and A. S. Landsman, Disease extinction in the presence of random vaccination, Phys. Rev. Lett., 101 (2008), p. 078101.
  * [27] W. E, W. Ren, and E. Vanden-Eijnden, Simplified and improved string method for computing the minimum energy paths in barrier-crossing events, The Journal of Chemical Physics, 126 (2007), p. 164103.
  * [28] V. Elgart and A. Kamenev, Rare event statistics in reaction-diffusion systems, Phys. Rev. E, 70 (2004), p. 041106.
  * [29] G. E. Falkovich, I. Kolokolov, V. Lebedev, and S. K. Turitsyn, Statistics of soliton-bearing systems with additive noise, Physical Review E, 63 (2001).
  * [30] M. Farazmand and T. P. Sapsis, A variational approach to probing extreme events in turbulent dynamical systems, Science Advances, 3.9 (2017), p. e1701533.
  * [31] R. P. Feynman and A. R. Hibbs, Quantum Mechanics and Path Integrals, McGraw-Hill, Inc., 1965.
  * [32] E. Forgoston, S. Bianco, L. B. Shaw, and I. B. Schwartz, Maximal sensitive dependence and the optimal path to epidemic extinction, B. Math. Biol., 73 (2011), pp. 495–514.
  * [33] E. Forgoston, L. Billings, P. Yecko, and I. B. Schwartz, Set-based corral control in stochastic dynamical systems: Making almost invariant sets more invariant, Chaos, 21 (2011), p. 013116.
  * [34] E. Forgoston and I. B. Schwartz, Escape rates in a stochastic environment with multiple scales, SIAM J. Appl. Dyn. Syst., 8 (2009), pp. 1190–1217.
  * [35] M. I. Freidlin and A. D. Wentzell, On small random perturbations of dynamical system, Russian Math. Surveys, 25 (1970), pp. 1–55.
  * [36] , Random Perturbations of Dynamical Systems, Springer-Verlag, 1984.
  * [37] T. A. Fulton and L. N. Dunkelberger, Lifetime of the zero-voltage state in Josephson tunnel junctions, Phys. Rev. B, 9 (1974), p. 4760.
  * [38] H. Gang, Stationary solution of master equations in the large-system-size limit, Phys. Rev. A, 36 (1987), pp. 5782–5790.
  * [39] C. W. Gardiner, Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences, Springer-Verlag, 2004.
  * [40] E. Gautier, Large deviations and support results for nonlinear Schrodinger equations with additive noise and applications, ESAIM: Probability and Statistics, 9 (2005), pp. 74–97.
  * [41] E. Gautier, Uniform large deviations for the nonlinear Schrodinger equation with multiplicative noise, Stochastic processes and their applications, 115 (2005), pp. 1904–1927.
  * [42] B. Gaveau, M. Moreau, and J. Toth, Decay of the metastable state in a chemical system: Different predictions between discrete and continuous models, Lett. Math. Phys., 37 (1996), pp. 285–292.
  * [43] D. T. Gillespie, Exact stochastic simulation of coupled chemical reactions, The Journal of Physical Chemistry, 81 (1977), pp. 2340–2361.
  * [44] P. Hänggi, P. Talkner, and M. Borkovec, Reaction-rate theory: fifty years after kramers, Reviews of modern physics, 62 (1990), p. 251.
  * [45] M. Hairer and H. Shen, A central limit theorem for the KPZ equation, The Annals of Probability, 45 (2017), pp. 4167–4221.
  * [46] M. Hairer and H. Weber, Large deviations for white-noise driven, nonlinear stochastic PDEs in two and three dimensions, Annales de la Faculte des Sciences de Toulouse Mathematiques, 24 (2015), pp. 55–92.
  * [47] C. Hartmann and C. Schütte, Efficient rare event simulation by optimal nonequilibrium forcing, Journal of Statistical Mechanics: Theory and Experiment, 2012 (2012), p. P11004.
  * [48] C. R. Heckman, M. A. Hsieh, and I. B. Schwartz, Going with the flow: enhancing switching rates in multi-gyre systems. accepted to the asme journal of dynamic systems, Measurement and Control, (2014).
  * [49] M. Heymann and E. Vanden-Eijnden, The geometric minimum action method: A least action principle on the space of curves, Communications on Pure and Applied Mathematics, 61 (2008), pp. 1052–1117.
  * [50] E. Iannone, F. Matera, A. Mecozzi, and M. Settembre, Nonlinear Optical Communications Networks, John Wiley & Sons, Inc., New York, 1998\.
  * [51] D. J. Jones, S. A. Diddams, J. K. Ranka, A. Stentz, R. S. Windeler, J. L. Hall, and S. T. Cundiff, Carrier-envelope phase control of femtosecond mode-locked lasers and direct optical frequency synthesis, Science, 288 (2000), pp. 635–639.
  * [52] H. Jónsson, G. Mills, and K. W. Jacobsen, Nudged elastic band method for finding minimum energy paths of transitions, 1998.
  * [53] D. G. Kendall, An artificial realization of a simple” birth-and-death” process, Journal of the Royal Statistical Society. Series B (Methodological), 12 (1950), pp. 116–119.
  * [54] D. A. Kessler and N. M. Shnerb, Extinction rates for fluctuation-induced metastabilities: A real space WKB approach, J. Stat. Phys., 127 (2007), pp. 861–886.
  * [55] E. Khain, M. Khasin, and L. M. Sander, Spontaneous formation of large clusters in a lattice gas above the critical point, Phys. Rev. E, 90 (2014), p. 062702.
  * [56] Y. S. Kivshar and B. A. Malomed, Dynamics of solitons in nearly integrable systems, Rev. Mod. Phys., 61 (1989), pp. 763–915.
  * [57] R. V. Kohn, M. G. Reznikoff, and E. Vanden-Eijnden, Magnetic elements at finite temperature and large deviation theory, J. Nonlinear Sci., 15 (2005), pp. 223–253.
  * [58] H. A. Kramers, Brownian motion in a field of force and the diffusion model of chemical reactions, Physica, 7 (1940), pp. 284–304.
  * [59] R. J. Kryscio and C. Lefévre, On the extinction of the s-i-s stochastic logistic epidemic, J. App. Prob., 26 (1989), pp. 685–694.
  * [60] R. Kubo, K. Matsuo, and K. Kitahara, Fluctuation and relaxation of macrovariables, J. Stat. Phys., 9 (1973), pp. 51–96.
  * [61] B. S. Lindley and I. B. Schwartz, An iterative action minimizing method for computing optimal paths in stochastic dynamical systems, Physica D, 255 (2013), pp. 22–30.
  * [62] D. G. Luchinsky, P. V. E. McClintock, and M. I. Dykman, Analogue studies of nonlinear systems, Rep. Prog. Phys., 61 (1998), pp. 889–997.
  * [63] B. J. Matkowsky and Z. Schuss, The exit problem for randomly perturbed dynamical systems, SIAM Journal on Applied Mathematics, 33 (1977), pp. 365–382.
  * [64] M. A. Mohamad, W. Cousins, and T. P. Sapsis, A probabilistic decomposition-synthesis method for the quantification of rare events due to internal instabilities, Journal of Computational Physics, 322 (2016), pp. 288–308.
  * [65] M. A. Mohamad and T. P. Sapsis, Probabilistic response and rare events in Mathieu’s equation under correlated parametric excitation, Ocean Engineering, 120 (2016), pp. 289–297.
  * [66] R. O. Moore, Trade-off between linewidth and slip rate in a mode-locked laser model, Optics Letters, 39 (2014), p. 3042.
  * [67] R. O. Moore, G. Biondini, and W. L. Kath, A Method to Compute Statistics of Large, Noise-Induced Perturbations of Nonlinear Schroedinger Solitons, SIAM Review, 50 (2008), p. 523.
  * [68] R. O. Moore, T. Schaefer, and C. K. Jones, Soliton broadening under random dispersion fluctuations: Importance sampling based on low-dimensional reductions, Optics Communications, 256 (2005), pp. 439–450.
  * [69] I. Nåsell, The quasi-stationary distribution of the closed endemic sis model, Advances App. Prob., 28 (1996), pp. 895–932.
  * [70] , On the time to extinction in recurrent epidemics, J. R. Statist. Soc. B, 61 (1999), pp. 309–330.
  * [71] G. Nieddu, L. Billings, and E. Forgoston, Analysis and control of pre-extinction dynamics in stochastic populations, B. Math. Biol., 76 (2014), pp. 3122–3137.
  * [72] G. T. Nieddu, L. Billings, J. H. Kaufman, E. Forgoston, and S. Bianco, Extinction pathways and outbreak vulnerability in a stochastic ebola model, Journal of The Royal Society Interface, 14 (2017), p. 20160847.
  * [73] B. $\O{}$ksendal, Stochastic differential equations: an introduction with applications, Universitext, Springer, Berlin, 6. ed., 6. corrected printing ed., 2013.  OCLC: 935584333.
  * [74] O. Ovaskainen and B. Meerson, Stochastic models of population extinction, Trends in Ecol. and Evol., 25 (2010), pp. 643–652.
  * [75] R. Paschotta, A. Schlatter, S. Zeller, H. Telle, and U. Keller, Optical phase noise and carrier-envelope offset noise of mode-locked lasers, Applied Physics B, 82 (2005), pp. 265–273.
  * [76] G. A. Pavliotis, Stochastic processes and applications, Springer, 2011\.
  * [77] H. Risken, The Fokker-Planck Equation, Springer-Verlag, 1996.
  * [78] E. Roberts, S. Be’er, C. Bohrer, R. Sharma, and M. Assaf, Dynamics of simple gene-network motifs subject to extrinsic fluctuations, Phys. Rev. E, 92 (2015), p. 062717.
  * [79] I. B. Schwartz, L. Billings, M. Dykman, and A. Landsman, Predicting extinction rates in stochastic epidemic models, J. Stat. Mech.-Theory E., (2009), p. P01005.
  * [80] I. B. Schwartz, E. Forgoston, S. Bianco, and L. B. Shaw, Converging towards the optimal path to extinction, J. R. Soc. Interface, 8 (2011), pp. 1699–1707.
  * [81] M. Shtaif and C. Menyuk, Carrier envelope phase locked modelocking in fiber lasers at ultra-high repetition rates, in CLEO: Science and Innovations, Optical Society of America, 2011, p. CMD6.
  * [82] R. C. Smith, Uncertainty quantification: theory, implementation, and applications, Computational science and engineering series, Society for Industrial and Applied Mathematics, Philadelphia, 2013.
  * [83] E. T. Spiller, W. L. Kath, R. O. Moore, and C. J. McKinstrie, Computing large signal distortions and bit-error ratios in DPSK transmission systems, IEEE Photonics Technology Letters, 17 (2005), pp. 1022 –1024.
  * [84] H. Touchette, The large deviation approach to statistical mechanics, Physics Reports, 478 (2009), pp. 1–69.
  * [85] L. S. Tsimring, Noise in biology, Reports on Progress in Physics, 77 (2014), p. 026601.
  * [86] N. G. van Kampen, Stochastic Processes in Physics and Chemistry, Elsevier, 2007.
  * [87] E. Vanden-Eijnden and J. Weare, Data assimilation in the low noise, accurate observation regime with application to the Kuroshio current, arXiv:1202.4952, (2012).
  * [88] , Rare event simulation of small noise diffusions, Communications on Pure and Applied Mathematics, 65 (2012), pp. 1770–1803.
  * [89] J. Ye, H. Schnatz, and L. Hollberg, Optical frequency combs: From frequency metrology to optical phase control, IEEE Journal of Selected Topics in Quantum Electronics, 9 (2003), pp. 1041–1058.
  * [90] W. Zhang, H. Wang, C. Hartmann, M. Weber, and C. Schuette, Applications of the cross-entropy method to importance sampling and optimal control of diffusions, SIAM Journal on Scientific Computing, 36 (2014), pp. A2654–A2672.