# Riemannian Geometry of Optimal Driving and Thermodynamic Length
and its Application to Chemical Reaction Networks

Dimitri Loutchko<EMAIL_ADDRESS>Yuki Sughiyama  Tetsuya J.
Kobayashi http://research.crmind.net Institute of Industrial Science, The
University of Tokyo, 4-6-1, Komaba, Meguro-ku, Tokyo 153-8505 Japan

###### Abstract

It is known that the trajectory of an endoreversibly driven system with
minimal dissipation is a geodesic on the equilibrium state space. Thereby, the
state space is equipped with the Riemannian metric given by the Hessian of the
free energy function, known as Fisher information metric. However, the
derivations given until now require both the system and the driving reservoir
to be in local equilibrium. In the present work, we rederive the framework for
chemical reaction networks and thereby enhance its scope of applicability to
the nonequilibrium situation. Moreover, because our results are derived
without restrictive assumptions, we are able to discuss phenomena that could
not been seen previously. We introduce a suitable weighted Fisher information
metric on the space of chemical concentrations and show that it characterizes
the dissipation caused by diffusive driving, with arbitrary diffusion rate
constants. This allows us to consider driving far from equilibrium. As the
main result, we show that the isometric embedding of a steady state manifold
into the concentration space yields a lower bound for the dissipation when the
system is driven along the manifold. We give an analytic expression for this
bound and for the corresponding geodesic, and thereby are able to dissect the
contributions from the driving kinetics and from thermodynamics. Finally, we
discuss in detail the application to quasi-thermostatic steady states.

††preprint: APS/123-QED

##  I Introduction

In general thermodynamics, the importance of the Hessian of the free energy
function as a Riemannian metric has been recognized as early as 1975 by
Weinhold [1, 2, 3] and by Ruppeiner [4]. In 1983 it was discovered by Salamon
and Berry that this metric controls the dissipation when the system is driven
endoreversibly between equilibrium states [5]. For this reason, the concept to
determine curves of minimal dissipation from the Riemannian geometry of the
parameter space was termed thermodynamic length. The derivation required that
the driving is given by the explicit form
$\mathrm{d}\eta^{i}/\mathrm{d}t=k(\eta^{i}_{e}-\eta)$ for the extensive
variables $\eta^{i}$ characterizing the equilibrium state of the system and
$\eta_{e}^{i}$ being the respective variables of the reservoir, with a
universal rate constant $k$. More recently, a derivation using a statistical
approach was given by Crooks [6, 7] based on a cycling driving protocol. In
2012, Zulkowski, Sivak, Crooks and DeWeese have presented an exact and
explicit computation of geodesics for a driven particle diffusing in a one-
dimensional harmonic potential [8].

These previous approaches were valid for endoreversible driving, i.e. for such
processes that both the system and the driving reservoir are in local chemical
equilibrium and dissipation takes place only at the boundary between the
system and the reservoir. This theory does not directly apply to systems in
nonequilibrium or those with algebraic constraints. Chemical reaction networks
(CRN) are an important class of such systems and we extend the concept of
thermodynamic to chemical reaction networks in equilibrium and nonequilibrium
steady states.

The importance of chemical reaction networks derives from the fact that they
form the basis to model and understand complex chemical processes on a
molecular level [9, 10]. They are central to the understanding of biological
function as well as in many industrial applications [11, 12, 13, 14, 15, 16],
where tight control and high efficiency of the reactions are mandatory.

For CRN, a physically meaningful driving dynamics should be given by
$\mathrm{d}x^{i}/\mathrm{d}t=k_{i}(x^{i}_{e}-x)$, where $x^{i}$ is the
concentration the $i$-th chemical, $x^{i}_{e}$ the concentration of the
respective chemical in the reservoir and $k_{i}$ the diffusion rate constant
across the boundary. This dynamics generalizes the one given by Salamon and
Berry and, more importantly, it can drive the CRN away from chemical
equilibrium or from a steady state. Furthermore, models of CRN found in
biochemical applications are often open CRN, and as such might not exhibit
equilibrium states. For such networks, the dissipation caused by the driving
between steady states is of interest, yet it cannot be determined within the
framework available thus far.

The mathematical foundation of CRN theory was laid by Horn, Jackson and
Feinberg [17, 18, 19, 20, 21], and the theory has been refined and developed
by employing various mathematical disciplines such as graph theory [22, 23,
24], homological algebra [25], and differential geometry [26]. In recent
years, it has become clear that CRN are most naturally analyzed using
geometrical methods: Algebraic geometry has led to various important
foundational results [27, 28, 29, 30] for CRN with mass action kinetics as
well as applications to biologically important networks and phenomena [31, 32,
33]. Meanwhile, information geometry has been successfully used to analyze the
geometry of thermodynamics within a stoichiometric compatibility class [34,
35]. Finally, both approaches could be merged within the framework of Hessian
geometry [36, 37, 38, 39], generalizing the global structure of the
concentration space as studied in algebraic geometry and equipping it with the
thermodynamic structure from information geometry as well as including the
more classical results [40, 41, 42, 43, 44, 45, 46, 47] on the thermodynamics
of CRN. In [34, 35, 37, 38, 39], the information geometric aspects of Hessian
geometry, based on the Bregman divergence, have been thoroughly explored and
linked to thermodynamics. In this work, we focus on the differential geometric
aspects of this Hessian geometry.

We show that the dissipation rate caused by the driving of a CRN is determined
by the square of the norm $\|j_{D}\|_{g_{X}^{K}}^{2}$ of the driving flux
$j_{D}$ with respect to the Riemannian metric
$\mathrm{diag}(x^{1}k_{1}^{-1},\dotsc,x^{n}k_{n}^{-1})$ on the concentration
space. This holds arbitrarily far from chemical equilibrium and even for CRN
without equilibrium states. The Riemannian distance
$L_{X}:=d_{X}(x_{0},x_{1})$ between two concentration vectors $x_{0}$ and
$x_{1}$ on the concentration space is called thermodynamic length. Its square
is proportional to the minimal dissipation due to the driving in the case that
the driving flux $j_{D}$ and the reaction flux $j_{R}$ obey the time scale
separation $\|j_{D}\|\gg\|j_{R}\|$, which allow to ignore the reaction
effects. The geodesic equations can be solved explicitly to yield the optimal
trajectory of the system $x(t)$ and the optimal driving protocol $x_{e}(t)$.
This is Result 1, which is derived in Sec. II. The total minimal dissipation
$\Sigma_{X}^{min}$ of this process can be explicitly computed as

$\displaystyle\Sigma_{X}^{min}=\frac{4}{T}\sum_{i}k_{i}^{-1}\left(\sqrt{x^{i}_{1}}-\sqrt{x^{i}_{0}}\right)^{2}.$

This is Result 2 of this work. Without any assumptions on the time scales of
$j_{D}$ and $j_{R}$, the minimal total dissipation
$\Sigma_{\mathcal{V}^{ss}}^{min}$ caused by the driving of the system along
the steady state manifold $\mathcal{V}^{ss}$ between $x_{0}$ anf $x_{1}$ on
$\mathcal{V}^{ss}$ is proportional to the square length of the Riemannian
distance between $x_{0}$ and $x_{1}$ on $\mathcal{V}^{ss}$ (Result 3). The
optimal system trajectory is thus a geodesic and the optimal driving protocol
can be determined from the respective geodesic equations. The main result
(Result 4) follows from the isometric embedding of $\mathcal{V}^{ss}$ into
concentration space and states that $\Sigma_{\mathcal{V}^{ss}}^{min}$ is
bounded from below by $\Sigma_{X}^{min}$, i.e.

$\displaystyle\Sigma_{\mathcal{V}^{ss}}^{min}\geq\frac{4}{T}\sum_{i}k_{i}^{-1}\left(\sqrt{x^{i}_{1}}-\sqrt{x^{i}_{0}}\right)^{2}.$

This lower bound is universal in the sense that it holds for all possible
steady state manifolds irrespective of the reaction kinetics of the CRN. The
Results 1-4 are derived in Sec. II.

In Sec. III, the general results from Sec. II are applied to the important
class of quasi-thermostatic CRN. By definition, these are all CRN whose steady
state manifold is given by

$\displaystyle\mathcal{V}^{ss}=\\{x\in\mathbb{R}^{n}_{>0}|\log x-\log
x_{ss}\in\mathrm{Ker}[S^{T}]\\},$

where $x_{ss}$ is any steady state and $S$ is the stoichiometric matrix. This
class includes all CRN with equilibrium states that obey the thermodynamics of
an ideal solution and all complex balanced nonequilibrium steady states. An
explicit parametrization akin to the exponential family in information
geometry is used to derive an explicit expression for the Riemannian metric on
the parameter space, which yields numerically solvable geodesic equations for
any optimal driving problem. This parametrization is complemented by a
parametrization via conserved quantities, and its physical importance is
discussed. In particular, the parametrization via conserved quantities is an
important tool to treat the problem of arbitrary driving protocols with the
time scale separation $\|j_{R}\|\gg\|j_{D}\|$, which we aim to investigate in
the future. Finally, in Sec. IV, we discuss possible extensions of our
approach to systems with non-diagonal Onsager relations.

### Setup and notation

#### Notation.

The logarithm and exponential functions of vectors are taken componentwise and
the resulting vectors are elements of the linear dual of the original space,
i.e.

$\displaystyle\exp\left(\sum_{i}x^{i}e_{i}\right)=\sum_{i}\exp(x^{i})e^{i},$

where $\\{e_{i}\\}$ is a basis for the vector space and $\\{e^{i}\\}$ the
respective dual basis. The reason for this convention is that the $\exp$ and
$\log$ maps appearing here are Legendre transformations.

By the symbol $\circ$ we denote the Hadamard product between vectors, i.e.

