# Transitions and Thermodynamics on Species Graphs of
Chemical Reaction Networks

Keisuke Sugie<EMAIL_ADDRESS>Department of Mathematical
Informatics, Graduate School of Information Science and Technology, The
University of Tokyo, Tokyo 113-8654, Japan  Dimitri Loutchko  Institute of
Industrial Science, The University of Tokyo, 4-6-1, Komaba, Meguro-ku, Tokyo
153-8505 Japan  Tetsuya J. Kobayashi  Department of Mathematical Informatics,
Graduate School of Information Science and Technology, The University of
Tokyo, Tokyo 113-8654, Japan  Institute of Industrial Science, The University
of Tokyo, 4-6-1, Komaba, Meguro-ku, Tokyo 153-8505 Japan

###### Abstract

Chemical reaction networks (CRNs) exhibit complex dynamics governed by their
underlying network structure. In this paper, we propose a novel approach to
study the dynamics of CRNs by representing them on species graphs (S-graphs).
By scaling concentrations by conservation laws, we obtain a graph
representation of transitions compatible with the S-graph, which allows us to
treat the dynamics in CRNs as transitions between chemicals. We also define
thermodynamic-like quantities on the S-graph from the introduced transitions
and investigate their properties, including the relationship between
specieswise forces, activities, and conventional thermodynamic quantities.
Remarkably, we demonstrate that this formulation can be developed for a class
of irreversible CRNs, while for reversible CRNs, it is related to conventional
thermodynamic quantities associated with reactions. The behavior of these
specieswise quantities is numerically validated using an oscillating system
(Brusselator). Our work provides a novel methodology for studying dynamics on
S-graphs, paving the way for a deeper understanding of the intricate interplay
between the structure and dynamics of chemical reaction networks.

††preprint: APS/123-QED

##  I Introduction.

The theoretical framework of chemical reaction network (CRNs for short) is
applied in various research fields, such as steady-state flux analysis in
metabolic systems [1] and modeling non-equilibrium processes like oscillatory
phenomena [2]. To investigate the properties that hold in a broad range of
processes in CRNs, including non-equilibrium states, thermodynamic approaches
have recently been used [3, 4]. Mathematically, CRNs possess a geometric
structure represented by hypergraphs with stoichiometric information [5],
whose structural constraints govern the behaviour of chemical reaction
dynamics, for example, the controllability of chemical reactions [6].

For large-scale CRNs, the species graph approach has been utilized. This is
due to the difficulty in intuitively understanding transitions between
chemical species over the hypergraph structure, as well as the fact that
graph-based characterisations such as small-world and scale-free properties
for graphs are relatively well established [7, 8]. With this background,
graph-based studies using species graphs (S-graph) or reaction graphs [9] are
often preferred over hypergraphs for large networks such as metabolic systems
[10, 11].

However, the dynamics corresponding to the S-graph have not been extensively
explored. Despite the potential importance of dynamical oscillations in actual
metabolic and glycolytic systems [12, 13], the dynamical aspects of species
transitions in chemical reaction systems have been neglected.

In this paper, through scaling concentrations by conservation laws, we
obtained a graph representation of transitions compatible with the S-graphs of
chemical reaction networks. This transition representation enabled us to
derive a novel thermodynamics-like formulation for CRNs. Remarkably, we
demonstrated that this thermodynamics-like formulation can be applied to a
class of irreversible CRNs, while to reversible CRNs, it is reduced to
conventional thermodynamic quantities associated with reactions. A key
contribution of our research is proposing a methodology to derive a master
equation on the S-graph governing the scaled concentrations from deterministic
CRN dynamics. Our work provides a novel approach to study of dynamics on
S-graphs, paving the way for a deeper understanding of the intricate interplay
between the structure and dynamics of chemical reaction networks.

The paper is structured as follows: The following section II introduces
fundamental terms related to the CRN dynamics and note thermodynamic
quantities defined in a reversible CRN. In the section III, we define a new
physical quantity called ”transition” and reinterpret the CRN dynamics as
dynamics on a species graph. The conservative nature of the dynamics on the
graph is realised through concentration scaling by the conservation laws of
the network. In the section IV, we define thermodynamic quantities on the
species graph from the transitions we have introduced and investigate their
properties. In particular, we discuss the relationship of the specieswise
forces and activities with conventional thermodynamic quantities. The
behaviour of the specieswise quantities is checked numerically using an
oscillating system (Brusselator) in the section V.

Figure 1:  Schematic of the method for obtaining the dynamics on the S-graph.
The physical quantity $T^{\ell}$ (right) compatible with the S-graph is
obtained from the CRN dynamics (left), given a conservation law $\ell$ of the
CRN.

##  II Foundation of chemical reaction network.

Fundamentals of chemical reaction networks and kinetics are reviewed in this
section. The detailed notations can be found in [14, 15].

###  II.1 Notation and Preliminaries

Consider a non empty finite set $\mathcal{I}$ of _species_. Let
$\mathcal{C}\subset\mathbb{Z}_{\geq 0}^{\mathcal{I}}$ be a finite set of
_complexes_. Suppose a non empty finite set of _reactions_
$\mathcal{R}\subseteq\mathbb{Z}_{\geq 0}^{\mathcal{I}}\times\mathbb{Z}_{\geq
0}^{\mathcal{I}}$ satisfying the following conditions,

$\displaystyle\emptyset$ $\displaystyle=$ $\displaystyle\quantity{(y,y)\mid
y\in\mathcal{C}}\cap\mathcal{R},$ (1) $\displaystyle\mathcal{C}$
$\displaystyle=$
$\displaystyle\quantity{y,y^{\prime}\mid(y,y^{\prime})\in\mathcal{R}}.$ (2)

We denote a reaction $(y,y^{\prime})\in\mathcal{R}$ as $y\to y^{\prime}$. For
a complex $y\in\mathcal{C}$, its support is defined as
$\mathrm{supp}(y)=\quantity{i\in\mathcal{I}\mid y_{i}>0}$. $y$ and
$y^{\prime}$ are respectively called reactant and product complexes of a
reaction $y\to y^{\prime}$. Each
$\mathrm{supp}(y),\mathrm{supp}(y^{\prime})\subseteq\mathcal{I}$ respectively
represents a set of reactants and products. The triple
$\mathcal{N}=(\mathcal{I},\mathcal{C},\mathcal{R})$ of species, complexes and
reactions is called a chemical reaction network (CRN) if $\mathcal{N}$
satisfies the above conditions.

We distinguish _proper reactions_ $y\to y^{\prime}$, which have $y\neq 0$ and
$y^{\prime}\neq 0$, _import reactions_ for which $y=0$, and _export reactions_
for which $y^{\prime}=0$. A CRN is _closed_ if all reactions $y\to
y^{\prime}\in\mathcal{R}$ are proper. We focus on a closed CRN
$(\mathcal{I},\mathcal{C},\mathcal{R})$, unless explicitly stated otherwise.

