# Hamilton-Jacobi-Bellman equations for Chemical Reaction Networks

Michael Snarski

###### Abstract

This is an expository note on large deviations, Hamilton-Jacobi-Bellman (HJB)
equations, and the role of the Freidlin-Wentzell quasipotential in chemical
reaction networks (CRNs). The note was motivated by observations in [1] which
identified Lyapunov functions for deterministic descriptions of CRNs by taking
appropriate scaling limits of invariant distributions for the corresponding
stochastic dynamics. We explain how this is a special case of a classical
theory due to Freidlin and Wentzell [2]. We also show that this Lyapunov
function is a solution to the HJB partial differential equation if and only if
the network is “complex-balanced”. The target audience are researchers in the
CRN community who are familiar with the Markov process description of CRNs,
but do not have the time to invest in learning the technical machinery of
large deviations. We conclude by exploring some possible relationships with
Kähler geometry which suggest an interesting and unexplored research
direction.

##  1 Introduction

Chemical reaction networks (CRNs) attempt to model the behavior of chemical
species interacting in a well-stirred homogeneous medium of large but finite
volume. Inflow and outflow of species is allowed, and we assume that the time-
scale over which particles meet is “sufficiently small” so that they are
guaranteed to interact “often”. Under such assumptions every particle of the
same species experiences identical conditions. As such, CRNs do not track the
spatial location of any individual particle but rather provide only
macroscopic information. In this note, the macroscopic description will be
either a discrete number of a certain type of species (if a _stochastic_
description of the system is employed) or a concentration value, which is a
non-negative real number (for a _deterministic_ description).

Our discrete, stochastic description of CRNs is a sequence of jump Markov
processes taking values in increasingly finer lattice subsets of
$\mathbb{R}^{n}_{+}$. The sequence is indexed by a scaling parameter that
represents the volume of the container. We will denote the scaling parameter
by $n\geqslant 1$, and the sequence of processes by
$\\{X^{n}(\cdot)\\}_{n\geqslant 1}$. Under appropriate conditions on the jump
intensities, this sequence has a Law of Large Numbers [3][4] limit which we
denote by $x(\cdot)$, a result we sometimes refer to as _Kurtz’s theorem_. The
limit $x(\cdot)$ is absolutely continuous and can be characterized as the
unique solution to a system of coupled ordinary differential equations; this
defines the continuous, deterministic description [4].

The focus of this paper is to explain some aspects of the relationship between
the stochastic and deterministic descriptions. An important class of networks
for us will be the _complex-balanced_ ones, introduced by Horn and Jackson
[5]. They showed that if $c\in\mathbb{R}^{d}_{>0}$ is a steady state for the
deterministic dynamics of a complex-balanced network, then the function of $x$
defined by

$V(c,x)\doteq\sum_{i=1}^{d}c_{i}\ell\left(\frac{x_{i}}{c_{i}}\right),\quad\ell(x)\doteq
x\log x-x+1$ (1.1)

is a Lyapunov function for the deterministic dynamics.

It is already known that there is a relationship between the Lyapunov function
$V$ and the stochastic dynamics. In [1], the authors explicitly write down the
stationary distribution associated to the Markov process $X^{n}$ using results
from [6], and then take an appropriate scaling limit of the stationary
distributions to obtain $V$.

It turns out that this procedure is a special case of a well-established
result in a subfield of probability known as Freidlin-Wentzell theory. The
trouble with this theory is that it can seem prohibitively technical to the
non-specialist. The goal of the present paper is to explain in the simplest
terms some implications of the theory, communicate the intuition behind it,
and provide ways of leveraging the results without knowing all the details.

Here is one contribution of the approach. We will associate to any CRN a
function $L:\mathbb{R}^{d}_{+}\times\mathbb{R}^{d}$, called the _local rate_
or _Lagrangian_ , and study the family of variational problems

$Q(c,x)=\inf_{\varphi}\left\\{\int_{0}^{T}L(\varphi(s),\dot{\varphi}(s))\text{ds}:\varphi(0)=c,\varphi(T)=x,T<\infty\right\\},$
(1.2)

where the infimum is over all absolutely continuous paths
$\varphi:[0,T]\rightarrow\mathbb{R}^{d}_{+}$. We will show that $Q(c,x)$ is
always a Lyapunov function for the deterministic dynamics, and $Q(c,x)=V(c,x)$
if and only if the network is complex-balanced.

A key part of our exposition is to emphasize a “Hamiltonian description” of a
chemical reaction network. In the case of mass-action kinetics, there is a
one-to-one correspondence between reaction networks and their naturally
associated Hamiltonians. In turn, each Hamiltonian defines an associated
partial differential equation, known as a Hamilton-Jacobi-Bellman (HJB)
equation. We will demonstrate how the steady-state HJB can be used to study
the stability of both the stochastic and deterministic descriptions.

The layout of the paper is as follows. Section 2 serves to define the most
basic objects in the paper, the CRNs and their associated Hamiltonian. In
Section 3 we sketch the basic framework connecting CRNs to large deviations
and classical mechanics. We do so in three steps: in Section 3.2 we show how
the Hamiltonian encodes the Law of Large Numbers (LLN) limit; in Section 3.3
we study the rate of convergence of the LLN and introduce the notion of a
Large Deviations Principle (LDP); and in Section 3.4 we show how the LDP is
related to the HJB equation through a family of variational problems.

Once we have established the main objects of large deviations theory, we focus
in Section 4 on a specific object known as the _quasipotential_ , defined
above in (1.2). In Section 4.2 we show that the quasipotential serves as a
Lyapunov function, and in Section 4.3 we show that the invariant measure
technique employed in [1] identifies the quasipotential. It is therefore
automatic that all the functions identified in [1] are Lyapunov functions. In
Section 4.4, we show that the quasipotential can be identified using a partial
differential equation known as the Hamilton-Jacobi-Bellman (HJB) equation.

The results described so far hold for general reaction networks with mass-
action kinetics. In Section 5 we specialize the results to complex-balanced
networks and show that the function $V$ solves the HJB equation if and only if
the network is complex-balanced. Finally, Section 6 shows some worked
examples. Section 6.1 aims to demonstrate how the HJB technique can be simpler
than the one employed in [1]. On the other hand, Section 6.2 shows a curious
coincidence between the form of $V$ and the form of a potential for a Kähler
metric obtained from the combinatorial data of a polytope in [7]. We suspect
that there is a deeper connection between the two, and that this is an
interesting area of research at the intersection of toric geometry and convex
analysis. In particular, the fan description of a toric variety may be closely
related to the fans used in work of Craciun [8].

We re-iterate that this paper serves as an explanatory guide, and contains
only a few novel observations. In particular, we do not go into detail about
the theory of large deviations or viscosity solutions to HJB equations. A good
rigorous reference on the geometry of CRNs within the context of Freidlin-
Wentzell theory is [9]. The definitive book on Freidlin-Wentzell theory is
[2]. The topic of solutions to HJB equations is a vast and subtle one and
relies heavily on the seminal work of Crandall and Lions on viscosity
solutions [10] (see also [11] and [12]).

##  2 CRNs and their Hamiltonians

In this section we define CRNs and their associated Hamiltonians for the case
of mass-action kinetics.

###  2.1 Reaction networks

######  Definition 2.1.

Let $\mathcal{A}=\\{A_{1},\ldots,A_{d}\\}$, let $\mathcal{C}$ be a finite
subset of $\mathbb{Z}^{d}_{+}$, and let
$\mathcal{R}\subset\mathcal{C}\times\mathcal{C}$ be a subset of ordered pairs
of complexes. To each ordered pair $(a,b)\in\mathcal{R}$ we associate the
vector $\zeta_{a,b}\doteq b-a\in\mathbb{Z}^{d}$. The elements
$A_{i}\in\mathcal{A}$ are called _species_ , the $a\in\mathcal{C}$ are called
_complexes_ , the vectors $\zeta_{a,b}$ are called _reaction vectors_ , and
the tuple $G=(\mathcal{A},\mathcal{C},\mathcal{R})$ is called a _chemical
reaction network (CRN)_.

######  Definition 2.2.

Let $G=(\mathcal{A},\mathcal{C},\mathcal{R})$ be a CRN. The subspace

$\mathcal{S}_{0}=\left\\{\sum_{(a,b)\in\mathcal{R}}c_{a,b}\zeta_{a,b}:c_{a,b}\in\mathbb{R},(a,b)\in\mathcal{R}\right\\}$

is called the _stoichiometric subspace_ of the network. For a given initial
condition $x(0)\in\mathbb{R}^{d}_{\geqslant 0}$, we let

$\mathcal{S}_{x(0)}\doteq\mathbb{R}^{d}_{\geqslant
0}\cap(\mathcal{S}_{0}+x(0))=\mathbb{R}^{d}_{\geqslant
0}\cap\\{s+x(0):s\in\mathcal{S}_{0}\\},$

and we refer to this as the stoichiometric compatibility class passing through
$x(0)$, or simply as the domain of the network. When the particular initial
condition $x(0)$ is irrelevant, or when $\mathcal{S}_{x}$ is the same for all
choices of $x$, we will denote the domain by $\mathcal{S}$.

We will use the following notation throughout.

######  Convention 2.3.

We use the convention of exponential notation for $x\in\mathbb{R}^{d}_{+}$ and
$a\in\mathbb{Z}^{d}_{+}$

$x^{a}=x_{1}^{a_{1}}\cdots x_{d}^{a_{d}}.$

######  Remark 2.4.

We will strive to use $x$ to denote a point in
$\mathcal{S}\subset\mathbb{R}^{d}$ and $x(\cdot)$ to denote a path
$x(\cdot)\doteq\\{x(t),t\in[0,T]\\}$. We do not call this a convention,
because we may not have been stringent enough in enforcing it, and hopefully
no confusion will arise.

To a CRN $G=(\mathcal{A},\mathcal{C},\mathcal{R})$ we associate _rate
constants_ $\kappa=\\{\kappa_{a,b}>0:(a,b)\in\mathcal{R}\\}$ which are
necessary to define the stochastic and deterministic dynamics. Given a
reaction network $G$ and rate constants $\kappa$, we associate an object we
will call the _Hamiltonian_ of the reaction network. This is an object which
is central in providing the link between the stochastic and deterministic
dynamics.

######  Definition 2.5.