$\displaystyle\left(\sum_{i}x^{i}e_{i}\right)\circ\left(\sum_{i}y^{i}e_{i}\right)=\sum_{i}x^{i}y^{i}e_{i}.$

For the differential geometric formalism, we refer to the textbook [48].

#### Chemical reaction networks (CRN).

In this work, we consider a chemical reaction network consisting of $n$
chemicals $X_{1},\dotsc,X_{n}$ and $r$ reactions $R_{1},\dotsc,R_{r}$. The
$j$th reaction is given by

$\displaystyle
R_{j}:\sum_{i=1}^{n}(S_{+})^{i}_{j}X_{i}\rightarrow\sum_{i=1}^{n}(S_{-})^{i}_{j}X_{i}$

with nonnegative integer coefficients $(S_{+})^{i}_{j}$ and $(S_{-})^{i}_{j}$.
These coefficients determine the reactants and products of the reaction
111Here we suppress the chemicals whose concentrations are kept constant by
the coupling to a reservoir, see [37] and [44] for more details on the role of
this coupling.. The structure of the network is thus encoded in the $n\times
r$ stoichiometric matrix $S=(S^{i}_{j})$ with matrix elements

$\displaystyle S^{i}_{j}=(S_{-})^{i}_{j}-(S_{+})^{i}_{j}.$

The state of the reaction network is characterized by a vector of nonnegative
concentration values $x=(x^{1},\dotsc,x^{n})\in\mathbb{R}^{n}_{>0}$, where
$x^{i}$ represents the concentration of the chemical $X_{i}$. We denote the
concentration space by $X:=\mathbb{R}^{n}_{>0}$. Finally, the dynamics of the
CRN is governed by the equation

$\displaystyle\frac{\mathrm{d}x}{\mathrm{d}t}=Sj_{R}+j_{D},$ (1)

where $j_{R}=(j_{R}^{1},\dotsc,j_{R}^{r})$ is the vector of reaction fluxes
and $j_{D}=(j_{D}^{1},\dotsc,j_{D}^{n})$ is the vector of the external driving
fluxes. Choosing a kinetic model is tantamount to specifying $j_{R}$ as a
function of $x$. For any given kinetic model $j_{R}(x)$, we denote the
respective steady state manifold by

$\displaystyle\mathcal{V}^{ss}=\\{x\in X|Sj_{R}(x)=0\\}.$

The results in this work have general validity as they are not based on the
choice of a particular kinetic model for $j_{R}$.

#### Chemical thermodynamics.

In chemical thermodynamics, the state of the system is characterized by
concentration vectors $x\in X\subset\mathbb{R}^{n}$ or, equivalently, by
chemical potential vectors $\mu$, which live in the dual space
$Y:=\left(\mathbb{R}^{n}\right)^{*}$. We denote the bilinear pairing between
the dual vector spaces by $\langle.,.\rangle$. The concentration and potential
vectors are connected by Legendre duality via the strongly convex free energy
function $\varphi(x)$ and its Legendre dual $\varphi^{*}(\mu)$ as
$\varphi^{*}(\mu)=\max_{x}\left\\{\langle x,\mu\rangle-\varphi(x)\right\\}$
and $\mu(x)=\operatorname*{argmax}_{\mu}\left\\{\langle
x,\mu\rangle-\varphi^{*}(\mu)\right\\}$ and analogous variational
characterizations of $\varphi(x)$ and $x(\mu)$. Additionally, the equality

$\displaystyle\varphi(x)+\varphi^{*}(\mu)=\langle x,\mu\rangle$ (2)

is satisfied for the Legendre dual pair of $x$ and $\mu$. Hereby, the
potential $\varphi(x)$ is the Gibbs free energy 222Here, we neglect the
constant term due to the coupling to external reservoirs. See [37] for a
thermodynamically thorough discussion of this potential function., which takes
the form

$\displaystyle\varphi(x)=\sum_{i=1}^{n}x^{i}\left(\mu^{0}_{i}+\log
x^{i}-1\right)$ (3)

for an ideal dilute solution (or, equivalently, an ideal gas). The vector
$\mu^{0}\in\left(\mathbb{R}^{n}\right)^{*}$ is the vector of standard chemical
potentials and we choose the energy scale such that $k_{B}T=1$. This explicit
form of $\varphi(x)$ will be used in calculations throughout the text. This
yields the $x$-dependence of $\mu$ as

$\displaystyle\mu=\mu^{0}+\log x$ (4)

and, using Eq. (2), the explicit form of the Legendre dual potential function
333As for $\varphi(x)$, the constant term arising due to the coupling to
external reservoirs is neglected. as

$\displaystyle\varphi^{*}(\mu)=\sum_{i=1}^{n}\exp(\mu_{i}-\mu^{0}_{i})=\sum_{i=1}^{n}x^{i}.$
(5)

The strictly convex function $\varphi(x)$ on the concentration space $X$ gives
rise to a Riemannian metric via its Hessian

$\displaystyle g_{X}\left(\frac{\partial}{\partial
x^{i}},\frac{\partial}{\partial
x^{j}}\right):=\frac{\partial^{2}\varphi(x)}{\partial x^{i}\partial x^{j}}.$

This is known as the Fisher information metric in information geometry [52]
and Weinhold or Ruppeiner metric in thermodynamics [1, 4]. Because the spaces
we work with have global coordinate systems, the metric can be globally
represented by a matrix. For the specific form of the convex function
$\varphi(x)$ given in Eq. (3), the metric $g_{X}$ is represented by the
diagonal matrix

$\displaystyle
g_{X}=\mathrm{diag}\left(\frac{1}{x^{1}},\dotsc,\frac{1}{x^{n}}\right).$ (6)

To account for the driving kinetics in the next section, we introduce a
diagonal weight matrix $K=\mathrm{diag}(k_{1},\dotsc,k_{n})$ with
$k_{i}\in\mathbb{R}_{>0}$ and the weighted Fisher information metric on $X$ as

$\displaystyle
g^{K}_{X}:=\mathrm{diag}\left(\frac{1}{k_{1}x^{1}},\dotsc,\frac{1}{k_{n}x^{n}}\right),$
(7)

which relates to the Hessian metric $g_{X}$ as
$g^{K}_{X}=K^{-1/2}g_{X}K^{-1/2}$.

The Legendre duality between $X$ and $Y$, together with the Hessian metric
$g_{X}$ and the Bregman divergence is the setup leading up the Hessian
geometry of CRN, which was thoroughly explored in [37] and [38]. Here, we
focus on the Riemannian geometry based on the metric $g_{X}^{K}$.

##  II Dissipation via thermodynamic length

Throughout this section, we analyze the case of an externally driven CRN
between two states $x_{0},x_{1}\in X$ in finite time $T$. The full dynamics of
the driven CRN is given by Eq. (1). We let $x(t)_{t\in[0,T]}$ denote the
resulting integral curve and also have $x_{0}=x(0)$ and $x_{1}=x(T)$. The
driving is due to the flux $j_{D}$, which is controlled by the concentrations
of chemicals in the reservoir, and the resulting dissipation is evaluated in
Sec. II.1. If the driving is faster than the chemical reactions, i.e.
$\|j_{R}\|\ll\|j_{D}\|$, then the dynamics reduces to

$\displaystyle\frac{\mathrm{d}x}{\mathrm{d}t}=j_{D}$ (8)

and the optimal trajectory $x(t)_{t\in[0,T]}$ is a geodesic in the
concentration space $X$, endowed with the weighted Fisher information metric
$g_{X}^{K}$. In this case, the optimal trajectory and the minimal dissipation
can be calculated analytically, as shown in Sec. II.2.

If such a time scale separation does not exist, then the problem of treating
an arbitrary driving field $j_{D}$ is rather intricate as the contributions to
the dissipation and to the dynamics generated by $j_{D}$ and by $j_{R}$ mix.
Even if the time scale separation $\|j_{R}\|\gg\|j_{D}\|$ holds, although this
constrains the trajectory $x(t)_{t\in[0,T]}$ to the steady state manifold, the
contribution of the dissipation due to chemical reactions cannot be neglected.
Because of the dependence on the details of the particular CRN and of the
driving protocol, we do not treat this general case here.

A possibility to circumvent the mixed contributions of $j_{D}$ and $j_{R}$
consists in restricting the driving field $j_{D}$ to the tangent space
$T_{x(t)}\mathcal{V}^{ss}$ of the steady state manifold and to impose the
initial condition $x_{0}\in\mathcal{V}^{ss}$. This implies that the dynamics
is again described by Eq. (8) and that the trajectory $x(t)_{t\in[0,T]}$ is
restricted to the steady state manifold. In this case, as shown in Sec. II.3,
the optimal trajectory is given by a geodesic on $\mathcal{V}^{ss}$ and can be
explicitly computed whenever a parameter manifold for $\mathcal{V}^{ss}$ is
known.

###  II.1 Dissipation

We assume that the system is in contact with an external reservoir and can
exchange all its chemicals $X_{1},\dotsc,X_{n}$ via diffusion with the
reservoir, with diffusion rate constants $k_{1},\dotsc,k_{n}$. The reservoir
is thereby characterized by the time-dependent concentration vector
$x_{e}(t)=(x_{e}^{1}(t),\dotsc,x_{e}^{n}(t))$ and the corresponding chemical
potential vector $\mu^{e}(t)$, which satisfy $\mu^{e}(t)=\mu^{0}+\log
x_{e}(t)$. The diffusion flux is given by Fick’s law, i.e. the flux vector
field on $X$ is

$\displaystyle j_{D}(x,t)=K\left(x_{e}(t)-x(t)\right),$ (9)

where $K=\mathrm{diag}(k_{1},\dotsc,k_{n})$ is the matrix of diffusion rate
coefficients. In the following, we suppress the time-dependence of the
variables for the sake of clarity of exposition. In addition, we assume that
the external reservoir potential vector is close to the system potential
vector, i.e. $\mu^{e}\approx\mu$. Then the relation