Consider $1$-norm $\norm{\cdot}$ on a complex set $\mathcal{C}$ of
$\mathcal{N}$; for a complex $y\in\mathcal{C}$,
$\norm{y}\coloneq\sum_{i\in\mathcal{I}}y_{i}>0$. The non-emptiness of a
complex’s support is equivalent to the positiveness of its norm;
$\mathrm{supp}(y)\neq\emptyset\Leftrightarrow\norm{y}\geq 0$. In the same way,
properness of a reaction is represented in terms of its support or norm,

$\displaystyle y\to y^{\prime}\text{ is a proper reaction}$ (3)
$\displaystyle\Leftrightarrow$
$\displaystyle\mathrm{supp}(y)\neq\emptyset\text{ and
}\mathrm{supp}(y^{\prime})\neq\emptyset$ (4) $\displaystyle\Leftrightarrow$
$\displaystyle\norm{y}>0\text{ and }\norm{y^{\prime}}>0.$ (5)

In this paper, $\norm{y}$ ($\norm{y^{\prime}}$) is called a _total reactant
(product) coefficient_ of a reaction $y\to y^{\prime}\in\mathcal{R}$. When a
reaction has the same value of its two total coefficients
$\norm{y}=\norm{y^{\prime}}$, the value is simply called a _total coefficient_
or _degree_ of $y\to y^{\prime}$.

In our notation, the directed graph $\mathcal{G}_{\mathcal{N}}$ induced by a
CRN $\mathcal{N}=(\mathcal{I},\mathcal{C},\mathcal{R})$ is defined as
$\mathcal{G}=(\mathcal{I},\mathcal{E})$ with the edge set

$\displaystyle\mathcal{E}=\quantity{i\leftarrow j\mid\exists y\to
y^{\prime}\in\mathcal{R}_{\mathcal{N}},y^{\prime}_{i}y_{j}>0}$ (6)

where the tuple $(i,j)$ is written as $i\leftarrow j$. $i$ ($j$) is called a
_head_ (_tail_) of the edge $e=i\leftarrow j$ and denoted as $i=h(e)$ or
$j=h(e)$. The induced graph $(\mathcal{I},\mathcal{E})$ is called a _species
graph_ (or _S-graph_ for short), which has a set of species as nodes and edges
whose head is a reactant and tail is a product of a certain reaction [9, 11].
The _incidence matrix_ $\mathbb{B}\in\quantity{0,\pm
1}^{\mathcal{I}\times\mathcal{E}}$ of the graph $(\mathcal{I},\mathcal{E})$ is
defined as

$\displaystyle\mathbb{B}_{ie}=\begin{cases}1&i\text{ is the tail of }e,\\\
-1&i\text{ is the head of }e,\\\ 0&\text{otherwise}.\end{cases}$ (7)

###  II.2 Conservation laws of CRN

$\mathcal{S}_{\mathcal{N}}=\mathrm{span}\quantity{y^{\prime}-y\mid y\to
y^{\prime}\in\mathcal{R}}$ is called a stoichiometric subspace of
$\mathcal{N}$. Define the orthogonal complementrary subspace of
$\mathcal{S}_{\mathcal{N}}$ for a canonical inner product of
$\mathbb{R}^{\mathcal{I}}$:
$\mathcal{S}^{\bot}_{\mathcal{N}}=\quantity{\ell\in\mathbb{R}^{\mathcal{I}}\mid\forall
z\in\mathcal{S}_{\mathcal{N}},\ell^{\top}z=0}$, which is called a set of
_conservation laws_. Each element of $S^{\bot}_{\mathcal{N}}$ is called a
conservation law of $\mathcal{N}$. The set of _semi-positive conservation
laws_ , or non-negative conservation laws, of $\mathcal{N}$ is denoted as
$\mathcal{L}_{\mathcal{N}}=\mathcal{S}^{\bot}_{\mathcal{N}}\cap\mathbb{R}_{\geq
0}^{\mathcal{I}}$ [16, 17]. We define _positive conservation laws_ , which is
namely an element of the set
$\mathcal{L}_{\mathcal{N}}^{+}=\mathcal{S}^{\bot}_{\mathcal{N}}\cap\mathbb{R}_{>0}^{\mathcal{I}}$.
Equivalently, positive laws are semi-positive conservation laws which have no
$0$ elements. A CRN $\mathcal{N}$ is called _conservative_ if there exists a
nontrivial positive law, i.e. $\mathcal{L}_{\mathcal{N}}^{+}\neq\emptyset$.
Denote a CRN $\mathcal{N}$ is _$\ell$ -conservative_ if $\mathcal{N}$ has a
positive conservation law $\ell\in\mathcal{L}_{\mathcal{N}}^{+}$. In
particular, we write a $1_{\mathcal{I}}$-conservative CRN simply as a
$1$-conservative CRN.

###  II.3 Kinetics on chemical reaction network.

Next we consider the kinetics on a network
$\mathcal{N}=(\mathcal{I},\mathcal{C},\mathcal{R})$. Suppose
$\mathbb{R}_{>0}^{\mathcal{I}}$ as a concentration space and assume that the
CRN concentration $x_{t}$ at any time $t$ is in the concentration space, which
means that the concentration of each species does not go to $0$. Consider a
reaction rate
$j:\mathbb{R}_{>0}^{\mathcal{I}}\to\mathbb{R}_{>0}^{\mathcal{R}}$ of
$\mathcal{N}$, which is also called kinetics on $\mathcal{N}$. The
deterministic dynamics on the CRN is written as

$\displaystyle\differential{}_{t}x_{t}=\sum_{y\to
y^{\prime}\in\mathcal{R}}(y^{\prime}-y)j_{y\to y^{\prime}}(x_{t}).$ (8)

The stoichiometric structure restricts the concentration dynamics of the
initial condition $x_{0}$ to a linear subspace
$(x_{0}+S)\cap\mathbb{R}_{>0}^{\mathcal{I}}$, which is called a stoichiometric
compatibility class. Because of this stoichiometric constraint, we see that
the inner product of a conservation law $\ell$ and concentration $x_{t}$ is
conserved throughout the dynamics, i.e. $\ell^{\top}x_{t}=\ell^{\top}x_{0}$.
The kinetics $j$ is a mass action kinetics if $j_{y\to y^{\prime}}(x)=k_{y\to
y^{\prime}}x^{y}=k_{y\to y^{\prime}}\cdot\prod_{i\in\mathcal{I}}x_{i}^{y_{i}}$
with time-independent reaction rate constants
$k\in\mathbb{R}_{>0}^{\mathcal{R}}$. We denote the value $j_{y\to
y^{\prime}}(x)$ as $j_{y\to y^{\prime}}$ for simplicity.

###  II.4 Thermodynamic quantities on reversible CRN.

