Also at ]Universal Biology Institute, The University of Tokyo, 7-3-1, Hongo,
Bunkyo-ku, 113-8654, Japan.

# Geometry of Nonequilibrium Chemical Reaction Networks
and Generalized Entropy Production Decompositions

Tetsuya J. Kobayashi<EMAIL_ADDRESS>[  Institute of Industrial
Science, The University of Tokyo, 4-6-1, Komaba, Meguro-ku, Tokyo 153-8505
Japan  Dimitri Loutchko  Institute of Industrial Science, The University of
Tokyo, 4-6-1, Komaba, Meguro-ku, Tokyo 153-8505 Japan  Atsushi Kamimura
Institute of Industrial Science, The University of Tokyo, 4-6-1, Komaba,
Meguro-ku, Tokyo 153-8505 Japan  Yuki Sughiyama  Institute of Industrial
Science, The University of Tokyo, 4-6-1, Komaba, Meguro-ku, Tokyo 153-8505
Japan

###### Abstract

We derive the Hessian geometric structure of nonequilibrium chemical reaction
networks (CRN) on the flux and force spaces induced by the Legendre duality of
convex dissipation functions and characterize their dynamics as a generalized
flow. With this structure, we can extend theories of nonequilibrium systems
with quadratic dissipation functions to more general ones with nonquadratic
ones, which are pivotal for studying chemical reaction networks. By applying
generalized notions of orthogonality in Hessian geometry to chemical reaction
networks, we obtain two generalized decompositions of the entropy production
rate, each of which captures gradient-flow and minimum-dissipation aspects in
nonequilibrium dynamics.

††preprint: APS/123-QED

##  I Introduction

Thermodynamics aims at establishing the general description of thermal
systems. Although such a description was obtained for equilibrium
thermodynamics [1, 2], its extension to nonequilibrium systems has been
limited to specific situations and models. For near-equilibrium situations,
Onsager and Machlup evaluated the entropy production rate using the linear
approximation known as the force-flux relation [3, 4, 5, 6], which corresponds
to a quadratic dissipation function. With the recent development of
macroscopic fluctuation theory and stochastic thermodynamics [7, 8, 9], this
result was extended to far-from equilibrium situations in fluid dynamics and
diffusion processes in a continuous space, but the dissipation functions are
still quadratic even though they are generalized to be state dependent. Those
systems are characterized geometrically via the inner product structure
induced by the quadratic functions and associated formal Riemannian geometry.
However, the knowledge obtained from such situations and models is not
directly applicable to other systems with a discrete state space or with
nonlinearities in the governing equations, where the dissipation functions are
no longer quadratic. In such cases, the inner product structure is no longer
an adequate mathematical basis. Although Wasserstein geometry has recently
been introduced into the thermodynamics of diffusion processes as a new
geometric approach [10, 11], it also employs the formal Riemannian geometric
structure associated with the Wasserstein distance via the Benamou-Brenier
formula and thereby is not directly applicable to systems with nonquadratic
dissipation functions [12, 13, 14].

In this work, we clarify that this problem can be resolved by using the
Hessian (information) geometric structure of the flux and force spaces induced
by their nonlinear Legendre duality via convex (not necessarily quadratic)
dissipation functions. For systems with quadratic convex functions, the force
and the corresponding flux are related by a linear transformation induced by
the Riemannian metric. Hessian geometry or information geometry [15, 16, 17]
enables us to relate them by a nonlinear Legendre transformation induced by
the convex dissipation functions. Because the Riemannian metric is regarded as
the Legendre transformation with quadratic dissipation functions, Hessian
geometry is a nonlinear extension of the inner product structure in Riemannian
geometry. Even with being a nonlinear extension, Hessian geometry still
preserves many important aspects of the inner product structure in the form of
generalized orthogonalities, Pythagorean theorem, among others [16, 17].
Therefore, we can naturally and consistently extend various previous results
for systems with quadratic dissipation functions to those with nonquadratic
ones.

We derive and demonstrate the structure by focusing mainly on chemical
reaction networks (CRNs) because they are representative thermodynamic systems
with both discrete state space and nonlinearity in the governing rate
equations [18, 19, 20, 21]. CRNs also include Markov jump processes (MJPs) on
a graph, which is an important class of systems in stochastic thermodynamics
with a discrete state space but with a linear governing equation111It should
be noted that a linear governing equation does not necessarily mean that its
dissipation functions are quadratic.. In addition, CRNs are also important in
light of their biological and engineering applications [22, 23, 24]. The
dissipation functions of these systems are not always quadratic, but are more
general convex functions.

The nonquadratic property makes it difficult to dissect the equilibrium-like
aspects from their nonequilibrium ones by the notion of orthogonality from the
usual inner product structure, as was done for overdamped diffusion processes
[25]. As a result, some entropy decomposition formulas have not yet been
generalized to MJP and CNR. By employing the generalized notions of
orthogonality in Hessian geometry [15, 16, 17], we demonstrate that different
aspects of nonequilibrium CRNs can be dissected as generalized decompositions
of the entropy production rate (EPR). In particular, we derive a CRN version
of the Maes-Netočnỳ relation by generalizing the Helmholtz decomposition. If
the dissipation functions are quadratic, these results are automatically and
formally reduced to the original results obtained for diffusion processes.
Thus, Hessian geometry provides a natural generalization for resolving the
problem. Finally, we also discuss how Hessian geometry can potentially fill
the gap between quadratic and nonquadratic cases by geometrically capturing
the nonlinear dual relation between force and flux, unify the description of
those systems, and thereby extend the applicability of nonequilibrium
thermodynamics. As such an example, we mention a relation of our results with
the thermodynamic uncertainty relation, and describe how our formulation can
contribute to network thermodynamics and variational characterizations of
nonequilibrium systems.

This work is organized as follows. In sec. II, we define MJP and CRN. In sec.
III, we introduce the Legendre duality between flux and force, and the
associated notion of generalized flow. In sec. IV, we clarify the Hessian
geometric structure in the flux-force space and generalized orthogonalities,
which are one of the main contributions of this work. In sec. V, we
demonstrate how the generalized notions of orthogonality lead to different
types of EPR decompositions and their geometric meaning. In sec. VI, we verify
the obtained decompositions using CRN that can have equilibrium, complex-
balanced, and noncomplex balanced steady states. In sec. VII, We summarize our
results and provide possible applications and contributions of the Hessian
structure in the flux-force space to other thermodynamic problems.

##  II Models

In this section, we define MJP and CRN and show how MJP can be regarded as a
special case of CRN.

###  II.1 Markov Jump Processes

A reversible Markov jump process describes random jumps of noninteracting
particles on a graph $\mathbb{G}$ consisting of $N_{\mathbb{v}}$ vertices,
$\\{\mathbb{v}_{i}\\}_{i\in[1,N_{\mathbb{v}}]}$, and $N_{\mathbb{e}}$ oriented
edges $\\{\mathbb{e}_{e}\\}_{e\in[1,N_{\mathbb{e}}]}$. $k_{e}^{+}\geq 0$ is
the forward jump rate from the head of the oriented edge $\mathbb{e}_{e}$ to
its tail. $k_{e}^{-}\geq 0$ is the reverse jump rate from the tail to the head
of $\mathbb{e}_{e}$. For infinitely many such particles, we consider
$p_{i}(t)\in[0,1]$, the fraction of particles on vertex $\mathbb{v}_{i}$ at
time $t$. Then, the forward and reverse one-way fluxes on the $e$th edge are

$\displaystyle j^{+}_{e}(\bm{p})$
$\displaystyle=k_{e}^{+}p_{\mathbb{v}^{+}(\mathbb{e}_{e})},$ $\displaystyle
j^{-}_{e}(\bm{p})$
$\displaystyle=k_{e}^{-}p_{\mathbb{v}^{-}(\mathbb{e}_{e})},$ (1)

where $\mathbb{v}^{+}(\mathbb{e}_{e})$ and $\mathbb{v}^{-}(\mathbb{e}_{e})$
are the head and tail vertices of edge $\mathbb{e}_{e}$222Here, we have abused
the notation $\mathbb{v}^{+}(\mathbb{e}_{e})$ to indicate the index of the
vertex $\mathbb{v}^{+}(\mathbb{e}_{e})$.. The total flux function is the
difference of one-way flux functions as
$\bm{j}(\bm{p})=\bm{j}^{+}(\bm{p})-\bm{j}^{-}(\bm{p})\in\mathbb{R}^{N_{\mathbb{e}}}$
where
$\bm{j}^{\pm}(\bm{p})=(j^{\pm}_{1}(\bm{p}),\cdots,j^{\pm}_{N_{\mathbb{e}}}(\bm{p}))^{T}$.
Then, the dynamics of the state vector
$\bm{p}(t):=(p_{1}(t),\cdots,p_{N_{\mathbb{v}}}(t))^{T}\in\mathbb{R}^{N_{\mathbb{v}}}_{\geq
0}$ is represented by the master equation:

$\displaystyle\dot{\bm{p}}=-B\bm{j}(\bm{p})=-\mathrm{div}\bm{j}(\bm{p}),$ (2)

where $B\in\\{0,\pm 1\\}^{N_{\mathbb{v}}\times N_{\mathbb{e}}}$ is the
incidence matrix of graph $\mathbb{G}$ and $\mathrm{div}:=B$. More
specifically, for $B=(b_{i,e})$,

$\displaystyle b_{i,e}$ $\displaystyle:=+1$ $\displaystyle\mbox{if
$\mathbb{v}_{i}$ is the head of the edge $\mathbb{e}_{e}$},$ $\displaystyle
b_{i,e}$ $\displaystyle:=-1$ $\displaystyle\mbox{if $\mathbb{v}_{i}$ is the
tail of the edge $\mathbb{e}_{e}$},$ $\displaystyle b_{i,e}$
$\displaystyle:=0$ $\displaystyle\mbox{otherwise}.$

Eq. (2) is the continuity equation for diffusion on a graph and $B$ can be
regarded as the discrete divergence operator on a graph [26]. We also define
the head and tail incidence matrices, respectively, as $B^{+}:=\max[B,0]$ and
$B^{-}:=\max[-B,0]$. Thus, $B=B^{+}-B^{-}$.

Then, the flux functions are compactly described in a vector form as

$\displaystyle\bm{j}^{\pm}(\bm{p})=\bm{k}^{\pm}\circ(B^{\pm})^{T}\bm{p},$ (3)

where $\circ$ is the component-wise product of two vectors. In this work, we
assume that all edges describe reversible transitions, i.e., $k_{e}^{\pm}>0$
for all $e$. While the representation of the master equation by Eq. (2) is
different from that used conventionally in stochastic thermodynamics, they are
equivalent and this representation suits our purpose of unveiling the
relationship between MJP and CRN as well as the underlying geometric
structure.

Figure 1: Diagrammatic illustration of CRN (A) Diagrammatic representation of
chemical reaction equations with reaction rate constants
$\\{k_{i}\\}_{i=1,2,3}$. (B) The graph-theoretic structure of the CRN in (A).
$\mathbb{v}_{1}$, $\mathbb{v}_{2}$, and $\mathbb{v}_{3}$ correspond to complex
$2\mathbb{X}_{1}$, $2\mathbb{X}_{2}$, and $\mathbb{X}_{1}+\mathbb{X}_{2}$,
respectively. Each directed edge represents a pair of forward and reverse
reactions, and the direction of the edge indicates the direction of the
forward reaction.

###  II.2 Chemical Reaction Networks

Deterministic chemical reaction networks (CRN) are an important class of
macroscopic thermodynamic models in light of their historical role played in
thermodynamics since Gibbs [27] and of their wide range of applications in
engineering and biology [22, 23, 28, 24]. A reversible CRN is composed of a
collection of forward and reverse reaction pairs(Fig. 1 (A)), the $e$th
forward reaction of which is described by the following chemical reaction
equation [22]:

$\displaystyle\alpha_{1,e}\mathbb{X}_{1}+\cdots+\alpha_{N_{\mathbb{X}},e}\mathbb{X}_{N_{\mathbb{X}}}\rightarrow\beta_{1,e}\mathbb{X}_{1}+\cdots+\beta_{N_{\mathbb{X}},e}\mathbb{X}_{N_{\mathbb{X}}},$
(4)

where $\mathbb{X}_{i}$ is the $i$th molecular species, $N_{\mathbb{X}}$ is the
number of different kinds of molecular species, and
$\alpha_{i,e}\in\mathbb{N}_{\geq 0}$ and $\beta_{i,e}\in\mathbb{N}_{\geq 0}$
are the numbers of molecule $\mathbb{X}_{i}$ involved as the reactants and
products of the $e$th reaction, respectively. The stoichiometric vector of the
$e$th forward reaction is defined as
$\bm{s}_{e}:=(\beta_{1,e}-\alpha_{1,e},\cdots,\beta_{N,e}-\alpha_{N,e})^{T}$.
The stoichiometric matrix is $S:=(\bm{s}_{1},\cdots,\bm{s}_{N_{\mathbb{e}}})$.
The stoichiometric vector of the $e$th reverse reaction is obtained by just
changing the sign of the forward one: $-\bm{s}_{m}$. Thus, the stoichiometric
matrix $S$ defined only for the forward reactions is sufficient to
characterize a reversible CRN. Let $j_{e}^{+}(\bm{x})$ and $j_{e}^{-}(\bm{x})$
be the one-way fluxes of the $e$th forward and reverse reactions in which
$\bm{x}=(x_{1},\cdots,x_{N_{\mathbb{X}}})^{T}\in\mathcal{X}:=\mathbb{R}_{>0}^{N_{\mathbb{X}}}$
represents the concentration of molecules. The total flux is
$\bm{j}(\bm{x})=\bm{j}^{+}(\bm{x})-\bm{j}^{-}(\bm{x})$ where
$\bm{j}^{\pm}(\bm{x})=(j^{\pm}_{1}(\bm{x}),\cdot,j^{\pm}_{N_{\mathbb{e}}}(\bm{x}))^{T}$.
Then, we have the deterministic chemical rate equation (CRE) [22, 21, 28] as

$\displaystyle\dot{\bm{x}}=S\bm{j}(\bm{x})=-\mathrm{div}\bm{j}(\bm{x}).$ (5)

Eq. (5) is the continuity equation for the CRN with the divergence operator
$\mathrm{div}=-S$.

The intrinsic graph structure of CRN can be manifested by considering the sets
of reactants and products as vertices of the graph $\mathbb{G}$ connected by
the reaction edges (Fig. 1 (B)). Specifically, the set of reactants
$\mathbb{v}_{e}^{+}:=(\alpha_{1,e}\mathbb{X}_{1}+\cdots+\alpha_{1,e}\mathbb{X}_{N_{\mathbb{v}}})$
and the set of products
$\mathbb{v}_{e}^{-}:=(\beta_{1,e}\mathbb{X}_{1}+\cdots+\beta_{1,e}\mathbb{X}_{N_{\mathbb{v}}})$
are regarded as the head and tail vertices of the $e$th reaction edge,
respectively. Such sets (vertices) are called complexes in the CRN theory
[28]. Because each complex is a set of molecular species, a CRN is a special
kind of hypergraph in which multiple molecular species are connected by
oriented edges (reactions) [29]. Reflecting this hypergraph nature, $S$ can be
decomposed as

$\displaystyle S=-\Gamma B=-\Gamma(B^{+}-B^{-}),$ (6)

where $N_{\mathbb{X}}$ is the number of different vertices (complexes), $B$ is
the incidence matrix of the complex graph $\mathbb{G}$, and
$\Gamma=(\bm{\gamma}_{1},\cdots,\bm{\gamma}_{N_{\mathbb{v}}})\in\mathbb{Z}^{N_{\mathbb{X}}\times
N_{\mathbb{v}}}$ where $\bm{\gamma}_{i}$ specifies the molecular species
involved in the $i$th vertex (complex) as
$\mathbb{v}_{i}:=(\gamma_{1,i}\mathbb{X}_{1}+\cdots+\gamma_{N_{\mathbb{X}},i}\mathbb{X}_{N_{\mathbb{X}}})$.

If we adopt the law of mass action kinetics, the $e$th forward and reverse
reaction fluxes can be represented as

$\displaystyle
j^{\pm}_{e}(\bm{x})=k^{\pm}_{e}\sum_{i=1}^{N_{\mathbb{v}}}b_{i,e}^{\pm}\prod_{j=1}^{N_{\mathbb{X}}}x_{j}^{\gamma_{j,i}},$
(7)

where $k_{e}^{+}\in\mathbb{R}_{\geq 0}$ and $k_{e}^{-}\in\mathbb{R}_{\geq 0}$
are the reaction rate constants of the $e$th forward and reverse reactions,
respectively. In vector form, we can compactly represent it as

$\displaystyle\bm{j}^{\pm}(\bm{x})=\bm{k}^{\pm}\circ(B^{\pm})^{T}\bm{x}^{\Gamma^{T}},$
(8)

where
$\bm{x}^{\bm{\gamma}}:=\prod_{j=1}^{N_{\mathbb{X}}}x_{j}^{\gamma_{j}}\in\mathbb{R}_{>0}$
and
$\bm{x}^{\Gamma^{T}}:=(\bm{x}^{\bm{\gamma}_{1}},\cdots,\bm{x}^{\bm{\gamma}_{N_{\mathbb{v}}}})^{T}$.
We should note an important relation,
$(B^{\pm})^{T}\bm{x}^{\Gamma^{T}}=\bm{x}^{(\Gamma B^{\pm})^{T}}$, which holds
because every column vector of $B^{\pm}$ contains only one $+1$ and the others
are $0$ (see also Appendix for notation).

By comparing Eq. (2) and Eq. (3) with Eq. (5) and Eq. (8), we can see that
CRNs contain MJPs as a special case. Specifically, if we identify $\bm{x}$
with $\bm{p}$ and set $\Gamma=I$, where $I$ is the identity matrix, Eq. (5)
and Eq. (8) are reduced to Eq. (2) and Eq. (3), respectively. In other words,
a MJP is a CRN, each complex of which contains only one molecular species.
Such a CRN is called a monomolecular reaction network and its thermodynamic
nature has been investigated in the context of equilibrium and nonequilibrium
chemical thermodynamics, especially by Hill [18]. Because CRNs include MJPs,
we work only on CRNs in the following sections333The state of MJP is the
probability vector $\bm{p}(t)$. Thus, the conservation of probability
$\sum_{i}p_{i}(t)=1$ should be satisfied. To regard MJP as a CNR, we here do
not impose such a conservation in advance because the conservation of
probability over time, i.e., $\sum_{i}p_{i}(t)=\mathrm{const.}$ is
automatically satisfied by the property of the incidence matrix $B$ such that
$\bm{1}^{T}B=\bm{0}$. Thus, by restricting the dynamics of the system with the
initial condition satisfying $\sum_{i}x_{i}(0)=1$, the constraint
$\sum_{i}p_{i}(t)=1$ is automatically obtained. .

##  III Legendre Duality of flux-force

Next, we introduce the Legendre duality between flux and force for CRNs444We
abbreviate the thermodynamic force as force here., and summarize their
relation to entropy production. The specific type of convex function
introduced here was recently derived in the Macroscopic Fluctuation Theorem
(MFT) [30, 31, 32, 33, 34].

For a given pair of one-way fluxes $\bm{j}^{\pm}$, the total flux and force
are obtained as

$\displaystyle\bm{j}:=\bm{j}^{+}-\bm{j}^{-}\in\mathcal{J},\qquad\bm{f}:=\frac{1}{2}\ln\frac{\bm{j}^{+}}{\bm{j}^{-}}\in\mathcal{F},$
(9)

where $\mathcal{J}=\mathbb{V}$ is an $N_{\mathbb{e}}$-dimensional vector space
and $\mathcal{F}=\mathbb{V}^{*}$ is its linear-algebraic dual space555The
definition of the force here includes $1/2$, which does not appear in the
conventional definition of thermodynamic force. This is because we adopt the
derivation of this form of force from the large deviation theory[34]. We can
remove it by including $1/2$ in the definition of the dissipation functions..
The force of this form comes from the local detailed balance condition [35],
and is also consistent with macroscopic chemical thermodynamics with the mass
action kinetics [19, 21, 2]. If $\bm{j}$ and $\bm{f}$ are defined as in Eq.
(9), the entropy production rate (EPR) is obtained as

$\displaystyle\dot{\Sigma}:=2\langle\bm{j},\bm{f}\rangle=(\bm{j}^{+}-\bm{j}^{-})^{T}\ln\frac{\bm{j}^{+}}{\bm{j}^{-}}\geq
0,$ (10)

where $\langle\bm{j},\bm{f}\rangle:=\sum_{e=1}^{N_{\mathbb{e}}}j_{e}f_{e}$ is
the usual bilinear pairing on $\mathbb{V}\times\mathbb{V}^{*}$.

As Eq. (9) implies, the pair of flux and force has a nonlinear relationship.
To show that the relation is a Legendre duality, we introduce the frenetic
activity [36]666The definition of the activity here includes $2$, which does
not appear in the definition of the activity in [36]. This is because we adopt
the derivation of this form of activity from the large deviation theory[34].
We can remove it by including $2$ in the definition of the dissipation
functions.:

$\displaystyle\bm{\omega}:=2\sqrt{\bm{j}^{+}\circ\bm{j}^{-}}\in\mathbb{R}_{>0}^{N_{\mathbb{e}}}.$
(11)

Then, we have
$\bm{j}=\frac{\bm{\omega}}{2}\circ\left[e^{\bm{f}}-e^{\bm{-f}}\right]$. Thus,
if $\bm{\omega}$ is given, the force $\bm{f}$ can be converted to the
corresponding flux $\bm{j}$ that satisfies Eq. (9). This relation between the
pair $(\bm{j},\bm{f})$ is a one-to-one Legendre duality (LD) induced by the
following strictly convex smooth functions:

$\displaystyle\Psi^{*}_{\bm{\omega}}(\bm{f})$
$\displaystyle:=\bm{\omega}^{T}\left[\cosh(\bm{f})-\bm{1}\right],$
$\displaystyle\Psi_{\bm{\omega}}(\bm{j})$
$\displaystyle:=\bm{j}^{T}\sinh^{-1}\left(\frac{\bm{j}}{\bm{\omega}}\right)-\bm{\omega}^{T}\left[\sqrt{\bm{1}+\left(\frac{\bm{j}}{\bm{\omega}}\right)^{2}}-\bm{1}\right],$
(12)

which lead to the Legendre transformations:

$\displaystyle\bm{j}$
$\displaystyle=\partial_{\bm{f}}\Psi^{*}_{\bm{\omega}}(\bm{f})=\bm{\omega}\circ\sinh(\bm{f}),$
(13) $\displaystyle\bm{f}$
$\displaystyle=\partial_{\bm{j}}\Psi_{\bm{\omega}}(\bm{j})=\sinh^{-1}\left(\frac{\bm{j}}{\bm{\omega}}\right).$
(14)

We can easily verify that Eq. (13) and Eq. (14) are equivalent to Eq. (9) by
direct computation (see also Appendix ). In the following, we abbreviate
$\partial_{\bm{f}}\Psi^{*}$ and $\partial_{\bm{j}}\Psi_{\bm{\omega}}$ with
$\partial\Psi^{*}$ and $\partial\Psi_{\bm{\omega}}$ because we do not use
differentiation of $\Psi^{*}$ and $\Psi$ with respect to other variables in
this work. Moreover, the pair $(\bm{j},\bm{f})$ satisfies the Legendre
identity:

$\displaystyle\Psi^{*}_{\bm{\omega}}(\bm{f})+\Psi_{\bm{\omega}}(\bm{j})-\langle\bm{j},\bm{f}\rangle=0.$
(15)

The potential functions $\Psi^{*}_{\bm{\omega}}(\bm{f})$ and
$\Psi_{\bm{\omega}}(\bm{j})$ of the form in Eq. (12) were recently derived via
large deviation functions of the corresponding microscopic models in MFT [30,
32, 33, 34], where they are called dissipation functions. Both dissipation
functions are nonnegative and symmetric:

$\displaystyle\Psi_{\bm{\omega}}(\bm{j})=\Psi_{\bm{\omega}}(-\bm{j})\geq
0,\qquad\Psi^{*}_{\bm{\omega}}(\bm{f})=\Psi^{*}_{\bm{\omega}}(-\bm{f})\geq 0,$
(16)

and satisfy
$\min_{\bm{f}}\Psi^{*}_{\bm{\omega}}(\bm{f})=\Psi^{*}_{\bm{\omega}}(\bm{0})=0$
and $\min_{\bm{j}}\Psi_{\bm{\omega}}(\bm{j})=\Psi_{\bm{\omega}}(\bm{0})=0$.
From these properties, we can verify that $\bm{f}=\bm{0}$ and $\bm{j}=\bm{0}$
are Legendre dual: $\bm{0}=\partial\Psi^{*}_{\bm{\omega}}(\bm{0})$ and
$\bm{0}=\partial\Psi_{\bm{\omega}}(\bm{0})$, This is consistent with the
thermodynamic requirement that if the force is zero, the corresponding flux
becomes zero, and vice versa.

In addition, from the Legendre identity (Eq. (15)), the nonnegativity of the
EPR can also be attributed to the nonnegativity of the dissipation functions:

$\displaystyle\dot{\Sigma}/2=\langle\bm{j},\bm{f}\rangle=\Psi^{*}_{\bm{\omega}}(\bm{f})+\Psi_{\bm{\omega}}(\bm{j})\geq
0.$ (17)

It should be noted that the dissipation functions naturally provide a
decomposition of the EPR into nonnegative terms.

In the following, a pair of flux and force with the same decoration, e.g.,
$(\bm{j},\bm{f})$, $(\bm{j}^{\prime},\bm{f}^{\prime})$, or
$(\bm{j}(\bm{x}),\bm{f}(\bm{x}))$, represents a Legendre dual pair linked by
Eq. (13) and Eq. (14). Because of the one-to-one LD, the CRE (Eq. (5)) can be
represented as :

$\displaystyle\dot{\bm{x}}=S\bm{j}(\bm{x})=$ $\displaystyle
S\partial\Psi^{*}_{\bm{\omega}(\bm{x})}[\bm{f}(\bm{x})]$
$\displaystyle=-\mathrm{div}\partial\Psi^{*}_{\bm{\omega}(\bm{x})}[\bm{f}(\bm{x})].$
(18)

The equation of this form is a generalized flow driven by the force
$\bm{f}(\bm{x})$ [14, 32, 34]. The representation of the dynamics in this form
is not specific to CRN or MJP. Thus, it can cover other systems such as
overdamped diffusion by appropriately defining $\mathrm{div}$, $\Psi(\bm{x})$,
$\Psi^{*}(\bm{x})$, and $\bm{f}(\bm{x})$.

It should be noted that the derivation of Legendre duality and associated
quantities and relations is not dependent on the specific functional form of
the one-way fluxes, i.e., Eq. (8) or Eq. (3), the former of which is from the
kinetic law of mass action assumption. Thus, the result here may be applied to
a wider class of kinetic laws.

##  IV Hessian Geometry and generalized orthogonalities

Because of the nonlinearity of Legendre duality with the nonquadratic
dissipation functions, we can no longer employ the inner product structure
between flux and force and the associated formal Riemannian geometric notions
(Fig. 2). We clarify that this problem is resolved by employing Hessian
geometry and its geometric notions777It should be noted that, in differential
geometry, Hessian geometry is a class of Riemannian geometry with a metric
given by a Hessian matrix, which additionally induces the Legendre dual
structure. . The elucidation of the fundamental roles of Hessian geometry in
nonequilibrium dynamics of MJP and CRN is one of the main contributions of
this work.

Figure 2: (A) Linear-algebraic dual spaces $(\mathcal{J},\mathcal{F})$ with
the inner product structure defined by a metric via $M$ and $M^{*}$. Linear
subspaces in $\mathcal{J}$, e.g., blue solid and red dashed lines, are mapped
to linear subspaces in $\mathcal{F}$. (B) Linear-algebraic dual spaces
$(\mathcal{J},\mathcal{F})$ with the Hessian geometric structure defined by
Legendre transformations $\partial\Psi$ and $\partial^{*}\Psi$. A linear
subspace in $\mathcal{J}$, e.g., red dashed line, is mapped to a curved
subspace in $\mathcal{F}$ whereas a linear subspace in $\mathcal{F}$, e.g.,
blue solid line, is mapped to a curved subspace in $\mathcal{J}$.

###  IV.1 Hessian geometry

Hessian geometry is the geometry induced by convex potential functions
$\Psi(\bm{j})$ and $\Psi^{*}(\bm{f})$ to a pair of dual spaces $\mathcal{J}$
and $\mathcal{F}$ [15, 16]888Hessian geometry is generally defined on affine
manifolds. However, because we relate the Hessian geometric structure with
inner product one, we here restrict $\mathcal{J}$ and $\mathcal{F}$ to be the
dual vector spaces for simplicity of presentation.. In this section,
$\Psi(\bm{j})$ and $\Psi^{*}(\bm{f})$ are not restricted to the specific form
of Eq. (12) to obtain general results. We consider the case for each fixed
$\bm{\omega}$999This also means fixed $\bm{x}$ if $\bm{\omega}$ is a function
of $\bm{x}$. This is similar to the inner product structure of the tangent and
cotangent spaces at each point $\bm{x}$ on a Reimannian manifold.. Thus, the
dependence of $\Psi(\bm{j})$ and $\Psi^{*}(\bm{f})$ on $\bm{\omega}$ is
omitted for the sake of notational simplicity. They are assumed more generally
to be just smooth and strictly convex functions satisfying the Legendre
duality:

$\displaystyle\Psi^{*}(\bm{f})$
$\displaystyle=\max_{\bm{j}}\left[\langle\bm{j},\bm{f}\rangle-\Psi(\bm{j})\right],$
(19) $\displaystyle\Psi(\bm{j})$
$\displaystyle=\max_{\bm{f}}\left[\langle\bm{j},\bm{f}\rangle-\Psi^{*}(\bm{f})\right].$
(20)

Instead of the association of $\bm{j}$ and $\bm{f}$ by a linear transformation
$\bm{j}=M^{*}\bm{f}$ on the inner product space (Fig. 2 A), $\bm{j}$ and
$\bm{f}$ are associated by the Legendre transformation(Fig. 2 B):

$\displaystyle\bm{j}$ $\displaystyle=\partial_{\bm{f}}\Psi^{*}(\bm{f}),$
$\displaystyle\bm{f}$ $\displaystyle=\partial_{\bm{j}}\Psi(\bm{j}).$ (21)

Hessian geometry is fundamental for capturing nonlinear geometry induced by
the convex functions. As an important application, it has played the essential
role in describing the geometry of statistical models in information geometry
[16, 17] and that of equilibrium thermodynamics[37, 2, 38].

To regard the Hessian geometric structure as a nonlinear generalization of the
inner product structure101010It should be noted that the Hessian geometric
structure is not necessarily restricted to this interpretation as an extension
of the inner product structure., we additionally assume that $\Psi(\bm{j})$
satisfies the symmetry condition $\Psi(\bm{j})=\Psi(-\bm{j})$. From this
symmetry condition for $\Psi(\bm{j})$, we also obtain the symmetry of
$\Psi^{*}(\bm{f})=\Psi^{*}(-\bm{f})$ via the Legendre duality. From the
symmetries, the minimum of $\Psi(\bm{j})$ and $\Psi^{*}(\bm{f})$ is attained
at $\bm{0}$: $\arg\min_{\bm{j}}\Psi(\bm{j})=\bm{0}$ and
$\arg\min_{\bm{f}}\Psi^{*}(\bm{f})=\bm{0}$. Because the Legendre
transformation Eq. (21) is independent of a constant in $\Psi(\bm{j})$ and
$\Psi^{*}(\bm{f})$, without losing generality, we can assume that
$\min_{\bm{j}}\Psi(\bm{j})=\Psi(\bm{0})=\min_{\bm{j}}\Psi(\bm{j})=\Psi(\bm{0})=0$.
Thus,

$\displaystyle\Psi(\bm{j})=\Psi(-\bm{j})\geq
0,\qquad\Psi^{*}(\bm{f})=\Psi^{*}(-\bm{f})\geq 0.$ (22)

In the following, we call $\Psi(\bm{j})$ and $\Psi^{*}(\bm{f})$ with this
property dissipation functions.

To demonstrate that the Hessian geometric structure includes the inner product
structure, suppose that $\Psi(\bm{j})$ is a quadratic function as

$\displaystyle\Psi(\bm{j})=\frac{1}{2}\langle\bm{j},M\bm{j}\rangle:=\frac{1}{2}\|\bm{j}\|^{2}_{M},$
(23)

where $M\in\mathbb{R}^{N_{\mathbb{e}}\times N_{\mathbb{e}}}$ is a positive
definite matrix. Then, its Legendre dual becomes

$\displaystyle\Psi^{*}(\bm{f})=\frac{1}{2}\langle
M^{-1}\bm{f},\bm{f}\rangle=\frac{1}{2}\langle
M^{*}\bm{f},\bm{f}\rangle:=\frac{1}{2}\|\bm{f}\|^{2}_{{M}^{*}},$ (24)

where ${M}^{*}:={M}^{-1}$. Thus, $\Psi(\bm{j})$ and $\Psi^{*}(\bm{f})$ are
reduced to the squared norms associated with the metric ${M}$ and ${M}^{*}$.
The Legendre transformations become

$\displaystyle\bm{j}$ $\displaystyle={M}\bm{f}=({M}^{*})^{-1}\bm{f},$
$\displaystyle\bm{f}$ $\displaystyle={M}^{*}\bm{j}={M}^{-1}\bm{j}.$ (25)

These are the linear pairing of $\bm{j}$ and $\bm{f}$ via the metric ${M}$.
Furthermore, for the paired $\bm{j}$ and $\bm{f}$ by ${M}$, $\Psi(\bm{j})$ and
$\Psi^{*}(\bm{f})$ are essentially identical to the entropy production as

$\displaystyle\frac{1}{2}\langle\bm{j},\bm{f}\rangle=\Psi(\bm{j})=\Psi^{*}(\bm{f})=\dot{\Sigma}/4.$
(26)

However, this identity of $\Psi(\bm{j})$ and $\Psi^{*}(\bm{f})$ no longer
holds for a nonquadratic $\Psi(\bm{j})$ and $\Psi^{*}(\bm{f})$ pair. In
addition, for quadratic cases, the generalized flow (Eq. (18)) reduces to the
flow on a Riemannian manifold:

$\displaystyle\dot{\bm{x}}=-\mathrm{div}M^{*}_{\bm{x}}[\bm{f}(\bm{x})].$ (27)

As a result, various notions of geometry in the inner product space are
generalized in Hessian geometry.

###  IV.2 Generalized distance and Bregman divergence

Owing to the nonlinearity, the notion of distance, i.e., is generalized at
least into two versions.

The first version is to define a generalized distance via the dissipation
functions as

$\displaystyle\mathcal{D}_{H}(\bm{j},\bm{j}^{\prime})$
$\displaystyle:=\Psi(\bm{j}-\bm{j}^{\prime}),$ (28)
$\displaystyle\mathcal{D}^{*}_{H}(\bm{f},\bm{f}^{\prime})$
$\displaystyle:=\Psi^{*}(\bm{f}-\bm{f}^{\prime}).$ (29)

From Eq. (22), $\mathcal{D}_{H}(\bm{j},\bm{j}^{\prime})$ is nonnegative,
symmetric
$\mathcal{D}_{H}(\bm{j},\bm{j}^{\prime})=\mathcal{D}_{H}(\bm{j}^{\prime},\bm{j})$
and satisfies $\mathcal{D}_{H}(\bm{j},\bm{j}^{\prime})=0$ if and only if
$\bm{j}=\bm{j}^{\prime}$. The same is true for
$\mathcal{D}^{*}_{H}(\bm{f},\bm{f}^{\prime})$. In addition, they are reduced
to the usual metric induced by the squared norm for quadratic cases, as shown
in Eq. (23) and Eq. (24).

The second version is the Bregman divergence, which is defined as

$\displaystyle\mathcal{D}[\bm{j}\|\bm{j}^{\prime}]$
$\displaystyle:=\Psi(\bm{j})-\Psi(\bm{j}^{\prime})-\langle\bm{j}-\bm{j}^{\prime},\partial_{\bm{j}}\Psi(\bm{j}^{\prime})\rangle$
$\displaystyle\mathcal{D}^{*}[\bm{f}^{\prime}\|\bm{f}]$
$\displaystyle:=\Psi^{*}(\bm{f})-\Psi^{*}(\bm{f}^{\prime})-\langle\partial_{\bm{f}}\Psi^{*}(\bm{f}^{\prime}),\bm{f}-\bm{f}^{\prime}\rangle.$

The Bregman divergence $\mathcal{D}[\bm{j}\|\bm{j}^{\prime}]$ is nonnegative
and strictly convex for $\bm{j}$111111Not necessarily for $\bm{j}^{\prime}$.
and also attains the minimum $0$ if and only if
$\bm{j}=\bm{j}^{\prime}$121212It should be noted that the symmetry property of
the dissipation function Eq. (22) is not required to derive these properties..
However, it is generally asymmetric
$\mathcal{D}[\bm{j};\bm{j}^{\prime}]\neq\mathcal{D}[\bm{j}^{\prime};\bm{j}]$.
If $\Psi(\bm{j})$ is quadratic (Eq. (23)),
$\mathcal{D}[\bm{j}\|\bm{j}^{\prime}]$ is reduced to
$\mathcal{D}[\bm{j}\|\bm{j}^{\prime}]=\frac{1}{2}\|\bm{j}-\bm{j}^{\prime}\|^{2}_{{M}}$
and $\mathcal{D}^{*}[\bm{f}\|\bm{f}^{\prime}]$ is reduced to
$\mathcal{D}^{*}[\bm{f}\|\bm{f}^{\prime}]=\frac{1}{2}\|\bm{f}-\bm{f}^{\prime}\|^{2}_{{M}^{*}}$.
It should be noted that for any flux-force pairs $(\bm{j},\bm{f})$ and
$(\bm{j}^{\prime},\bm{f}^{\prime})$, the divergences
$\mathcal{D}[\bm{j}\|\bm{j}^{\prime}]$ and
$\mathcal{D}^{*}[\bm{f}^{\prime}\|\bm{f}]$ are the same:

$\displaystyle\mathcal{D}[\bm{j}\|\bm{j}^{\prime}]=\mathcal{D}^{*}[\bm{f}^{\prime}\|\bm{f}].$
(30)

Thus, $\mathcal{D}[\bm{j}\|\bm{j}^{\prime}]$ and
$\mathcal{D}^{*}[\bm{f}^{\prime}\|\bm{f}]$ are different representations of
the same quantity in the flux space $\mathcal{J}$ and the force space
$\mathcal{F}$, respectively. As a result,
$\mathcal{D}[\bm{j}\|\bm{j}^{\prime}]$ and
$\mathcal{D}^{*}[\bm{f}^{\prime}\|\bm{f}]$ are also the same as the following
mixed representation:

$\displaystyle\mathcal{D}[\bm{j};\bm{f}^{\prime}]:=\Psi(\bm{j})+\Psi^{*}(\bm{f}^{\prime})-\langle\bm{j},\bm{f}^{\prime}\rangle,$
(31)

i.e.,
$\mathcal{D}[\bm{j}\|\bm{j}^{\prime}]=\mathcal{D}[\bm{j};\bm{f}^{\prime}]=\mathcal{D}^{*}[\bm{f}^{\prime}\|\bm{f}]$
(see also Appendix ).

###  IV.3 Generalized Hilbert orthogonality

The two generalized distances naturally lead to two generalized notions of
orthogonality.

For any two force vectors $\bm{f}_{S},\bm{f}_{A}\in\mathcal{F}$, the first
orthogonality is defined to satisfy

$\displaystyle\Psi^{*}(\bm{f}_{S}+\bm{f}_{A})=\Psi^{*}(\bm{f}_{S}-\bm{f}_{A}).$
(32)

This definition comes from the fact that, if $\bm{f}_{S}$ and $\bm{f}_{A}$ are
orthogonal in an inner product space, then
$\|\bm{f}_{S}+\bm{f}_{A}\|_{{M}^{*}}=\|\bm{f}_{S}-\bm{f}_{A}\|_{{M}^{*}}$
holds. We call this orthogonality a generalized Hilbert orthogonality or
simply Hilbert orthogonality [39]. Then an orthogonal decomposition of a given
force $\bm{f}$ into $\bm{f}=\bm{f}_{S}+\bm{f}_{A}$ is obtained by finding an
iso-dissipation force $\bm{f}_{iso}$ satisfying
$\Psi^{*}(\bm{f})=\Psi^{*}(\bm{f}_{iso})$ and by computing $\bm{f}_{S}$ and
$\bm{f}_{A}$ as

$\displaystyle\bm{f}_{S}$ $\displaystyle=\frac{\bm{f}+\bm{f}_{iso}}{2},$
$\displaystyle\bm{f}_{A}=\frac{\bm{f}-\bm{f}_{iso}}{2}.$ (33)

Note that $\bm{f}_{iso}=\bm{f}_{S}-\bm{f}_{A}$ and also that we can have
different decompositions of $\bm{f}$ by choosing different iso-dissipation
forces $\bm{f}^{\prime}_{iso}$. In addition, from the symmetry of $\Psi^{*}$,
we have
$\Psi^{*}(\bm{f})=\Psi^{*}(\bm{f}_{iso})=\Psi^{*}(-\bm{f}_{iso})=\Psi^{*}(-\bm{f})$.
As a consequence of the orthogonality, we obtain

$\displaystyle\mathcal{D}[\bm{j}\|\bm{j}_{iso}]=$ $\displaystyle
2\langle\bm{j},\bm{f}_{A}\rangle\geq 0,$ (34)

where we use $\mathcal{D}[\bm{j}\|\bm{j}_{iso}]\geq 0$ (see Appendix for a
proof). Similarly, by considering the symmetry of $\Psi^{*}(\bm{f})$, we also
have

$\displaystyle\mathcal{D}_{\bm{\omega}}[\bm{j}\|-\bm{j}_{iso}]=2\langle\bm{j},\bm{f}_{S}\rangle\geq
0.$ (35)

From these relations, we have a decomposition of $\langle\bm{j},\bm{f}\rangle$
into two nonnegative terms:

$\displaystyle\langle\bm{j},\bm{f}\rangle$
$\displaystyle=\langle\bm{j},\bm{f}_{S}\rangle+\langle\bm{j},\bm{f}_{A}\rangle$
(36)
$\displaystyle=\frac{\mathcal{D}[\bm{j}\|-\bm{j}_{iso}]}{2}+\frac{\mathcal{D}[\bm{j}\|\bm{j}_{iso}]}{2}.$
(37)

Because $\langle\bm{j},\bm{f}\rangle$ is proportional to the EPR (Eq. (10)),
this decomposition provides a way to obtain nonnegative decompositions of the
EPR. It should be noted that we have infinitely many decompositions of this
type. In addition, by choosing $\bm{j}_{iso}=\bm{j}$, the decomposition
reduces to the fluctuation relation
$\dot{\Sigma}=\mathcal{D}[\bm{j}\|-\bm{j}]$, where we used
$\bm{j}=\bm{j}_{iso}$. This orthogonality has been employed in the MFT [32,
33, 34] to characterize the quasi-potential and the gradient flow aspects of
equilibrium and nonequilibrium systems as we see in the next section. The iso-
dissipation hypersurface is known as a level surface of the convex function in
Hessian geometry [15] and as a central affine surface in affine differential
geometry [40]. In both cases, it works as a central geometric object.
Moreover, this central affine surface also plays an important role in isobaric
thermodynamics, where volume can change in conjunction with reactions [38].
The Hilbert orthogonality is also defined in the flux space using
$\Psi(\bm{j})$ similarly to
$\Psi(\bm{j}_{A}+\bm{j}_{S})=\Psi(\bm{j}_{A}-\bm{j}_{S})$ while we do not use
it in this work. It should be noted that the orthogonality of $\bm{f}_{A}$ and
$\bm{f}_{S}$ generally does not mean the orthogonality of the corresponding
$\bm{j}_{A}$ and $\bm{j}_{S}$.

###  IV.4 Information geometric orthogonality

The second orthogonality comes from the information geometry [15, 16, 17]. For
three flux vectors
$\bm{j},\bm{j}^{\prime},\bm{j}^{\prime\prime}\in\mathcal{J}$ satisfying
$\bm{j}=\bm{j}^{\prime}+\bm{j}^{\prime\prime}$, the Bregman divergence
$\mathcal{D}[\bm{j}\|\bm{0}]$ can be decomposed as

$\displaystyle\mathcal{D}[\bm{j}\|\bm{0}]=\mathcal{D}[\bm{j}\|\bm{j}^{\prime\prime}]+\mathcal{D}[\bm{j}^{\prime\prime}\|\bm{0}]+\langle\bm{j},\bm{f}^{\prime\prime}\rangle,$
(38)

where $\bm{f}^{\prime\prime}$ is the Legendre dual of $\bm{j}^{\prime\prime}$.
The second orthogonality of $\bm{j}^{\prime}$ and $\bm{j}^{\prime\prime}$ is
defined by $\langle\bm{j}^{\prime},\bm{f}^{\prime\prime}\rangle=0$. Then the
following generalized Pythagorean theorem (GPT) holds:

$\displaystyle\mathcal{D}[\bm{j}\|\bm{0}]=\mathcal{D}[\bm{j}\|\bm{j}^{\prime\prime}]+\mathcal{D}[\bm{j}^{\prime\prime}\|\bm{0}].$
(39)

For quadratic dissipation functions, the GPT is reduced to the conventional
Pythagorean theorem:

$\displaystyle\|\bm{j}\|^{2}_{{M}}=\|\bm{j}^{\prime}\|^{2}_{{M}}+\|\bm{j}^{\prime\prime}\|^{2}_{{M}}.$
(40)

Thus, Eq. (39) is a generalization of the usual orthogonality. We call this
orthogonality the information geometric orthogonality. Similarly, for
$\bm{f}=\bm{f}^{\dagger}+\bm{f}^{\ddagger}$, the dual orthogonality between
$\bm{f}^{\dagger},\bm{f}^{\ddagger}\in\mathcal{F}$ is also defined by
$\langle\bm{j}^{\dagger},\bm{f}^{\ddagger}\rangle=0$, leading to the dual
version of GPT:

$\displaystyle\mathcal{D}^{*}[\bm{f}\|\bm{0}]=\mathcal{D}^{*}[\bm{f}\|\bm{f}^{\dagger}]+\mathcal{D}^{*}[\bm{f}^{\dagger}\|\bm{0}].$
(41)

The Bregman divergence and the information geometric orthogonality are the
central geometric tools in information geometry for analyzing manifolds of
statistical models and for conducting statistically meaningful projections
onto submanifolds [16, 17]. In addition, they have also been recently employed
in the thermodynamics of diffusion processes, MJP, and CRN [41, 42, 43, 37,
38, 38]. However, these works investigate the information and Hessian
geometric structures of state space, i.e., the space of the probability vector
$\bm{p}$ or the concentration vector $\bm{x}$, in which the thermodynamic
potential functions such as the Gibbs free energy work as the convex potential
function inducing the Hessian structure. As a result, the Hessian structure of
the state space captures the energetic or equilibrium aspects of the systems.
In contrast, the Hessian structure of the flux-force space in this work
captures the kinetic and nonequilibrium aspects of the systems.

It should be noted that the Hilbert and information geometric orthogonalities
are reduced to usual orthogonality when the dissipation functions are
quadratic (see also Appendix ).

##  V Chemical Reaction Dynamics and Entropy Production Decompositions

In this section, we demonstrate how the two generalizations of orthogonality
enable us to extend the different entropy production decompositions to CRN
(and MJP). From the explicit functional forms of the one-way fluxes under mass
action kinetics (Eq. (8)), the total flux, force, and activity defined in Eq.
(9) and Eq. (11) are expressed as

$\displaystyle\begin{aligned}
\bm{j}_{0}(\bm{x};\bm{k}^{\pm})&=(\bm{k}^{+}\circ B_{+}^{T}-\bm{k}^{-}\circ
B_{-}^{T})\bm{x}^{\Gamma^{T}},\\\
\bm{f}_{0}(\bm{x};\bm{k}^{\pm})&=\frac{1}{2}\left[\ln\frac{\bm{k}^{+}}{\bm{k}^{-}}-S^{T}\ln\bm{x}\right],\\\
\bm{\omega}_{0}(\bm{x};\bm{k}^{\pm})&=2\sqrt{\bm{k}^{+}\circ\bm{k}^{-}}\circ\bm{x}^{[\Gamma(B^{+}+B^{-})]^{T}/2}.\end{aligned}$
(42)

We use the subscript $(.)_{0}$ to designate the flux, force, and activity
given by these particular forms. With these specific forms of force and
activity, the CRN dynamics is described as a generalized flow (Eq. (18) and
Fig. 3).

We further transform the kinetic parameters $\bm{k}^{\pm}$ into force parts
($\bm{K}$) and activity part ($\bm{\kappa}$) as

$\displaystyle\bm{k}^{+}=\bm{\kappa}\circ\bm{K}^{1/2},\qquad\bm{k}^{-}=\bm{\kappa}\circ\bm{K}^{-1/2},$
(43)

where $\bm{\kappa}:=\sqrt{\bm{k}^{+}\circ\bm{k}^{-}}$ and
$\bm{K}:=\bm{k}^{+}/\bm{k}^{-}$. Thus, $(\bm{\kappa},\bm{K})$ has the same
information as $\bm{k}^{\pm}$. Moreover, we can verify that the force and
activity are dependent only on $\bm{K}$ and $\bm{\kappa}$, respectively, i.e.,
$\bm{f}_{0}(\bm{x};\bm{k}^{\pm})=\bm{f}_{0}(\bm{x};\bm{K})$ and
$\bm{\omega}_{0}(\bm{x};\bm{k}^{\pm})=\bm{\omega}_{0}(\bm{x};\bm{\kappa})$.

Figure 3: Schematic representation of the flux-force relationship for CRN with
mass action kinetics and induced CRE (Eq. (18)). Depending on the current
state $\bm{x}(t)$, (a), the force $\bm{f}_{0}(\bm{x}_{t})$ is induced (b). The
force is mapped to the corresponding flux $\bm{j}_{0}(\bm{x}_{t})$ by the
Legendre transformation
$\partial\Psi^{*}_{\bm{\omega}_{0}(\bm{x}_{t})}[\bm{f}_{0}(\bm{x}_{t})]$ (c).
The flux induces the change in $\bm{x}_{t}$ via the continuity equation (d).

###  V.1 Equilibrium dynamics

First, we consider the case in which the dynamics is equilibrium. For CNR, the
equilibrium states are defined as the set
$\mathcal{M}_{eq}(\bm{K},\bm{\kappa})$ satisfying the detailed balance
condition (DBC)

$\displaystyle\mathcal{M}_{eq}(\bm{K},\bm{\kappa}):=\\{\bm{x}|\bm{j}_{0}(\bm{x};\bm{K},\bm{\kappa})=0\\}.$
(44)

For a parameter set $(\bm{K},\bm{\kappa})$ that admits the existence of
equilibrium states, i.e., $\mathcal{M}_{eq}(\bm{K},\bm{\kappa})\neq\emptyset$,
the CRN becomes an equilibrium CRN. The condition
$\mathcal{M}_{eq}(\bm{K},\bm{\kappa})\neq\emptyset$ is satisfied when the
parameter $\bm{K}$ satisfies the Wegscheider equilibrium (EQ) condition [21,
2, 37]:

$\displaystyle\ln\bm{K}\in\mathrm{Im}S^{T}.$ (45)

We denote such $\bm{K}$ as $\bm{K}_{eq}$. Note that the EQ condition is
independent of the activity parameter $\bm{\kappa}$. Then, for each initial
condition $\bm{x}_{0}>\bm{0}$, a unique equilibrium state $\bm{x}_{eq}$
exists, which is determined by the intersection of
$\mathcal{M}_{eq}(\bm{K}_{eq})$ and the stoichiometric compatibility class
[44, 37]

$\displaystyle\mathcal{P}(\bm{x}_{0}):=\\{\bm{x}|(\bm{x}-\bm{x}_{0})\in\mathrm{Im}[S]\\}$
(46)

as $\bm{x}_{eq}=\mathcal{M}_{eq}(\bm{K}_{eq})\cap\mathcal{P}(\bm{x}_{0})$
131313The uniqueness of the intersection is a consequence of Hessian geometric
structure on the concentration-chemical potential space.  [44, 37]. Moreover,
from the DBC, we have $\ln\bm{K}_{eq}=S^{T}\ln\bm{x}_{eq}$. Then, the force
$\bm{f}_{0}(\bm{x})$ of the equilibrium dynamics is represented as

$\displaystyle\bm{f}_{0}(\bm{x})$
$\displaystyle=-\frac{1}{2}S^{T}\ln\left(\frac{\bm{x}}{\bm{x}_{eq}}\right)=-\frac{1}{2}S^{T}\partial_{\bm{x}}\varphi_{eq}(\bm{x}),$
(47)

where
$\varphi_{eq}(\bm{x})=\mathcal{D}_{KL}[\bm{x}\|\bm{x}_{eq}]:=\bm{x}^{T}\ln\frac{\bm{x}}{\bm{x}_{eq}}-\bm{1}^{T}(\bm{x}-\bm{x}_{eq})$
is the Gibbs free energy of the equilibrium CRN. The matrix $S^{T}$ can be
considered as the discrete version of the gradient $\mathrm{grad}=S^{T}$ for
the chemical hypergraph because $S^{T}$ is the adjoint operator of
$S=-\mathrm{div}$. Thus, the equilibrium force is a gradient of
$\varphi_{eq}(\bm{x})$:

$\displaystyle\bm{f}_{0}(\bm{x})$
$\displaystyle=-\frac{1}{2}\mathrm{grad}\left[\partial_{\bm{x}}\varphi_{eq}(\bm{x})\right].$
(48)

Additionally, the Wegscheider condition (Eq. (45)) is interpreted as
$\ln\bm{K}_{eq}\in\mathrm{Im}[\mathrm{grad}]$. Then the dynamics (Eq. (18))
under the EQ parameter condition becomes a generalized gradient flow (GF) of
$\varphi_{eq}(\bm{x})$ [45]:

$\displaystyle\dot{\bm{x}}=S\bm{j}_{0}(\bm{x})$
$\displaystyle=S\partial\Psi^{*}_{\bm{x}}\left[-\frac{1}{2}S^{T}\partial_{\bm{x}}\varphi_{eq}(\bm{x})\right]$
(49)
$\displaystyle=-\mathrm{div}\,\partial\Psi^{*}_{\bm{x}}\left[-\frac{1}{2}\mathrm{grad}\,\left[\partial_{\bm{x}}\varphi_{eq}(\bm{x})\right]\right],$
(50)

where $\Psi^{*}_{\bm{x}}:=\Psi^{*}_{\bm{\omega}_{0}(\bm{x})}$. Because
$\bm{x}(t)$ is the gradient flow of $\varphi_{eq}(\bm{x}(t))$,
$\mathrm{d}\varphi_{eq}(\bm{x}(t))/\mathrm{d}t$ is always non-positive as

$\displaystyle-\frac{\mathrm{d}\varphi_{eq}(\bm{x})}{\mathrm{d}t}$
$\displaystyle=\langle\dot{\bm{x}},-\partial_{\bm{x}}\varphi_{eq}(\bm{x})\rangle=\langle
S\bm{j}_{0}(\bm{x}),-\partial_{\bm{x}}\varphi_{eq}(\bm{x})\rangle$
$\displaystyle=\langle\bm{j}_{0}(\bm{x}),-S^{T}\partial_{\bm{x}}\varphi_{eq}(\bm{x})\rangle=2\langle\bm{j}_{0}(\bm{x}),\bm{f}_{0}(\bm{x})\rangle$
$\displaystyle=2\left[\Psi_{\bm{x}}(\bm{j}_{0}(\bm{x}))+\Psi_{\bm{x}}^{*}(\bm{f}_{0}(\bm{x}))\right]=\dot{\Sigma}\geq
0,$ (51)

where we use Eq. (15) and Eq. (10). By integration, we have the relation
between the change of free energy and the entropy production:

$\displaystyle\varphi_{eq}(\bm{x}(t))-\varphi_{eq}(\bm{x}(0))=\int_{t^{\prime}=0}^{t}\dot{\Sigma}_{t^{\prime}}\mathrm{d}t^{\prime}.$
(52)

###  V.2 Complex balanced dynamics, Hilbert orthogonality, and Hatano-Sasa
decomposition

Figure 4: Schematic illustration of the Hilbert orthogonality between
$\bm{f}_{S}(\bm{x}_{t})$ (blue arrows) and $\bm{f}_{A}$(magenta arrows). The
green solid curve represents the iso-dissipation hypersurface:
$\Psi^{*}_{\bm{x}}(\bm{f})=\Psi^{*}_{\bm{x}}(\bm{f}_{0}(\bm{x}_{t}))=\mathrm{const.}$.
The black circle, red square, and green star in $\mathcal{F}$ are
$\bm{f}_{0}(\bm{x}_{t})$, $\bm{f}_{S}(\bm{x}_{t})$, and
$\bm{f}_{iso}(\bm{x}_{t})$, respectively. The black circle, red square, and
green star in $\mathcal{J}$ are their Legendre transform, i.e.,
$\bm{j}_{0}(\bm{x}_{t})$, $\bm{j}_{S}(\bm{x}_{t})$, and
$\bm{j}_{iso}(\bm{x}_{t})$, respectively. The light blue lines in
$\mathcal{F}$ are the NEF subspaces, $\mathcal{P}^{f}(\bm{f}_{A})$(solid line)
and $\mathcal{P}^{f}(-\bm{f}_{A})$ (dashed line), and the EF subspace
$\mathcal{P}^{f}_{eq}$ (dotted line). In $\mathcal{J}$, the corresponding NEF
manifolds, $\mathcal{M}^{f}_{\bm{x}}(\bm{f}_{A})$ (the light solid blue curve)
and $\mathcal{M}^{f}_{\bm{x}}(-\bm{f}_{A})$ (the light dashed blue curve), and
the EF manifold $\mathcal{M}^{f}_{\bm{x},eq}$ (the light dotted blue curve)
are also depicted.

Next, we consider the nonequilibrium complex balanced (CB) steady state, and
the EPR decomposition induced by the Hilbert orthogonality. The result here is
a reinterpretation of the result in [33, 32, 34] from the viewpoint of Hessian
geometry. Nevertheless, it is worth including because the decomposition is
contrasted with the information geometric decomposition in the next
subsection.

CB states of CRN are defined as the set satisfying

$\displaystyle\mathcal{M}_{cb}(\bm{K},\bm{\kappa}):=\\{\bm{x}|B\bm{j}_{0}(\bm{x};\bm{K},\bm{\kappa})=\bm{0}\\}.$
(53)

A parameter set $(\bm{K}_{cb},\bm{\kappa}_{cb})$ induces CB dynamics if the CB
condition $\mathcal{M}_{cb}(\bm{K}_{cb},\bm{\kappa}_{cb})\neq\emptyset$ holds.
Then, for each initial condition $\bm{x}_{0}$, the CB steady state
$\bm{x}_{cb}$ is uniquely determined by the intersection of
$\mathcal{M}_{cb}(\bm{K}_{cb},\bm{\kappa}_{cb})$ and the stoichiometric
compatibility class $\mathcal{P}(\bm{x}_{0})$ as
$\bm{x}_{cb}=\mathcal{M}_{cb}(\bm{K}_{cb},\bm{\kappa}_{cb})\cap\mathcal{P}(\bm{x}_{0})$
[44, 37].

The CB state inherits important properties of the equilibrium state such as
the uniqueness, global stability, and the gradient-flow-like aspect [44, 20,
21, 37]. By definition, an EQ state is a CB state, but not vice versa.
Additionally, the steady states of MJP always satisfy the CB condition because
$B\bm{j}_{0}(\bm{x};\bm{K},\bm{\kappa})=\bm{0}$ is nothing but the steady
state condition for MJP (Eq. (2)). Thus, MJP are unconditionally complex
balanced. In contrast, CRN are not always complex balanced depending on the
parameter value.

Using $\bm{x}_{cb}$, the force $\bm{f}_{0}(\bm{x})$ is represented as

$\displaystyle\bm{f}_{0}(\bm{x})$
$\displaystyle=-\frac{1}{2}S^{T}\partial_{\bm{x}}\varphi_{cb}(\bm{x})+\frac{1}{2}\ln\frac{\bm{K}_{cb}}{\tilde{\bm{K}}_{eq}},$
(54)

where $\tilde{\bm{K}}_{eq}:=S^{T}\ln\bm{x}_{cb}$ and
$\varphi_{cb}(\bm{x}):=\mathcal{D}_{KL}[\bm{x}\|\bm{x}_{cb}]$. Let
$\bm{f}_{S}(\bm{x}):=-\frac{1}{2}S^{T}\partial_{\bm{x}}\varphi_{cb}(\bm{x})$
and $\bm{f}_{A}:=\frac{1}{2}\ln\frac{\bm{K}_{cb}}{\tilde{\bm{K}}_{eq}}$. Then,
for all $\bm{x}\in\mathcal{X}$, the generalized Hilbert orthogonality holds
(Fig. 4):

$\displaystyle\Psi^{*}_{\bm{x}}(\bm{f}_{S}(\bm{x})+\bm{f}_{A})=\Psi^{*}_{\bm{x}}(\bm{f}_{S}(\bm{x})-\bm{f}_{A}),$
(55)

because of the CB condition $B\bm{j}_{0}(\bm{x}_{cb})=\bm{0}$ [32, 33, 34]
(see Appendix for a proof). Thus, Eq. (36) admits a nonnegative decomposition
of the EPR:

$\displaystyle\dot{\Sigma}(\bm{x})=\dot{\Sigma}_{ex}^{GF}(\bm{x})+\dot{\Sigma}_{hk}^{GF}(\bm{x}),$
(56)

where
$\dot{\Sigma}_{ex}^{GF}(\bm{x}):=2\langle\bm{j}_{0}(\bm{x}),\bm{f}_{S}(\bm{x})\rangle$
and
$\dot{\Sigma}_{hk}^{GF}(\bm{x}):=2\langle\bm{j}_{0}(\bm{x}),\bm{f}_{A}\rangle$.
Owing to this nonnegative decomposition, $\varphi_{cb}(\bm{x})$ works as a
Lyapunov function (quasi-potential):

$\displaystyle-\frac{\mathrm{d}\varphi_{cb}(\bm{x}_{t})}{\mathrm{d}t}$
$\displaystyle=\dot{\Sigma}_{ex}^{GF}(\bm{x}_{t})=\mathcal{D}[\bm{j}_{0}(\bm{x}_{t})\|\bm{f}_{iso}(\bm{x}_{t})]\geq
0,$

and $\dot{\Sigma}_{ex}^{GF}(\bm{x})$ behaves like an excess EPR, i.e.,
$\lim_{t\to\infty}\dot{\Sigma}_{ex}^{GF}(\bm{x}_{t})=0$. Thus,
$\varphi_{cb}(\bm{x})$ as a Lyapunov function guarantees the stability of the
CB state $\bm{x}_{cb}$ and the gradient-flow-like property of the CB dynamics.

This decomposition is equivalent to those proposed for CRN in [20, 21] as a
generalization of the Hatano-Sasa decomposition [46]. Therefore, this result
is just a reinterpretation of the previous ones. However, owing to the
geometric perspective endowed by Hessian geometry, the non-negativity of
$\dot{\Sigma}_{ex}^{GF}(\bm{x})$ and $\dot{\Sigma}_{hk}^{GF}(\bm{x})$ becomes
transparent, the proof of which in previous works required more technical and
less transparent computations[20, 21].

###  V.3 Information geometric orthogonality and Maes-Netočnỳ decomposition

Figure 5: Schematic illustration of the information geometric orthogonalities.
(Left) The orthogonality between $\bm{j}_{st}(\bm{x}_{t})$ and
$\bm{j}_{0}(\bm{x}_{t})-\bm{j}_{st}(\bm{x}_{t})$ (blue arrows) and the dual
orthogonality between $\bm{j}_{eq}(\bm{x}_{t})$ and
$\bm{j}_{0}(\bm{x}_{t})-\bm{j}_{eq}(\bm{x}_{t})$ (red arrows) in
$\mathcal{J}$. The blue planes are the iso-velocity subspace
$\mathcal{P}^{v}(\bm{j}_{0}(\bm{x}_{t}))$ (the upper plane) and the zero-
velocity subspace $\mathcal{P}^{v}_{st}$ (the lower plane). The grey curves
are the NEF manifold $\mathcal{M}^{f}_{\bm{x}}(\bm{f}_{0}(\bm{x}_{t}))$ and
the EF manifold $\mathcal{M}^{f}_{\bm{x},eq}$. (Right) The same
orthogonalities shown in $\mathcal{F}$ space. The orthogonality between
$\bm{f}_{st}(\bm{x}_{t})$ and $\bm{f}_{0}(\bm{x}_{t})-\bm{f}_{st}(\bm{x}_{t})$
(blue arrows) and the dual orthogonality between $\bm{f}_{eq}(\bm{x}_{t})$ and
$\bm{f}_{0}(\bm{x}_{t})-\bm{f}_{eq}(\bm{x}_{t})$ (red arrows) in
$\mathcal{F}$. The blue surfaces are the iso-velocity manifold
$\mathcal{M}^{v}_{\bm{x}_{t}}(\bm{j}_{0}(\bm{x}_{t}))$ (the upper surface) and
the zero-velocity manifold $\mathcal{M}^{v}_{\bm{x}_{t},st}$ (the lower
surface). The grey lines are the NEF subspaces
$\mathcal{P}^{f}(\bm{f}_{0}(\bm{x}_{t}))$ and the EF subspace
$\mathcal{P}^{f}_{\bm{x},eq}$.

While the gradient flow aspect and associated EPR decomposition by the quasi-
potential have been investigated in CRN theory, MFT, and applied mathematics
[20, 21, 33, 44], the decomposition does not capture other important
thermodynamic aspects such as minimum dissipation or minimum entropy
production principles. In addition, the Hilbert orthogonality for
$\varphi_{cb}(\bm{x})$ (Eq. (55)) holds only when the CB condition is
satisfied. In CRN, we can have non-complex-balanced steady states depending on
the parameter values.

As another type of the EPR decomposition, we have Maes-Netočnỳ decomposition,
which generalizes the Komatsu-Nakagawa-Sasa-Tasaki approach [47, 48].
Moreover, its geometric interpretation for overdamped diffusion processes was
recently provided in terms of the formal Riemannian geometry of the flux-force
space [25]. Its extension to MJP and CRN has yet to be achieved because of
their nonquadratic dissipation functions. We show that the information
geometric orthogonality is central to this extension.

In the Maes-Netočnỳ decomposition, for a given $\bm{x}_{t}$141414$\bm{x}_{t}$
is not necessarily a steady state here., the non-equilibrium flux
$\bm{j}_{0}(\bm{x})$ is decomposed into the quasi-steady flux
$\bm{j}_{st}(\bm{x})$ and the remaining part as
$\bm{j}_{0}(\bm{x})=\bm{j}_{st}(\bm{x})+(\bm{j}_{0}(\bm{x})-\bm{j}_{st}(\bm{x}))$.
The quasi-steady flux $\bm{j}_{st}(\bm{x})$ is the flux that makes a given
state $\bm{x}$ steady obtained by modulating the conserved force, i.e., being
generated by the gradient of potential force in the case of overdamped
diffusion. This means that $\bm{j}_{st}(\bm{x})\in\mathrm{Ker}[\mathrm{div}]$
and $\bm{f}_{0}(\bm{x})-\bm{f}_{st}(\bm{x})\in\mathrm{Im}[\mathrm{grad}]$
hold. Furthermore, $\bm{j}_{st}(\bm{x})$ is also characterized as the minimum
EPR flux under the constraint of
$\bm{f}_{0}(\bm{x})-\bm{f}_{st}(\bm{x})\in\mathrm{Im}[\mathrm{grad}]$ [47].
From the viewpoint of vector calculus, the decomposition is a special case of
the Helmholtz-Hodge decomposition[49]. With the the Riemannian inner product
structure between $\mathcal{J}$ and $\mathcal{F}$ for overdamped diffusion, we
can also regard $\bm{j}_{st}(\bm{x})$ as the orthogonal projection of
$\bm{j}_{0}(\bm{x})$ onto $\mathrm{Ker}[\mathrm{div}]$ along
$\mathrm{Im}[\mathrm{grad}]$. Because the Maes-Netočnỳ decomposition is
tightly linked to the orthogonality in the inner product space, it is not
trivial how to extend the decomposition to other systems such as CRN (and MJP)
with nonquadratic dissipation functions. We show how information-geometric
orthogonality can resolve this problem.

First, we define a linear affine subspace $\mathcal{P}^{v}(\bm{j})$ as

$\displaystyle\mathcal{P}^{v}(\bm{j})=\left\\{\bm{j}^{\prime}|\bm{j}-\bm{j}^{\prime}\in\mathrm{Ker}S\right\\}\subset\mathcal{J}.$
(57)

Because $\bm{j}^{\prime}\in\mathcal{P}^{v}(\bm{j})$ satisfies
$S\bm{j}=S\bm{j}^{\prime}$, $\mathcal{P}^{v}(\bm{j})$ is the set of the fluxes
that induces the same instantaneous velocity $\dot{\bm{x}}$ as $\bm{j}$. Thus,
we call $\mathcal{P}^{v}(\bm{j})$ an iso-velocity subspace in $\mathcal{J}$.
Because $\dot{\bm{x}}=0$ for the steady state flux, we also denote
$\mathcal{P}^{v}_{st}:=\mathcal{P}^{v}(\bm{0})$ as the steady state (zero-
velocity) subspace.

As the complementary subspace of $\mathcal{P}^{v}(\bm{j})$, we also define a
linear affine subspace $\mathcal{P}^{f}(\bm{f})\subset\mathcal{F}$ as

$\displaystyle\mathcal{P}^{f}(\bm{f}):=\left\\{\bm{f}^{\prime}|\bm{f}-\bm{f}^{\prime}\in\mathrm{Im}S^{T}\right\\}.$
(58)

Because $S^{T}$ is the gradient of a CRN, i.e., $\mathrm{grad}=S^{T}$,
$\mathcal{P}^{f}(\bm{f})$ is the subspace obtained by shifting
$\mathrm{Im}[\mathrm{grad}]=\mathrm{Im}[S^{T}]$ by $\bm{f}$. We call
$\mathcal{P}^{f}(\bm{f})$ an nonequilibrium force (NEF) subspace in
$\mathcal{F}$ if $\bm{f}\notin\mathrm{Im}S^{T}$ and the equilibrium force (EF)
subspace if $\bm{f}\in\mathrm{Im}S^{T}$ because, if the equilibrium condition
is satisfied, $\bm{f}_{0}(\bm{x})$ is always in $\mathrm{Im}[S^{T}]$ as in Eq.
(48). We denote the EF subspace by
$\mathcal{P}^{f}_{eq}:=\mathcal{P}^{f}(\bm{0})$.

To obtain the generalized version of $\bm{j}_{st}(\bm{x})$, we transform
$\mathcal{P}^{f}(\bm{f})$ and $\mathcal{P}^{f}_{eq}$ to $\mathcal{J}$ via the
Legendre transformation:

$\displaystyle\mathcal{M}^{f}_{\bm{x}}(\bm{f})$
$\displaystyle:=\partial\Psi^{*}_{\bm{x}}\left[\mathcal{P}^{f}(\bm{f})\right]\subset\mathcal{J},$
(59) $\displaystyle\mathcal{M}^{f}_{\bm{x},eq}$
$\displaystyle:=\partial\Psi^{*}_{\bm{x}}\left[\mathcal{P}^{f}_{eq}\right]\subset\mathcal{J}.$
(60)

It should be noted that $\mathcal{M}^{f}_{\bm{x}}(\bm{f})$ depends on the
current state $\bm{x}$ of the system even though $\mathcal{P}^{f}(\bm{f})$
does not because the Legendre transformation is dependent on $\bm{x}$ via the
activity $\bm{\omega}(\bm{x})$. From the definition,
$\mathcal{M}^{f}_{\bm{x}}(\bm{f})$ and $\mathcal{M}^{f}_{\bm{x},eq}$ are the
set of nonequilibrium and equilibrium fluxes obtained by modulating the
gradient force.

Then the generalized quasi-steady flux $\bm{j}_{st}(\bm{x})$ is obtained by
the intersection of $\mathcal{P}^{v}_{st}$ and
$\mathcal{M}^{f}_{\bm{x}}(\bm{f}_{0}(\bm{x}))$:

$\displaystyle\bm{j}_{st}(\bm{x})$
$\displaystyle:=\mathcal{P}^{v}_{st}\cap\mathcal{M}^{f}_{\bm{x}}(\bm{f}_{0}(\bm{x})).$
(61)

For the case of overdamped diffusion processes, the Legendre transformation of
the corresponding quadratic dissipation function is linear and, therefore,
$\mathcal{M}^{f}_{\bm{x}}(\bm{f}_{0}(\bm{x}))$ becomes a flat subspace. The
existence and uniqueness of $\bm{j}_{st}(\bm{x})$ follow from the linear
algebra. However, for the generalized $\bm{j}_{st}(\bm{x})$, the existence and
uniqueness of $\bm{j}_{st}(\bm{x})$ are not guaranteed a priori, because
$\mathcal{M}^{f}_{\bm{x}}(\bm{f}_{0}(\bm{x}))$ is a curved manifold as shown
in Fig. 5 (A). Nonetheless, $\bm{j}_{st}(\bm{x})$ in Eq. (61) exists uniquely.
More generally, for any $\bm{j}^{\prime}$ and $\bm{f}^{\prime\prime}$, the
following intersection is unique and transversal (see Appendix for proof):

$\displaystyle\bm{j}_{int}(\bm{j}^{\prime},\bm{f}^{\prime\prime})=\mathcal{P}^{v}(\bm{j}^{\prime})\cap\mathcal{M}^{f}_{\bm{x}}(\bm{f}^{\prime\prime}).$
(62)

This result generalizes the uniqueness of the intersection of complementary
subspaces under an inner product structure. This is one of the notable
properties of Hessian geometry. In information geometry and Hessian geometry,
$\mathcal{P}^{v}(\bm{j})$ and $\mathcal{M}^{f}_{\bm{x}}(\bm{f})$ are called
dual flat spaces [16, 17, 15].

By virtue of the uniqueness, we obtain a unique decomposition of the flux
$\bm{j}_{0}(\bm{x})$:

$\displaystyle\bm{j}_{0}(\bm{x})=\bm{j}_{st}(\bm{x})+(\bm{j}_{0}(\bm{x})-\bm{j}_{st}(\bm{x})).$
(63)

Because $\bm{j}_{st}(\bm{x})\in\mathrm{Ker}S$ and
$\bm{j}_{st}(\bm{x})\in\mathcal{M}^{f}_{\bm{x}}(\bm{f}_{0}(\bm{x}))$ where
$\mathcal{M}^{f}_{\bm{x}}(\bm{f}_{0}(\bm{x}))=\partial_{\bm{j}}\Psi^{*}_{\bm{x}}[\bm{f}_{0}(\bm{x})+\mathrm{Im}[S^{T}]]$,
$\bm{j}_{st}(\bm{x})$ is the steady flux that is obtained by modulating the
nonequilibrium force $\bm{f}_{0}(\bm{x})$ by adding the gradient equilibrium
force. We also have
$\bm{f}_{0}(\bm{x})-\bm{f}_{st}(\bm{x})\in\mathrm{Im}S^{T}$. Thus, Eq. (63) is
a generalization of Maes-Netočnỳ decomposition and also of Helmholtz
decomposition. Moreover, from
$\langle\bm{j}_{st}(\bm{x}),\bm{f}_{0}(\bm{x})-\bm{f}_{st}(\bm{x})\rangle=0$,
we have the generalized Pythagorean relation:

$\displaystyle\mathcal{D}_{\bm{x}}^{*}[\bm{f}_{0}(\bm{x})\|\bm{0}]$
$\displaystyle=\mathcal{D}_{\bm{x}}^{*}[\bm{f}_{0}(\bm{x})\|\bm{f}_{st}(\bm{x})]+\mathcal{D}_{\bm{x}}^{*}[\bm{f}_{st}(\bm{x})\|\bm{0}],$
(64)

which is further reduced to a decomposition of the dual dissipation function:

$\displaystyle\Psi_{\bm{x}}^{*}(\bm{f}_{0}(\bm{x}))$
$\displaystyle=\mathcal{D}_{\bm{x}}^{*}[\bm{f}_{0}(\bm{x})\|\bm{f}_{st}(\bm{x})]+\Psi_{\bm{x}}^{*}(\bm{f}_{st}(\bm{x})).$
(65)

Then, $\bm{f}_{st}$ is characterized variationally as (see Appendix ):

$\displaystyle\bm{f}_{st}(\bm{x})$
$\displaystyle=\arg\min_{\bm{f}\in\mathcal{P}^{f}(\bm{f}_{0}(\bm{x}))}\Psi_{\bm{x}}^{*}(\bm{f}).$
(66)

This variational formula indicates that $\bm{f}_{st}(\bm{x})$ is the minimum
dissipation force that can be obtained by tilting the nonequilibrium force
$\bm{f}_{0}(\bm{x})$ with the equilibrium gradient force. Thus,
$\bm{j}_{st}(\bm{x})$ is the minimum dissipation (MD) steady state flux. For
quadratic cases where the equivalence of the dissipation functions and the EPR
holds (Eq. (63)), the result reduces to the minimum entropy production
principle.

However, for nonquadratic cases, Eq. (65) alone does not provide a relevant
EPR decomposition because EPR is not determined only by
$\Psi_{\bm{x}}^{*}(\bm{f}_{0}(\bm{x}))$. To resolve the problem, we consider
the dual decomposition of $\bm{j}_{0}$ as shown in Fig. 5 (A):

$\displaystyle\bm{j}_{0}(\bm{x})=\bm{j}_{eq}(\bm{x})+(\bm{j}(\bm{x})-\bm{j}_{eq}(\bm{x})),$
(67)

where we define $\bm{j}_{eq}(\bm{x})$ as

$\displaystyle\bm{j}_{eq}(\bm{x}):=\mathcal{P}^{v}(\bm{j}_{0}(\bm{x}))\cap\mathcal{M}^{f}_{\bm{x},eq}.$
(68)

Because $\bm{j}_{eq}(\bm{x})\in\mathcal{P}^{v}(\bm{j}_{0}(\bm{x}))$, the flux
$\bm{j}_{eq}(\bm{x})$ induces the same instantaneous velocity $\dot{\bm{x}}$
as $\bm{j}_{0}(\bm{x})$. Moreover, because
$\bm{j}_{eq}(\bm{x})\in\mathcal{M}^{f}_{\bm{x},eq}$, it is the flux induced by
an equilibrium force, i.e., a pure gradient force. Thus, $\bm{j}_{eq}(\bm{x})$
is the equilibrium flux that induces the same dynamics as
$\bm{j}_{0}(\bm{x})$. It should be noted that $\bm{j}_{0}(\bm{x})$ is
generally induced by a nonequilibrium force. In addition, we have
$\bm{f}_{eq}(\bm{x})\in\mathrm{Im}S^{T}$ and
$\bm{j}(\bm{x})-\bm{j}_{eq}(\bm{x})\in\mathrm{Ker}S$ from
$\bm{j}_{eq}(\bm{x})\in\mathcal{M}^{f}_{\bm{x},eq}$ and
$\bm{j}_{eq}(\bm{x})\in\mathcal{P}^{v}(\bm{j}_{0}(\bm{x}))$, respectively.
Thus, Eq. (67) is another generalization of Helmholtz-Hodge decomposition of
the flux $\bm{j}_{0}(\bm{x})$. Moreover, owing to
$\langle\bm{j}(\bm{x})-\bm{j}_{eq}(\bm{x}),\bm{f}_{eq}(\bm{x})\rangle=0$, we
have the Pythagorean relation:

$\displaystyle\mathcal{D}_{\bm{x}}[\bm{j}_{0}(\bm{x})\|\bm{0}]$
$\displaystyle=\mathcal{D}_{\bm{x}}[\bm{j}_{0}(\bm{x})\|\bm{j}_{eq}(\bm{x})]+\mathcal{D}_{\bm{x}}[\bm{j}_{eq}(\bm{x});\bm{0}],$
(69)

which is reduced to a decomposition of the primal dissipation function:

$\displaystyle\Psi_{\bm{x}}(\bm{j}_{0}(\bm{x}))$
$\displaystyle=\mathcal{D}_{\bm{x}}[\bm{j}_{0}(\bm{x})\|\bm{j}_{eq}(\bm{x})]+\Psi_{\bm{x}}(\bm{j}_{eq}(\bm{x})).$
(70)

From this, $\bm{j}_{eq}$ can be characterized variationally (see Appendix for
proof):

$\displaystyle\bm{j}_{eq}(\bm{x})$
$\displaystyle=\arg\min_{\bm{j}\in\mathcal{P}^{v}(\bm{j}_{0}(\bm{x}))}\Psi_{\bm{x}}(\bm{j}).$
(71)

This variational formula means that, among all fluxes that induce the same
velocity as $\bm{j}_{0}(\bm{x})$, the equilibrium flux $\bm{j}_{eq}(\bm{x})$
is the flux that minimizes the primal dissipation function. Thus,
$\bm{j}_{eq}(\bm{x})$ is the MD flux.

By combining Eq. (70) and Eq. (65), we have an EPR decomposition (see Appendix
for proof):

$\displaystyle\dot{\Sigma}(\bm{x})=\dot{\Sigma}_{hk}^{MD}(\bm{x})+\dot{\Sigma}_{ex}^{MD}(\bm{x})$
(72)

where we define housekeeping and excess EPR as

$\displaystyle\dot{\Sigma}_{hk}^{MD}(\bm{x})$
$\displaystyle:=2\left[\mathcal{D}_{\bm{x}}[\bm{j}_{0}(\bm{x})\|\bm{j}_{eq}(\bm{x})]+\Psi^{*}_{\bm{x}}(\bm{f}_{st}(\bm{x}))\right]\geq
0,$ $\displaystyle\dot{\Sigma}_{ex}^{MD}(\bm{x})$
$\displaystyle:=2\left[\mathcal{D}_{\bm{x}}^{*}[\bm{f}_{0}(\bm{x})\|\bm{f}_{st}(\bm{x})]+\Psi_{\bm{x}}(\bm{j}_{eq}(\bm{x}))\right]\geq
0.$

If the state $\bm{x}(t)$ converges to a steady state
$\bm{x}_{st}$151515$\bm{x}_{st}$ is not necessarily a CB state but can be a
general steady state, then

$\displaystyle\lim_{t\to\infty}\dot{\Sigma}_{hk}(\bm{x}(t))$
$\displaystyle=\dot{\Sigma}_{hk}(\bm{x}_{st})=2\langle\bm{j}_{0}(\bm{x}_{st}),\bm{f}_{0}(\bm{x}_{st})\rangle$
(73) $\displaystyle\lim_{t\to\infty}\dot{\Sigma}_{ex}(\bm{x}(t))$
$\displaystyle=\dot{\Sigma}_{ex}(\bm{x}_{st})=0,$ (74)

where we used $\bm{j}_{0}(\bm{x}_{st})=\bm{j}_{st}(\bm{x}_{st})$ and
$\bm{j}_{eq}(\bm{x}_{st})=\bm{0}$. Thus, this EPR decomposition geometrically
generalizes the Maes-Netočnỳ one to MJP and CRN and clarifies its dualistic
minimum dissipation principles. If the state $\bm{x}(t)$ converges to a CB
state $\bm{x}_{cb}$, then
$\lim_{t\to\infty}\dot{\Sigma}_{hk}(\bm{x}(t))=2\langle\bm{j}_{A},\bm{f}_{A}\rangle$
holds.

##  VI Numerical Demonstration

We numerically demonstrate the two decompositions in Eq. (56) and Eq. (72) and
the geometric relations among $\bm{j}_{0}(\bm{x})$, $\bm{j}_{iso}(\bm{x})$,
$\bm{j}_{S}(\bm{x})$, and $\bm{j}_{A}$ and those among $\bm{j}_{0}(\bm{x})$,
$\bm{j}_{st}(\bm{x})$, and $\bm{j}_{eq}(\bm{x})$ by using a CRN in Fig 1
(A)[50, 51].

Table 1: The parameter values used for simulations.  Type | $k_{1}^{+}$ | $k_{1}^{-}$ | $k_{2}^{+}$ | $k_{2}^{-}$ | $k_{3}^{+}$ | $k_{3}^{-}$
---|---|---|---|---|---|---
Equilibrium | $4$ | $1$ | $\frac{3}{\sqrt{2}}$ | $3\sqrt{2}$ | $2\sqrt{2}$ | $4\sqrt{2}$
Complex balanced | $\frac{1}{2}$ | $2$ | $4$ | $\frac{47}{4}$ | $\sqrt{2}$ | $\frac{15}{2}+2\sqrt{2}$
Non-complex balanced | $\frac{1}{2}$ | $2$ | $\frac{1}{17}$ | $\frac{85}{8}$ | $\frac{273}{68}$ | $\frac{137}{68}$
Figure 6: Heatmap plots of EPR $\dot{\Sigma}$ (left panels) and the two
dissipation functions, $\Psi_{\bm{x}}(\bm{j}_{0}(\bm{x})$ (center panels) and
$\Psi_{\bm{x}}^{*}(\bm{f}_{0}(\bm{x})$ (right panels) as functions of $\bm{x}$
for the equilibrium parameter set (A), the complex balanced parameter set (B),
and the non-complex balanced parameter set(C). The parameter sets are listed
in Tab. 1. In each panel, gray and black points are the initial and the steady
states, respectively. The white arrow is the trajectory of $\bm{x}(t)$. The
black line in each panel is the set of equilibrium states (A), complex-
balanced states (B), and non-complex balanced steady states (C), respectively.
Figure 7: The EPR and its decompositions along the trajectory $\bm{x}_{t}$ for
complex balanced (A,B) and non-complex balanced (C,D) parameter sets. The
decomposition of EPR into two dissipation functions (A, C). The decompositions
of EPR by Hilbert ($\dot{\Sigma}_{ex}^{GF}(\bm{x})$,
$\dot{\Sigma}_{hk}^{GF}(\bm{x})$) and information geometric
($\dot{\Sigma}_{ex}^{MD}(\bm{x})$, $\dot{\Sigma}_{hk}^{MD}(\bm{x})$)
orthogonalities (B, D).

The CRN depicted in Fig. 1 (A) is defined by the following set of chemical
reaction equations:

$\displaystyle 2\mathbb{X}_{1}$
$\displaystyle\xleftrightharpoons[k^{+}_{1}]{k^{-}_{1}}2\mathbb{X}_{2},$
$\displaystyle 2\mathbb{X}_{2}$
$\displaystyle\xleftrightharpoons[k^{+}_{2}]{k^{-}_{2}}\mathbb{X}_{1}+\mathbb{X}_{2},$
$\displaystyle\mathbb{X}_{1}+\mathbb{X}_{2}$
$\displaystyle\xleftrightharpoons[k^{+}_{3}]{k^{-}_{3}}2\mathbb{X}_{1}.$ (75)

The corresponding structural quantities are

$\displaystyle S$ $\displaystyle=\begin{pmatrix}-2&+1&+1\\\
+2&-1&-1\end{pmatrix},$ $\displaystyle B$
$\displaystyle=\begin{pmatrix}+1&0&-1\\\ -1&+1&0\\\ 0&-1&+1\end{pmatrix},$
(76) $\displaystyle\Gamma$ $\displaystyle=\begin{pmatrix}2&0&1\\\
0&2&1\end{pmatrix}.$ (77)

One can show that this simple nonlinear CRN has any of the equilibrium, CB,
and non-CB steady states depending on the kinetic parameter values[50, 51]. In
addition, the flux and force spaces are three-dimensional, and, thereby, the
relevant geometric objects and quantities can be visualized computationally.
We used the parameter values in Tab. 1 for the simulations in Fig. 6 and in
Fig. 7. Because the state space is two-dimensional and the stoichiometric
compatibility class is one-dimensional, i.e.,
$\mathcal{P}(\bm{x}_{0}):=\\{\bm{x}|\bm{x}-\bm{x}_{0}\in\mathrm{Im}[S]\\}=\\{\bm{x}|\bm{x}=\bm{x}_{0}+\xi(1,-1)^{T},\xi\in\mathbb{R}\\}$,
the trajectory is restricted in the one-dimensional $\mathcal{P}(\bm{x}_{0})$
independently of the parameter values. We also selected the parameter values
in Tab. 1 so that all the sets of equilibrium, CB and non-CB steady states are
equivalent to the set, $\\{\bm{x}|x_{2}=2x_{1}\\}$. Thus, the topological
properties of the state space are the same for the three sets of parameters.

Figure 6 shows the values of the EPR and the associated dissipation functions
in the state space. For the equilibrium parameter (Fig. 6 (A)), the EPR and
dissipation functions attain their global minimum value of zero on the
equilibrium states (black line). For the CB parameter (Fig. 6 (B)), the global
minimum is not necessarily attained, which results in the non-zero EPR at the
CB steady state $\bm{x}_{cb}$ (Fig. 7 (A)). For the non-CB parameter (Fig. 6
(C)), the landscapes of the EPR and dissipation functions become more
complicated, and their values along the trajectory $\bm{x}_{t}$ are not
monotonous (Fig. 7 (C)).

For the CB and non-CB cases, the two decompositions, Eq. (56) and Eq. (72),
are computed (Fig. 7 (B) and (D)). For the CB case, we verify that both
decompositions provide the expected behaviors: both housekeeping and excess
components stay non-negative; the excess one converges to $0$; the
housekeeping one does to a finite non-negative value(Fig. 7 (B)).

For the non-CB case, the expected behaviors are also produced (Fig. 7 (D)) for
the information-geometric decomposition, Eq. (72) (a generalization of Maes-
Netočnỳ decomposition. On the contrary, the housekeeping component of the
Hilbert decomposition, Eq. (56), becomes negative within a certain time window
because the Hilbert orthogonality, Eq. (55), does not hold in this case. It
should be noted that this result does not mean that the Hilbert decomposition
is not applicable to the non-CB case. The result only indicates that the force
decomposition by
$\bm{f}_{S}(\bm{x}):=-\frac{1}{2}S^{T}\partial_{\bm{x}}\varphi_{cb}(\bm{x})$
and $\bm{f}_{A}:=\frac{1}{2}\ln\frac{\bm{K}_{cb}}{\bm{K}_{eq}}$ with the
specific form of $\varphi_{cb}(\bm{x}):=\mathcal{D}_{KL}[\bm{x}\|\bm{x}_{cb}]$
no longer satisfies the Hilbert orthogonality. We may be able to recover the
orthogonality and the associated decomposition[34] by choosing other
functional form of $\varphi_{cb}(\bm{x})$. Finding such $\varphi_{cb}(\bm{x})$
is linked to the computation of quasi-potential, but this problem is still
challenging for generic non-CB cases.

Finally, we computationally visualized the various geometric objects
introduced in this work in Figs. 8 and 9 for the CB parameter set in Tab. 1.

Figure 8: Computational visualization of the Hilbert orthogonality between
$\bm{f}_{S}(\bm{x}_{t})$ and $\bm{f}_{A}$ for the CRN in Fig 1 (A). The green
solid surface in $\mathcal{F}$ represents the iso-dissipation hypersurface:
$\Psi^{*}_{\bm{x}_{t}}(\bm{f})=\Psi^{*}_{\bm{x}_{t}}(\bm{f}_{0}(\bm{x}_{t}))=\mathrm{const.}$.
The black circle with a white border, red square with a black border, green
star, and magenta circle with a white border in $\mathcal{F}$ are respectively
$\bm{f}_{0}(\bm{x}_{t})$, $\bm{f}_{S}(\bm{x}_{t})$,
$\bm{f}_{iso}(\bm{x}_{t})$, and $\bm{f}_{A}$ evaluated at
$\bm{x}_{t}=(5/2,1/2)^{T}$. The black circle, red square, green star, and
magenta circle with borders in $\mathcal{J}$ are their Legendre transform,
i.e., $\bm{j}_{0}(\bm{x}_{t})$, $\bm{j}_{S}(\bm{x}_{t})$,
$\bm{j}_{iso}(\bm{x}_{t})$, and $\bm{j}_{A}(\bm{x}_{t})$, respectively. The
black, red, and green circles without a border in $\mathcal{F}$ are the
trajectories of $\\{\bm{f}_{0}(\bm{x}_{t})\\}$,
$\\{\bm{f}_{S}(\bm{x}_{t})\\}$, and $\\{\bm{f}_{iso}(\bm{x}_{t})\\}$. The
black, red, green, and magenta circles without a border in $\mathcal{J}$ are
the trajectories of $\\{\bm{j}_{0}(\bm{x}_{t})\\}$,
$\\{\bm{j}_{S}(\bm{x}_{t})\\}$, $\\{\bm{j}_{iso}(\bm{x}_{t})\\}$, and
$\\{\bm{j}_{A}(\bm{x}_{t})\\}$, respectively. The light blue lines in
$\mathcal{F}$ are the NEF subspaces, $\mathcal{P}^{f}(\bm{f}_{A})$(solid line)
and $\mathcal{P}^{f}(-\bm{f}_{A})$ (dashed line), and the EF subspace
$\mathcal{P}^{f}_{eq}$ (dotted line). In $\mathcal{J}$, the corresponding NEF
manifolds, $\mathcal{M}^{f}_{\bm{x}}(\bm{f}_{A})$ (the light blue curve with
dots), is also depicted. The blue plane in $\mathcal{J}$ is the zero-velocity
subspace $\mathcal{P}^{v}_{st}$, and the blue surface in $\mathcal{F}$ is its
Legendre transformation: $\mathcal{M}^{v}_{\bm{x}_{t},st}$.  Figure 9:
Computational visualization of the information geometric orthogonality for the
CRN in Fig 1 (A). (Left) The orthogonality between $\bm{j}_{st}(\bm{x}_{t})$
(blue dashed arrow) and $\bm{j}_{0}(\bm{x}_{t})-\bm{j}_{st}(\bm{x}_{t})$ and
the dual orthogonality between $\bm{j}_{eq}(\bm{x}_{t})$ and
$\bm{j}_{0}(\bm{x}_{t})-\bm{j}_{eq}(\bm{x}_{t})$ (red solid arrow) in
$\mathcal{J}$. The blue planes are the iso-velocity subspace
$\mathcal{P}^{v}(\bm{j}_{0}(\bm{x}_{t}))$ (the upper plane) and the zero-
velocity subspace $\mathcal{P}^{v}_{st}$ (the lower plane). The grey curve
with dots is the NEF manifold
$\mathcal{M}^{f}_{\bm{x}}(\bm{f}_{0}(\bm{x}_{t}))$. Black, red, and blue
circles are the trajectories of $\\{\bm{j}_{0}(\bm{x}_{t})\\}$,
$\\{\bm{j}_{eq}(\bm{x}_{t})\\}$, and $\\{\bm{j}_{st}(\bm{x}_{t})\\}$,
respectively. (Right) The same orthogonalities shown in $\mathcal{F}$ space.
The orthogonality between $\bm{f}_{st}(\bm{x}_{t})$ (blue curved arrow) and
$\bm{f}_{0}(\bm{x}_{t})-\bm{f}_{st}(\bm{x}_{t})$ and the dual orthogonality
between $\bm{f}_{eq}(\bm{x}_{t})$ and
$\bm{f}_{0}(\bm{x}_{t})-\bm{f}_{eq}(\bm{x}_{t})$ (red curved arrow) in
$\mathcal{F}$. Black, red circles are the trajectories of
$\\{\bm{f}_{0}(\bm{x}_{t})\\}$, $\\{\bm{f}_{eq}(\bm{x}_{t})\\}$, respectively.
The blue surfaces are the iso-velocity manifold
$\mathcal{M}^{v}_{\bm{x}_{t}}(\bm{j}_{0}(\bm{x}_{t}))$ (the upper surface) and
the zero-velocity manifold $\mathcal{M}^{v}_{\bm{x}_{t},st}$ (the lower
surface). The grey lines are the NEF subspaces
$\mathcal{P}^{f}(\bm{f}_{0}(\bm{x}_{t}))$ and the EF subspace
$\mathcal{P}^{f}_{\bm{x},eq}$.

First, we show the behaviors of the forces, $\bm{f}_{0}(\bm{x}_{t})$,
$\bm{f}_{iso}(\bm{x}_{t})$, $\bm{f}_{S}(\bm{x}_{t})$, and $\bm{f}_{A}$, which
satisfy the Hilbert orthogonality, and the corresponding fluxes. In the force
space $\mathcal{F}$ (the right panel of Fig. 8), we observe that the
trajectories of $\bm{f}_{0}(\bm{x}_{t})$, $\bm{f}_{iso}(\bm{x}_{t})$, and
$\bm{f}_{S}(\bm{x}_{t})$ are actually restricted on one dimensional subspaces
(lines), $\mathcal{P}^{f}(\bm{f}_{A})$, $\mathcal{P}^{f}(-\bm{f}_{A})$, and
$\mathcal{P}^{f}_{eq}=\mathcal{P}^{f}(\bm{0})$, respectively.
$\bm{f}_{0}(\bm{x}_{t})$ and $\bm{f}_{iso}(\bm{x}_{t})$ are also on the iso-
dissipation hypersuface (the green surface) satisfying
$\Psi^{*}_{\bm{x}_{t}}(\bm{f}_{0}(\bm{x}_{t}))=\Psi^{*}_{\bm{x}_{t}}(\bm{f}_{iso}(\bm{x}_{t}))$.
In the flux space $\mathcal{J}$ (the left panel of Fig. 8),
$\mathcal{P}^{f}(\bm{f}_{A})$ is transformed to the one-dimensional curve,
$\mathcal{M}^{f}_{\bm{x}_{t}}(\bm{f}_{A})$. The trajectories of
$\bm{j}_{0}(\bm{x}_{t})$, $\bm{j}_{iso}(\bm{x}_{t})$, and
$\bm{j}_{S}(\bm{x}_{t})$ are also curved in $\mathcal{J}$. All the
trajectories converge onto the zero-velocity subspace $\mathcal{P}^{v}_{st}$.

Next, we verify the information-geometric orthogonality in Fig. 9. In the flux
space $\mathcal{J}$ (the left panel in Fig. 9), we observe that
$\bm{j}_{0}(\bm{x}_{t})$ and the corresponding equilibrium flux
$\bm{j}_{eq}(\bm{x}_{t})$ are on the iso-velocity subspace
$\mathcal{P}^{v}(\bm{j}_{0}(\bm{x}_{t}))$. The steady state flux
$\bm{j}_{st}(\bm{x}_{t})$ are on the intersection of the zero-velocity
subspace $\mathcal{P}^{v}_{st}$ and the NEF manifold
$\mathcal{M}^{f}_{\bm{x}_{t}}(\bm{f}_{0}(\bm{x}_{t}))$. By mapping these
objects into $\mathcal{F}$ by the Legendre transform, we observe that
$\bm{f}_{0}(\bm{x}_{t})$ is on the intersection of the NEF subspace
$\mathcal{P}^{f}(\bm{f}_{0}(\bm{x}_{t}))$ and the iso-velocity manifold
$\mathcal{M}^{v}_{\bm{x}_{t}}(\bm{j}_{0}(\bm{x}_{t}))$. Similarly,
$\bm{f}_{eq}(\bm{x}_{t})$ is on the intersection of the EF subspace
$\mathcal{P}^{f}_{eq}$ and the iso-velocity manifold
$\mathcal{M}^{v}_{\bm{x}_{t}}(\bm{j}_{0}(\bm{x}_{t}))$. By the Legendre
transformation, the iso-velocity manifold
$\mathcal{M}^{v}_{\bm{x}_{t}}(\bm{j}_{0}(\bm{x}_{t}))$, and the zero-velocity
manifold $\mathcal{M}^{v}_{\bm{x}_{t},st}$ are curved in $\mathcal{F}$. These
graphical representations demonstrate how the nonlinear Legendre transform
relates the subspaces and manifolds in $\mathcal{J}$ and $\mathcal{F}$ and
also how the generalized notions of orthogonality of Hessian geometry are
realized in these spaces.

##  VII Summary and Discussion

We have clarified that Hessian geometry is the natural geometric structure of
nonequilibrium and nonlinear CRNs as well as MJPs. More generally, it can
capture the geometry induced by nonquadratic convex dissipation functions. By
employing generalized notions of orthogonality, various aspects of
nonequilibrium dynamics are dissected as decompositions of EPRs, which
generalize the well-established ones [47, 48, 25] to CRNs and MJPs.

###  VII.1 Thermodynamic uncertainty relation and Fisher information

The nonquadratic property of the dissipation functions of CRN and MJP also
appears implicitly in different problems of thermodynamics. One such example
is thermodynamic uncertainty relation (TUR) and its extension to MJP. TUR is a
relation that bounds the fluctuation of a generalized current in a
nonequilibrium steady state by the entropy production[52]:

$\displaystyle\frac{\mathbb{E}[\mathcal{J}_{d}^{ss}]^{2}}{\mathbb{Var}[\mathcal{J}_{d}^{ss}]}\leq\frac{1}{2}\int_{0}^{\tau}\dot{\Sigma}(t^{\prime})\mathrm{d}t^{\prime}$
(78)

where $\mathcal{J}_{d}^{ss}$ is a time-integrated generalized current at the
steady state, and $\mathbb{E}[\mathcal{J}_{d}^{ss}]$ and
$\mathbb{Var}[\mathcal{J}_{d}^{ss}]$ are the mean and variance of
$\mathcal{J}_{d}^{ss}$. Since the proposal of the TUR conjecture and its first
proof [53, 54], the TUR has been extended to various models and situations,
including MJP [52] and CRN [55]. However, even though Eq. (78) is tight for
overdamped diffusion processes, it is not tight for MJP. Instead, for MJP, the
fluctuation of the current is bounded tightly by the pseudo entropy
production(pEPR) [55]:

$\displaystyle\frac{\mathbb{E}[\mathcal{J}_{d}^{ss}]^{2}}{\mathbb{Var}[\mathcal{J}_{d}^{ss}]}\leq\frac{1}{2}\int_{0}^{\tau}\dot{\Pi}(t^{\prime})\mathrm{d}t^{\prime}$
(79)

where the pEPR is defined as

$\displaystyle\dot{\Pi}(t):=\frac{(\bm{j}^{+}(\bm{x}_{t})-\bm{j}^{-}(\bm{x}_{t}))^{2}}{\bm{j}^{+}(\bm{x}_{t})+\bm{j}^{+}(\bm{x}_{t})}.$
(80)

and $\dot{\Pi}(t)\leq\dot{\Sigma}(t)$ holds [52]. We here mention that this
gap between the EPR and pEPR is linked to the Fisher information of the
dissipation function. The Fisher information matrix (or metric) for a strictly
convex function $\phi^{*}(\bm{f})$ is defined by its Hessian as

$\displaystyle
G_{e,e^{\prime}}^{*}(\bm{f}):=\frac{\partial^{2}\phi^{*}(\bm{f})}{\partial
f_{e}\partial f_{e^{\prime}}}.$ (81)

For $\phi^{*}(\bm{f})=\Psi_{\bm{\omega}}^{*}(\bm{f})$ defined by Eq. (12), we
have

$\displaystyle
G^{*}(\bm{\omega},\bm{f})=\mathrm{diag}[\bm{\omega}\circ\cosh(\bm{f})]=(\bm{j}^{+}+\bm{j}^{-}),$
(82)

where we used Eq. (9) and Eq. (11). For
$\phi(\bm{j})=\Psi_{\bm{\omega}}(\bm{j})$, we have the Fisher information
matrix:

$\displaystyle
G(\bm{\omega},\bm{j})=G^{*}(\bm{\omega},\bm{f})^{-1}=\frac{1}{\bm{j}^{+}+\bm{j}^{-}}.$
(83)

Then, the pEPR is represented as

$\displaystyle\dot{\Pi}(\bm{j})=\langle\bm{j},G(\bm{\omega},\bm{j})\bm{j}\rangle,$
(84)

which allows us to regard $\dot{\Pi}(t)$ as an approximation of the EPR by
replacing the actual force $\bm{f}$ associated with $\bm{j}$ with a pseudo-
force $\bm{f}^{p}:=G(\bm{\omega},\bm{j})\bm{j}$. If the dissipation function
is quadratic, the actual force $\bm{f}$ and the pseudo-force $\bm{f}^{p}$
coincide. Thus, the non-tightness of the TUR for MJP is a manifestation of the
nonquadratic nature of the dissipation functions.

###  VII.2 Network thermodynamics

The information geometric decomposition that we introduced is tightly related
to the network thermodynamics [56, 18, 19, 57, 58]. The network thermodynamics
is an attempt to extend the methodology of network theory developed mainly for
linear electric circuits to other physical systems [59, 60]. While the linear
network theory successfully works for linear electric circuits, which have a
linear force and flux relation (Ohm’s law) and resulting quadratic dissipation
functions, its application to MJP and CRN has encountered difficulties due to
the nonlinear relation between force and flux, especially when we evaluate the
EPR decomposition for a transient and far from equilibrium state. The Hessian
geometric structure and the generalized decompositions in this work are
extensions of those in network thermodynamics. Specifically, the information
geometric decomposition is a nonlinear generalization of the cycle-cocycle
decomposition [61]. However, network thermodynamics and network theory
accommodate a wide variety of methods that are not yet exploited in this work
[62, 61]. For example, the algebraic structure underlying the network is
explored by using the integral basis defined by cycles, cords, and spanning
trees in the network [19, 63]. As an application, the steady EPR is related to
cyclic fluxes and affinities in the network. We may incorporate and exploit
the algebraic structure of MJP and CRN more explicitly into our framework.

###  VII.3 Duality and variational characterizations in thermodynamics

The notion of duality is the core of thermodynamics. In equilibrium
thermodynamics, the duality of extensive and intensive variables induced by
thermodynamic potential functions characterizes the energetic aspect of
physical systems. The roles and differences of dual pairs of thermodynamic
potential functions are well recognized. In addition, Hessian geometry is
central to comprehend the geometric properties of equilibrium thermodynamics
[37, 2]. For the duality between the force and flux explored in this work, the
dual dissipation functions have been recognized for a while in MFT [7] and
metric gradient flow theory [14]. Nevertheless, nonquadratic dissipation
functions have been investigated only very recently even in MFT [31, 32, 33,
34] and also in other physics communities [64]. As in the case of equilibrium
thermodynamics, grasping the roles played by the dual functions should be
essential for understanding nonequilibrium and kinetic aspects of
thermodynamics. Hence, Hessian geometry becomes an indispensable tool for
investigating the geometry induced by the duality and also the variational
aspects of nonequilibrium phenomena. In thermodynamics, there exist continued
attempts to characterize nonequilibrium states and relations variationally,
which go by the names of minimum entropy production principle [65, 66, 67],
maximum entropy production principle [68, 69], the least dissipation principle
[5, 70], and others. However, all these principles still have limitations in
their applicability for nonlinear and far from equilibrium situations [71].
Hessian geometry and also information geometry, which has an ability to handle
nonlinearity induced by convex functions, may contribute to resolving a part
of problems in such principles. Moreover, they can serve as the natural
language to integrate the equilibrium (energetic) and nonequilibrium (kinetic)
descriptions, and thereby provide us a more universal understanding of
thermodynamic systems and extend the applicability of nonequilibrium
thermodynamics.

##  VIII Acknowledgement

This research is supported by JST (JPMJCR2011,JPMJCR1927) and JSPS (19H05799).
We thank Shuhei Horiguchi for valuable suggestions.

## References

  * Callen and Callen [1985] H. B. Callen and H. B. Callen,  _Thermodynamics and an Introduction to Thermostatistics_ (Wiley, 1985).
  * 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, ArXiv211212403 Cond-Mat Physicsphysics (2021), arXiv:2112.12403 [cond-mat, physics:physics] .
  * Onsager [1931a] L. Onsager, Reciprocal Relations in Irreversible Processes. I., Phys. Rev. 37, 405 (1931a).
  * Onsager [1931b] L. Onsager, Reciprocal Relations in Irreversible Processes. II., Phys. Rev. 38, 2265 (1931b).
  * Onsager and Machlup [1953] L. Onsager and S. Machlup, Fluctuations and Irreversible Processes, Phys. Rev. 91, 1505 (1953).
  * Machlup and Onsager [1953] S. Machlup and L. Onsager, Fluctuations and Irreversible Process. II. Systems with Kinetic Energy, Phys. Rev. 91, 1512 (1953).
  * Bertini  _et al._ [2015] L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio, and C. Landim, Macroscopic fluctuation theory, Rev. Mod. Phys. 87, 593 (2015).
  * Taniguchi and Cohen [2007] T. Taniguchi and E. G. D. Cohen, Onsager-Machlup Theory for Nonequilibrium Steady States and Fluctuation Theorems, J Stat Phys 126, 1 (2007).
  * Schmiedl and Seifert [2007] T. Schmiedl and U. Seifert, Stochastic thermodynamics of chemical reaction networks, J. Chem. Phys. 126, 044101 (2007).
  * Aurell  _et al._ [2011] E. Aurell, C. Mejía-Monasterio, and P. Muratore-Ginanneschi, Optimal Protocols and Optimal Transport in Stochastic Thermodynamics, Phys. Rev. Lett. 106, 250601 (2011).
  * Nakazato and Ito [2021] M. Nakazato and S. Ito, Geometrical aspects of entropy production in stochastic thermodynamics based on Wasserstein distance, Phys. Rev. Research 3, 043093 (2021).
  * Jordan  _et al._ [1997] R. Jordan, D. Kinderlehrer, and F. Otto, Free energy and the Fokker-Planck equation, Physica D: Nonlinear Phenomena 16th Annual International Conference of the Center for Nonlinear Studies, 107, 265 (1997).
  * Benamou and Brenier [2000] J.-D. Benamou and Y. Brenier, A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem, Numer. Math. 84, 375 (2000).
  * Ambrosio  _et al._ [2006] L. Ambrosio, N. Gigli, and G. Savare,  _Gradient Flows: In Metric Spaces and in the Space of Probability Measures_ (Springer Science & Business Media, 2006).
  * Shima [2007] H. Shima,  _The Geometry of Hessian Structures_ (World Scientific, 2007).
  * Amari [2016] S.-i. Amari,  _Information Geometry and Its Applications_ (Springer, 2016).
  * Nielsen [2020] F. Nielsen, An Elementary Introduction to Information Geometry, Entropy 22, 1100 (2020).
  * Hill [2005] T. L. Hill,  _Free Energy Transduction and Biochemical Cycle Kinetics_ (Courier Corporation, 2005).
  * Schnakenberg [1976] J. Schnakenberg, Network theory of microscopic and macroscopic behavior of master equation systems, Rev. Mod. Phys. 48, 571 (1976).
  * 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, Phys. Rev. X 6, 041064 (2016).
  * Beard and Qian [2008] D. A. Beard and H. Qian,  _Chemical Biophysics: Quantitative Analysis of Cellular Systems_, Cambridge Texts in Biomedical Engineering (Cambridge University Press, Cambridge, 2008).
  * 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).
  * Dechant  _et al._ [2022] A. Dechant, S.-i. Sasa, and S. Ito, Geometric decomposition of entropy production in out-of-equilibrium systems, ArXiv210912817 Cond-Mat (2022), arXiv:2109.12817 [cond-mat] .
  * Grady and Polimeni [2010] L. J. Grady and J. R. Polimeni,  _Discrete Calculus: Applied Analysis on Graphs for Computational Science_ (Springer Science & Business Media, 2010).
  * Gibbs [1878] J. W. Gibbs, On the equilibrium of heterogeneous substances, Am. J. Sci. s3-16, 441 (1878).
  * Feinberg [2019] M. Feinberg,  _Foundations of Chemical Reaction Network Theory_ (Springer, 2019).
  * Klamt  _et al._ [2009] S. Klamt, U.-U. Haus, and F. Theis, Hypergraphs and Cellular Networks, PLOS Computational Biology 5, e1000385 (2009).
  * Maes [2017a] C. Maes, Frenetic Bounds on the Entropy Production, Phys. Rev. Lett. 119, 160601 (2017a).
  * Kaiser  _et al._ [2018] M. Kaiser, R. L. Jack, and J. Zimmer, Canonical Structure and Orthogonality of Forces and Currents in Irreversible Markov Chains, J Stat Phys 170, 1019 (2018).
  * Renger [2018] D. R. M. Renger, Gradient and GENERIC Systems in the Space of Fluxes, Applied to Reacting Particle Systems, Entropy 20, 596 (2018).
  * Renger and Zimmer [2021] D. R. M. Renger and J. Zimmer, Orthogonality of fluxes in general nonlinear reaction networks, Discrete Contin. Dyn. Syst. - S 14, 205 (2021).
  * Patterson  _et al._ [2021] R. I. A. Patterson, D. R. M. Renger, and U. Sharma, Variational structures beyond gradient flows: A macroscopic fluctuation-theory perspective, ArXiv210314384 Math-Ph (2021), arXiv:2103.14384 [math-ph] .
  * Maes [2021] C. Maes, Local detailed balance, SciPost Phys. Lect. Notes , 032 (2021).
  * Maes [2017b] C. Maes,  _Non-Dissipative Effects in Nonequilibrium Systems_ (Springer, 2017).
  * 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, ArXiv211214910 Phys. (2021), arXiv:2112.14910 [physics] .
  * Sughiyama  _et al._ [2022] Y. Sughiyama, A. Kamimura, D. Loutchko, and T. J. Kobayashi, The Chemical Thermodynamics for Growing Systems, ArXiv220109417 Cond-Mat Physicsphysics (2022), arXiv:2201.09417 [cond-mat, physics:physics] .
  * Renger [2021] D. R. M. Renger, Anisothermal Chemical Reactions: Onsager-Machlup and Macroscopic Fluctuation Theory, ArXiv211112164 Cond-Mat Physicsmath-Ph (2021), arXiv:2111.12164 [cond-mat, physics:math-ph] .
  * Nomizu [1994] K. Nomizu,  _Affine Differential Geometry: Geometry of Affine Immersions_ (Cambridge University Press, 1994).
  * Ito [2018] S. Ito, Stochastic Thermodynamic Interpretation of Information Geometry, Phys. Rev. Lett. 121, 030605 (2018).
  * Kolchinsky and Wolpert [2021] A. Kolchinsky and D. H. Wolpert, Work, Entropy Production, and Thermodynamics of Information under Protocol Constraints, Phys. Rev. X 11, 041024 (2021).
  * Ohga and Ito [2021] N. Ohga and S. Ito, Information-geometric Legendre duality in stochastic thermodynamics, ArXiv211211008 Cond-Mat (2021), arXiv:2112.11008 [cond-mat] .
  * Craciun  _et al._ [2009] G. Craciun, A. Dickenstein, A. Shiu, and B. Sturmfels, Toric dynamical systems, Journal of Symbolic Computation In Memoriam Karin Gatermann, 44, 1551 (2009).
  * Mielke  _et al._ [2014] A. Mielke, M. A. Peletier, and D. R. M. Renger, On the Relation between Gradient Flows and the Large-Deviation Principle, with Applications to Markov Chains and Diffusion, Potential Anal 41, 1293 (2014).
  * Hatano and Sasa [2001] T. Hatano and S.-i. Sasa, Steady-State Thermodynamics of Langevin Systems, Phys. Rev. Lett. 86, 3463 (2001).
  * Maes and Netočný [2014] C. Maes and K. Netočný, A Nonequilibrium Extension of the Clausius Heat Theorem, J Stat Phys 154, 188 (2014).
  * Komatsu  _et al._ [2008] T. S. Komatsu, N. Nakagawa, S.-i. Sasa, and H. Tasaki, Steady-State Thermodynamics for Heat Conduction: Microscopic Derivation, Phys. Rev. Lett. 100, 230602 (June 13, 2008).
  * Bhatia  _et al._ [2013] H. Bhatia, G. Norgard, V. Pascucci, and P.-T. Bremer, The Helmholtz-Hodge Decomposition—A Survey, IEEE Trans. Vis. Comput. Graph. 19, 1386 (2013).
  * Horn and Rowlinson [1973] F. Horn and J. S. Rowlinson, Stability and complex balancing in mass-action systems with three short complexes, Proc. R. Soc. Lond. Math. Phys. Sci. 334, 331 (1973).
  * Pérez Millán  _et al._ [2012] M. Pérez Millán, A. Dickenstein, A. Shiu, and C. Conradi, Chemical Reaction Systems with Toric Steady States, Bull Math Biol 74, 1027 (2012).
  * Shiraishi [2021] N. Shiraishi, Optimal Thermodynamic Uncertainty Relation in Markov Jump Processes, J Stat Phys 185, 19 (2021).
  * Gingrich  _et al._ [2016] T. R. Gingrich, J. M. Horowitz, N. Perunov, and J. L. England, Dissipation Bounds All Steady-State Current Fluctuations, Phys. Rev. Lett. 116, 120601 (2016).
  * Barato and Seifert [2015] A. C. Barato and U. Seifert, Thermodynamic Uncertainty Relation for Biomolecular Processes, Phys. Rev. Lett. 114, 158101 (2015).
  * Yoshimura and Ito [2021] K. Yoshimura and S. Ito, Thermodynamic Uncertainty Relation and Thermodynamic Speed Limit in Deterministic Chemical Reaction Networks, Phys. Rev. Lett. 127, 160601 (2021).
  * Oster  _et al._ [1971] G. Oster, A. Perelson, and A. Katchalsky, Network Thermodynamics, Nature 234, 393 (1971).
  * Schnakenberg [2012] J. Schnakenberg,  _Thermodynamic Network Analysis of Biological Systems_ (Springer Science & Business Media, 2012).
  * van der Schaft  _et al._ [2013] A. van der Schaft, S. Rao, and B. Jayawardhana, On the Mathematical Structure of Balanced Chemical Reaction Networks Governed by Mass Action Kinetics, SIAM J. Appl. Math. 73, 953 (2013).
  * Belevitch [1962] V. Belevitch, Summary of the History of Circuit Theory, Proc. IRE 50, 848 (1962).
  * Chen [2012] W.-K. Chen,  _Applied Graph Theory_ (Elsevier, 2012).
  * Wachtel  _et al._ [2015] A. Wachtel, J. Vollmer, and B. Altaner, Fluctuating currents in stochastic thermodynamics. I. Gauge invariance of asymptotic statistics, Phys. Rev. E 92, 042132 (2015).
  * Altaner  _et al._ [2012] B. Altaner, S. Grosskinsky, S. Herminghaus, L. Katthän, M. Timme, and J. Vollmer, Network representations of nonequilibrium steady states: Cycle decompositions, symmetries, and dominant paths, Phys. Rev. E 85, 041133 (2012).
  * Polettini [2015] M. Polettini, System/Environment Duality of Nonequilibrium Network Observables, in  _Mathematical Technology of Networks_, Springer Proceedings in Mathematics & Statistics, edited by D. Mugnolo (Springer International Publishing, Cham, 2015) pp. 191–205.
  * Suzuki [2012] M. Suzuki, Irreversibility and entropy production in transport phenomena, II: Statistical–mechanical theory on steady states including thermal disturbance and energy supply, Physica A: Statistical Mechanics and its Applications 391, 1074 (2012).
  * Prigogine [1968] I. Prigogine,  _Introduction to Thermodynamics of Irreversible Processes_ (Wiley, 1968).
  * Klein and Meijer [1954] M. J. Klein and P. H. E. Meijer, Principle of Minimum Entropy Production, Phys. Rev. 96, 250 (1954).
  * Polettini [2011] M. Polettini, Macroscopic constraints for the minimum entropy production principle, Phys. Rev. E 84, 051117 (2011).
  * Županović  _et al._ [2004] P. Županović, D. Juretić, and S. Botrić, Kirchhoff’s loop law and the maximum entropy production principle, Phys. Rev. E 70, 056108 (2004).
  * Martyushev [2021] L. M. Martyushev, Maximum entropy production principle: History and current status, Phys.-Usp. 64, 558 (2021).
  * Gyarmati [2013] I. Gyarmati,  _Non-Equilibrium Thermodynamics: Field Theory and Variational Principles_ (Springer Science & Business Media, 2013).
  * Polettini [2013] M. Polettini, Fact-Checking Ziegler’s Maximum Entropy Production Principle beyond the Linear Regime and towards Steady States, Entropy 15, 2570 (2013).