$\displaystyle\mu^{e}-\mu$ $\displaystyle=\frac{\partial\mu}{\partial
x}\left(x^{e}-x\right)=g_{X}\left(x^{e}-x\right)$ (10)

holds 444This can be directly seen from the Taylor expansion
$x_{e}^{i}-x^{i}=\exp(\mu^{e}_{i}-\mu^{0}_{i})-\exp(\mu_{i}-\mu^{0}_{i})=\exp(\mu_{i}-\mu^{0}_{i})\left(\exp(\mu^{e}_{i}-\mu_{i})-1\right)=x^{i}\left(\mu^{e}_{i}-\mu_{i}+\mathcal{O}\left((\mu^{e}_{i}-\mu_{i})^{2}\right)\right)$,
where $g_{X}$ is the Fisher information metric introduced in Eq. (6). This
leads, using Eqs. (9) and (10), to the dissipation rate at the boundary as

$\displaystyle\sigma_{D}$ $\displaystyle=\langle
j_{D},\mu^{e}-\mu\rangle=\langle j_{D},g_{X}\left(x^{e}(t)-x\right)\rangle$
$\displaystyle=\langle
j_{D},K^{-\frac{1}{2}}g_{X}K^{-\frac{1}{2}}j_{D}\rangle.$
$\displaystyle=\|j_{D}\|_{g^{K}_{X}}^{2},$ (11)

where $g_{X}^{K}$ is the weighted Fisher information metric introduced in Eq.
(7). This expression is the reason for introducing the weights to the standard
Fisher information metric: The dissipation rate at the boundary between the
system and the reservoir is given by the squared norm of the driving vector
$j_{D}(x,t)$. Together with Eq. (8), this allows to relate the dissipation
rate to the squared speed of the integral curve $x(t)$ and the total
dissipation to its energy (here, energy is used in the context of differential
geometry, which has no relation to the thermodynamical energy).

###  II.2 Optimal fast driving

If the time scale separation $\|j_{R}\|\ll\|j_{D}\|$ holds, which we call the
fast driving regime, the optimal system trajectory $x(t)$, the corresponding
driving protocol $x_{e}(t)$ and the minimal dissipation can be determined
analytically. In the fast driving regime, the dynamics of the driven system
obeys Eq. (8) and therefore the dissipation rate is given by Eq. (11):

$\displaystyle\sigma_{D}=\left\|\frac{\mathrm{d}x}{\mathrm{d}t}\right\|^{2}_{g_{X}^{K}}$

The total dissipation $\Sigma_{X}$ is therefore

$\displaystyle\Sigma_{X}=\int_{0}^{T}\left\|\frac{\mathrm{d}x}{\mathrm{d}t}\right\|^{2}_{g_{X}^{K}}\mathrm{d}t.$
(12)

Applying the Cauchy-Schwarz inequality to the integral yields

$\displaystyle\Sigma_{X}\geq\frac{1}{T}\left(\int_{0}^{T}\left\|\frac{\mathrm{d}x(t)}{\mathrm{d}t}\right\|_{g_{X}^{K}}\mathrm{d}t\right)^{2}=\frac{1}{T}L_{X}^{2},$
(13)

where $L_{X}$ is the length of the curve $x(t)$. For the trivial weight matrix
$K=I$, this is known as thermodynamic length [5] and we also use this term for
the lengths of curves on spaces endowed with the weighted Fisher metric. The
equality in Eq. (13) holds if and only if the speed of the curve,
$\left\|\frac{\mathrm{d}x}{\mathrm{d}t}\right\|_{g_{X}^{K}}$, is constant. In
this case, $x(t)$ is a geodesic. If it is a minimal length geodesic, then both
the length $L_{X}$ of the curve and its energy $\Sigma_{X}$ are minimized. In
this case, $L_{X}$ is equal to the distance between $x_{0}$ and $x_{1}$ on the
space $X$, i.e.

$\displaystyle
L_{X}^{min}=d_{X}(x_{0},x_{1})=\min_{\begin{subarray}{c}\gamma:[0,T]\rightarrow
X\\\ \gamma(0)=x_{0}\\\
\gamma(T)=x_{1}\end{subarray}}\int_{0}^{T}\left\|\frac{\mathrm{d}\gamma(t)}{\mathrm{d}t}\right\|_{g_{X}^{K}}\mathrm{d}t.$

For the driving field given by diffusion, i.e. $j_{D}=K(x_{e}-x)$, and the
weighted Fisher information metric $g_{X}^{K}$, the geodesic equations are

$\displaystyle\frac{\mathrm{d}^{2}x^{i}}{\mathrm{d}t^{2}}-\frac{1}{2x^{i}}\left(\frac{\mathrm{d}x^{i}}{\mathrm{d}t}\right)^{2}=0$

for $i=1,\dotsc,n$, cf. Appendix A. They have the explicit solution

$x^{i}(t)=-\left(\sqrt{x^{i}_{1}}-\sqrt{x^{i}_{0}}\right)^{2}\frac{t}{T}\left(1-\frac{t}{T}\right)+\\\
+x^{i}_{1}\frac{t}{T}+x^{i}_{0}\left(1-\frac{t}{T}\right)$ (14)

for $t\in[0,T]$. The second term of the expression describes a linear
trajectory in concentration space, which would yield the shortest path in the
standard Euclidean metric, and the first term shows the correction due to the
non-Euclidean geometry of the concentration space. Interestingly, the optimal
curve $x(t)$ does not depend on the details of the driving kinetics, i.e. it
is identical for all possible values of the diffusion rate constants $k_{i}$.
In other words, the optimal curve is determined purely by the thermodynamics
of the system via the Hessian $g_{X}$. This is the first main result of this
section:

######  Result 1.

The optimally driven curve $x(t)$ in the fast driving limit is given by Eq.
(14). It does not depend on the diffusion rate constants but only on the
thermodynamics of the system. The optimal driving protocol $x_{e}(t)$ is given
by

$\displaystyle x_{e}(t)=K^{-1}\frac{\mathrm{d}x}{\mathrm{d}t}+x(t),$ (15)

which follows from Eqs. (8) and (9).

With the explicit solution at hand, Eq. (11) yields the dissipation rate and
Eq. (12) the total dissipation in the driving process and thus the second main
result:

######  Result 2.

The optimally driven system under fast driving has the constant dissipation
rate

$\displaystyle\sigma_{D}^{min}=\frac{4}{T^{2}}\sum_{i}k_{i}^{-1}\left(\sqrt{x^{i}_{1}}-\sqrt{x^{i}_{0}}\right)^{2}$

and the minimal total dissipation is

$\displaystyle\Sigma_{X}^{min}=\frac{4}{T}\sum_{i}k_{i}^{-1}\left(\sqrt{x^{i}_{1}}-\sqrt{x^{i}_{0}}\right)^{2}.$
(16)

The results are not only significant in the fast driving regime but they
provide lower bounds for dissipation when the trajectory $x(t)$ is restricted
to the steady state manifold, as is shown in the following section.

###  II.3 Driving between arbitrary steady states

In this section we consider the case that the driving field $j_{D}$ is
constrained to the tangent space $T_{x(t)}\mathcal{V}^{ss}$ of the steady
state manifold with the initial condition $x_{0}\in\mathcal{V}^{ss}$. This
implies that the integral curve $x(t)$ is restricted to $\mathcal{V}^{ss}$ and
therefore the term $Sj_{R}(x(t))$ in Eq. (1) vanishes for all points on the
curve by definition of $\mathcal{V}^{ss}$. Thus, the dynamics is given by

$\displaystyle\frac{\mathrm{d}x}{\mathrm{d}t}=j_{D}$

and the dissipation rate is given by
$\sigma_{D}=\|\frac{\mathrm{d}x}{\mathrm{d}t}\|_{g_{X}^{K}}$ as before. Now,
the integral curve is restricted to $\mathcal{V}^{ss}$ and by the same
argument as before, the minimal length curve on $\mathcal{V}^{ss}$ between
$x_{0}$ and $x_{1}$ will minimize the dissipation. This can be made precise
(and accessible to direct computation as demonstrated in Sec. III) as follows:
Without loss of generality we assume that $\mathcal{V}^{ss}$ is connected and
thus parametrized by a connected manifold $M$, i.e. there is an embedding

$\displaystyle f:M\rightarrow X$

with $\mathrm{Im}[f]=\mathcal{V}^{ss}$. Defining the metric on $M$ as the
pullback

$\displaystyle g_{M}^{K}:=f^{*}g^{K}_{X}$

makes the embedding $f$ isometric and thus the dissipation rate can be written
as

$\displaystyle\sigma_{D}=\left\|\frac{\mathrm{d}x}{\mathrm{d}t}\right\|^{2}_{g_{X}^{K}}=\left\|\frac{\mathrm{d}m}{\mathrm{d}t}\right\|^{2}_{g_{M}^{K}},$

where $f(m)=x$. The total entropy production is bounded from below by

$\displaystyle\Sigma_{\mathcal{V}^{ss}}=\int_{0}^{T}\left\|\frac{\mathrm{d}m}{\mathrm{d}t}\right\|^{2}_{g_{M}^{K}}\mathrm{d}t\geq\frac{1}{T}L_{\mathcal{V}^{ss}}^{2},$

where $L_{\mathcal{V}^{ss}}$ is the length of the curve $m(t)$ on $M$, which
is bounded from below by the distance
$d_{\mathcal{V}^{ss}}\left(m_{0},m_{1}\right)$ between $m_{0}:=f^{-1}(x_{0})$
and $m_{1}:=f^{-1}(x_{1})$ on $M$, i.e.