We introduce the thermodynamics on a reversible chemical kinetic system. A CRN
$(\mathcal{I},\mathcal{C},\mathcal{R})$ is called _reversible_ iff for any
reaction $y\to y^{\prime}\in\mathcal{R}$, the reverse also exists as a
reaction: $y^{\prime}\to y\in\mathcal{R}$. The reaction set of a reversible
CRN $\mathcal{R}$ can be represented as a disjoint sum of two sets
$\mathcal{R}^{+},\mathcal{R}^{-}$:
$\mathcal{R}=\mathcal{R}^{+}\sqcup\mathcal{R}^{-}$, which are satisfied

$\displaystyle y\to y^{\prime}\in\mathcal{R}^{\pm}\implies y^{\prime}\to
y\in\mathcal{R}^{\mp}.$ (9)

Reactions of $\mathcal{R}^{\pm}$ are called forward / backward reactions,
respectively. Some thermodynamic quantities are defined in the reversible CRN
system. $J_{y\to y^{\prime}}=j_{y\to y^{\prime}}-j_{y^{\prime}\to y}$ is
called _net current_ or _flow_ of reaction $y\to
y^{\prime}\in\mathcal{R}^{+}$. $F_{y\to y^{\prime}}=\ln j_{y\to
y^{\prime}}-\ln j_{y^{\prime}\to y}$ is _thermodynamic force_ or _affinity_ of
$y\to y^{\prime}\in\mathcal{R}^{+}$, which indicates the total entropy change
in the network and its surroundings. $A_{y\to y^{\prime}}=j_{y\to
y^{\prime}}+j_{y^{\prime}\to y}$ is the quantity called _dynamical activity_ ,
which is the total flow in forward and backward reactions. The dynamics on the
reversible CRN can be rewritten as a linear summation of stoichiometric
vectors $y^{\prime}-y$ for $\to y^{\prime}\in\mathcal{R}^{+}$ with
coefficients of net currents $J$,

$\displaystyle\differential{}_{t}x$ $\displaystyle=\sum_{y\to
y^{\prime}\in\mathcal{R}^{+}}(y^{\prime}-y)J_{y\to y^{\prime}}.$ (10)

##  III species-transitional representation of stoichiometric dynamics.

In this section, we provide an insight into CRN dynamics with newly proposed
physical quantities. The quantities represent concentration transitions
between chemical species in discrete states. A CRN uniquely induces a directed
graph with chemical species as nodes and transitions between species as edges.

###  III.1 Chemical Transitions on Directed Graphs.

Our first result is the probability-like representation of the CRN dynamics
that is followed by the scaled concentration. Suppose $\ell$ to be a positive
conservation law of the CRN. For the conservation law $\ell$, we call
$\ell\circ x$ the $\ell$-scaled concentration of the ordinary concentration
vector $x\in X$ of the CRN dynamics and denote as $x^{\ell}=\ell\circ
x\in\mathbb{R}_{>0}^{\mathcal{I}}$. Since the sum over $i\in\mathcal{I}$ is
conserved, $x^{\ell}$ can be regarded as a probability disrtibution on
$\mathcal{I}$, which follows a deterministic time evolution. For the
concentration following the dynamics of the equation 8, the time evolution of
$\ell$-scaled concentration $x^{\ell}$ satisfies the following equation:

$\displaystyle\differential{}_{t}x^{\ell}_{i}=\sum_{j\in\mathcal{I}_{\ell}\setminus\quantity{i}}\mathcal{T}_{ij}^{\ell}-\sum_{j\in\mathcal{I}_{\ell}\setminus\quantity{i}}\mathcal{T}_{ji}^{\ell}.$
(11)

$\mathcal{T}_{ji}^{\ell}$ is the concentration-dependent _transition_ from $i$
to $j$ with respect to the law $\ell$ defined as

$\displaystyle\mathcal{T}_{ij}^{\ell}=\sum_{\begin{subarray}{c}y\to
y^{\prime}\in\mathcal{R}\end{subarray}}\norm{y^{\ell}}^{-1}{y^{\prime}}_{i}^{\ell}y_{j}^{\ell}j_{y\to
y^{\prime}}(x_{t}),$ (12)

where $\norm{\cdot}$ is $1$-norm and $y^{\ell}=\ell\circ y$
($y^{\ell}=\ell\circ y^{\prime}$) is the scaled reactant (product)
coefficients. Note that $\norm{y^{\ell}}=\norm{y^{\prime\ell}}$ holds for the
two total coefficients of any reaction $y\to y^{\prime}\in\mathcal{R}$ because
$\ell$ is a conservation law. The equation 11 can explicitly be rewritten as a
continuity-equation-like form [5, 18] with the connection matrix $\mathbb{B}$
of the S-graph $(\mathcal{I},\mathcal{E})$ as

$\displaystyle\differential_{t}x^{\ell}=-\mathbb{B}\mathcal{T}^{\ell}.$ (13)

In the equation 13, $\mathcal{T}^{\ell}$ is identified with the vector in
$\mathcal{R}_{\geq 0}^{\mathcal{E}}$; $(\mathcal{T}_{ij}^{\ell})_{i\leftarrow
j\in\mathcal{E}}$, which is justified by the conformity of S-graph edge sets
and transitions 29.

This equation 11 is similar to the master equation, but generally this
equation alone cannot describe the time evolution of $x^{\ell}$ because the
equation is not closed with respect to $x^{\ell}$ and not equivalent to the
original one 8 in general. If the law
$\ell\not\in\mathcal{L}_{\mathcal{N}}^{+}$, where $\ell$’s support
$\mathcal{I}_{\ell}=\mathrm{supp}(\ell)$ is not equal to the original set
$\mathcal{I}$, the transition matrix $\mathcal{T}^{\ell}$ depends not only on
the $x^{\ell}$ (i.e. $(x_{i})_{i\in\mathcal{I}_{\ell}}$) but also on the whole
concentration $x=(x_{i})_{i\in\mathcal{I}}$, and avoiding this kind of
argument is the reason why $\ell$ is assumed to be positive (or
$\mathcal{I}_{\ell}=\mathcal{I}$) in the previous paragraph.

It can be verified by a simple calculation that the formula 11 is correct. The
sum of $\mathcal{T}_{ij}^{\ell}$ and $\mathcal{T}_{ji}^{\ell}$ over
$j\in\mathcal{I}$ are respectively,

$\displaystyle\sum_{j\in\mathcal{I}}T_{ij}^{\ell}$ $\displaystyle=\sum_{y\to
y^{\prime}\in\mathcal{R}}{y_{i}^{\prime}}^{\ell}j_{y\to y^{\prime}}(x_{t}),$
(14) $\displaystyle\sum_{j\in\mathcal{I}}T_{ji}^{\ell}$
$\displaystyle=\sum_{y\to y^{\prime}\in\mathcal{R}}y_{i}^{\ell}j_{y\to
y^{\prime}}(x_{t}),$ (15)

because of
$\norm{y^{\ell}}^{-1}\sum_{j}y^{\ell}_{j}=\norm{y^{\prime\ell}}^{-1}\sum_{j}{y^{\prime}}^{\ell}_{j}=1$.
by taking the difference of the above two, 14 and 15, 11 is proved as follows.