Let $G=(\mathcal{A},\mathcal{C},\mathcal{R})$ be a CRN with associated rate
constants $\kappa=\\{\kappa_{a,b}>0:(a,b)\in\mathcal{R}\\}$. We define the
_Hamiltonian_ associated to $G$ by

$H_{G}(x,p)\doteq\sum_{(a,b)\in\mathcal{R}}\kappa_{a,b}x^{a}(e^{\langle\zeta_{a,b},p\rangle}-1).$
(2.1)

The set of all functions
$H:\mathbb{R}^{d}_{+}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ of the form
(2.1) with $\kappa_{a,b}>0$, $a\in\mathbb{Z}^{d}_{+}$ and
$\zeta_{a,b}\in\mathbb{Z}^{d}$ is in one-to-one correspondence with the set of
CRNs with mass-action kinetics.

######  Example 2.1.

Consider the one-dimensional reaction network with a single species
$\mathcal{A}=\\{A\\}$, two complexes $\mathcal{C}=\\{0,1\\}$ and two reactions
$\mathcal{R}=\\{(0,1),(1,0)\\}$. To this CRN we associate the rate constants
$\kappa_{0,1}=\kappa_{1,0}=1$. This corresponds to the reaction network

$A\rightleftarrows\emptyset$

where both rate constants are equal to 1. The associated Hamiltonian is

$H(x,p)=x(e^{-p}-1)+e^{p}-1.$

We will refer to Example 2.1 to illustrate various objects introduced in this
note.

##  3 From CRNs to Hamilton-Jacobi-Bellman

In this section we sketch the connection between CRNs and the Hamilton-Jacobi-
Bellman equation by appealing to large deviations theory. We begin in Section
3.1 by providing a roadmap for the reader which identifies the important
definitions and the equations relating them. The key idea is that in the rare
event regime, the qualitative prediction of an event is governed
asymptotically by its cost. The difficulty is that evaluating the cost of an
event is in general a non-convex infinite-dimensional optimization problem.

The remaining subsections provide greater detail on parts of the roadmap.
Section 3.2 focuses on the law of large numbers / deterministic dynamics in
terms of the Hamiltonian; Section 3.3 provides a brief overview of large
deviations theory, and 3.4 explores the calculus of variations and HJB
approaches to the optimization problem.

In the end we cover only the very basic concepts and omit many important
technical details, such as the existence and uniqueness of minimizers of
variational problems and the notion of spatial derivatives on the boundary.
The interested reader can consult [13] or [14] for an introduction to large
deviations.

###  3.1 Roadmap

We associate to any CRN $G$ a sequence of pure jump Markov processes
$\\{X^{n},n\geqslant 1\\}$ defined on a fixed time interval $[0,T]$. Each
process $\\{X^{n}(t),t\in[0,T]\\}$ is an element of the space of right-
continuous trajectories with left limits known as Skorokhod space, which we
denote $D_{T}\doteq D([0,T]:\mathcal{S})$. It will turn out that under
appropriate conditions, the sequence $\\{X^{n},n\geqslant 1\\}$ satisfies a
Large Deviations Principle on $D_{T}$. Roughly speaking, such a “path-space”
LDP for $\\{X^{n}\\}$ implies that there is a function
$I_{T}:D_{T}\rightarrow\infty$ such that

$\mathbb{P}(X^{n}\in A)\approx\inf_{\varphi\in A}e^{-nI_{T}(\varphi)}.$ (3.1)

for “nice enough” $A\subset D_{T}$. In other words, large deviations provides
a relationship between the _probability_ $\mathbb{P}(X^{n}\in A)$ of an event
$A$ and its _cost_ $\inf_{\varphi\in A}I_{T}(\varphi)$. The more “expensive”
an event, the less likely it is. In the case of CRNs, the function $I_{T}$
will be of the form of an _action functional_ ,

$I_{T}(\varphi)=\int_{0}^{T}L(\varphi(t),\dot{\varphi}(t))\text{dt},$

where $L(x,\beta)$ is real-valued and non-negative, defined for
$(x,\beta)\in\mathbb{R}^{d}_{+}\times\mathbb{R}^{d}$, possibly taking the
value $+\infty$, and $L(x,\cdot)$ is convex for each
$x\in\mathbb{R}^{d}_{>0}$. The function $L$ will be called the _local rate_ or
_Lagrangian_ , and the integral of $L$ over a trajectory will be called the
$I_{T}(\varphi)$ the _cost_ of a trajectory. Under general conditions (which
hold for CRNs), if $I_{T}(\varphi)=0$ then $\varphi$ satisfies the
deterministic dynamics.

It is in general difficult to evaluate $\inf_{\varphi\in A}I_{T}(\varphi)$ for
general $A\subset D_{T}$. A well-studied class of events $A\subset D_{T}$ is
the set of trajectories satisfying _finite time endpoint constraints_ : for
any $x,y\in\mathcal{S}$ and $T>0$, we consider

$A_{T}(x,y)=\\{\varphi\in D_{T}:\varphi(0)=x,\varphi(T)=y\\}.$ (3.2)

The problem of estimating the exponential decay rate of the probability
$\mathbb{P}(X^{n}\in A_{T}(x,y))$ is reduced by the large deviations relation
in (3.1) to solving the _calculus of variations_ problem

$\inf_{\varphi}\left\\{\int_{0}^{T}L(\varphi(t),\dot{\varphi}(t))\text{dt}:\varphi(0)=x,\varphi(T)=y\right\\}.$

A standard way to identify the minimizers of a calculus of variations problem
is to consider the _Euler-Lagrange_ equation

$\frac{d}{dt}\left(\frac{\partial
L}{\partial\beta}(\varphi(t),\dot{\varphi}(t))\right)=\frac{\partial
L}{\partial x},$

with boundary conditions $\varphi(0)=x,\varphi(T)=y$. One can express the
Euler-Lagrange equation in the “Hamiltonian formalism” by introducing the
_Hamiltonian_

$H(x,p)=\sup_{\beta\in\mathbb{R}^{d}}\\{\langle p,\beta\rangle-L(x,\beta)\\}.$
(3.3)

Formally, by introducing the variable $p(x,\beta)=\partial
L(x,\beta)/\partial\beta$ and separating the second order equation into a pair
of two first order equations, we obtain _Hamilton’s equations of motion_

$\displaystyle\dot{x}(t)$ $\displaystyle=$ $\displaystyle D_{p}H(x(t),p(t)),$
(3.4) $\displaystyle\dot{p}(t)$ $\displaystyle=$ $\displaystyle-
D_{x}H(x(t),p(t)),$ (3.5)

with boundary conditions given by $x(0)=x$ and $x(T)=y$. If
$(x(\cdot),p(\cdot))$ is a pair of solutions to the equations of motion
(3.4)-(3.5), then $\\{x(t):t\in[0,T]\\}$ is (under suitable conditions) a
trajectory which achieves the minimum in the variational problem
$\mathcal{J}_{T}(x,y)$. (For further reading, consider p. 255, Section D: _The
Hamilton-Jacobi equation_ in [15].)

Observe that while a solution $(x(\cdot),p(\cdot))$ to the equation of motion
can in principle be used to evaluate the cost
$\inf_{\varphi}I_{T}(A_{T}(x,y))$ by evaluating $I_{T}(x(\cdot))$, it provides
more information than needed for the task of evaluating the decay rate of the
event, as in (3.1). For that, one only needs _the_ _minimizing value_ , not
_a_ _minimizer_.

Consider the idea of parametrizing the cost $\mathcal{J}_{T}(x_{0},x)$ for
some fixed initial position $x_{0}$ in terms of the terminal condition
$x\in\mathbb{R}^{d}_{\geqslant 0}$ and the terminal time $T>0$. We obtain a
function $V(T,x)=\mathcal{J}_{T}(x_{0},x)$ which formally solves the Hamilton-
Jacobi-Bellman (HJB) partial differential equation (PDE)

$\partial_{T}V(T,x)+H(y,D_{x}V(T,x))=0,$ (3.6)

with boundary condition $V(0,x)=0$ if $x=x_{0}$ and $V(0,x)=+\infty$
otherwise. The function $V(T,x)$ represents “the cost to move from $x_{0}$ to
$x$ in time $T$”. Under additional conditions on $x_{0}$ which will be
discussed in Section 4, if one instead solves the time-independent (or _steady
state_) HJB,

$H(x,D_{x}V(x))=0,$

with boundary condition $V(x)=0$ for $x=x_{0}$, then the solution $V(x)$
represents “the cost to move from $x_{0}$ to $x$ in an arbitrarily long but
finite amount of time”.

To summarize,

  1. 1.

the exponential decay rate of the probability of an event with LLN scaling is
governed by the cost of the event;

  2. 2.

evaluating the cost of an event can be seen as a calculus of variations
problem with constraints;

  3. 3.

for the event of moving between one point and another in a fixed amount of
time there is a PDE for the cost;

  4. 4.

the steady-state form of the PDE corresponds to the cost of an event over an
“arbitrarily long but finite” time horizon.

###  3.2 Stochastic dynamics and its Law of Large Lumbers

We now provide the details for the steps outlined at the beginning of the
section, starting with a more precise definition of the stochastic dynamics
and the associated law of large numbers.

Let $G=(\mathcal{A},\mathcal{C},\mathcal{R})$ be a $d$-dimensional CRN with
associated rate constants $\kappa=\\{\kappa_{a,b}:(a,b)\in\mathcal{R}\\}$, and
choose any $n\geqslant 1$. This will be the scaling parameter used for a law
of large numbers. It is intended to represent the volume of the container. For
$x\in n^{-1}\mathbb{Z}^{d}_{+}$, we define

$\lambda^{n}_{a,b}(x)=\frac{\kappa_{a,b}}{n^{|a|}}\prod_{i=1}^{d}\frac{(nx_{i})!}{(nx_{i}-a_{i})!}1\\{nx_{i}\geqslant
a_{i}\\}$

and

$\lambda_{a,b}(x)=\kappa_{a,b}x^{a}.$

Note that if $x\in\mathbb{R}^{d}_{\geqslant 0}$ and $x_{n}\in
n^{-1}\mathbb{Z}^{d}_{+}$, $x_{n}\rightarrow x$, then
$\lambda^{n}_{a,b}(x_{n})\rightarrow\lambda_{a,b}(x)$.

For any fixed $T>0$, we associate to $G$ a sequence of jump Markov processes
$\\{X^{n}(t),t\in[0,T]\\}_{n\geqslant 1}$ with generator