$\displaystyle L_{\mathcal{V}^{ss}}\geq
d_{\mathcal{V}^{ss}}\left(m_{0},m_{1}\right)=\min_{\begin{subarray}{c}\gamma:[0,T]\rightarrow
M\\\ \gamma(0)=m_{0}\\\
\gamma(T)=m_{1}\end{subarray}}\int_{0}^{T}\left\|\frac{\mathrm{d}\gamma(t)}{\mathrm{d}t}\right\|_{g_{\eta}^{K}}\mathrm{d}t.$

This leads to the following result:

######  Result 3.

In the case that the system is driven on the steady state manifold, the
dissipation $\Sigma_{\mathcal{V}^{ss}}$, which is caused by the driving, is
minimized if and only if $m(t)$ is a geodesic on $M$. In this case, the
minimal dissipation reaches its lower bound, which is given by

$\displaystyle\Sigma_{\mathcal{V}^{ss}}^{min}=\frac{1}{T}d_{\mathcal{V}^{ss}}\left(m_{0},m_{1}\right)^{2}.$

The integral curve on the concentration space is $x(t)=f(m(t))$ and the
corresponding driving protocol is given explicitly by

$\displaystyle x_{e}(t)$
$\displaystyle=K^{-1}J_{f}\frac{\mathrm{d}m}{\mathrm{d}t}+f(m),$

which follows from Eq. (15) and where $J_{f}$ is the Jacobian of $f$ at
$m(t)$. This result holds independently of any time scale separation between
$j_{D}$ and $j_{R}$.

Finally, the isometry of the embedding $f$ implies that

$\displaystyle d_{X}(x_{0},x_{1})\leq
d_{\mathcal{V}^{ss}}\left(m_{0},m_{1}\right)$

and thus $\Sigma_{X}^{min}\leq\Sigma_{\mathcal{V}^{ss}}^{min}$. This is
illustrated in Fig. 1. We are led to the main result of this work:

######  Result 4.

The total dissipation $\Sigma_{\mathcal{V}^{ss}}$ for a driving process on the
steady state manifold is bounded from below by

$\displaystyle\Sigma_{\mathcal{V}^{ss}}\geq\frac{4}{T}\sum_{i}k_{i}^{-1}\left(\sqrt{x^{i}_{1}}-\sqrt{x^{i}_{0}}\right)^{2}.$

This lower bound is achieved if and only of the system is driven on the
optimal curve for the fast driving regime, given by Eq. (14).

Figure 1: Illustration of the approach to the lower bound for the total
dissipation $\Sigma_{\mathcal{V}^{ss}}^{min}$. The isometric embedding of
$\mathcal{V}^{ss}$ (represented by the green surface) into the concentration
space $X$ allows to compare the lengths of the geodesics
$d_{\mathcal{V}^{ss}}$ and $d_{X}$ on $\mathcal{V}^{ss}$ and $X$ between the
two points $x_{0},x_{1}\in\mathcal{V}^{ss}$ as $d_{\mathcal{V}^{ss}}\geq
d_{X}$. The relations $\Sigma_{X}^{min}=d_{X}^{2}/T$ and
$\Sigma_{\mathcal{V}^{ss}}^{min}=d_{\mathcal{V}^{ss}}^{2}/T$ show that
$\Sigma_{\mathcal{V}^{ss}}^{min}\geq\Sigma_{X}^{min}$. The latter can be
explicitly computed, cf. Eq. (16), and the former is a lower bound for the
dissipation when the system is driven on $\mathcal{V}^{ss}$ with any protocol,
thus yielding Result 4.

Note that nonequilibrium steady states have a nonzero dissipation rate caused
by continuously ongoing chemical reactions. Our approach does not take this
dissipation into account, but only gives a recipe for minimizing the
dissipation caused through the driving.

For an arbitrary CRN, the parametrization of the steady state manifold
$f:M\rightarrow X$ is, in general, not available. In fact, even when $j_{R}$
is given by mass action kinetics, the determination of steady states involves
highly nontrivial algebraic geometry and must be carried out for each CRN
individually [54]. Yet, the Result 2 holds irrespectively of any of the
complicated details of $\mathcal{V}^{ss}$.

For the class of quasi-thermostatic CRN, an explicit parametrization of
$\mathcal{V}^{ss}$ is available and the Result 3 can be made precise and
accessible for direct calculation. This is presented in the next section.

##  III Driving in quasi-thermostatic CRN

Quasi-thermostatic CRN, as defined by Horn and Jackson in [17], are
characterized by the following particular form of the steady state manifold

$\displaystyle\mathcal{V}^{ss}=\\{x\in X|\log x-\log
x_{ss}\in\mathrm{Ker}[S^{T}]\\},$ (17)

where $x_{ss}$ is a particular solution of the equation $Sj_{R}(x_{ss})=0$.
The quasi-thermostatic property derives its importance from the fact that all
equilibrium states of ideal solutions and complex balanced steady states of
CRN with mass action kinetics are quasi-thermostatic. The class of quasi-
thermostatic CRN is, however, even broader because there exist quasi-
thermostatic CRN that belong to neither of the two classes mentioned above
[55]. In the following, we give two possible parametrizations of quasi-
thermostatic steady states and illustrate each with an explicit calculation of
the geodesic and minimal dissipation based on Result 3.

###  III.1 Exponential Parametrization of Quasi-thermostatic CRN

By choosing a basis $u_{1}^{*},\dotsc,u_{q}^{*}$ of $\mathrm{Ker}[S^{T}]$ and
writing ${U^{*}}=(u^{*}_{1},\dotsc,u^{*}_{q})$ for the $n\times q$ matrix of
basis vectors, the characterization in Eq. (17) can be rewritten as

$\displaystyle\mathcal{V}^{ss}=\\{x_{ss}\circ\exp({U^{*}}\eta^{*})|\eta^{*}\in\mathbb{R}^{q}\\}.$
(18)

Thus, we obtain the parametrization of $\mathcal{V}^{ss}$ via the embedding

$\displaystyle f:\mathbb{R}^{q}$ $\displaystyle\rightarrow X$ (19)
$\displaystyle\eta^{*}$ $\displaystyle\mapsto
x_{ss}\circ\exp({U^{*}}\eta^{*}).$

This structure has been thoroughly studied in algebraic geometry [27, 30],
where $\mathcal{V}^{ss}$ was shown to be a toric variety and, using the
matrix-tree theorem explicit formulae for $x_{ss}$ were obtained. In the
particular case that $\mathcal{V}^{ss}$ is the equilibrium manifold of a
closed CRN, $x_{ss}$ is given by $\exp(-\mu^{0})$. In this work, the
analytification $\mathcal{V}^{ss}$ of the toric variety from algebraic
geometry is used. We call $\mathcal{V}^{ss}$ a toric manifold.

We now treat the Riemannian geometry of the quasi-thermostatic steady state
manifold $\mathcal{V}^{ss}$. The pullback of the metric $g_{X}$ to
$\mathbb{R}^{q}$ can be computed as

$\displaystyle g_{\mathbb{R}^{q}}^{K}$
$\displaystyle=f^{*}g_{X}^{K}=\left(\frac{\partial
x}{\partial\eta^{*}}\right)^{T}g_{X}^{K}\left(\frac{\partial
x}{\partial\eta^{*}}\right)$ $\displaystyle=J_{f}^{T}g_{X}^{K}J_{f}.$

With the Jacobian explicitly given by $J_{f}=D{U^{*}}$, where
$D:=\mathrm{diag}(x^{1},\dotsc,x^{n})$. Note that $D$ can be explicitly
expressed in the $\eta^{*}$ coordinates by using Eq. (19). We obtain the
desired metric on the parameter space $\mathbb{R}^{q}$ as

$\displaystyle
g_{\mathbb{R}^{q}}^{K}=({U^{*}})^{T}K^{-1/2}DK^{-1/2}({U^{*}})^{T}.$ (20)

Although we could not solve the geodesic equations on $\mathbb{R}^{q}$ in an
analytically closed form, they can be solved numerically for any given CRN.
The following example illustrates this.

#### Example.

Consider the following nonlinear chemical reaction network

A${4\textrm{B},}$

with stoichiometric matrix $S=(-1~{}4)^{T}$. The space $\mathrm{Ker}[S^{T}]$
is spanned by the vector $u^{*}=(4~{}1)$, the space of $\eta^{*}$ variables is
one-dimensional, i.e. it is given by $\mathbb{R}^{1}$, and the metric
$g_{\mathbb{R}^{1}}^{K}$ is given by

$\displaystyle
g_{\mathbb{R}^{1}}^{K}=16k_{\textrm{A}}^{-1}x^{\textrm{A}}_{ss}\exp(4\eta^{*})+k_{\textrm{B}}^{-1}x^{\textrm{B}}_{ss}\exp(\eta^{*})$

according to Eq. (20). The geodesic equation for $\eta^{*}$ is

$\displaystyle\frac{\mathrm{d}^{2}\eta^{*}}{\mathrm{d}t^{2}}+\left(\frac{64k\exp(3\eta^{*})+1}{32k\exp(3\eta^{*})+2}\right)\left(\frac{\mathrm{d}\eta^{*}}{\mathrm{d}t}\right)^{2}=0,$
(21)