$\displaystyle\sum_{j\in\mathcal{I}}\mathcal{T}_{ij}^{\ell}-\sum_{j\in\mathcal{I}}\mathcal{T}_{ij}^{\ell}$
$\displaystyle=\sum_{y\to
y^{\prime}\in\mathcal{R}}({y_{i}^{\prime}}^{\ell}-y_{i}^{\ell})j_{y\to
y^{\prime}}(x_{t})$ (16) $\displaystyle=\ell\circ\sum_{y\to
y^{\prime}\in\mathcal{R}}({y_{i}^{\prime}}-y_{i})j_{y\to y^{\prime}}(x_{t})$
(17) $\displaystyle=\ell\circ\differential_{t}x_{i}$ (18)
$\displaystyle=\differential_{t}x^{\ell}_{i}.$ (19)

The representation 13 is also proved as

$\displaystyle\sum_{j\in\mathcal{I}}\mathcal{T}_{ij}^{\ell}-\sum_{j\in\mathcal{I}}\mathcal{T}_{ij}^{\ell}$
$\displaystyle=\sum_{e\in\mathcal{E}:h(e)=i}\mathcal{T}_{e}^{\ell}-\sum_{e\in\mathcal{E}:t(e)=i}\mathcal{T}_{e}^{\ell}$
(20)
$\displaystyle=-\quantity(\sum_{e\in\mathcal{E}:t(e)=i}\mathcal{T}_{e}^{\ell}-\sum_{e\in\mathcal{E}:h(e)=i}\mathcal{T}_{e}^{\ell})$
(21) $\displaystyle=-\mathbb{B}\mathcal{T}^{\ell}.$ (22)

The transformation is derived from the compatibility of $\mathcal{E}$ and
$\mathcal{T}^{\ell}$, see 29.

###  III.2 Description of transition by sampling probability.

We describe the stochastic aspects of this transition representation of
dynamics. In this subsection, the conservation law is assumed to be a one
vector for simplicity $\ell=1_{\mathcal{I}}$, where
$\norm{y}=\norm{y^{\prime}}$ holds for any reaction $y\to y^{\prime}$.
Replacing $y$ by $y^{\ell}$ and $y^{\prime}$ by ${y^{\prime}}^{\ell}$, a
similar argument holds for a general conservation law $\ell$. The transition
from $j$ to $i$ is written as the following form,

$\displaystyle\mathcal{T}_{ij}$ $\displaystyle=\sum_{\begin{subarray}{c}y\to
y^{\prime}\in\mathcal{R}\end{subarray}}\mathbb{P}(i\leftarrow j|y\to
y^{\prime})\cdot\norm{y}j_{y\to y^{\prime}},$ (23)

where $\mathcal{T}^{1_{\mathcal{I}}}$ is denoted as $\mathcal{T}$ and

$\displaystyle\mathbb{P}(i\leftarrow j|y\to
y^{\prime})=\frac{y^{\prime}_{i}}{\norm{y^{\prime}}}\cdot\frac{y_{j}}{\norm{y}}.$
(24)

$j_{y\to y^{\prime}}(x_{t})\differential t$ can be interpreted as the
probability that the reaction $y\to y^{\prime}$ occurs in the interval
$[t,t+\differential t]$ under the assumption that reactions occur at most once
in the short-time $\differential t$. $\mathbb{P}(i\leftarrow j|y\to
y^{\prime})$ is proportional to $y^{\prime}_{i}\cdot y_{j}$ where
$y_{i}^{\prime}$ is $i$’s product coefficient and $y_{j}$ is $i$’s reactant
coefficient of $y\to y^{\prime}$, and indicates the chemical species sampling
with the conditional probability that the state transitions from $j$ (at $t$)
to $i$ (at $t+\differential t$) under the condition that the reaction $y\to
y^{\prime}$ occurs in $[t,t+\differential t]$. Thus the joint probability

$\displaystyle\begin{split}&\mathbb{P}(t+\differential t,i;t,j;y\to
y^{\prime})\\\ =&\mathbb{P}(i\leftarrow j|y\to y^{\prime})\cdot j_{y\to
y^{\prime}}(x_{t})\differential t.\end{split}$ (25)

is defined such that a molecule is in state $j$ at $t$ & in state $i$ at
$t+\differential t$ and $y\to y^{\prime}$ occurs in $[t,t+\differential t]$.
Then the transitions is rewritten as

$\displaystyle\mathcal{T}_{ij}\differential{t}$
$\displaystyle=\mathbb{E}_{y\to y^{\prime}\sim\mathbb{P}(t+\differential
t,i;t,j;y\to y^{\prime})}(\norm{y}),$ (26)

which is the expected value of the total product/reactant coefficient for
reactions $y\to y^{\prime}$. The dynamics of $x$ is

$\displaystyle\differential x_{i}=\sum_{y\to
y^{\prime}\in\mathcal{R}}\norm{y}\sum_{j\in\mathcal{I}}I_{ij},$ (27)

where $\differential x_{i}$ is the small concentration change from $t$ to
$t+\differential t$ and $I_{ij}$ indicates the net probablity flow from $j$ to
$i$ or the assymetry in the two-time joint probablities,

$\displaystyle\begin{split}I_{ij}=&\mathbb{P}(t+\differential t,i;t,j;y\to
y^{\prime})\\\ -&\mathbb{P}(t+\differential t,j;t,i;y\to
y^{\prime}).\end{split}$ (28)

The equation 27 can be interpreted that the change in $x$ is driven by the
coefficient-weighted asymmetries in the two-time joint probabilities of each
reaction.

###  III.3 Reversibility of CRN and species transitions.

We discuss the reversibility of the CRN structure (or its species graph) and
its relation to our transition quantity in this subsection. First, we see that
the existence of a reaction producing $j$ from $i$ is equivalent to the $ij$
component of the transition matrix $\mathcal{T}^{\ell}$ being not vanishing.
Derived from this, we study the relationship between the reversibility of CRNs
and species graphs and the transitions.

For a positive law $\ell$ and its transition $\mathcal{T}^{\ell}$, the
following two condition is equivalent.

  1. 1.

$\exists y\to y^{\prime}\in\mathcal{R}$ such that $y_{i}^{\prime}y_{j}>0$.

  2. 2.

$\mathcal{T}_{ij}^{\ell}>0$.

The proof is given as follows. Assuming the condition 1 and denoting such
reaction as $\tilde{y}\to\tilde{y}^{\prime}$, the positiveness can be found
from the expression of $\mathcal{T}^{\ell}_{ij}$ as
$\mathcal{T}^{\ell}_{ij}\geq\norm{\tilde{y}^{\ell}}^{-1}\tilde{y}^{\prime\ell}_{i}\tilde{y}_{j}^{\ell}j_{\tilde{y}\to\tilde{y}^{\prime}}(x_{t})>0$,
using $\tilde{y}_{i}^{\prime}\tilde{y}_{j}>0$ and
$j_{\tilde{y}\to\tilde{y}^{\prime}}(x)>0\,(\forall
x\in\mathbb{R}_{>0}^{\mathcal{I}})$. Assuming that there are no reactions that
satisfies $y_{i}y_{j}>0$, then the transition is always zero:
$\mathcal{T}_{ij}^{\ell}=0$, and the proof is completed by its contraposition.