$\mathcal{L}^{n}f(x)=\sum_{(a,b)\in\mathcal{R}}n\lambda^{n}_{a,b}(x)(f(x+\zeta_{a,b}/n)-f(x)).$

We assume that $X^{n}(0)\in n^{-1}\mathbb{Z}^{d}$, and consequently $X^{n}$ is
a $D([0,T]:n^{-1}\mathbb{Z}^{d})$-valued random variable. The family of
generators associated to the CRN in Example 2.1 is

$\mathcal{L}^{n}f(x)=nx1_{\\{nx\geqslant
1\\}}\left(f\left(x-\frac{1}{n}\right)-f(x)\right)+n\left[f\left(x+\frac{1}{n}\right)-f(x)\right].$

It turns out that $X^{n}$ satisfies a Law of Large Numbers (LLN), that is,
there is a function $x(\cdot):[0,T]\rightarrow\mathcal{S}$ such that, for any
$\delta>0$,

$\mathbb{P}\left(\sup_{0\leqslant t\leqslant
T}|X^{n}(t)-x(t)|>\delta\right)\xrightarrow{n\rightarrow\infty}0.$ (3.7)

We refer to the limit in (3.7) as Kurtz’s theorem [3][4]. One can formally
guess the LLN limit $x(\cdot)$ by the heuristic calculation

$\lambda_{a,b}^{n}(x)\left(\frac{f(x+\zeta_{a,b}/n)-f(x)}{1/n}\right)\approx\lambda_{a,b}(x)\langle
Df(x),\zeta_{a,b}\rangle,$

for $f\in C^{1}(\mathbb{R}_{+}).$ It turns out that $x(\cdot)$ is absolutely
continuous and satisfies $\dot{x}(t)=D_{p}H(x(t),0)$, where

$D_{p}H(x(t),0)\doteq\sum_{(a,b)\in\mathcal{R}}\lambda_{a,b}(x(t))\zeta_{a,b}.$
(3.8)

We call (3.8) the _deterministic dynamics_ of a CRN. For the network in
Example 2.1, we find

$D_{p}H(x,0)=1-x.$

###  3.3 Large devations: LLN’s rate of convergence

Whenever one has a sequence of random variable satisfying a law of large
numbers, one can ask whether the sequence satisfies a large deviations
principle. Roughly speaking, a large deviations principle identifies the _rate
of decay_ in $n$ of the probabilities (3.7).

In this section we define rate functions and a large deviations principle in
abstract terms. We let $\mathcal{X}$ denote a complete, seperable metric
space, also known as a Polish space.

######  Definition 3.1.

A function $I:\mathcal{X}\rightarrow[0,\infty)$ is a _rate function_ if for
each $M<\infty$, the sublevel set $\\{x\in\mathcal{X}:I(x)\leqslant M\\}$ is
compact.

We recall that a having _closed_ sublevel sets is equivalent to being lower
semicontinuous, and that any lower semicontinuous function achieves its
minimum on a compact set. It follows that a rate function achieves its minimum
on any sublevel set.

######  Definition 3.2.

A sequence $Y^{n}$ satisfies a _Large Deviations Principle (LDP)_ with rate
function $I$ and speeed $n^{-1}$ if, for any open set $E\subset\mathcal{X}$,
we have the large deviations lower bound

$\liminf_{n\rightarrow\infty}\frac{1}{n}\log\mathbb{P}(Y^{n}\in
E)\geqslant-\inf_{y\in E}I(y)$

and for any closed set $F\subset\mathcal{X}$, we have the large deviations
upper bound

$\limsup_{n\rightarrow\infty}\frac{1}{n}\log\mathbb{P}(Y^{n}\in
F)\leqslant-\inf_{y\in F}I(y).$

We define, for any set $A\subset\mathcal{X}$,

$I(A)\doteq\inf_{y\in A}I(y),$

and we will say that an event is _rare_ if $I(A)>0$.

####  3.3.1 Large deviations for CRNs

There are at least two obstacles when trying to establish a large deviations
principle for chemical reaction networks. First, one must place appropriate
assumptions on the jump intensities to ensure that they do not diverge “too
quickly” and that there is an appropriate controllability condition; see [16]
for such a result. Another important issue is boundary behavior: one needs to
account for the possibility of species vanishing. A proof of the LDP which
accounts for this kind of behavior has recently been proposed by [17].

The result is likely to hold on the space $D([0,T])$ with the usual Skorokhod
topology under milder conditions than what has been established as of yet. If
the result does hold, then one can be confident the rate function will be