with $k:=\frac{k_{B}x^{\textrm{A}}_{ss}}{k_{A}x^{\textrm{B}}_{ss}}$. Fig. 2A
shows the geodesic on the steady state manifold determined by Eq. (21) and the
geodesic on the space $X$, which is determined by Eq. (14), with the numerical
parameter values given in the figure caption. The curves are color coded by
the speed of the driving, which is given by the Euclidean norms of
$\frac{\mathrm{d}x}{\mathrm{d}t}$ and
$\frac{\mathrm{d}f(\eta^{*})}{\mathrm{d}t}$, respectively. The speed of the
curves is also shown in Fig. 2B and as a function of time. The speed increases
with increasing norm of the coordinate $x$ for both curves, whereby the system
driven on $\mathcal{V}^{ss}$ has the greater acceleration. As both curves are
geodesics, their speed with respect to the metric $g_{X}^{K}$ must be
constant, as is confirmed by the numerical results shown in Fig. 2(B). The
total dissipation for the optimal driving along $\mathcal{V}^{ss}$ is
$\Sigma^{min}_{\mathcal{V}^{ss}}=T\left\|\frac{\mathrm{d}f(\eta^{*})}{\mathrm{d}t}\right\|_{g^{K}_{X}}^{2}\approx
0.20$, with the lower bound
$\Sigma_{X}^{min}=T\left\|\frac{\mathrm{d}x}{\mathrm{d}t}\right\|_{g^{K}_{X}}^{2}\approx
0.17$.

For linear CRN, i.e. CRN where all the reaction complexes consist of exactly
one chemical, the geodesic on $X$, given by Eq. (14), and the geodesic on
$\mathcal{V}^{ss}$ coincide (see Appendix C). In this case the equality
$\Sigma^{min}_{\mathcal{V}^{ss}}=\Sigma_{X}^{min}$ holds and
$\Sigma^{min}_{\mathcal{V}^{ss}}$ can be computed by Eq. (16). In order to
illustrate the difference between $\Sigma^{min}_{\mathcal{V}^{ss}}$ and
$\Sigma_{X}^{min}$ we have thus chosen an example with a strong nonlinearity.
But even in the nonlinear case, $\Sigma_{X}^{min}$ provides a good estimate
for $\Sigma^{min}_{\mathcal{V}^{ss}}$. This is not only the case for the given
example but for various other CRN that we have analyzed numerically. In future
work, we aim to analytically quantify the error in estimating
$\Sigma^{min}_{\mathcal{V}^{ss}}$ by $\Sigma_{X}^{min}$ by relating it to the
curvature of $\mathcal{V}^{ss}$.

Figure 2:  Illustration of optimal driving of the
$\textrm{A}\rightleftharpoons 4\textrm{B}$ CRN, using the numerical values
$x^{\textrm{A}}_{ss}=0.1$, $x^{\textrm{B}}_{ss}=0.5$, $k_{\textrm{A}}=100$,
$k_{\textrm{B}}=20$, $\eta^{*}_{0}=-2$, $\eta^{*}_{1}=\log 2$ and $T=1$; with
$x_{0}=f(\eta^{*}_{0})$ and $x_{1}=f(\eta^{*}_{1})$. A. Trajectories of the
system being driven from $x_{0}$ to $x_{1}$. The lower curve $x(t)$ is the
optimal curve on the whole concentration space, given by Eq. (14) and the
curve $x_{e}(t)$ (yellow) is the corresponding optimal driving protocol, given
by Eq. (15). The curve $f(\eta^{*}(t))$ is the optimal curve on the steady
state manifold $\mathcal{V}^{ss}$ (light blue). The curves $x(t)$ and
$f(\eta^{*}(t))$ are color coded according to the speed of the driving, given
by the Euclidean norm of the tangent vector $\frac{\mathrm{d}x}{\mathrm{d}t}$.
In B, the Euclidean norms, which give the speed of the driving, and the
$\|.\|_{g_{X}^{K}}$ norms (their 10-fold values are plotted here for better
visibility), which govern the dissipation, of the tangent vectors $x(t)$ and
$f(\eta^{*}(t))$ are plotted for the course of the driving.

###  III.2 Parametrization via conserved quantities

In addition to the parametrization of the steady state manifold
$\mathcal{V}^{ss}$ of quasi-thermostatic CRN as a toric manifold, given by Eq.
(19), the manifold $\mathcal{V}^{ss}$ can be parametrized by the vector of
conserved quantities $\eta=({U^{*}})^{T}x$, with the matrix ${U^{*}}$
introduced in the previous section and $x\in\mathcal{V}^{ss}$ 555Because
${U^{*}}$ spans $\mathrm{Ker}[S^{T}]$ and thus
$\frac{\mathrm{d}({U^{*}})^{T}x}{\mathrm{d}t}=({U^{*}})^{T}Sj_{R}=0$, the
quantities $({U^{*}})^{T}x$ are conserved by the reaction dynamics of the
CRN.. This is based on Birch’s theorem [57, 27], which states that the
intersection of $\mathcal{V}^{ss}$ with the stoichiometric polytope
$P(\eta):=\left\\{x\in X|({U^{*}})^{T}x=\eta\right\\}$ exists and is unique,
i.e. the map $({U^{*}})^{T}:\mathcal{V}^{ss}\rightarrow\mathbb{R}^{q}$,
defined by $x\mapsto({U^{*}})^{T}x$, has a unique inverse given by

$\displaystyle h:\mathcal{E}$ $\displaystyle\rightarrow X$ (22)
$\displaystyle\eta$ $\displaystyle\mapsto\mathcal{V}^{ss}\cap P(\eta),$

where $\mathcal{E}:=({U^{*}})^{T}\mathcal{V}^{ss}$ is the image of
$\mathcal{V}^{ss}$. This is the parametrization of quasi-thermostatic steady
states by the vector of conserved quantities. The geometrical reason for the
uniqueness of this intersection is the generalized orthogonality between
$\mathcal{V}^{ss}$ and $P(\eta)$, as discussed in [38, 37]. Fig. 3 illustrates
this.

Figure 3: Illustration of the parametrization of points on a manifold of
quasi-thermostatic steady states $\mathcal{V}^{ss}$ via the vector of
conserved quantities $\eta$. The codimension of $\mathcal{V}^{ss}$ in $X$ is
the dimension of the stiochiometric polytope $P(\eta)$ and therefore a zero
dimensional intersection set is expected in general. The dual orthogonality
from Hessian geometry [37, 38] or, alternatively, Birch’s theorem from
algebraic geometry [27], ensure that the intersection consists of exactly one
point for each $\eta\in\mathcal{E}$. As such, this structure is a folitation
of $X$, with the base manifold being $\mathcal{V}^{ss}$ and the leaves being
$P(\eta)$. In other words, $\mathcal{V}^{ss}$ is the moduli space of
stiochiometric compatibility classes of the CRN. In the figure, it is
illustrated how the $P(\eta)$ and $P(\eta^{\prime})$ are parallel affine
spaces in $X$, that cover all of $X$ and are indexed by points on
$\mathcal{V}^{ss}$.

Although the map $h$ cannot be given in analytically closed form, its Jacobian
can be computed. First note that the map $h^{-1}f$ is the coordinate change
from $\eta^{*}$ to $\eta$ variables and its Jacobian is given by

$\displaystyle
J_{h^{-1}f}=\left(\frac{\partial\eta^{*}}{\partial\eta}\right)=({U^{*}})^{T}D{U^{*}},$

where $f$ is the map given in Eq. (19). Using $J_{f}=D{U^{*}}$ leads to

$\displaystyle
J_{h}=J_{f}J_{f^{-1}h}=J_{f}J_{h^{-1}f}^{-1}=D{U^{*}}\left(({U^{*}})^{T}D{U^{*}}\right)^{-1}.$
(23)

Thus, defining the weighted Fisher information metric $g_{\mathcal{E}}^{K}$ on
$\mathcal{E}$ by

$g_{\mathcal{E}}^{K}=\left(\frac{\partial\eta^{*}}{\partial\eta}\right)^{T}g_{\mathbb{R}^{q}}^{K}\left(\frac{\partial\eta^{*}}{\partial\eta}\right)=J_{f^{-1}h}^{T}g_{\mathbb{R}^{q}}^{K}J_{f^{-1}h}\\\
=\left(({U^{*}})^{T}D{U^{*}}\right)^{-1}({U^{*}})^{T}K^{-1}D{U^{*}}\left(({U^{*}})^{T}D{U^{*}}\right)^{-1}$
(24)

makes the parametrization $h:\mathcal{E}\rightarrow\mathcal{V}^{ss}$ isometric
666Note that an analytical expression for $g_{\mathcal{E}}^{K}$ in the $\eta$
coordinates is not available in general, because the matrix $D$ requires an
analytically closed form of the map $h$.. This is summarized in the following
diagram

${\left(\mathcal{E},g_{\mathcal{E}}^{K}\right)}$${\left(\mathcal{V}^{ss},\left.g^{K}_{X}\right|_{\mathcal{V}^{ss}}\right)}$${\left(\mathbb{R}^{q},g_{\mathbb{R}^{q}}^{K}\right)}$${\left(X,g^{K}_{X}\right).}$$\scriptstyle{h}$$\scriptstyle{\iota}$$\scriptstyle{f^{-1}}$$\scriptstyle{({U^{*}})^{T}}$$\scriptstyle{f}$$\scriptstyle{({U^{*}})^{T}}$
(25)

With the respective metrics $g^{K}_{\mathcal{E}}$, $g^{K}_{\mathbb{R}^{q}}$
and $\left.g^{K}_{X}\right|_{\mathcal{V}^{ss}}$, the maps in the upper row of
the diagram are isometries and thus the dissipation through finite time
driving can be computed on either of the three spaces. The embedding
$\iota:\mathcal{V}^{ss}\rightarrow X$ is isometric and allows for the explicit
computation of lower bounds for the total dissipation due to the driving as
stated in Result 4.

Finally, we remark that the coordinates $\eta$ and $\eta^{*}$ are Legendre
dual, see [38] for details.

###  III.3 Applicability of the two parametrizations