Our transition matrix is compatible with the species graph in the following
manner. The edge set of species graph is characterized by the transition
$\mathcal{T}^{\ell}$ as

$\displaystyle\mathcal{E}=\quantity{i\leftarrow
j\mid\mathcal{T}_{ij}^{\ell}>0},$ (29)

which is checked by the definition of the edge set. As a corollary, we note
that the reversibility of edges in the species graph can be expressed in terms
of reversibility with respect to the transition T, that is, the equivalence of
the following two conditions: For a positive law $\ell$ and its transition
$\mathcal{T}^{\ell}$, the following two condition is equivalent.

  1. 1.

$i\leftarrow j\in\mathcal{E}$ is reversible.

  2. 2.

$\mathcal{T}_{ij}^{\ell}>0\Leftrightarrow\mathcal{T}_{ji}^{\ell}>0$.

We comment on the reversibility of the original CRN and its species graph. For
a reversible conservative CRN, its S-graph is reversible, which is because
$y_{i}^{\prime}y_{j}>0$ for $y\to y^{\prime}\in\mathcal{R}$ implies
$y_{j}y_{i}^{\prime}>0$ for $y^{\prime}\to y\in\mathcal{R}$ in reversible
CRNs. However, the converse does not hold in general, i.e. a CRN with the
reversible S-graph is not necessarily reversible, which can be checked by the
CRN with the following reactions,

$\displaystyle A$ $\displaystyle\rightharpoonup C,$ (30) $\displaystyle B$
$\displaystyle\rightharpoonup C,$ (31) $\displaystyle 2C$
$\displaystyle\rightharpoonup A+B.$ (32)

This CRN has $1$ as a positive conservation law and its transition is positive
on the edges

$\displaystyle\mathcal{E}=\quantity{C\leftarrow A,A\leftarrow C,B\leftarrow
A,A\leftarrow B},$ (33)

which constitutes a reversible S-graph.

##  IV Thermodynamical quantities on species graphs.

The distribution on nodes of an reversible graph and its dynamics has been
studied in the stochastic thermodynamics [19]. Our formulation enables us to
treat the CRN dynamics in the style of stochastic thermodynamics. For
simplicity, consider a reversible CRN whose conservation law graphs are
undirected, which means that the graphs contain any reverse edges.

We can analogously define some themodynamic quantities for the dynamics on the
species graph $\mathcal{G}=(\mathcal{I},\mathcal{E})$: specieswise net
currents, forces and activities,

$\displaystyle\mathcal{J}_{ij}^{\ell}$ $\displaystyle=$
$\displaystyle\mathcal{T}_{ij}^{\ell}-\mathcal{T}_{ji}^{\ell},$ (34)
$\displaystyle\mathcal{F}_{ij}^{\ell}$ $\displaystyle=$
$\displaystyle\ln\mathcal{T}_{ij}^{\ell}-\ln\mathcal{T}_{ji}^{\ell},$ (35)
$\displaystyle\mathcal{A}_{ij}^{\ell}$ $\displaystyle=$
$\displaystyle\mathcal{T}_{ij}^{\ell}+\mathcal{T}_{ji}^{\ell}.$ (36)

We call these physical quantities _specieswise thermodynamic quantities_ of
the conservative CRN, to be distinguished from original reaction-related flows
$J$, forces $F$ and activities $A$. The three types of quantities are
expressed with the summations of the conventional thermodynamical quantities
over the reactions ${y\to y^{\prime}}$ as

$\displaystyle\mathcal{J}_{ij}^{\ell}$ $\displaystyle=$
$\displaystyle\sum_{y\to
y^{\prime}\in\mathcal{R}^{+}}\ell_{i}\ell_{j}\norm{y^{\ell}}^{-1}(y^{\prime}_{i}y_{j}-y^{\prime}_{j}y_{i})J_{y\to
y^{\prime}},$ (37) $\displaystyle\mathcal{F}_{ij}^{\ell}$ $\displaystyle=$
$\displaystyle\ln\frac{\sum_{y\to
y^{\prime}\in\mathcal{R}}\norm{y^{\ell}}^{-1}y_{i}^{\prime}y_{j}j_{y\to
y^{\prime}}}{\sum_{y\to
y^{\prime}\in\mathcal{R}}\norm{y^{\ell}}^{-1}y_{i}^{\prime}y_{j}j_{y^{\prime}\to
y}},$ (38) $\displaystyle\mathcal{A}_{ij}^{\ell}$ $\displaystyle=$
$\displaystyle\sum_{y\to
y^{\prime}\in\mathcal{R}^{+}}\ell_{i}\ell_{j}\norm{y^{\ell}}^{-1}(y^{\prime}_{i}y_{j}+y^{\prime}_{j}y_{i})A_{y\to
y^{\prime}}.$ (39)

In the case of the reversible S-graph, using the net current
$\mathcal{J}^{\ell}$, the dynamics can be written as a form of the continuity
equation [5, 18], that is,

$\displaystyle\differential_{t}x^{\ell}=-\mathbb{B}_{\mathcal{E}^{+}}\mathcal{J}^{\ell}.$
(40)

$\mathcal{E}^{+}\subseteq\mathcal{E}$ is a set of forward edge set which
satisfies $i\leftarrow j\in\mathcal{E}^{+}\Leftrightarrow j\leftarrow
i\not\in\mathcal{E}^{+}$ and $\mathbb{B}_{\mathcal{E}^{+}}$ is the restriction
of the S-graph’s incidence matrix $\mathbb{B}$ to $\quantity{0,\pm
1}^{\mathcal{I}\times\mathcal{E}^{+}}$. $\mathcal{J}^{\ell}$ is identified
with the vector $(\mathcal{J}_{ij}^{\ell})_{i\leftarrow j\in\mathcal{E}^{+}}$.

Note that some qualities or inequalities hold between thermodynamic quantities
on the graph and ones on the CRN, which give physical meanings of
$\mathcal{J}_{ij}^{\ell},\mathcal{F}_{ij}^{\ell},\mathcal{A}_{ij}^{\ell}$.

###  IV.1 Specieswise Forces

Assume that there exists at least one reaction which has $i$ as a product and
$j$ as a reactant. $\mathcal{F}_{ij}^{\ell}$ is regarded as a quantity
representative of thermodynamic forces of reactions which have $j$ as a
reactant and $i$ as a product because the inequality below holds for any $y\to
y^{\prime}\in\mathcal{R}_{ij}$,

$\displaystyle\min_{y\to y^{\prime}\in\mathcal{R}_{ij}}{F_{y\to
y^{\prime}}}\leq\mathcal{F}_{ij}^{\ell}\leq\max_{y\to
y^{\prime}\in\mathcal{R}_{ij}}{F_{y\to y^{\prime}}},$ (41)