$I_{T}(\varphi)=\left\\{\begin{array}[]{ll}\int_{0}^{T}L(\varphi(s),\dot{\varphi}(s))\text{
ds }&\text{ if }\varphi\text{ is }\text{AC}\\\
+\infty&\text{else}\end{array}\right.,$

where for $(x,\beta)\in\mathbb{R}^{d}_{+}\times\mathbb{R}^{d}$,

$L(x,\beta)=\inf_{u_{a,b}}\left\\{\sum_{(a,b)\in\mathcal{R}}\lambda_{a,b}(x)\ell\left(\frac{u_{a,b}}{\lambda_{a,b}(x)}\right):u_{a,b}\geqslant
0\text{ for
}(a,b)\in\mathcal{R},\beta=\sum_{(a,b)\in\mathcal{R}}\zeta_{a,b}u_{a,b}\right\\},$

and

$\ell(x)\doteq x\log x-x+1,$

with $\ell(0)=0$ by convention. We call $I_{T}(\varphi)$ the _cost_ of the
trajectory $\varphi$. An important regularity property of $L$ is that for any
$x\in\mathcal{S}$, there is a unique $b(x)$ such that $L(x,b(x))=0$, and
moreover, the function $b(x)$ is Lipschitz continuous. Indeed, $L(x,\beta)=0$
if and only if $u_{a,b}=\lambda_{a,b}(x)$ for all $(a,b)\in\mathcal{R}$, or
equivalently
$\beta=\sum_{(a,b)\in\mathcal{R}}\zeta_{a,b}\lambda_{a,b}(x)=D_{p}H(x,0)$.

Moreover, since $L$ is nonnegative, it must be that any $\varphi$ which
satisfies $I_{T}(\varphi)=0$ satisfies the deterministic dynamics, which are
sometimes referred to as the _zero-cost_ dynamics. We remark that the rate
function $I_{T}(\cdot)$ depends on the time interval $[0,T]$ on which the
trajectories are defined.

####  3.3.2 Heuristics for large deviations

Rigorously justifying the form of the large deviations rate function is beyond
the scope of this paper, and we refer the reader interested in the details to
an introductory textbook such as [18]. In the interest in providing some
motivation for the form of $L$ and $H$, we provide two examples with heuristic
calculations.

######  Example 3.1.

Here is a basic calculation that motivates the definition of the Hamiltonian
and the Legendre transform. Suppose $\\{\xi_{i}\\}_{i=1}^{\infty}$ are i.i.d.
random variables with “light tails”. Then, for $\beta>0$ and any $p\geqslant
0$, we have by Chebyshev’s inequality that

$\mathbb{P}\left(\frac{1}{n}\sum_{i=1}^{n}\xi_{i}\geqslant\beta\right)=\mathbb{P}\left(e^{p\sum_{i}\xi_{i}}\geqslant
e^{pn\beta}\right)\leqslant\frac{\mathbb{E}\left[e^{p\sum_{i}\xi_{i}}\right]}{e^{pn\beta}}=e^{-n(p\beta-H(p))},$

where

$H(p)=\log\mathbb{E}[e^{p\xi_{1}}].$

We can optimize over the parameter $p$ by taking the infimum on the right-hand
side with the inequality remaining valid. It can be shown under appropriate
conditions that, in fact,

$\lim_{n\rightarrow\infty}\frac{1}{n}\log\mathbb{P}\left(\frac{1}{n}\sum_{i=1}^{n}\xi_{i}\geqslant\beta\right)=L(\beta)$

where $L(\beta)=\sup_{p\in\mathbb{R}}\\{p\beta-H(p)\\}.$ This is a result
covered by what is known as Cramér’s theorem; see [19], p. 26.

The preceding example allows us to guess the rate function on path space.

######  Example 3.2.

Consider a stochastic process of the form

$Y^{n}\left(\frac{i+1}{n}\right)=Y^{n}\left(\frac{i}{n}\right)+\frac{1}{n}\xi_{i},$
(3.9)

where the $\xi_{i}$ are i.i.d. random variables with “light tails”. The
scaling in (3.9) is a LLN-type scaling, and it is reasonable to expect that

$Y^{n}(t)\xrightarrow{n\rightarrow\infty}t\mathbb{E}X_{1},\quad t\in[0,1].$

The key here is a _time-scale separation_ due to the $1/n$ scaling. Over a
time interval $[t,t+\delta]$, $Y^{n}(t)$ will change at most $O(\delta)$ with
high probability, since the increments $\xi_{i}$ are light-tailed (say, sub-
Gaussian). For $\delta>0$,

$\frac{Y^{n}(t+\delta)-Y^{n}(t)}{\delta}\approx\frac{1}{n\delta}\sum_{i=\lfloor
nt\rfloor}^{\lfloor nt+n\delta\rfloor}\xi_{i}.$ (3.10)

For large $n$ and $\delta>0$ small, we can think of the left-hand side of
(3.10) as a sum of $\lfloor n\delta\rfloor$ i.i.d. random variables. By
Cramér’s theorem,

$\mathbb{P}\left(\frac{Y^{n}(t+\delta)-Y^{n}(t)}{\delta}\approx\beta\right)\approx
e^{-n\delta L(\beta)},$ (3.11)

where $L$ is as in Example 3.1. The random variables in this example are
trajectories, and $\beta$ can be interpreted as the _velocity_ of $Y^{n}$.
Heuristically, the probability of $Y^{n}$ being close to an (absolutely
continuous) path $\\{\varphi(t),t\in[0,1]\\}$ is, for small $\delta>0$ an,
roughly

$\displaystyle\mathbb{P}(Y^{n}\approx\varphi)$ $\displaystyle=$
$\displaystyle\mathbb{P}(\|Y^{n}-\varphi\|_{\infty}\leqslant\varepsilon)$
$\displaystyle\approx$
$\displaystyle\mathbb{P}(Y^{n}(j\delta)\approx\varphi(j\delta),j=0,\ldots,\lfloor
1/\delta\rfloor)$ $\displaystyle\approx$
$\displaystyle\mathbb{P}\left(\frac{Y^{n}(j\delta+\delta)-Y^{n}(j\delta)}{\delta}\approx\frac{\varphi(j\delta+\delta)-\varphi(j\delta)}{\delta},j=0,\ldots,\lfloor
1/\delta\rfloor\right)$ $\displaystyle\approx$
$\displaystyle\prod_{j=0}^{\lfloor 1/\delta\rfloor}e^{-n\delta
L\left(\frac{\varphi(j\delta+\delta)-\varphi(j\delta)}{\delta}\right)}$
$\displaystyle\approx$ $\displaystyle
e^{-n\int_{0}^{1}L(\dot{\varphi}(s))\text{ds}}.$

In the above we used (3.11) with (essentially) $\beta=\dot{\varphi}$, and we
crucially used independence of the increments to break up the intersection of
events into a product. This example shows one way of guessing the form of the
rate function for processes.

######  Remark 3.3.

If the distribution of the increment $\xi_{n}$ in the preceding example is
taken to be dependent on the previous state $Y_{n-1}$, then under appropriate
assumptions the local rate $L$ will be state-dependent as well, and
$L(x,\beta)$ can be defined for fixed $x$ by the duality formula (3.3),

$L(x,\beta)=\sup_{p\in\mathbb{R}^{d}}\\{\langle p,\beta\rangle-H(x,p)\\}.$
(3.12)

###  3.4 From large deviations to Hamilton-Jacobi-Bellman

While a large deviations principle allows one to identify the asymptotic decay
rate of any sufficiently nice event $A\subset\mathcal{X}$, it is in general
difficult to solve the optimization problem $\inf_{\varphi\in A}I_{T}(A)$. In
this section we elaborate on the details of taking a PDE perspective on the
optimization problem for a particular class of events. The underlying PDE is
called the Hamilton-Jacobi-Bellman (HJB) equation, and is sometimes referred
to as the _dynamic programming_ equation, especially in discrete space and
time. The techniques in continuous space and time trace back to Lagrangian and
Hamiltonian classical mechanics, as well as optimal control and the calculus
of variations.

We first outline the main ideas. As in Section 3.1, consider the family of
events parametrized by the initial and terminal conditions,

$A_{T}(x,y)=\\{\varphi\in D([0,T]):\varphi(0)=x,\varphi(T)=y\\}.$

We let

$\mathcal{J}_{T}(x,y)=\inf_{\varphi\in A_{T}(x,y)}I_{T}(\varphi).$ (3.13)

Thus, the probability of going from $x\in\mathbb{R}^{d}$ to
$y\in\mathbb{R}^{d}$ in time $T>0$ is roughly

$\mathbb{P}_{x}(X^{n}(T)\in B_{\delta}(y))\approx e^{-n\mathcal{J}_{T}(x,y)}.$

The problem of evaluating $\mathcal{J}_{T}(x,y)$ is a finite time problem with
endpoint constraints. Roughly speaking, $\mathcal{J}_{T}(x,y)$ measures the
cost of the stochastic process $X^{n}$ moving from $x$ to $y$ in time $T$. The
most likely trajectory for $X^{n}$ to follow starting at $X^{n}(0)=x$ is the
law of large numbers dynamics, $\dot{x}(t)=D_{p}H(x(t),0)$ with $x(0)=x$, and
in this case we expect $X^{n}(T)\approx x(T)$. The cost of this trajectory is
identically zero: $\mathcal{J}_{T}(x(0),x(T))=0$.

We now sketch the connection of $\mathcal{J}_{T}(x,y)$ to the HJB PDE. For
some fixed $y\in\mathbb{R}^{d}_{>0}$ and all $x\in\mathcal{S}_{y}$,
$t\in[0,T]$, define

$\displaystyle V(t,x)$
$\displaystyle=\inf_{\varphi}\left\\{\int_{t}^{T}L(\varphi(s),\dot{\varphi}(s))\text{ds}:\varphi(t)=x,\varphi(T)=y\right\\}$
(3.14) $\displaystyle=\mathcal{J}_{T-t}(x,y),$

and define $V(t,x)=+\infty$ if the set over which the infimum is taken is
empty.

The connection to the PDE is made by observing that if $\varphi$ is an optimal
trajectory starting at $(t,x)$ and ending at $(T,y)$, and $\varphi$ passes
through $(t+\delta,\varphi(t+\delta))$, then $\varphi$ should also be optimal
to move from $(t+\delta,\varphi(t+\delta))$ to $(T,y)$. This is known as the
“principle of optimality”. In other words, we expect

$\displaystyle V(t,x)$
$\displaystyle=\min_{\varphi}\bigg{\\{}V(t+\delta,\varphi(t+\delta))$
$\displaystyle\quad+\int_{t}^{t+\delta}L(\varphi(s),\dot{\varphi}(s))ds:\varphi(t)=x,\varphi(T)=y\bigg{\\}}.$

By subtracting $V(t,x)$ from both sides, dividing by $\delta$, sending
$\delta\rightarrow 0$ and using the Fenchel-Legendre duality (3.12) with
$\dot{\varphi}(t)$ in place of $\beta$ and $-D_{x}V(t,x)$ in place of $p$, we
find that $V$ formally satisfies the HJB

$\partial_{t}V(t,x)-H(x,-D_{x}V(t,x))=0$ (3.15)

with boundary condition $V(T,x)=0$ if $x=y$ and $V(T,x)=+\infty$ otherwise.
Note that the PDE (3.15) differs by two negative signs from the equation in
(3.6), due to the fact that in the present case the terminal time and position
are fixed, while in the other case the terminal time and position are the
variables. The two choices represent fundamentally different, but related
problems.

A precise statement relating the variational definition of $V$ in (3.14) to
the solution of (3.15) is difficult since the very notion of solution for HJB
equations is itself quite subtle. For instance, the function $V$ need not be
differentiable. Crandall and Lions established the notion of _viscosity_
solutions to resolve this issue in their seminal work [10]. The theory is
technical and outside the scope (or purpose) of this note, but see: [12],
[11], Theorem 5, p. 128 in [20] for a simplified situation, or [21] for more
general extensions in the “infinite horizon” case (which will be discussed
below). We will primarily use the results in [2] Chapters 3-5.

The HJB PDE provides another approach to establishing a LDP (this is outlined
in [12]). It is also useful for the design of rare event simulation schemes
such as importance sampling or splitting (see [22, 23, 24, 25]). Certain kinds
of rare events, such as situations where a stochastic process “jumps” between
two metastable states, may occur over very long timescales. In fact, one can
prove that under appropriate assumptions on a domain $D$ containing a unique
stable steady state, and an appropriately scaled stochastic process $X^{n}$
satisfying an appropriate “uniform” LDP, we have

$\lim_{n\rightarrow\infty}-\frac{1}{n}\log\mathbb{P}_{x}(X^{n}(t)\notin
D,\text{ for some }t\leq n^{k})=V(x),\quad x\in D,$

for any $k\geq 1$, where $V(x)$ is given by the variational problem

$V(x)=\inf_{T<\infty}\inf_{y\notin D}\mathcal{J}_{T}(x,y).$

See [26] for a proof and an application to splitting methods. The function
$V(x)$ formally satisfies the _steady-state_ HJB equation,

$H(x,-D_{x}V(x))=0,$

with boundary condition $V(x)=0$ for $x\in\partial D$. The variational
formulation $V(x)=\inf_{T<\infty}\inf_{y\notin D}\mathcal{J}_{T}(x,y)$ allows
us to interpret $V(x)$ as being the cost to exit $D$ starting from $x\in D$.

There is a related problem which considers the cost of hitting a point $x\in
D$ starting from the steady state in an arbitrarily long but finite amount of
time. Given an appropriate domain $D\subset\mathbb{R}^{d}$ with a unique
stable steady state $c\in D$, the solution $V(x)$ of

$H(x,DQ(x))=0$

with $Q(c)=0$ identifies the cost of moving from $c$ to
$x\in\\{y:Q(y)\leq\min_{z\in\partial D}Q(z)\\}$; see Theorem 4.3 in Chapter 5
of [2]. The solution $Q$ is a classical object known as the Freidlin-Wentzell
_quasipotential_ and is the central object of study in this note; see the next
section.

We end the section with a conjecture that is related to the ideas discussed in
this section but has no relation to the rest of the note.

##### A conjecture.

This conjecture is related to necessary conditions for a trajectory to be a
minimizer of the optimization problem $\mathcal{J}_{T}(x,y)$. claim that a
“Lagrange multiplier” approach to the optimization problem
$\mathcal{J}_{T}(x,y)$. The infinite-dimensional analogue of the Karush-Kuhn-
Tucker conditions which provide necessary conditions for Lagrange multipliers
is the celebrated _Pontryagin Maximum Principle_ ; see [27], also [28] for the
same result under weaker conditions.

We conjecture that if the pair $(x(t),p(t))$ satisfies

$\displaystyle\dot{x}(t)$ $\displaystyle=$
$\displaystyle\sum_{(a,b)\in\mathcal{R}}\zeta_{a,b}\lambda_{a,b}(x(t))e^{\langle
p(t),\zeta_{a,b}\rangle}$ $\displaystyle\dot{p}(t)$ $\displaystyle=$
$\displaystyle\sum_{(a,b)\in\mathcal{R}}D_{x}\lambda_{a,b}(x(t))(1-e^{\langle
p(t),\zeta_{a,b}\rangle}),$

with boundary conditions $x(0)=x$, $x(T)=y$, then
$I_{T}(x(\cdot))=\mathcal{J}_{T}(x,y)$. One cannot apply even the weakest
forms of the Pontryagin Maximum Principle since the underlying class of
controls does not satisfy the integrability assumptions. This lack of
integrability is related to the blowup of $L(x,\beta)$ as
$x\rightarrow\partial\mathcal{S}$ for certain directions $\beta$. We also note
that for $x\in\partial S$ the spatial gradient must be given an appropriate
interpretation.

##  4 Quasipotential, Lyapunov functions, and HJB

We consider a family of variational problems which define an object known as
the _quasipotential_. The quasipotential was first introduced by Freidlin and
Wentzell and is often used to study the long-term dynamics of stochastic
systems.

There are three important results we wish to communicate: the relationship of
the quasipotential to Lyapunov functions, invariant measures, and the
Hamilton-Jacobi-Bellman equation. In Section 4.1, we define the
quasipotential, as well as the notion of a Lyapunov function, and we show in
Section 4.2 that the quasipotential serves as a Lyapunov function for the
noiseless dynamics. In Section 4.3 we demonstrate how the quasipotential can
be deduced by an appropiately scaled limit of invariant measures for the
process $X^{n}$, and in Section 4.4 we show how it can instead by obtained
from the HJB PDE.

###  4.1 Definitions

Recall that $\mathcal{J}_{T}(x,y)$ denotes the optimal value of a variational
problem, as defined in (3.13).

######  Definition 4.1.

The quasipotential between two points $x,y\in\mathbb{R}^{d}$ is defined as

$Q(x,y)=\inf_{T<\infty}\mathcal{J}_{T}(x,y),$

with $Q(x,y)=+\infty$ if $y\notin\mathcal{S}_{x}$ (recall Definition 2.2).

Thus the quasipotential in some sense measures the “optimal cost” of moving
from $x$ to $y$ when allowing for an arbitrarily long but finite amount of
time. Note that the infimum need not be achieved, or if it is achieved the
minimizer need not be unique.

Next, we define more precisely a term that has been used throughout the note.

######  Definition 4.2 (Stable steady state).

A steady state $c$ of the dynamical system $\dot{x}(t)=\Phi(x(t))$ is said to
be _stable_ if for every neighbourhood $U$ of $c$ there is a smaller
neighbourhood $V\subset U$ such that the trajectories of
$\dot{x}(t)=\Phi(x(t))$ starting in $V$ converge to $c$ without leaving $U$.
We say that a domain $D$ containing $c$ is attracted to $c$ if the
trajectories $\dot{x}(t)=\Phi(x(t))$ converge to $c$ without leaving $D$.

We will assume throughout without further mention that $D$ will denote a
domain that is attracted to a stable steady state $c$ under the deterministic
dynamics $\Phi(x)=D_{p}H(x,0)$, and we refer to the maximal open set of points
that are attracted to $c$ as the _domain of attraction_.

The study of stability is aided by the notion of a Lyapunov function.

######  Definition 4.3 (Lyapunov function).

Let $D\subset\mathbb{R}^{d}_{+}$ be an open subset of $\mathbb{R}^{d}_{+}$ and
let $\Phi:\mathbb{R}^{d}_{+}\rightarrow\mathbb{R}^{d}$. Suppose that
$c\in\mathbb{R}^{d}_{+}$ is a stable steady state for the system

$\dot{x}(t)=\Phi(x(t)).$

A continuously differentiable function $U:D\rightarrow\mathbb{R}$ is called a
_strict Lyapunov function_ for the system $\dot{x}(t)=\Phi(x(t))$ at $c$ if
$U(c)=0$, $U(x)>0$ for all $x\neq c$, and

$\langle DU(x),\Phi(x)\rangle\leqslant 0$

for all $x\in D$, with equality if and only if $x=c$.

###  4.2 Quasipotential and Lyapunov functions

In this section we show that the quasipotential serves as a Lyapunov function
for the deterministic dynamics. This can be understood heuristically as
follows. The quasipotential $Q(x)$ represents the minimal “cost” to deviate
from a stable steady state $c$ over arbitrarily long but finite time
intervals. The deterministic dynamics have $c$ as a stable steady state, so
any deviation from $c$ necessarily carries a positive cost. The cost increases
the “further away” one deviates, so the cost is a natural quantity which
increases as one moves away from $c$, which is precisely a Lyapunov function
type property.

Though this property is implicit in the work of Freidlin and Wentzell, we
provide a proof for completeness.

######  Proposition 4.4.

Let $c$ be a stable steady state for the system $\dot{x}(t)=D_{p}H(x(t),0)$.
The quasipotential $Q(x)\doteq Q(c,x)$ is a Lyapunov function for this system
in the sense of Definition 4.3 on the domain of attraction $D$ of $c$.

###### Proof.

Let $x\in D$. We first show that $Q(x)\geqslant 0$ with equality if and only
if $x=c$. Since $c$ is assumed to be a steady state, we have $D_{p}H(c,0)=0$.
By properties of the Legendre transform, $L(c,0)=0$, so any trajectory which
stays at $c$ for a positive amount of time achieves the minimum in the
variational definition of the quasipotential and incurs zero cost. Thus
$Q(c)=0$.

Next, for each $x\in\mathbb{R}^{d}$ the unique $b(x)\in\mathbb{R}^{d}$ such
that $L(x,b(x))=0$ is $b(x)=D_{p}H(x,0)$. Since $c$ is a stable steady state,
it can be shown that any absolutely continuous trajectory $\varphi$ which
leads from $c$ to $x$ in time $T$ must necessarily satisfy
$\dot{\varphi}(t)\neq D_{p}H(\varphi(t),0)$ for all $t$ in some interval. In
particular, $\varphi$ must incur a positive cost over some time interval,
hence $Q(x)>0$ for all $x\neq c$.

It remains to show that $Q(x(t+h))-Q(x(t))\leqslant 0$ for all $t\geqslant 0$
and $h>0$, with equality if and only if $x(t)=c$. Let $\varepsilon>0$ be
arbitrary and $T<\infty$, and let $\varphi\in D([0,T])$ be such that
$\varphi(0)=c$, $\varphi(T)=x(t)$, and $I_{T}(\varphi)<Q(x(t))+\varepsilon$.
Define a new trajectory $\psi\in D([0,T+h])$ by concatenating it with
$\varphi$ on $[0,T]$ and with the trajectory $\\{x(s):s\in[t,t+h]\\}$ on
$(T,T+h]$. Since $x(\cdot)$ follows the noiseless dynamics, it has zero cost
and $I_{T}(\psi)=I_{T}(\varphi)$. Since $\psi(0)=c$ and $\psi(T+h)=x(t+h)$, it
follows that $Q(x(t+h))\leqslant Q(x(t))+\varepsilon$. Since $\varepsilon>0$
was arbitrary, $Q(x(t+h))\leqslant Q(x(t))$. See Figure 1 for a mental
picture.

Figure 1: The trajectory $\psi(\cdot)$ defined on $[0,T+h]$ is the
concatenation of $\varphi$ on $[0,T]$ with $x(\cdot)$ on $(T,T+h]$.

To rule out the case of equality, note that if there are $t\geqslant 0$ and
$h>0$ such that $Q(x(t+h))=Q(x(t))$, then by the mean value theorem there is
$\xi\in(t,t+h)$ such that $\langle DQ(x(\xi)),\dot{x}(\xi)\rangle h=0$. Since
$c$ is the unique stable steady state in $D$, $\dot{x}(\xi)\neq 0$ and
$DQ(x(\xi))\neq 0$ unless $x(\xi)=c$. Thus equality cannot hold unless
$x(\xi)=c$. This completes the proof. ∎

To illustrate the role of the quasipotential as a Lyapunov function, consider
the following example.

######  Example 4.1.

Consider the reaction network $A\rightleftarrows\emptyset$ from Example 2.1.
We will see in Examples 4.2 and 4.3 that the quasipotential for this CRN is
$Q(x)=x\log x-x+1$. This is a Lyapunov function for the deterministic dynamics
$\dot{x}(t)=1-x(t)$ relative to the stable steady state at $c=1$. Since $c$ is
globally attracting, $Q$ is a Lyapunov function on all of $\mathbb{R}_{+}$.
Indeed, $Q(1)=0$, $Q(x)>0$ for $x\neq 1$, and

$\displaystyle\frac{d}{dt}Q(x(t))$ $\displaystyle=$
$\displaystyle\dot{x}(t)\log x(t)$ $\displaystyle=$ $\displaystyle(1-x(t))\log
x(t).$

If $x(t)\in[0,1)$ then $\log x(t)<0$, whereas if $x(t)\in(1,\infty)$ we have
$(1-x(t))<0$. In either case, the above expression is negative, with equality
if and only if $x(t)=1$.

###  4.3 Quasipotential and invariant measure

The following is a classical result of Freidlin-Wentzell theory; see for
instance Theorem 4.3 in Chapter 4 of [2] for the case of diffusions. The
modification of the proof for jump processes is relatively straightforward, as
the estimates depend on “uniform” properties of the rate function; see p. 74
of [2]. This relationship between the invariant distribution and the
quasipotential is what was used to determine the Lyapunov functions in [1].

######  Theorem 4.5.

Let $c$ be the unique stable steady state of the system
$\dot{x}(t)=D_{p}H(x(t),0)$, and suppose that its domain of attraction is all
of $\mathbb{R}^{d}$. Assume that for each $n\geqslant 1$, $X^{n}$ has a unique
invariant measure $\pi^{n}$. Then

$\lim_{n\rightarrow\infty}\frac{1}{n}\log\pi^{n}(A)=-\inf_{x\in A}Q(c,x),$

where $A\subset\mathbb{R}^{d}$ is any domain with compact boundary $\partial
A$ common for both $A$ and the closure of $A$.

######  Example 4.2.

The sequence of processes $\\{X^{n},n\geqslant 1\\}$ associated to the
reaction network $A\rightleftarrows\emptyset$ are birth-death process with
state-dependent jump rates as defined by the generator in Example 3.2, and the
density of its unique stationary distribution is given by

$\pi^{n}(x)=e^{-n}\frac{n^{nx}}{(nx)!},\quad x\in n^{-1}\mathbb{Z}^{d}_{+}.$

Applying the large deviations scaling and using Stirling’s formula, $\log
k!=k\log k-k+O(\log k)$ for $k\geqslant 1$, we find

$\lim_{n\rightarrow\infty}-\frac{1}{n}\log\pi^{n}(x)=x\log x-x+1.$

We note that Example 4.2 illustrates a technique which was used to determine
Lyapunov functions for the deterministic dynamics in [1].

###  4.4 Quasipotential and HJB

The quasipotential is related to the HJB in the following way.

######  Theorem 4.6.

Let $c$ be a stable steady state for the noiseless dynamics
$\dot{x}(t)=D_{p}H(x(t),0)$, and let $D$ denote the domain of attraction of
$c$. Let $U:D\rightarrow[0,\infty)$ be continuously differentiable on $D$ and
continuous on its closure $\bar{D}$, with $U(x)>0$ and $DU(x)\neq 0$ for
$x\neq c$, and $U(c)=0,DU(c)=0$. Moreover, suppose that for all $x\in D$,

$H(x,DU(x))=0.$

Then $U(x)=Q(c,x)$ for all $x\in\\{z:U(z)\leqslant\min_{y\in\partial
D}U(y)\\}$.

The proof of this proposition can be found in Theorem 4.3, Chapter 5 of [2].
Since the HJB PDE allows us to identify the quasipotential and the
quasipotential is always a Lyapunov function, the HJB can be used to identify
Lyapunov functions for the noiseless dynamics. The example below shows how the
HJB can be used to identify the Lyapunov function $x\log x-x+1$ introduced in
Example 4.1.

######  Example 4.3.

Consider the time-independent HJB for the CRN $A\rightleftarrows\emptyset$,

$H(x,DU(x))=x(e^{-DU(x)}-1)+e^{DU(x)}-1=0.$

One can check by direct calculation that $U(x)=x\log x-x+1$ satisfies
$H(x,DU(x))=0$ for all $x\in\mathbb{R}_{+}$. Moreover, $U(x)\geqslant 0$ with
equality if and only if $x=1$, and the same is true for $DU(x)=\log x$. $U$ is
also smooth and the domain of attraction is all of $\mathbb{R}_{+}$, so $U$
satisfies the conditions of Theorem 4.6 for all $x\in\mathbb{R}_{+}$ and is
equal to the quasipotential everywhere.

##  5 Quasipotential for complex-balanced systems

In this section we consider the quasipotential for the class of “complex-
balanced” reaction networks. We will show that the time-independent Hamilton-
Jacobi-Bellman equation can be used to characterize such networks, as it will
admit a solution of a specific form if and only if the network is complex-
balanced.

###  5.1 Complex-balanced networks

We first define a structural property of reaction networks.

######  Definition 5.1.

A CRN $G=(\mathcal{A},\mathcal{C},\mathcal{R})$ is said to be _weakly
reversible_ if for any $(a,b)\in\mathcal{R}$, there is a sequence of reactions
$(b_{0},b_{1}),(b_{1},b_{2}),\ldots,(b_{k-1},b_{k})\in\mathcal{R}$ such that
$b_{0}=b$ and $b_{k}=a$.

We now define complex-balanced networks. The notion was introduced in [5] and
has been extensively studied since [29, 30, 31, 6, 1].

######  Definition 5.2.

A CRN $G=(\mathcal{A},\mathcal{C},\mathcal{R})$ is said to be complex-balanced
for a choice of rate constants $\\{\kappa_{a,b}:(a,b)\in\mathcal{R}\\}$ if,
for any $z\in\mathcal{C}$,

$\sum_{a:(a,z)\in\mathcal{R}}\kappa_{a,z}c^{a}=\sum_{b:(z,b)\in\mathcal{R}}\kappa_{z,b}c^{z}.$
(5.1)

In other words, complex-balanced systems are ones which admit a steady state
where the net flow in each _complex_ is zero. This is a stronger requirement
than simply being a steady state, which requires only that the net flow in
each _species_ is zero.

There are two important consequences of a network being complex-balancing. The
first is that a complex-balanced network must be weakly reversible in the
sense of Definition 5.1. The second is that complex-balanced networks admit a
unique steady state in each stoichiometric compatibility class, and that this
steady state is stable. The stability of the steady state is typically
established using the well-known Lyapunov function

$V(c,x)=\sum_{i=1}^{d}c_{i}\ell(x_{i}/c_{i}),$ (5.2)

where $\ell(x)=x\log x-x+1$. We will show in the next section that a stronger
characterization of $V$ is that it is a solution to the HJB PDE corresponding
to a complex-balanced reaction network, which by Proposition 4.4 automatically
guarantees it is a Lyapunov function.

###  5.2 Main theorems

We first show that (5.2) is a solution to the HJB if and only if the network
is complex-balanced. We note that this has already been observed in [32],
Theorem 11. There, the approach taken is the typical approximation of the
Fokker-Planck equation with a LLN scaling111This is typically referred to as a
“WKB approximation”: one passes from the Schrödinger equation to the HJB
equation in the “small noise”, i.e. classical limit.

######  Theorem 5.3.

Let $G=(\mathcal{A},\mathcal{C},\mathcal{R})$ be a reaction network with rate
constants $\kappa=\\{\kappa_{a,b}>0:(a,b)\in\mathcal{R}\\}$, and let $c$ be a
steady state in any stoichiometric compatibility class. Then $V(c,x)$ as
defined in (5.2) satisfies $H(x,D_{x}V(c,x))=0$ for all $x\in\mathcal{S}$ if
and only if the state $c$ is complex-balancing.

###### Proof.

The proof is a straightforward calculation. Suppose first that
$c\in\mathbb{R}^{d}_{>0}$ satisfies, for each $z\in\mathcal{C}$, the complex-
balanced relation

$\sum_{a:(a,z)\in\mathcal{R}}\kappa_{a,z}c^{a}=\sum_{b:(z,b)\in\mathcal{R}}\kappa_{z,b}c^{z}.$

We verify that $V$ satisfies the steady-state HJB by direct calculation.

$\displaystyle H(x,DV(x))$ $\displaystyle=$
$\displaystyle\sum_{(a,b)\in\mathcal{R}}\kappa_{a,b}x^{a}(e^{\langle\zeta_{a,b},DV(x)\rangle}-1)$
$\displaystyle=$
$\displaystyle\sum_{(a,b)\in\mathcal{R}}\kappa_{a,b}x^{a}\left(\frac{\left(\frac{x}{c}\right)^{b}}{\left(\frac{x}{c}\right)^{a}}-1\right)$
$\displaystyle=$
$\displaystyle\sum_{(a,b)\in\mathcal{R}}\kappa_{a,b}x^{b}\frac{c^{a}}{c^{b}}-\sum_{(a,b)\in\mathcal{R}}\kappa_{a,b}x^{a}$
$\displaystyle=$
$\displaystyle\sum_{(a,b)\in\mathcal{R}}\left(\left[\sum_{z\in\mathcal{C}}1_{\\{b=z\\}}\right]\kappa_{a,b}x^{b}\frac{c^{a}}{c^{b}}\right)-\sum_{(a,b)\in\mathcal{R}}\left(\left[\sum_{z\in\mathcal{C}}1_{\\{a=z\\}}\right]\right)\kappa_{a,b}x^{a}$
$\displaystyle=$
$\displaystyle\sum_{z\in\mathcal{C}}\left[\sum_{a:(a,z)\in\mathcal{R}}\kappa_{a,z}x^{z}\frac{c^{a}}{c^{z}}-\sum_{b:(z,b)\in\mathcal{R}}\kappa_{z,b}x^{z}\right]$
$\displaystyle=$
$\displaystyle\sum_{z\in\mathcal{C}}\left[\frac{x^{z}}{c^{z}}\sum_{a:(a,z)\in\mathcal{R}}\kappa_{a,z}c^{a}-\sum_{b:(z,b)\in\mathcal{R}}\kappa_{z,b}x^{z}\right]$
$\displaystyle=$
$\displaystyle\sum_{z\in\mathcal{C}}\left[\frac{x^{z}}{c^{z}}\sum_{b:(z,b)\in\mathcal{R}}\kappa_{z,b}c^{z}-\sum_{b:(z,b)\in\mathcal{R}}\kappa_{z,b}x^{z}\right]$
$\displaystyle=$ $\displaystyle 0.$

In the penultimate equality we used the complex-balanced relation.

For the reverse direction, observe that until the penultimate equality we did
not use the complex-balanced relation. The relation $H(x,DV(x))=0$ for all
$x\in\mathcal{S}$ implies

$\sum_{z\in\mathcal{C}}\left[x^{z}\left(\sum_{a:(a,z)\in\mathcal{R}}\kappa_{a,z}\frac{c^{a}}{c^{z}}-\sum_{b:(z,b)\in\mathcal{R}}\kappa_{z,b}\right)\right]=0$

for all $x\in\mathcal{S}$. Since each $z\in\mathcal{C}$ corresponds to a
different monomial $x^{z}$, for the above expression to be identically zero on
all of $\mathcal{S}$, we must have each term in the sum over $z$ being zero.
Thus

$\sum_{a:(a,z)\in\mathcal{R}}\kappa_{a,z}\frac{c^{a}}{c^{z}}-\sum_{b:(z,b)\in\mathcal{R}}\kappa_{z,b}=0,$

which, after multiplying by $c^{z}$, is the complex balanced relation (5.1). ∎

######  Corollary 5.4.

The network is complex-balanced if and only if $Q(c,x)=V(c,x)$ on a
neighbourhood of $c$.

###### Sketch of proof..

Suppose that $Q(c,x)=V(c,x)$ for all $x$ in some neighbourhood $N$ of $c$. By
restricting to a smaller neighbourhood if necessary, we may assume that $N$ is
contained in the domain of attraction of $c$. Since $Q=V$ is the
quasipotential, it solves the steady-state HJB $H(x,DV(x))=0$ for $x\in N$
(see e.g. [21]). By Theorem 5.3, $c$ satisfies the complex-balanced relations,
and so the network must be complex-balanced.

Conversely, suppose that the network is complex-balanced. Then by Theorem 5.3,
$V$ solves the HJB PDE for all $x\in\mathcal{S}_{c}$. Moreover, $V(c,x)$ is
non-negative, and $V(c,x)$ and $DV(c,x)$ equal respectively the scalar and
vector 0 if and only if $x=c$. Thus $V$ satisfies the conditions of Theorem
4.6, and so $V$ equals the quasipotential on
$N=\\{x:V(c,x)<\min_{y\notin\text{int}(\mathcal{S}_{c})}V(c,y)\\}$. (Here,
$\text{int}(\mathcal{S}_{c})$ refers to the interior of $\mathcal{S}_{c}$.) ∎

##  6 Examples

###  6.1 Calculations in [1]

We demonstrate how Theorem 9 in [1] can be easily derived from the theory
developed in Section 4.4. We formulate the statement in a way that is more
natural for the quasipotential.

######  Proposition 6.1.

Let $G=(\mathcal{S},\mathcal{C},\mathcal{R})$ be a one-dimensional birth-death
CRN with $\mathcal{S}=\\{A\\}$ and rate constants
$\\{\kappa_{a,b}\\}_{(a,b)\in\mathcal{R}}$. Let $D_{c}$ denote the domain of
attraction of $c>0$. Then, for $x\in D_{c}$, the gradient of the
quasipotential relative to $c$ is

$DQ(c,x)=\int_{c}^{x}\log\left(\frac{\sum_{(n,n-1)\in\mathcal{R}}\kappa_{n,n-1}s^{n}}{\sum_{(n,n+1)\in\mathcal{R}}\kappa_{n,n+1}s^{n}}\right)\text{ds}.$
(6.1)

###### Proof.

Fix $c>0$. The equation $H(x,p(x))=0$ can be made algebraic by using the
substitution $y(x)\doteq e^{p(x)}$ and clearing the denominator of $y$’s. For
$x\in D_{c}$, consider

$\displaystyle K(x,y)$ $\displaystyle\doteq$ $\displaystyle yH(x,\log y)$
$\displaystyle=$
$\displaystyle\sum_{(n,n+1)\in\mathcal{R}}\kappa_{n,n+1}x^{n}(y^{2}-y)+\sum_{(n,n-1)\in
R}\kappa_{n,n-1}x^{n}(1-y).$

Write $A=A(x)=\sum_{(n,n+1)\in\mathcal{R}}\kappa_{n,n+1}x^{n}$ and
$B=B(x)=\sum_{(n,n-1)}\kappa_{n,n-1}x^{n}$. Then

$K(x,y)=A(x)y^{2}-(A(x)+B(x))y+B(x).$

Viewing this as a quadratic equation in $y$, we find the pair of solutions to
be

$y_{\pm}=\frac{1}{2A}\left(A+B\pm\sqrt{(A-B)^{2}}\right).$

There are four cases and two possible solutions, depending on whether we take
the $+$ or $-$, and whether $A>B$ or $A\leqslant B$. We have

| $A>B$ | $A\leqslant B$
---|---|---
$+$ | $1$ | $B/A$
$-$ | $B/A$ | $1$

.

Since $p(x)=\log(y(x))$, we find that a solution to the original equation is

$p(x)=\log\left(\frac{B(x)}{A(x)}\right),$

which is the same as (6.1). It remains only to verify that this solution
satisfies the conditions of Theorem 4.6. Observe that $\Phi(x)=0$ if and only
if $A(x)=B(x)$. If $c$ is stable, then $\Phi(c)>0$ for $x<c$, $x\in D_{c}$,
and $\Phi(c)<0$ for $x>c,x\in D_{c}$, which implies respectively that
$A(x)<B(x)$ and $A(x)>B(x)$. Consequently, $p(x)<0$ for $x<c,x\in D_{c}$ and
$p(x)>0$ for $x>c,x\in D_{c}$. This is all that needs to be checked.

∎

As another illustration, consider Example 13 in [1].

######  Example 6.1.

Consider the reaction network

$A\rightarrow\emptyset\rightarrow 2A,$

with associated Hamiltonian

$H(x,p)=x(e^{-p}-1)+e^{2p}-1.$

This is a cubic equation in $y(x)=e^{p(x)}$, a solution of which is

$y(x)=\frac{1}{2}\left(\sqrt{1+4x}-1\right).$

The other solutions are $y(x)=-\left(\sqrt{1+4x}+1\right)/2$ and $y(x)=1$, the
first of which is not admissible because $y(x)=e^{p(x)}\geqslant 0$, while the
second corresponds to the deterministic system $\dot{x}=\Phi(x)$.

One can conclude that the gradient of the quasipotential $Q$ satisfies

$DQ(x)=\log\left(\sqrt{1+4x}-1\right)-\log 2.$

Indeed,

$\displaystyle H(x,DQ)$ $\displaystyle=$ $\displaystyle
x\left(\frac{1}{\frac{1}{2}\left(\sqrt{1+4x}-1\right)}-1\right)+\left(\frac{1}{2}\right)^{2}\left(\sqrt{1+4x}-1\right)^{2}-1$
$\displaystyle=$ $\displaystyle 0.$

Our procedure is much simpler than the one outlined in [1]. The same style of
calculation can be repeated for each example in [1], and in each case we find
that the steady-state Hamilton-Jacobi-Bellman equation is satisfied.

###  6.2 Examples from symplectic and Kähler geometry

While the results developed so far clarify why an appropriately scaled limit
of the invariant distribution identifies the quasipotential, they are simply a
recasting of known results in probabilistic terms. We suspect, however, that
there is some underlying geometric structure which describes the level curves
which determine the quasipotential, i.e. $(x,p):H(x,p)=0$.

One of the clues we have found most compelling is a series of papers by Miguel
Abreu, which builds on work of Delzant and Guillemin. The main idea is to view
the stoichiometric compatibility class as a convex polytope, which can then be
seen as the image of a moment map of an effective Hamiltonian action . Theorem
2.8 in [7] then says that there is a “canonical” toric complex structure
determined by the combinatorial data which defines the polytope. The
associated Riemannian Kähler metric is then determined by a potential
function, which turns out to be closely related to the quasipotential.

It is important to note that the quasipotential carries information only about
certain aspects of the CRN, since it solves the HJB for only a particular
class of infinite-horizon problems. We mention this because of related work on
toric varieties in CRNs by [33, 34]. The authors demonstrate that the
algebraic relationships between the variables $x$ and rate constants $\kappa$
imposed by the equation $D_{p}H(x,0)=0$ (see (3.8)) define an ideal that
corresponds to a toric variety if the network is complex-balanced. Though the
observations here might be related, we note that two different CRNs can have
the same quasipotential, but whose algebraic varieties $D_{p}H(x,0)=0$ can
look very different.

We present three examples. The first two illustrate the issue of viewing the
stoichiometric compatibility classes as convex polytopes. We formally change
coordinates to show that the quasipotential for complex balanced networks
agrees with the Kähler potential on an appropriately reparametrized polytope.
The last example is somewhat unrelated, but provides a “classical mechanics”
perspective on the network $A\rightleftarrows\emptyset$. We are compelled to
include this example because there is an elegant transformation which maps the
corresponding Hamiltonian into that of an “antiharmonic oscillator”.

######  Definition 6.2.

The potential introduced in Example 6.2 is defined in terms of combinatorial
data by

$g_{P}(x)=\sum_{i=1}^{k}(\langle\mu_{i},x\rangle-\lambda_{i})\log(\langle\mu_{i},x\rangle-\lambda_{i}),$
(6.2)

where $\mu_{i}\in\mathbb{Z}^{d}$ and $\lambda_{i}\in\mathbb{R}$. The polytope
is the set of points $x$ where $\langle\mu_{i},x\rangle\geq\lambda_{i}$,

$P\doteq\\{a_{i}(x)\geq 0,i=1,\dots,k\\},\quad
a_{i}(x)\doteq\langle\mu_{i},x\rangle-\lambda_{i},$

and for points in the interior of $P$, $g_{P}(x)$ is well-defined and smooth.
There are additional combinatorial assumptions on this setup and the reader is
referred to [7] for details.

######  Remark 6.3.

The potential $g_{P}$ is strictly convex on the interior of its domain of
definition, and in particular, $G_{P}(x)\doteq\text{Hess}(g_{P}(x))$ is
positive definite there. In particular, $G_{P}$ defines a Riemannian metric.
It is this metric which we suspect is related to the notion of “cost” or
“distance” induced by the quasipotential. Note that the Hessian with respect
to $x$ of $V(c,x)$ as defined in (5.2) is independent of $c$.

######  Example 6.2.

Consider the reaction network

$A_{1}\rightleftarrows A_{2},$

with both rate constants set equal to 1. It is easy to verify that this
network is complex balanced, and that $c=(1,1)$ is the unique stable steady
state in its stoichiometric compatibility class
$\mathcal{S}_{c}=\\{x\in\mathbb{R}^{2}_{+}:x_{1}+x_{2}=2\\}$. The
quasipotential relative to $c$ is

$Q(x)=\ell(x_{1})+\ell(x_{2})=x_{1}\log x_{1}+x_{2}\log x_{2},$

where $\ell(z)=z\log z-z+1$ for $z\geqslant 0$. Since $x_{1}+x_{2}=2$ and
$x_{i}\in[0,2]$ for $i=1,2$, we can parametrize $x=(x_{1},x_{2})$ in terms of
a single variable, say $y\in[-1,1]$. Let $x_{1}=1-y$ and $x_{2}=1+y$. Then

$Q(x(y))=(1-y)\log(1-y)+(1+y)\log(1+y).$

We will compare this network with Example 2.6 in [7], which we transcribe here
for convenience. Let $M=S^{2}\subset\mathbb{R}^{3}$ be the standard unit two-
sphere, equipped with its standard symplectic form, which is the standard area
form with total area $4\pi$. Consider the $S^{1}$ action on $M$ given by
rotation around any axis. The moment map $\phi:M\rightarrow\mathbb{R}$ is the
projection map to this axis of rotation, so the moment polytope is $P=[-1,1]$.
This polytope is determined by two affine functions,

$a_{-1}(x)=1+x,\quad a_{1}(x)=1-x,$

and the potential function $g_{P}$ defined on the interior of $P$ is

$g_{P}(y)=\frac{1}{2}[a_{1}(y)\log(a_{1}(y))+a_{2}(y)\log(a_{2}(y))]=\frac{1}{2}Q(y).$

######  Example 6.3.

Consider the network

$A_{1}+A_{2}\rightleftarrows 2A_{2}\rightleftarrows A_{2}+A_{3}.$

It is easy to verify that this network is also complex-balanced, and that
$c=(1,1,1)$ is the unique stable steady state in it stoichiometric
compatibility class. The quasipotential on
$\mathcal{S}_{c}=\left\\{x:\sum_{i=1}^{3}x_{i}=3\right\\}$ is therefore

$Q(x)=x_{1}\log x_{1}+x_{2}\log x_{2}+x_{3}\log x_{3}.$

Using the relation $x_{1}+x_{2}+x_{3}=3$, we can re-express the quasipotential
in terms of two variables $y_{1},y_{2}$, where $x_{3}=3-y_{1}-y_{2}$. Then

$\displaystyle Q(x(y))$ $\displaystyle=$ $\displaystyle y_{1}\log
y_{1}+y_{2}\log y_{2}+(3-y_{1}-y_{2})\log(3-y_{1}-y_{2}),$

which is exactly what one would obtain from the Kähler potential (6.2) using
$\mu_{1}=(1,0)$, $\lambda_{1}=0$, $\mu_{2}=(0,1)$, $\lambda_{2}=0$, and
$\mu_{3}=(-1,-1)$ and $\lambda_{3}=-3$. These three pairs define the convex
polytope $\\{y\in\mathbb{R}^{2}:y_{1}\geqslant 0,y_{2}\geqslant
0,y_{1}+y_{2}\leqslant 3\\}$, which is the projection of the simplex
$\mathcal{S}_{c}$ onto the $(x_{1},x_{2})$ plane.

######  Example 6.4.

In this example we use basic techniques from classical mechanics. A standard
reference is [15]. Consider again the reaction network from Example 2.1,

$A\rightleftarrows\emptyset,$

with both rate constants set equal to 1. The Hamiltonian is

$H(x,p)=x(e^{-p}-1)+e^{p}-1.$

We initially sought to identify the minimizers in $A_{T}(x,y)$, the space of
trajectories $\varphi$ which go from $x$ to $y$ in time $T$ (see (3.2)), and
we conjectured that if $(x(\cdot),p(\cdot))$ satisfies Hamilton’s equations of
motion (3.4)-(3.5) then $x(\cdot)$ would be a minimizing trajectory in
$A_{T}(x,y)$. Since the Hamiltonian $H$ is “time-independent”, we have the
_conservation of energy_ relation along curves $(x(\cdot),p(\cdot))$,

$\frac{d}{dt}H(x(t),p(t))=D_{x}H(x(t),p(t))\dot{x}(t)+D_{p}H(x(t),p(t))\dot{p}(t)=0,$

which can be seen by using the equations for $\dot{x}(t)$ and $\dot{p}(t)$.
Thus, we expect the pairs $(x(\cdot),p(\cdot))$ to live on level sets
$L_{E}\doteq\\{(x,p):H(x,p)=E\\}$. If we can determine the relationship
between the level sets and points $(x(t),p(t))$, $t\in[0,T]$, then we can have
a better understanding of the energy landscape of the network.

##### Harmonic oscillator.

To illustrate the approach we will take with the CRN, consider perhaps the
most well-studied Hamiltonian system is the _harmonic oscillator_ ,

$H_{\text{osc}}(x,p)=\frac{x^{2}}{2}+\frac{p^{2}}{2},$

with equations of motion

$\displaystyle\dot{x}(t)$ $\displaystyle=$ $\displaystyle p$ (6.3)
$\displaystyle\dot{p}(t)$ $\displaystyle=$ $\displaystyle-x.$ (6.4)

The level sets of $H_{\text{osc}}$ are ellipses (in this case, circles). While
the pair (6.3)-(6.4) can be solved directly, part of the success of classical
mechanics relies on parametrizing systems in terms of conserved quantities.
Consider the change of variables

$\displaystyle x(E,Q)$ $\displaystyle=$ $\displaystyle\sqrt{2E}\cos(Q)$
$\displaystyle p(E,Q)$ $\displaystyle=$ $\displaystyle\sqrt{2E}\sin(Q).$

Let $J=\frac{\partial(x,p)}{\partial(E,Q)}$ denote the Jacobian of this
transformation. It turns out that the above change of variables is
_symplectic_ , in the sense that

$J^{T}\Omega J=\Omega,$

where

$\Omega=\left(\begin{array}[]{cc}0&1\\\ -1&0\end{array}\right).$

Symplectic transformations preserve Hamilton’s equations, and in this case the
new Hamiltonian is $H_{\text{new}}(E,Q)=E$.

##### Antiharmonic oscillator.

Consider now the _antiharmonic oscillator_ ,

$H_{\text{anti}}(x,p)=\frac{p^{2}}{2}-\frac{x^{2}}{2}.$

Similarly to the previous example, we can change variables to obtain a new
Hamiltonian which depends only the “energy level” $E$. In this case we would
use hyperbolic trigonometric functions and see that the level sets are
hyperbolas.

##### Back to $A\rightleftarrows\emptyset$.

The level sets of the Hamiltonian are displayed in Figure 2.

Figure 2: Level sets of the Hamiltonian $H(x,p)=x(e^{-p}-1)+e^{p}-1$.

The level curves corresponding to $H=0$ look like the asymptotes of a
hyperbola. We will make this evident. Consider the change of variables

$\displaystyle x(X,P)$ $\displaystyle=$ $\displaystyle(X+1)(P+1)$
$\displaystyle p(X,P)$ $\displaystyle=$ $\displaystyle\log(P+1).$

It turns out this change of variables is also symplectic, and the new
Hamiltonian is given by

$H_{\text{new}}(X,P)=-XP.$

If $(x,p)\in\mathbb{R}_{+}\times\mathbb{R}$, then
$(X,P)\in[-1,\infty)\times(-1,\infty)$. It is straightforward to find a
symplectic transformation which turns this Hamiltonian into the antiharmonic
oscillator, though not necessarily with real coordinates. For more details,
see for instance [15].

## References

  * [1] D Anderson, G Craciun, M Gopalkrishnan, and C Wiuf.  Lyapunov functions, stationary distributions, and non-equilibrium potential for reaction networks.  Bulletin of Mathematical Biology, 77(9):1744–1767, 2015.
  * [2] Mark Freidlin and Alexander Wentzell.  Random Perturbations of Dynamical Systems.  Springer, 3rd edition, 2012.
  * [3] Thomas Kurtz.  Strong approximation theorems for density dependent markov chains.  Stochastic Processes Applications, 6:223–240, 1977/78.
  * [4] Thomas Kurtz.  The relationship between stochastic and deterministic models for chemical reactions.  Journal of Chemical Physics, 57(7):2976–2978, 1972.
  * [5] R. Jackson and F. Horn.  General mass action kinetics.  Archive for Rational Mechanics and Analysis, 47(2), 1972.
  * [6] David Anderson, Gheorghe Craciun, and Thomas Kurtz.  Product-form stationary distributions for deficiency zero chemical reaction networks.  Bulletin of Mathematical Biology, 72(8):1947–1970, 2010.
  * [7] Miguel Abreu.  Kahler geometry of toric manifolds in symplectic coordinates.  Arxiv, 0004122v1, 2008.
  * [8] G. Craciun.  Toric differential inclusions and a proof of the global attractor conjecture.  Arxiv preprint, 2015.
  * [9] Andrea Agazzi, Amir Dembo, and Jean-Pierre Eckmann.  On the geometry of chemical reaction networks: Lyapunov function and large deviations.  Arxiv preprint, 2018.
  * [10] Michael G. Crandall and Pierre-Louis Lions.  Viscosity solutions of hamilton-jacobi equations.  Trans. Amer. Math. Soc., 277:1–42, 1983.
  * [11] M. Bardi and I Capuzzo-Dolcette.  Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman equations.  Birkhauser Basel, 1997.
  * [12] Jin Feng and Thomas Kurtz.  Large Deviations for Stochastic Processes.  Mathematical Surveys and Monographs, 131., 2006.
  * [13] S. R. S. Varadhan.  Large deviations (special invited paper).  The Annals of Probability, 36(2):397–419, 2008.
  * [14] Paul Dupuis.  Representations and weak convergence methods for the analysis and approximation of rare events.  http://www.dam.brown.edu/people/documents/Padova.pdf, Last checked: February 23rd, 2019.
  * [15] Vladimir Arnold.  Mathematical Methods of Classical Mechanics.  Springer, 2nd edition, 1997.
  * [16] Andrea Agazzi, Amir Dembo, and Jean-Pierre Eckmann.  Large deviations theory for markov jump models of chemical reaction networks.  The Annals of Applied Probability, 28(3):1821–1855, 2018.
  * [17] Robert Patterson and Michiel Renger.  Large deviations of reaction fluxes.  Arxiv Preprint, 2018.
  * [18] Adam Schwartz and Alan Weiss.  Large deviations for performance analysis.  Chapman and Hall, 1995.
  * [19] A. Dembo and O. Zeitouni.  Large Deviations Techniques and Applications.  Springer, 2009.
  * [20] Lawrence C. Evans.  Partial Differential Equations.  American Mathematical Society, 1998.
  * [21] P. Nyquist, H. Hult, and B. Djehiche.  Min-max representations of viscosity solutions of hamilton-jacobi equations and applications in rare-event simulation.  Arxiv, 2014.
  * [22] Paul Dupuis and Hui Wang.  Subsolutions of an isaacs equation and efficient schemes for importance sampling.  Mathematics of Operations Research, 32(3):723–757, 2007.
  * [23] P. Dupuis, K. Spiliopoulos, and X. Zhou.  Escaping from an attractor: importance sampling and rest points i.  The Annals of Applied Probability, 25(5):2909–2598, 2015.
  * [24] Paul Dupuis and Thomas Dean.  Splitting for rare event simulation: a large deviation approach to design and analysis.  Stochastic Processes and their Applications, 119(2):562–587, 2009\.
  * [25] Paul Dupuis and Thomas Dean.  The design and analysis of a generalized restart/dpr algorithm for rare event simulation.  Annals of Operations Research, 189:63–102, 2011.
  * [26] P. Dupuis, M. Snarski, and A. Buijsrogge.  Splitting algorithms for rare event simulation over long time intervals.  (Submitted), (2019).
  * [27] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko.  The Mathematical Theory of Optimal Processes.  Interscience. English translation, 1962.
  * [28] Richard Vinter and Aram Arutyunov.  A simple ’finite approximations’ proof of the pontryagin maximum principle, under reduced differentiability hypotheses.  J. of Set Valued Analysis, 12:5–24, 2004.
  * [29] F. Horn.  Necessary and sufficient conditions for complex balancing in chemical kinetics.  Arch. Rational Mech. Anal., 49, 172-186, 1972.
  * [30] Martin Feinberg.  Complex balancing in general kinetic systems.  Arch. Rational Mech. Anal, 49:172–186, 1972.
  * [31] Martin Feinberg.  The existence and uniqueness of steady states for a class of chemical reaction networks.  Archive for Rational Mechanics and Analysis, 132:311–370, 1995\.
  * [32] Zhou Fang and Chuanhou Gao.  Lyapunov function partial differential equations for chemical reaction networks: Some special cases.  SIAM J. Appl. Dyn. Syst., 2(18):1163–1199, 2019.
  * [33] B. Sturmfels, G. Craciun, A. Dickenstein, and A. Shiu.  Toric dynamical systems.  Journal of Symbolic Computation, 44, 2009.
  * [34] Anne Shiu.  Algebraic methods for biochemical reaction network theory.  PhD Thesis, 2010.