The metric $g^{K}_{\mathbb{R}^{q}}$ does have an analytically closed
expression in the $\eta^{*}$-coordinates via Eqs. (19) and (20). Therefore,
the space $\mathbb{R}^{q}$ of $\eta^{*}$-parameters, equipped with the metric
$g^{K}_{\mathbb{R}^{q}}$, is suitable for explicit computations. The example
in Sec. III.1 can be generalized to arbitrary quasi-thermostatic CRN, i.e. the
geodesic equations in $\eta^{*}$-coordinates on $\mathbb{R}^{q}$ can be
explicitly written and numerically solved for a given set of parameters to
obtain the optimal system trajectory, the optimal driving protocol and the
minimal total dissipation.

In contrast, there is no analytically closed expression for the metric
$g_{\mathcal{E}}^{K}$ on the space $\mathcal{E}$ of $\eta$-coordinates.
However, the vector of conserved quantities $\eta$ has an intuitive physical
meaning and is a natural tool to analyze the case of the time scale separation
$\|j_{R}\|\gg\|j_{D}\|$. In the case of arbitrary driving fields $j_{D}$, the
time scale separation $\|j_{R}\|\gg\|j_{D}\|$ results in a slow time scale
dynamics on $\mathcal{V}^{ss}$, which is determined by the pushforward
$({U^{*}})^{T}x(t)=\eta(t)$ of the trajectory $x(t)$ to the space of conserved
quantities. While the detailed analysis of this scenario is dependent on the
details of the CRN, the slow dynamics on $\mathcal{V}^{ss}$ can be reproduced
by the trajectory $h(({U^{*}})^{T}x(t))$. The corresponding driving protocol
is given by

$\displaystyle
x_{e}^{slow}(t)=K^{-1}J_{h}({U^{*}})^{T}\frac{\mathrm{d}x(t)}{\mathrm{d}t}+h(({U^{*}})^{T}x(t)).$
(26)

However, this formula requires the evaluation of the function $h$. This can be
circumvented by transforming the original driving field $j_{D}$ into a driving
field $j_{D}^{*}$ on the tangent space $T\mathcal{V}^{ss}$ of the steady state
manifold. For each $x\in\mathcal{V}^{ss}$ and $t\in[0,T]$, the change in
steady state, which is the state that the system will eventually relax to
after the disturbance caused by $j_{D}(x,t)$, is described by the change of
the vector of conserved quantities. The latter quantity is given by the
pushforward $J_{({U^{*}})^{T}}j_{D}(x,t)$ to the tangent space
$T_{({U^{*}})^{T}x}\mathcal{E}$. There is one and only one vector field
$j_{D}^{*}(x,t)$ on $T_{x}\mathcal{V}^{ss}$, which causes the same change in
the vector of conserved quantities as $j_{D}(x,t)$. It must be given by the
pushforward of $J_{({U^{*}})^{T}}j_{D}(x,t)$ to $\mathcal{V}^{ss}$, i.e.

$\displaystyle j_{D}^{*}(x,t):=J_{h}J_{({U^{*}})^{T}}j_{D}(x,t).$

Using the explicit formula for the Jacobian $J_{h}$ given in Eq. (23), one
obtains the desired vector field as

$\displaystyle
j_{D}^{*}(x,t)=D{U^{*}}\left(({U^{*}})^{T}D{U^{*}}\right)^{-1}({U^{*}})^{T}j_{D}(x,t).$

This approach is illustrated in Fig. 4. One verifies by a direct calculation
that $J_{({U^{*}})^{T}}j^{*}_{D}=J_{({U^{*}})^{T}}j_{D}$ and therefore $j_{D}$
and $j_{D}^{*}$ indeed induce the same slow dynamics on $\mathcal{V}^{ss}$.
Moreover, in the case that $j_{D}$ is already a vector field on the tangent
space of $\mathcal{V}^{ss}$, the equality $j^{*}_{D}=j_{D}$ holds. Therefore,
when considering arbitrary driving fields and the time scale separation
$\|j_{R}\|\gg\|j_{D}\|$, this construction plays an important role. In future
work, it has to be supplemented by an analysis of the dissipation caused by
the relaxation of the system to $\mathcal{V}^{ss}$ caused by $j_{R}$.

Figure 4: The construction of a driving vector field $j^{*}_{D}$ tangent to
$\mathcal{V}^{ss}$ from an arbitrary driving field $j_{D}$ such that both
fields generate the same slow dynamics on $\mathcal{V}^{ss}$. At each point
$x_{ss}\in\mathcal{V}^{ss}$, the vector field $j_{D}$ will drive the system
away from $\mathcal{V}^{ss}$ in general. However, if the time scale separation
$\|j_{R}\|\gg\|j_{D}\|$ holds, the steady state reached after such a
perturbation will be determined by the change in the vector of conserved
quantities $\Delta\eta=\eta^{\prime}-\eta$, i.e. by the pushforward
$J_{({U^{*}})^{T}}j_{D}$ to the tangent space
$T_{({U^{*}})^{T}x_{ss}}\mathcal{E}$. All tangent vectors at $x_{ss}$, which
lie in $P(\eta^{\prime})$ have the same pushforward. Among them, there is
exactly one, which is tangent to $\mathcal{V}^{ss}$. This is the desired
driving field $j_{D}^{*}$. Note that in this construction, both stoichiometric
polytopes $P(\eta)$ and $P(\eta^{\prime})$ lie in $T_{x_{ss}}X$ and not in
$X$.

##  IV Discussion and Outlook

In this work, we have analyzed the application of the Riemannian geometry with
a weighted Fisher information metric to the dissipation in driven chemical
reaction networks. This geometry should be thought of as an infinitesimal and
weighted version of the Hessian geometry established in [37, 38], where
information geometric aspects based on the Bregman divergence were analyzed.
Inspired by the classical results on thermodynamic length, we have shown that
the dissipation rate is given by the speed of the integral curve if either the
diffusion is fast or the driving generates curves that lie on the steady state
manifold. Thereby, we have presented a mathematically and physically precise
derivation of this formula based on the physical model of diffusion of
chemicals between the system and the reservoir. Moreover, we were able to
explicitly include the diffusion rate constants into the calculation. This was
the reason to introduce weights in the various Fisher information metrics.

The careful mathematical setup allowed us, almost effortlessly, to obtain
several new results for the dissipation in driven chemical reaction networks.
If the diffusion of the driving is faster than the chemical reactions, a
complete analytical solution was given: The optimal curve in concentration
space, the optimal driving protocol and the minimal total dissipation were
determined. Interestingly, it was found that the optimal curve does not depend
on the diffusion rate constants, i.e. it is determined purely by the
thermodynamics of the system, whereas the optimal driving protocol depends on
the kinetics.

Then, without any assumption on the time scales, if the driving proceeds along
the steady state manifold, the geodesic equations were formulated on the
respective parameter space. For the case of quasi-thermostatic steady states,
the equations were derived explicitly and thus they can be numerically solved
for any given CRN. More importantly, the isometric embedding of the steady
state manifold into the concentration space was used to obtain a lower bound
for the total dissipation. This bound is universal as it is independent of the
particular kinetic model $j_{R}(x)$ and is valid for any steady state
manifold. In all numerical examples we investigated thus far, this bound
turned out to provide a good estimate for the actual total dissipation.
Thereby the estimation error seems to grow with the curvature of the steady
state manifold and it will be a future challenge to make this precise.

Seen through the lens of general thermodynamics, this work shows that the
concept of thermodynamic length can be extended to the nonequilibrium
situation. This is necessitated by the physically meaningful driving of the
concentration variables $x$ instead of the extensive parameters $\eta$, which
characterize the equilibrium system. In regard to the original derivation in
[5], where only a specific driving field on the space of $\eta$ parameters is
considered, we have shown how the isometric embedding into a larger state
space can yield a more detailed understanding of the driving kinetics and
enable the calculation of explicit bounds.

The explicit analytical results obtained in this work shed new light on time-
reversed driving: Results 1 and 3 show that the optimal system trajectory for
the driving from $x_{1}$ to $x_{0}$ will be the time-reversed trajectory
caused by the optimal driving from $x_{0}$ to $x_{1}$. However, the driving
protocols will follows different paths (see Appendix C for a discussion of
time reversal illustrated with an example). This leads to the conclusion that
the approach taken in [6] to compute the total dissipation for the driving
between $x_{0}$ and $x_{1}$ is strictly valid only in the $T\rightarrow\infty$
limit for the trivial weight matrix $K=kI$ with a uniform rate constant
$k\in\mathbb{R}_{>0}$.

In future work, it will be interesting to address the general case of an
arbitrary driving protocol and without any assumptions on time scale
separation or restrictions of the direction of the driving field. Thereby, the
results from [37], in particular the Pythagorean theorems, will be useful to
evaluate the contribution to dissipation through chemical reactions.

However, there are more conceptual questions that have been raised by the
presented analysis. The weighted metric $g^{K}_{X}$ should be thought of as an
Onsager matrix for the fluxes $j_{D}$ and the corresponding forces
$\mu^{e}-\mu$. This is in line with the interpretation of the metric as the
generalized friction tensor given in [8]. It will be a rewarding future
challenge to work out the geometry of optimal driving for arbitrary Onsager
matrices. Moreover, the fact that the calculated optimal curve did not depend
on the particular details of the driving kinetics and was of thermodynamical
origin requires a thorough physical explanation. It might be necessary to
employ a mathematically more advanced framework to disentangle the
contributions from kinetics and thermodynamics, which are mixed within the
matrix $g^{K}_{X}$ in our current approach.

Apart from the theoretical insights and questions raised in this work, we have
provided another tool to better understand biochemical reaction networks and
to aid in the design and optimization of operational protocols in industrial
chemical applications.

###### Acknowledgements.

We thank Atsushi Kamimura and all other members of our lab for fruitful
discussions. This research is supported by JSPS KAKENHI Grant Numbers 19H05799
and 21K21308, and by JST CREST JPMJCR2011 and JPMJCR1927.

##  Appendix A Geodesic equations