where $\mathcal{R}^{ij}$ is a set of reactions which have $j$ as a reactant
and $i$ as a product,

$\displaystyle\mathcal{R}_{ij}=\quantity{y\to y^{\prime}\in\mathcal{R}\mid
y^{\prime}_{j}>0\text{ and }y_{i}>0}.$ (42)

From the first assumption, this reaction subset is non-empty.

The inequality derives from the expression of $\mathcal{F}_{ij}^{\ell}$ by the
reaction set $\mathcal{R}_{ij}$,

$\displaystyle\mathcal{F}_{ij}^{\ell}=\ln\frac{\sum_{y\to
y^{\prime}\in\mathcal{R}_{ij}}\norm{y^{\ell}}^{-1}y_{i}^{\prime}y_{j}j_{y\to
y^{\prime}}}{\sum_{y\to
y^{\prime}\in\mathcal{R}_{ij}}\norm{y^{\ell}}^{-1}y_{i}^{\prime}y_{j}j_{y^{\prime}\to
y}}.$ (43)

Each term $\norm{y\to y^{\prime}}_{\ell}^{-1}y_{i}^{\prime}y_{j}F_{y\to
y^{\prime}}$ is positive (non-zero), thus the bound of
$\mathcal{F}_{ij}^{\ell}$ 41 is obtained from the inequality
$\min_{i}a_{i}/b_{i}\leq\sum_{i}a_{i}/\sum_{i}b_{i}\leq\max_{i}a_{i}/b_{i}\,(a_{i},b_{i}>0)$
and the monotonicity of the logarithmic function $\ln(\cdot)$. From the
derivation, we know that the equality of 41 holds if and only if the original
forces $F_{y\to y^{\prime}}$ take the same value for each $y\to
y^{\prime}\in\mathcal{R}_{ij}$.

Finally, we mention the asymptotic behaviour of $F_{ij}^{\ell}$. Consider a
situation where a certain reaction kinetics $j_{y\to y^{\prime}}$ is dominant;
$j_{y\to y^{\prime}}\gg j_{\tilde{y}\to\tilde{y}^{\prime}}$ for any
$\tilde{y}\to\tilde{y}^{\prime}\in\mathcal{R}_{ij}\setminus\quantity{y\to
y^{\prime}}$. In this case, the specieswise force will be approximately the
thermodynamic force of the reaction $y\to y^{\prime}$, because

$\displaystyle\begin{split}\mathcal{F}_{ij}^{\ell}&\simeq\ln\frac{\norm{y^{\ell}}^{-1}y_{i}^{\prime}y_{j}j_{y\to
y^{\prime}}}{\norm{y^{\ell}}^{-1}y_{i}^{\prime}y_{j}j_{y^{\prime}\to y}}\\\
&=\ln\frac{j_{y\to y^{\prime}}}{j_{y^{\prime}\to y}}=F_{y\to
y^{\prime}}.\end{split}$ (44)

###  IV.2 Specieswise Activities

Figure 2:  The time evolutions of (a) concentrations of species $X$ and $Y$,
(b) kinetics $j_{m}\,(m=\pm 1,\pm 2,\pm 3)$ and (c) thermodynamic forces
$F_{m}\,(m=1,2,3)$. The parameters and initial conditions are
$k_{1}^{+}=1\times 10^{-3},k_{1}^{-}=k_{2}^{+}=k_{2}^{-}=1,k_{3}^{+}=1\times
10^{-2},k_{3}^{-}=1\times 10^{-4},[X]_{0}=1,[Y]_{0}=6,[A]_{0}=[B]_{0}=1\times
10^{3}$ , which is the same as in [20].

Next we see the relationship between the specieswise dynamical activities on
the induced graph and the covariance of the chemical dynamics. The quantity

$\displaystyle\tilde{D}_{ij}^{\ell}=\sum_{y\to
y^{\prime}\in\mathcal{R}^{+}}((y^{\prime}-y)^{\ell})_{i}((y^{\prime}-y)^{\ell})_{j}A_{y\to
y^{\prime}}$ (45)

indicates the (scaled) covariance between $i$ and $j$ of the chemical Fokker-
Plank equation [21].

$\mathcal{A}_{ij}^{\ell}$ gives a lower bound for the covariance
$\tilde{D}_{ij}^{\ell}$:

$\displaystyle\tilde{D}_{ij}^{\ell}\geq-\quantity(\max_{y\to
y^{\prime}\in\mathcal{R}^{+}}\norm{y^{\ell}})\cdot\mathcal{A}_{ij}^{\ell}.$
(46)

The proof is completed by a simple transformation of the equation 45 and the
process is shown below.

$\displaystyle\tilde{D}_{ij}^{\ell}$ $\displaystyle=\sum_{y\to
y^{\prime}\in\mathcal{R}^{+}}\quantity{(y^{\prime}_{i}y_{i}+y^{\prime}_{j}y_{j})-(y^{\prime}_{i}y_{j}+y^{\prime}_{j}y_{i})}A_{y\to
y^{\prime}}$ (47) $\displaystyle\geq-\sum_{y\to
y^{\prime}\in\mathcal{R}^{+}}(y^{\prime}_{i}y_{j}+y^{\prime}_{j}y_{i})A_{y\to
y^{\prime}}$ (48) $\displaystyle=-\quantity(\max_{y\to
y^{\prime}\in\mathcal{R}^{+}}\norm{y^{\ell}})\cdot\mathcal{A}_{ij}^{\ell}.$
(49)

The equality holds if and only if the total coefficient
$\norm{y^{\ell}}=\norm{y^{\prime\ell}}$ takes the same value for any reaction
$y\to y^{\prime}\in\mathcal{R}^{+}$ and no reaction exists where $i,j$ are
both reactants or products.

The asymptotic analysis differs slightly from the force case. In order for the
reaction $y\to y^{\prime}$ to asymptotically coincide with the dominant case
$j_{y\to y^{\prime}}\gg j_{\tilde{y}\to\tilde{y}^{\prime}}$, an additional
condition is required such that this reaction must be of maximum order
throughout the entire reaction system:
$\norm{y^{\ell}}=\max_{\tilde{y}\to\tilde{y}^{\prime}\in\mathcal{R}_{ij}}\norm{\tilde{y}^{\ell}}$.

##  V numerical experiment

Through a conservative CRN example, we confirm the behaviour of our defined
physical quantities. The Brusselator is used for numerical experiment, which
is a well-known autocatalytic reaction model showing chemical oscillation
[22]. We use the damped version of Brusselator system in [20] which has $4$
species $\mathcal{I}=\quantity{A,B,X,Y}$ and $6$ reactions whose reaction
equations are represented as

$\displaystyle A$ $\displaystyle\rightleftharpoons$ $\displaystyle X,$ (50)
$\displaystyle 2X+Y$ $\displaystyle\rightleftharpoons$ $\displaystyle 3X,$
(51) $\displaystyle X+B$ $\displaystyle\rightleftharpoons$ $\displaystyle
Y+A.$ (52)

Formally, the reaction set is defined as

$\displaystyle\mathcal{R}=\quantity{\begin{aligned}
r_{1}&=(0,0,1,0)^{\top}\to(1,0,0,0)^{\top},\\\
r_{2}&=(2,1,0,0)^{\top}\to(3,0,0,0)^{\top},\\\
r_{3}&=(1,0,0,1)^{\top}\to(0,1,1,0)^{\top},\\\
r_{-1}&=(1,0,0,0)^{\top}\to(0,0,1,0)^{\top},\\\
r_{-2}&=(3,0,0,0)^{\top}\to(2,1,0,0)^{\top},\\\
r_{-3}&=(0,1,1,0)^{\top}\to(1,0,0,1)^{\top}\end{aligned}}$ (53)

where $\mathbb{Z}_{\geq 0}^{\mathcal{I}}$ is identified with $\mathbb{Z}_{\geq
0}^{\absolutevalue{\mathcal{I}}}=\mathbb{Z}_{\geq 0}^{4}$. The kinetics $j$ is
assumed to be mass-actional as follows.

$\displaystyle j_{1}$ $\displaystyle=k_{1}^{+}[A],$ (54) $\displaystyle j_{2}$
$\displaystyle=k_{2}^{+}[X]^{2}[Y],$ (55) $\displaystyle j_{1}$
$\displaystyle=k_{3}^{+}[X][B],$ (56) $\displaystyle j_{-1}$
$\displaystyle=k_{1}^{-}[X],$ (57) $\displaystyle j_{-2}$
$\displaystyle=k_{2}^{-}[X]^{3},$ (58) $\displaystyle j_{-3}$
$\displaystyle=k_{3}^{-}[Y][A],$ (59)

where the notation $[Z]$ means the concentration of the chemical species
$Z\in\mathcal{I}$. The conventional flux and driving forces are defined as
$J_{m}=j_{m}-j_{-m}$ and $F_{m}=\ln j_{m}-\ln j_{-m}$ for $m=1,2,3$. The
system exhibits damped oscillations in concentration, kinetics and
thermodynamic driving forces in 2.

The S-graph $(\mathcal{I},\mathcal{E})$ of this system has an edge set

$\displaystyle\mathcal{E}=\quantity{\begin{aligned} &X\leftarrow Y,Y\leftarrow
X,\\\ &X\leftarrow A,A\leftarrow X,\\\ &Y\leftarrow B,B\leftarrow Y,\\\
&B\leftarrow A,A\leftarrow B.\end{aligned}},$ (60)

whose reversibility is inherited from the original CRN. The transition
$\mathcal{T}$ for a conservation law $1$ has elements with

$\displaystyle\begin{split}\mathcal{T}_{XY}&=\quantity(k_{2}^{+}[X]^{2}+k_{3}^{-}[A]/2)[Y],\\\
\mathcal{T}_{YX}&=\quantity(k_{2}^{-}[X]^{2}+k_{3}^{+}[B]/2)[X],\\\
\mathcal{T}_{XA}&=\quantity(k_{1}^{-}+k_{3}^{-}[Y]/2)[A],\\\
\mathcal{T}_{AX}&=\quantity(k_{1}^{+}+k_{3}^{+}[B]/2)[X],\\\
\mathcal{T}_{XB}&=0,\\\ \mathcal{T}_{BX}&=0,\\\ \mathcal{T}_{YA}&=0,\\\
\mathcal{T}_{AY}&=0,\\\ \mathcal{T}_{YB}&=\quantity(k_{3}^{+}[X]/2)[B],\\\
\mathcal{T}_{BY}&=\quantity(k_{3}^{-}[A]/2)[Y],\\\
\mathcal{T}_{AB}&=\quantity(k_{3}^{+}[X]/2)[B],\\\
\mathcal{T}_{BA}&=\quantity(k_{3}^{-}[Y]/2)[A].\end{split}$ (61)

The behaviour of the specieswise force $\mathcal{F}_{XY}$ and activity
$\mathcal{A}_{XY}$ from $Y$ to $X$ and their bounds are observed. As seen in
the subsection IV.1, $\mathcal{F}_{XY}$ could take values in the interval

$\displaystyle\Lambda_{XY}=\quantity[\min_{r\in\quantity{r_{2},r_{-3}}}F_{r},\max_{r\in\quantity{r_{2},r_{-3}}}F_{r}]$
(62)

at each time. The time evolutions of $\mathcal{F}_{XY}$ and the possible
interval $\Lambda_{XY}$ are showed in Fig 3. In this simulation, we can
comfirm that the area is actually $\Lambda_{XY}=[F_{r_{-2}},F_{r_{3}}]$ for
any time $t$. The force and its upper bound is quite close at certain times,
which is interpreted by the asymtotic approximation discussed in the last of
the subsection IV.1, together with the time series of the kinetics Fig 2 (b).
The time evolutions of the covariance $D_{XY}$ between $X,Y$ and the lower
bound $-3\mathcal{A}_{XY}$ is seen in the Fig 4. Again, it can be seen that
the two physical quantities are approximately coincident in situations where
the largest reaction order, reaction 2, is dominant.

Figure 3: The time evolutions of the specieswise force $\mathcal{F}_{XY}$ from
$Y$ to $X$ and its possible interval $\Lambda_{XY}$. The upper and lower line
in the figure represent $F_{r_{2}}$ and $F_{r_{-3}}$, respectively. Figure 4:
The time evolutions of the covariance $D_{XY}$ and the bound by the
specieswise activities $\mathcal{A}_{XY}$.

By observing the numerical results, two properties can be identified for our
defined quantities below. First, we argue that the physical quantities we have
defined are properly interpretable. Indeed, we have so far confirmed both
theoretically and numerically that specieswise forces and activities are
related to conventional driving forces and covariance. In particular, for the
former, Figure 2 (a) shows that $X$ increases and $Y$ decreases in the long
term, which is consistent with a positive $Y$ to $X$ force for most of the
time $0\leq t\leq 60$. Secondly, the specieswise force $\mathcal{F}_{XY}$ and
activity $\mathcal{A}_{XY}$ is found to inherit the oscillatory nature of the
system in the simulation.

##  VI conclusion and discussion

This paper introduced a novel framework for analyzing CRN dynamics by defining
”transitions” as a new physical quantity and reinterpreting the dynamics on a
species graph. The dynamics on this graph were proven conservative by scaling
concentrations according to the network’s conservation laws. Thermodynamic
quantities like forces and activities were then formulated at the species
level based on these transitions and the relation between the specieswise
thermodynamic quantities and conventional thermodynamic quantities were
investigated. When applied numerically to the Brusselator oscillating reaction
system, the specieswise quantities exhibited intuitive behavior aligned with
theoretical predictions. This transition-based graphical approach provides an
alternative perspective for understanding the thermodynamics of complex CRNs
with potential for enhanced modeling and analysis. Further research is merited
into the broader implications and applications of this specieswise
thermodynamic formulation.