On a Riemannian manifold $X$ with local coordinates $(x^{1},\dotsc,x^{n})$ and
Riemannian metric

$g^{ij}=g\left(\frac{\partial}{\partial x^{i}},\frac{\partial}{\partial
x^{j}}\right)$

the geodesic equations read

$\displaystyle\frac{\mathrm{d}^{2}x^{k}}{\mathrm{d}t^{2}}+\sum_{i,j=1}^{n}\Gamma^{k}_{ij}\frac{\mathrm{d}x^{i}}{\mathrm{d}t}\frac{\mathrm{d}x^{j}}{\mathrm{d}t}=0.$

Hereby, the Christoffel symbols are given by

$\displaystyle\Gamma^{m}_{ij}=\frac{1}{2}\sum_{k=1}^{n}\left[\frac{\partial}{\partial
x^{i}}g^{jk}+\frac{\partial}{\partial x^{j}}g^{ki}-\frac{\partial}{\partial
x^{k}}g^{ij}\right]g_{km}$

and $g_{km}$ are the matrix elements of the inverse of $g=(g^{ij})$, which is
given by a matrix in local coordinates. See, for example, [48] for more
details.

##  Appendix B Optimal driving in linear CRN

A CRN is called linear if each reaction has exactly one reactant and one
product. Linear CRN with mass action kinetics are mathematically equivalent to
Markov chain models on a graph, which are widely used in stochastic
thermodynamics. Whereas the latter strictly obey the conservation of
probability, the conserved quantities in linear CRN can be driven. In this
appendix it is shown that the optimal driving on the steady state manifold of
a linear quasi-thermostatic CRN, as described in Result 3, coincides with the
optimal driving on $X$, given in Result 1.

Without loss of generality, consider a linear CRN with one linkage class [59].
The stoichiometric matrix $S$ of the CRN is identical to the incidence matrix
of the digraph with $n$ vertices given by the chemicals $X_{1},\dotsc,X_{n}$
and directed edges between $X_{i}$ and $X_{j}$ iff there is a reaction
$X_{i}\rightarrow X_{j}$. Therefore $S^{T}$ has a one-dimensional kernel
spanned by $u^{*}=(1~{}1~{}\dots~{}1)$ and there is one conserved quantity
given by

$\displaystyle\eta=\sum_{i=1}^{n}x^{i}.$

The steady state manifold is given by
$\mathcal{V}^{ss}=\\{x_{ss}\circ\exp(\eta^{*})|\eta^{*}\in\mathbb{R}\\}$,
which yields the parametrization of all $x\in\mathcal{V}^{ss}$ via the
conserved quantity $\eta\in\mathbb{R}_{>0}$ as

$\displaystyle x=\frac{\eta}{\eta_{ss}}x_{ss},$

where $\eta_{ss}=\langle u^{*},x_{ss}\rangle=\sum_{i=1}^{n}x^{i}_{ss}$. Thus,
the metric on $\mathcal{E}=\mathbb{R}_{>0}$ is given by Eq. (24) as

$\displaystyle g_{\mathcal{E}}^{K}$
$\displaystyle=\left(({U^{*}})^{T}D{U^{*}}\right)^{-1}({U^{*}})^{T}DK^{-1}{U^{*}}\left(({U^{*}})^{T}D{U^{*}}\right)^{-1}$
$\displaystyle=\frac{1}{\eta\eta_{ss}}\sum_{i=1}^{n}x^{i}_{ss}k_{i}^{-1},$

which gives the Christoffel symbol $\Gamma=-\frac{1}{2\eta}$ and the geodesic
equation for $\eta$

$\displaystyle\frac{\mathrm{d}^{2}\eta}{\mathrm{d}t^{2}}-\frac{1}{2\eta}\left(\frac{\mathrm{d}\eta}{\mathrm{d}t}\right)^{2}=0.$

This equation has the explicit solution

$\eta(t)=-\left(\sqrt{\eta_{1}}-\sqrt{\eta_{0}}\right)^{2}\frac{t}{T}\left(1-\frac{t}{T}\right)+\\\
+\eta_{1}\frac{t}{T}+\eta_{0}\left(1-\frac{t}{T}\right),$ (27)

where $\eta_{0}=\eta(0)$ and $\eta_{1}=\eta(T)$. The ideal driving on
$\mathcal{V}^{ss}$ between the points $x_{0}=\frac{\eta_{0}}{\eta_{ss}}x_{ss}$
and $x_{1}=\frac{\eta_{1}}{\eta_{ss}}x_{ss}$ thus follows the trajectory

$\displaystyle x(t)=\frac{\eta(t)}{\eta_{ss}}x_{ss}.$ (28)

The total dissipation is given by

$\displaystyle\Sigma_{\mathcal{V}^{ss}}^{min}$
$\displaystyle=T\sigma^{\min}=T\frac{\mathrm{d}\eta(t)}{\mathrm{d}t}g_{\mathcal{E}}^{K}\frac{\mathrm{d}\eta(t)}{\mathrm{d}t}$
$\displaystyle=\frac{4}{\eta_{ss}T}\left(\sqrt{\eta_{1}}-\sqrt{\eta_{0}}\right)^{2}\sum_{i=1}^{n}x^{i}_{ss}k_{i}^{-1}$
$\displaystyle=\frac{4}{T}\sum_{i=1}^{n}k_{i}^{-1}\left(\sqrt{\frac{\eta_{1}}{\eta_{ss}}x^{i}_{ss}}-\sqrt{\frac{\eta_{0}}{\eta_{ss}}x^{i}_{ss}}\right)^{2}$
$\displaystyle=\frac{4}{T}\sum_{i=1}^{n}k_{i}^{-1}\left(\sqrt{x^{i}_{1}}-\sqrt{x^{i}_{0}}\right)^{2}.$

This dissipation $\Sigma_{\mathcal{V}^{ss}}^{min}$ is equal to the lower bound
$\Sigma_{X}^{min}$ stated in Result 2 and thus the driving on the equilibrium
manifold coincides with the optimal driving for the fast driving regime given
by Eq. (14). This can also be verified by directly substituting Eq. (27) into
Eq. (28).

##  Appendix C Reversal of driving

The Results 1 and 3 imply that a system which is optimally driven from $x_{0}$
to $x_{1}$ will follow the same trajectory as the optimally driven system from
$x_{1}$ to $x_{0}$ with time-reversal: Denoting the former trajectory by
$x^{f}(t)$ and the latter by $x^{b}(t)$, the relation $x^{b}(t)=x^{f}(T-t)$
for all $t\in[0,T]$ holds. However, the same results show that this is not
true for the driving protocols $x_{e}^{f}(t)$ and $x_{e}^{b}(t)$. They lie on
different curves and cannot be identified under time reversal, even if one
would introduce a lag time. This is illustrated in Fig. 5. Therein, the
example from Sec. III.1 is considered with the additional assumption that
there are no chemical reactions between $A$ and $B$ and thus that all of $X$
is the equilibrium manifold. The driving protocols $x_{e}^{f}(t)$ and
$x_{e}^{b}(t)$ are different non-intersecting curves.

Figure 5: Reversal of driving: an example. The example from Fig. 2 with the
assumption that no reactions take place between A and B is used here. The
optimal protocol $x^{f}_{e}(t)$ for the driving from $x_{0}$ to $x_{1}$ is
shown in yellow and the optimal protocol $x^{b}_{e}(t)$ for the driving from
$x_{1}$ to $x_{0}$ is shown in blue.

For $x_{e}^{f}(t)$ and $x_{e}^{b}(t)$ to become the same curves, it is
necessary and sufficient that the matrix $K$ is proportional to the identity
matrix and that the time limit $T\rightarrow\infty$ is taken. In this case,
they are identified with the geodesic $x(t)$.

## References

  * Weinhold [1975a] F. Weinhold, Metric geometry of equilibrium thermodynamics, The Journal of Chemical Physics 63, 2479 (1975a).
  * Weinhold [1975b] F. Weinhold, Metric geometry of equilibrium thermodynamics. II. scaling, homogeneity, and generalized gibbs–duhem relations, The Journal of Chemical Physics 63, 2484 (1975b).
  * Salamon  _et al._ [1984] P. Salamon, J. Nulton, and E. Ihrig, On the relation between entropy and energy versions of thermodynamic length, The Journal of chemical physics 80, 436 (1984).
  * Ruppeiner [1979] G. Ruppeiner, Thermodynamics: A riemannian geometric model, Physical Review A 20, 1608 (1979).
  * Salamon and Berry [1983] P. Salamon and R. S. Berry, Thermodynamic length and dissipated availability, Physical Review Letters 51, 1127 (1983).
  * Crooks [2007] G. E. Crooks, Measuring thermodynamic length, Physical Review Letters 99, 100602 (2007).
  * Sivak and Crooks [2012] D. A. Sivak and G. E. Crooks, Thermodynamic metrics and optimal paths, Phys. Rev. Lett. 108, 190602 (2012).
  * Zulkowski  _et al._ [2012] P. R. Zulkowski, D. A. Sivak, G. E. Crooks, and M. R. DeWeese, Geometry of thermodynamic control, Physical Review E 86, 041148 (2012).
  * Mikhailov and Ertl [2017] A. S. Mikhailov and G. Ertl,  _Chemical Complexity: Self-Organization Processes in Molecular Systems_ (Springer, 2017).
  * Alon [2019] U. Alon,  _An Introduction to Systems Biology: Design Principles of Biological Circuits_ (CRC Press LLC, 2019).
  * Eungdamrong and Iyengar [2004] N. J. Eungdamrong and R. Iyengar, Modeling cell signaling networks, Biology of the Cell 96, 355 (2004).
  * Rabinowitz and Vastag [2012] J. D. Rabinowitz and L. Vastag, Teaching the design principles of metabolism, Nature chemical biology 8, 497 (2012).
  * Shinar and Feinberg [2010] G. Shinar and M. Feinberg, Structural sources of robustness in biochemical reaction networks, Science 327, 1389 (2010).
  * Murray [2011] J. D. Murray,  _Mathematical Biology: I. An Introduction_ (Springer Science & Business Media, 2011).
  * Epstein and Pojman [1998] I. R. Epstein and J. A. Pojman,  _An Introduction to Nonlinear Chemical Dynamics: Oscillations, Waves, Patterns, and Chaos_ (Oxford University Press, New York, 1998).
  * Beard and Qian [2008] D. A. Beard and H. Qian,  _Chemical Biophysics: Quantitative Analysis of Cellular Systems_ (Cambridge University Press, Cambridge, 2008).
  * Horn and Jackson [1972] F. Horn and R. Jackson, General mass action kinetics, Archive for rational mechanics and analysis 47, 81 (1972).
  * Horn [1972] F. Horn, Necessary and sufficient conditions for complex balancing in chemical kinetics, Archive for Rational Mechanics and Analysis 49, 172 (1972).
  * Feinberg and Horn [1974] M. Feinberg and F. J. Horn, Dynamics of open chemical systems and the algebraic structure of the underlying reaction network, Chemical Engineering Science 29, 775 (1974).
  * Feinberg [1987] M. Feinberg, Chemical reaction network structure and the stability of complex isothermal reactors—I. the deficiency zero and deficiency one theorems, Chemical engineering science 42, 2229 (1987).
  * Feinberg [1988] M. Feinberg, Chemical reaction network structure and the stability of complex isothermal reactors—II. multiple steady states for networks of deficiency one, Chemical Engineering Science 43, 1 (1988).
  * Okada and Mochizuki [2016] T. Okada and A. Mochizuki, Law of localization in chemical reaction networks, Physical review letters 117, 048101 (2016).
  * Mochizuki and Fiedler [2015] A. Mochizuki and B. Fiedler, Sensitivity of chemical reaction networks: a structural approach. 1. examples and the carbon metabolic network, Journal of theoretical biology 367, 189 (2015).
  * Fiedler and Mochizuki [2015] B. Fiedler and A. Mochizuki, Sensitivity of chemical reaction networks: a structural approach. 2. regular monomolecular systems, Mathematical methods in the applied sciences 38, 3519 (2015).
  * Hirono  _et al._ [2021] Y. Hirono, T. Okada, H. Miyazaki, and Y. Hidaka, Structural reduction of chemical reaction networks based on topology, Physical Review Research 3, 043123 (2021).
  * Shinar  _et al._ [2009] G. Shinar, U. Alon, and M. Feinberg, Sensitivity and robustness in chemical reaction networks, SIAM Journal on Applied Mathematics 69, 977 (2009).
  * Craciun  _et al._ [2009] G. Craciun, A. Dickenstein, A. Shiu, and B. Sturmfels, Toric dynamical systems, Journal of Symbolic Computation 44, 1551 (2009).
  * Gopalkrishnan  _et al._ [2014] M. Gopalkrishnan, E. Miller, and A. Shiu, A geometric approach to the global attractor conjecture, SIAM Journal on Applied Dynamical Systems 13, 758 (2014).
  * Craciun  _et al._ [2019] G. Craciun, S. Muller, C. Pantea, P. Y. Yu, G. Craciun, S. Muller, C. Pantea, and P. Y. Yu, A generalization of Birchs theorem and vertex-balanced steady states for generalized mass-action systems, Mathematical Biosciences and Engineering 16, 8243 (2019).
  * Craciun and Sorea [2020] G. Craciun and M.-S. Sorea, The structure of the moduli spaces of toric dynamical systems, arXiv preprint arXiv:2008.11468 (2020).
  * Dickenstein  _et al._ [2019] A. Dickenstein, M. P. Millan, A. Shiu, and X. Tang, Multistationarity in structured reaction networks, Bulletin of Mathematical Biology 81, 1527 (2019).
  * Joshi and Shiu [2017] B. Joshi and A. Shiu, Which small reaction networks are multistationary?, SIAM Journal on Applied Dynamical Systems 16, 802 (2017).
  * Gross  _et al._ [2016] E. Gross, B. Davis, K. L. Ho, D. J. Bates, and H. A. Harrington, Numerical algebraic geometry for model selection and its application to the life sciences, Journal of The Royal Society Interface 13, 20160256 (2016).
  * Yoshimura and Ito [2021a] K. Yoshimura and S. Ito, Thermodynamic Uncertainty Relation and Thermodynamic Speed Limit in Deterministic Chemical Reaction Networks, Physical Review Letters 127, 160601 (2021a).
  * Yoshimura and Ito [2021b] K. Yoshimura and S. Ito, Information geometric inequalities of chemical thermodynamics, Physical Review Research 3, 013175 (2021b).
  * Shima [2007] H. Shima,  _The geometry of Hessian structures_ (World Scientific, 2007).
  * Sughiyama  _et al._ [2021] Y. Sughiyama, D. Loutchko, A. Kamimura, and T. J. Kobayashi, A hessian geometric structure of chemical thermodynamic systems with stoichiometric constraints, arXiv preprint arXiv:2112.12403 (2021).
  * Kobayashi  _et al._ [2021] T. J. Kobayashi, D. Loutchko, A. Kamimura, and Y. Sughiyama, Kinetic derivation of the hessian geometric structure in chemical reaction systems, arXiv preprint arXiv:2112.14910 (2021).
  * Sughiyama  _et al._ [2022] Y. Sughiyama, A. Kamimura, D. Loutchko, and T. J. Kobayashi, The chemical thermodynamics for growing systems, arXiv preprint arXiv:2201.09417 (2022).
  * Schnakenberg [1976] J. Schnakenberg, Network theory of microscopic and macroscopic behavior of master equation systems, Reviews of Modern Physics 48, 571 (1976).
  * Hill [1966] T. L. Hill, Studies in irreversible thermodynamics IV. diagrammatic representation of steady state fluxes for unimolecular systems, Journal of Theoretical Biology 10, 442 (1966).
  * Hill [2005] T. L. Hill,  _Free Energy Transduction and Biochemical Cycle Kinetics_ (Courier Corporation, 2005).
  * Schmiedl and Seifert [2007] T. Schmiedl and U. Seifert, Stochastic thermodynamics of chemical reaction networks, The Journal of Chemical Physics 126, 044101 (2007).
  * Polettini and Esposito [2014] M. Polettini and M. Esposito, Irreversible thermodynamics of open chemical networks. I. emergent cycles and broken conservation laws, The Journal of chemical physics 141, 07B610_1 (2014).
  * Ge and Qian [2016] H. Ge and H. Qian, Nonequilibrium thermodynamic formalism of nonlinear chemical reaction systems with Waage–Guldberg’s law of mass action, Chemical Physics 472, 241 (2016).
  * Rao and Esposito [2016] R. Rao and M. Esposito, Nonequilibrium Thermodynamics of Chemical Reaction Networks: Wisdom from Stochastic Thermodynamics, Physical Review X 6, 041064 (2016).
  * Rao and Esposito [2018] R. Rao and M. Esposito, Conservation laws shape dissipation, New Journal of Physics 20, 023007 (2018).
  * Do Carmo and Flaherty Francis [1992] M. P. Do Carmo and J. Flaherty Francis,  _Riemannian geometry_ , Vol. 6 (Springer, 1992).
  * Note [1] Here we suppress the chemicals whose concentrations are kept constant by the coupling to a reservoir, see [37] and [44] for more details on the role of this coupling.
  * Note [2] Here, we neglect the constant term due to the coupling to external reservoirs. See [37] for a thermodynamically thorough discussion of this potential function.
  * Note [3] As for $\varphi(x)$, the constant term arising due to the coupling to external reservoirs is neglected.
  * Amari [2016] S.-I. Amari,  _Information geometry and its applications_ , Vol. 194 (Springer, 2016).
  * Note [4] This can be directly seen from the Taylor expansion $x_{e}^{i}-x^{i}=\mathop{exp}\nolimits(\mu^{e}_{i}-\mu^{0}_{i})-\mathop{exp}\nolimits(\mu_{i}-\mu^{0}_{i})=\mathop{exp}\nolimits(\mu_{i}-\mu^{0}_{i})\left(\mathop{exp}\nolimits(\mu^{e}_{i}-\mu_{i})-1\right)=x^{i}\left(\mu^{e}_{i}-\mu_{i}+\mathcal{O}\left((\mu^{e}_{i}-\mu_{i})^{2}\right)\right)$.
  * Dickenstein [2016] A. Dickenstein, Biochemical reaction networks: An invitation for algebraic geometers, in  _Mathematical congress of the Americas_ , Vol. 656 (Contemp. Math, 2016) pp. 65–83.
  * Horn [1973] F. Horn, Stability and complex balancing in mass-action systems with three short complexes, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 334, 331 (1973).
  * Note [5] Because ${U^{*}}$ spans $\mathrm{Ker}[S^{T}]$ and thus $\frac{\mathrm{d}({U^{*}})^{T}x}{\mathrm{d}t}=({U^{*}})^{T}Sj_{R}=0$, the quantities $({U^{*}})^{T}x$ are conserved by the reaction dynamics of the CRN.
  * Pachter and Sturmfels [2005] L. Pachter and B. Sturmfels, eds.,  _Algebraic Statistics for Computational Biology_ (Cambridge University Press, Cambridge, 2005).
  * Note [6] Note that an analytical expression for $g_{\mathcal{E}}^{K}$ in the $\eta$ coordinates is not available in general, because the matrix $D$ requires an analytically closed form of the map $h$.
  * Feinberg [2019] M. Feinberg,  _Foundations of chemical reaction network theory_ (Springer, 2019).