One potential limitation of our framework is that the newly introduced concept
of ”transitions” may not precisely mirror the actual physical transformations
occurring at the molecular level. To illustrate, consider a second-order
chemical reaction $A+B\to A^{\prime}+B^{\prime}$ through which the physical
chemistry explains the molecules $A$ and $B$ undergo structural changes to
become $A^{\prime}$ and $B^{\prime}$ respectively. Following this depiction,
the reaction should solely contribute to the transitions $A\to A^{\prime}$ and
$B\to B^{\prime}$. However, our methodology additionally accounts for the
transitions $A\to B^{\prime}$ and $B\to A^{\prime}$, which do not directly
correspond to the molecular rearrangements. Despite this potential disconnect
from microscopic mechanisms, our transition-based approach offers the
significant advantage of being systematically computable from the system’s
stoichiometric and kinetic information, lending itself to practical analysis
workflows.

Looking ahead, we anticipate our framework enabling deeper dynamical analysis
of reaction networks, particularly metabolic systemsand others central to
cellular processes. Theoretically, exploring how stochastic thermodynamics
laws manifest in CRNs under our specieswise formulation is intriguing. For
instance, it is an interesting topic that how the relations between
correlations and driving forces at a steady state for Markov jump system [23]
potentially transform in the CRN context. As this specieswise approach
develops, prospects include analyzing complex network dynamics and bridging
disparate fields through a unified reaction network lens, with potential for
significant discoveries.

## References

  * Orth  _et al._ [2010] J. Orth, I. Thiele, and B. Palsson, What is flux balance analysis?, Nature biotechnology 28, 245 (2010).
  * van Zon  _et al._ [2007] J. S. van Zon, D. K. Lubensky, P. R. H. Altena, and P. R. ten Wolde, An allosteric model of circadian kaic phosphorylation, Proceedings of the National Academy of Sciences 104, 7420 (2007), https://www.pnas.org/doi/pdf/10.1073/pnas.0608665104 .
  * Avanzini  _et al._ [2024] F. Avanzini, M. Bilancioni, V. Cavina, S. D. Cengio, M. Esposito, G. Falasco, D. Forastiere, N. Freitas, A. Garilli, P. E. Harunari, V. Lecomte, A. Lazarescu, S. G. M. Srinivas, C. Moslonka, I. Neri, E. Penocchio, W. D. Piñeros, M. Polettini, A. Raghu, P. Raux, K. Sekimoto, and A. Soret, Methods and conversations in (post)modern thermodynamics, SciPost Phys. Lect. Notes , 80 (2024).
  * Yoshimura and Ito [2021a] K. Yoshimura and S. Ito, Thermodynamic uncertainty relation and thermodynamic speed limit in deterministic chemical reaction networks, Phys. Rev. Lett. 127, 160601 (2021a).
  * Kobayashi  _et al._ [2023] T. J. Kobayashi, D. Loutchko, A. Kamimura, S. A. Horiguchi, and Y. Sughiyama, Information geometry of dynamics on graphs and hypergraphs, Information Geometry 10.1007/s41884-023-00125-w (2023).
  * Himeoka  _et al._ [2024] Y. Himeoka, S. A. Horiguchi, and T. J. Kobayashi, Stoichiometric rays: A simple method to compute the controllable set of enzymatic reaction systems (2024), arXiv:2403.02169 [physics.bio-ph] .
  * Watts and Strogatz [1998] D. J. Watts and S. H. Strogatz, Collective dynamics of ‘small-world’ networks, Nature 393, 440–442 (1998).
  * Newman [2018] M. Newman,  _Networks_ (Oxford University Press, 2018).
  * Angeli  _et al._ [2009] D. Angeli, P. De Leenheer, and E. Sontag, Graph-theoretic characterizations of monotonicity of chemical networks in reaction coordinates, Journal of Mathematical Biology 61, 581–616 (2009).
  * Wagner and Fell [2001] A. Wagner and D. A. Fell, The small world inside large metabolic networks, Proceedings of the Royal Society of London. Series B: Biological Sciences 268, 1803–1810 (2001).
  * Klamt  _et al._ [2009] S. Klamt, U.-U. Haus, and F. Theis, Hypergraphs and cellular networks, PLOS Computational Biology 5, 1 (2009).
  * Himeoka and Kaneko [2016] Y. Himeoka and K. Kaneko, Enzyme oscillation can enhance the thermodynamic efficiency of cellular metabolism: consequence of anti-phase coupling between reaction flux and affinity, Physical Biology 13, 026002 (2016).
  * Kim and Hyeon [2021] P. Kim and C. Hyeon, Thermodynamic optimality of glycolytic oscillations, The Journal of Physical Chemistry B 125, 5740 (2021), pMID: 34038120, https://doi.org/10.1021/acs.jpcb.1c01325 .
  * Fontanil and Mendoza [2022] L. Fontanil and E. Mendoza, Common complexes of decompositions and complex balanced equilibria of chemical reaction networks, MATCH Communications in Mathematical and in Computer Chemistry 87, 329 (2022).
  * Feinberg [2019] M. Feinberg,  _Foundations of Chemical Reaction Network Theory_ (Springer International Publishing, 2019).
  * Müller  _et al._ [2022] S. Müller, C. Flamm, and P. F. Stadler, What makes a reaction network ”chemical”?, J Cheminform 14, 63 (2022).
  * De Martino  _et al._ [2014] A. De Martino, D. De Martino, R. Mulet, and A. Pagnani, Identifying all moiety conservation laws in genome-scale metabolic networks, PLOS ONE 9, 1 (2014).
  * Grady and Polimeni [2010] L. J. Grady and J. R. Polimeni,  _Discrete Calculus_ (Springer London, 2010).
  * Seifert [2012] U. Seifert, Stochastic thermodynamics, fluctuation theorems and molecular machines, Reports on Progress in Physics 75, 126001 (2012).
  * Yoshimura and Ito [2021b] K. Yoshimura and S. Ito, Information geometric inequalities of chemical thermodynamics, Phys. Rev. Res. 3, 013175 (2021b).
  * Gillespie [2000] D. T. Gillespie, The chemical Langevin equation, The Journal of Chemical Physics 113, 297 (2000), https://pubs.aip.org/aip/jcp/article-pdf/113/1/297/19042305/297_1_online.pdf .
  * Lefever  _et al._ [1988] R. Lefever, G. Nicolis, and P. Borckmans, The brusselator: it does oscillate all the same, Journal of the Chemical Society, Faraday Transactions 1: Physical Chemistry in Condensed Phases 84, 1013 (1988).
  * Ohga  _et al._ [2023] N. Ohga, S. Ito, and A. Kolchinsky, Thermodynamic bound on the asymmetry of cross-correlations, Phys. Rev. Lett. 131, 077101 (2023).