# Trade-offs between number fluctuations and response in nonequilibrium
chemical reaction networks

Hyun-Myung Chun<EMAIL_ADDRESS>School of Physics, Korea Institute for
Advanced Study, Seoul, 02455, Korea  Jordan M. Horowitz<EMAIL_ADDRESS>Department of Biophysics, University of Michigan, Ann Arbor, Michigan, 48109,
USA  Center for the Study of Complex Systems, University of Michigan, Ann
Arbor, Michigan 48104, USA  Department of Physics, University of Michigan, Ann
Arbor, Michigan 48109, USA

###### Abstract

We study the response of chemical reaction networks driven far from
equilibrium to logarithmic perturbations of reaction rates. The response of
the mean number of a chemical species is observed to be quantitively limited
by number fluctuations as well as the maximum thermodynamic driving force. We
prove these trade-offs for linear chemical reaction networks and a class of
nonlinear chemical reaction networks with a single chemical species. Numerical
results for several model systems support the conclusion that these trade-offs
continue to hold for a broad class of chemical reaction networks, though their
precise form appears to depend sensitively on the deficiency of the network.

##  I Introduction

Driven networks of chemically reacting species display a rich phenomenology
that allows them to operate as a scaffold for an array of functions, both
within engineered setups as well as in living organisms. In many of these
situations, a fruitful method for assessing performance is to measure how
sensitively the network responds to external stimuli, and in particular how
the supply of energy shapes that response. Energy transduction can amplify
sensitivity in biological signaling cascades [1, 2, 3, 4, 5] and in gene
transcription networks [6, 7], improving their efficacies. It can induce
negative differential response that enhances robustness in biochemical
networks [8]. Moreover, response can serve as an order parameter to
interrogate a dynamical nonequilibrium phase transition [9]. While the
nonequilibrium nature of the driving is integral to these scenarios, most of
our understanding to date about the general principles that shape response is
limited to near-equilibrium situations.

Near thermodynamic equilibrium, response is dictated by fluctuations as
encoded in the fluctuation-dissipation theorem (FDT) [10, 11, 12, 13]. In the
present context, it relates the static (or zero-frequency) response of the
mean number of a chemical species $\langle n\rangle$ due to a change in its
chemical potential $\mu$ with the number variance ${\rm Var}\\{n\\}$:

$\frac{\partial\langle n\rangle}{\partial\mu}=\beta{\rm Var}\\{n\\},$ (1)

where $\beta$ is the inverse temperature. Sensitive networks must also be
noisy; therefore, by driving the dynamics we can create functional networks
that are low-noise and responsive.

The fundamental nature of the FDT suggests that response and fluctuations
remain linked in some form even for systems driven far from equilibrium. To
date, a number of such relationships have been identified that link response
to some measure of fluctuations [14, 15, 16, 17, 18, 19, 20, 21]. While the
predictions offer fundamental insights, putting them to use requires system-
specific information. For some proposals, identifying the correct measure of
fluctuations requires either knowledge of the steady-state distribution or
information about the dynamics [15, 22, 23, 24] (see, e.g., [21] for a
review). Others require an appropriate choice of perturbation [25, 26, 27,
28]. An alternative approach is offered by relaxing the desire for equalities
and instead investigating potential trade-offs (or inequalities). For example,
it was recently shown that the maximum response to any perturbation is
generally limited by an information-theoretic measure of dynamical
fluctuations [29]. The inspiration for the present work is the recently
established thermodynamic-response inequalities for discrete and continuous
Markov processes [30, 31]. Together these works offer hope that simple and
general inequalities can be identified that can serve as guideposts for
rationalizing response in nonequilibrium dynamics.

In this paper, we study fluctuation-response relations for chemical reaction
networks (CRNs) in nonequilibrium steady states driven by thermodynamic
forces. Inspired by [30], we consider specific classes of perturbations that
have proven to lead to concrete predictions in the context of small stochastic
systems, and analyze potential inequalities between the easily-measured
quantities that traditionally appear in the FDT (1): number response and
fluctuations. In particular, we have studied how the CRN responds to changes
in the rate constants $k_{+\rho}$ and $k_{-\rho}$ of a chemical reaction
channel, hypothesizing that there are simple trade-offs between number
response, fluctuations, and the maximum thermodynamic driving force
${\mathcal{F}}$:

$\left|k_{\pm\rho}\frac{\partial\langle n\rangle}{\partial
k_{\pm\rho}}\right|\leq\alpha_{1}{\rm Var}\\{n\\},$ (2)
$\left|k_{+\rho}\frac{\partial\langle n\rangle}{\partial
k_{+\rho}}+k_{-\rho}\frac{\partial\langle n\rangle}{\partial
k_{-\rho}}\right|\leq\alpha_{2}{\rm Var}\\{n\\}g(\mathcal{F}),$ (3)

where $\alpha_{1,2}$ are system-dependent prefactors, $g(\mathcal{F})$ is an
increasing function with $g(0)=0$, and the shorthand notation $k_{\pm\rho}$
means ‘$k_{+\rho}$ or $k_{-\rho}$’ throughout. The values of $\alpha_{1,2}$
and the form of $g(\mathcal{F})$, in general, depend on the topological
structure of the CRN, the chemical species of interest, and the perturbed
reaction. We prove the validity of (2) and (3) with
$g(\mathcal{F})=\tanh(\mathcal{F}/4)$ for linear CRNs and for CRNs with birth-
death dynamics. For CRNs with zero deficiency, we further provide numerical
evidence that supports the accuracy of (2) and (3) with
$g(\mathcal{F})=\tanh(\mathcal{F}/4)$. Numerical studies of CRNs with nonzero
deficiency also suggest trade-offs between number response and fluctuations,
though we have not been able to identify a functional form for
$g({\mathcal{F}})$. Our analysis builds on a rich literature regarding the
topological structure of CRNs [32, 33] and how that interfaces with
nonequilibrium thermodynamics, both in the context of macroscopic
deterministic dynamics [34, 35, 36, 37, 38, 39, 40] and within microscopic
stochastic descriptions [41, 42]. These theoretical developments have already
led to the identification of universal thermodynamic properties, such as
fluctuation theorems [43, 44, 45], as well as information-geometric trade-offs
[46, 47].

The fluctuation-response relations explored here are intimately related to the
inequalities predicted in [30] for the response of nonequilibrium Markov jump
processes; however, the trade-offs (2) and (3) are not a trivial consequence.
The theory developed in [30] applies only to a system with a finite number of
microscopic states, and when applied to rate constant perturbations in CRNs
leads to bounds that scale with system size. However, the structure of the
microscopic state space of a CRN is restricted by symmetry due to repeated
subunits, not taken into account in previous work. Consequently, the
inequalities (2) and (3) can provide finite bounds on the response even in
CRNs with an infinite number of microscopic states, and in macroscopic
systems. Thus, this work complements the previous analyses, demonstrating that
inequalities predicted in [30] can be modified to a larger class of scenarios.

This paper is organized as follows. In Sec. II, we set the stage by
introducing notation. In Sec. III, we begin with a class of CRNs with a single
chemical species to illustrate how the trade-offs arise from the topological
structure of the microscopic state space. The trade-offs for linear CRNs are
proved in Sec. IV. We also derive weaker inequalities between fluctuation and
response for deficiency zero nonlinear CRNs. In Sec. IV.3, we demonstrate the
compact forms (2) and (3) may also hold for deficiency zero nonlinear CRNs by
providing numerical results for various models. In Sec. V, we numerically
examine the validity of the trade-offs (2) and (3) for deficiency nonzero
CRNs. We close the paper in Sec. VI with a summary of the main results and
concluding remarks.

##  II Chemical Reaction Networks

We begin by laying out the definitions and notations we will need to analyze
response and noise in CRNs. As a guide, we provide a concrete example of a CRN
in Fig. 1.

Figure 1: Example of a chemical reaction network: It is characterized by three
chemical species $\\{X_{1}=A,X_{2}=B,X_{3}=C\\}$, five complexes
$\\{\bm{y}_{1}=(1,1,0)^{\rm T},\bm{y}_{2}=(2,0,0)^{\rm
T},\bm{y}_{3}=(0,1,1)^{\rm T},\bm{y}_{4}=(1,0,0)^{\rm
T},\bm{y}_{5}=(0,0,1)^{\rm T}\\}$, and three reaction channels
{$\bm{\nu}_{1}^{+}(=\bm{y}_{1})\rightleftharpoons\bm{\nu}_{1}^{-}(=\bm{y}_{2})$,
$\bm{\nu}_{2}^{+}(=\bm{y}_{2})\rightleftharpoons\bm{\nu}_{2}^{-}(=\bm{y}_{3})$,
$\bm{\nu}_{3}^{+}(=\bm{y}_{4})\rightleftharpoons\bm{\nu}_{3}^{-}(=\bm{y}_{5})$}.
The total number of chemical species $\bm{l}\cdot\bm{n}=n_{1}+n_{2}+n_{3}$ is
conserved over time.

###  II.1 Constituents and Network Structure

A CRN consists of $N_{S}$ chemical species,
${\bm{X}}=\\{X_{i}\\}_{i\in\mathcal{S}}$ with
${\mathcal{S}}=\\{1,\dots,N_{S}\\}$, reacting via $N_{R}$ reaction channels,
$\rho\in\mathcal{R}=\\{1,\cdots,N_{R}\\}$,

${\bm{\nu}}_{\rho}^{+}\cdot{\bm{X}}\rightleftharpoons{\bm{\nu}}_{\rho}^{-}\cdot{\bm{X}}.$
(4)

Here, the elements of the stoichiometric vectors
${\bm{\nu}}_{\rho}^{+}=\\{\nu_{i\rho}^{+}\\}_{i\in\mathcal{S}}$ and
${\bm{\nu}}_{\rho}^{-}=\\{\nu_{i\rho}^{-}\\}_{i\in\mathcal{S}}$ are the
stoichiometric coefficients that denote the number of species $X_{i}$
participating in the forward and reverse chemical reactions of each channel.
Note that the designation of forward or reverse (and reactant or product) is
arbitrary, but fixed from the onset. Within this framework we can also include
annihilation or creation due to exchanges with the surroundings by setting
either ${\bm{\nu}}_{\rho}^{+}={\bm{\emptyset}}$ or
${\bm{\nu}}_{\rho}^{-}={\bm{\emptyset}}$, where
${\bm{\emptyset}}=(0,\dots,0)^{\rm T}$. Lastly, to ensure thermodynamic
consistency, we will assume all reaction channels are microscopically
reversible, i.e., if a chemical reaction is present then so is its reverse.

Many of the features of the reaction dynamics will turn out to depend on the
topological structure encoded in the CRN. We can exhibit this structure by
collecting the changes caused by each reaction channel
$\Delta\bm{\nu}_{\rho}={\bm{\nu}}_{\rho}^{-}-{\bm{\nu}}_{\rho}^{+}$ into the
$N_{S}\times N_{R}$ stoichiometric matrix

$\mathsf{S}=\begin{pmatrix}|&|&&|\\\
\Delta\bm{\nu}_{1}&\Delta\bm{\nu}_{2}&\cdots&\Delta\bm{\nu}_{N_{R}}\\\
|&|&&|\end{pmatrix}.$ (5)

Each element $S_{i\rho}=\Delta\nu_{i\rho}$ represents the change in the number
of species $X_{i}$ when reaction $\rho$ occurs. Denoting the kernel (or right
null space) of $\mathsf{S}$ as $\ker(\mathsf{S})$, we then have by the rank-
nullity theorem that only $s=N_{R}-\dim(\ker(\mathsf{S}))$ columns of
$\mathsf{S}$ are linearly independent. Without loss of generality, we can then
reserve the first $s$ indices for linearly independent reactions: in other
words, a vector $\Delta\bm{\nu}_{\rho}$ with $\rho>s$ is a linear combination
of vectors $\Delta\bm{\nu}_{\rho}$ with $\rho\leq s$ [36].

The vectors in the left and right null space of $\mathsf{S}$ then have
physical meanings that help to characterize the reaction dynamics [36, 38,
42]. Vectors ${\bm{l}}$ in the left null space, ${\bm{l}}\cdot\mathsf{S}=0$,
represent conservation laws: weighted combinations of chemical species that do
not change when any reaction occurs. Indeed, when reaction $\rho$ occurs the
change in the quantity ${\bm{l}}\cdot{\bm{X}}$ is zero,
${\bm{l}}\cdot\Delta{\bm{\nu}}_{\rho}=0$. We can then choose a set of linearly
independent basis vectors
$\\{\bm{l}_{\alpha}|\alpha=1,\cdots,\dim(\ker(\mathsf{S}^{\rm T}))\\}$,
representing a set of independent conservation laws. Similarly, a column
vector $\bm{r}$ in the right null space, $\mathsf{S}\cdot\bm{r}=\bm{0}$,
represents a chemical cycle: a series of reactions that leaves the numbers of
all of chemical species unchanged. Indeed, since the elements of $\mathsf{S}$
are all integers, the solution of the linear equation
$\mathsf{S}\cdot\bm{r}=0$ are rational numbers in general, which in turn can
be rescaled to be integers. In this way, a vector $\bm{r}$ can identify a
series of reactions that do not change the numbers of each species, where each
element ${\bm{r}}_{\rho}$ represents the number of occurrences of chemical
reaction $\rho$ in the cycle and whose sign denotes the direction (see Fig.
1). Again we can choose a set of linearly independent basis vectors
$\\{\bm{r}_{\alpha}|\alpha=1,2,\cdots,\dim(\ker(\mathsf{S}))\\}$.

Chemical reaction network theory has identified an alternative decomposition
of the network structure that turns out to be important for characterizing the
dynamics [32, 33]. Instead of analyzing the relationships between chemical
species using the stoichiometric matrix, as we have done, we look at the
relationships between the complexes: the collection of chemical species on
either side of a reaction channel together with their stoichiometric
coefficients. For example, in the reaction channel
$2X_{1}+X_{2}\rightleftharpoons 3X_{3}$, the complexes are
$\\{2X_{1}+X_{2},3X_{3}\\}$ while the chemical species are
$\\{X_{1},X_{2},X_{3}\\}$. Each complex can then be identified with the
corresponding stoichiometric vector ${\bm{\nu}}_{\rho}^{\pm}$. For example, in
the reaction $2X_{1}+X_{2}\rightleftharpoons 3X_{3}$, the complexes can be
represented by the column vectors $(2,1,0)^{\rm T}$ and $(0,0,3)^{\rm T}$.
Now, not all $\bm{\nu}_{\rho}^{\pm}$ are distinct, so the number of complexes
$N_{C}$ may be smaller than or equal to the number of chemical reactions
($2N_{R}$). To emphasize this, we introduce a separate notation for the
complexes $\bm{y}_{l}$ with $l\in\mathcal{C}=\\{1,2,\cdots,N_{C}\\}$,
collecting them into a $N_{S}\times N_{C}$ matrix ${\mathsf{Y}}$ whose
elements contain information about which chemical species constitute each
complex,

$\mathsf{Y}=\begin{pmatrix}|&|&&|\\\
\bm{y}_{1}&\bm{y}_{2}&\cdots&\bm{y}_{N_{C}}\\\ |&|&&|\end{pmatrix}.$ (6)

The complexes further form a complex-reaction graph, where the complexes are
the vertices and the reaction channel are the edges (see Fig. 1). The topology
of this graph is captured by the $N_{C}\times N_{R}$ incidence matrix
$\mathsf{C}$, whose elements represent which reactions link which complexes,
with signs denoting incoming or outgoing:

$C_{l\rho}=\begin{cases}1&{\rm if}~{}\bm{y}_{l}=\bm{\nu}_{\rho}^{-},\\\
-1&{\rm if}~{}\bm{y}_{l}=\bm{\nu}_{\rho}^{+},\\\ 0&{\rm
otherwise}.\end{cases}$ (7)

Now, one of the central observations of the mathematical theory of CRNs is
that the stoichiometric matrix can be decomposed as
$\mathsf{S}=\mathsf{Y}\cdot\mathsf{C}$, linking how reactions change chemical
species to how they change complexes. An immediate consequence is the fact
that the $\ker(\mathsf{S})$ may not coincide with $\ker(\mathsf{C})$, implying
that there could be chemical cycles that are not visible in the complex-
reaction graph. This discrepancy can be captured by introducing the deficiency
of the CRN $\delta=\dim(\ker(\mathsf{S}))-\dim(\ker(\mathsf{C}))$. Its value
has profound implications for what kinds of dynamical phenomenology the CRN
can exhibit [33]. When $\delta=0$, the dynamics has many features similar to
equilibrium systems, whereas $\delta\neq 0$ has been linked to bistability
[48], oscillations [49], and has been observed to affect the stochastic
thermodynamics [50]. We will also see that the CRN’s deficiency has
implications for nonequilibrium response as well. The formula for deficiency
we have introduced has a clear interpretation, but the original definition can
be more useful for determining the value of $\delta$. To this end, we note
that a complex-reaction graph may be dissected into $\ell$ isolated connected
components called linkage classes. Then, the deficiency of a CRN can be
expressed as $\delta=N_{C}-\ell-s$ [33], which clearly shows it is a
nonnegative integer. For example, the CRN in Fig. 1 has a deficiency of
$\delta=5-2-2=1$. It has no visible cycle in the complex-reaction graph but a
single chemical cycle visible in the microscopic state space.

###  II.2 Dynamics and Fluctuations

To specify the dynamics, we consider a dilute, well-stirred mixture in a fixed
volume $\Omega$. In this situation, an accurate dynamical model is provided by
assuming mass-action kinetics, i.e., the rate of a chemical reaction is
proportional to the concentration of the reactants. In other situations,
distinct types of effective kinetics can emerge. For example, Michaelis-Menten
kinetics are widely used to model enzymatic reactions [51]. However, such
emergent effective kinetics may be restored by taking appropriate limits from
mass action kinetics [52, 39], and thus we do not explicitly consider them
here.

Under our assumptions, an accurate description of the fluctuations is to model
the chemical reactions as stochastic Poisson processes that randomly change
the population of the chemical species $\bm{n}=(n_{1},\cdots,n_{N_{S}})$,
where $n_{i}$ denotes the number of species $X_{i}$. Then by the law of mass
action, the reaction rates are given by [53, 54]

$w_{\pm\rho}(\bm{n})=\frac{k_{\pm\rho}}{\Omega^{\sum_{i\in\mathcal{S}}\nu_{i\rho}^{\pm}-1}}\frac{\bm{n}!}{(\bm{n}-\bm{\nu}_{\rho}^{\pm})!},$
(8)

where the $k_{\pm\rho}$ are the rate constants, and $\bm{a}!$ is a shorthand
for $\prod_{i\in\mathcal{S}}a_{i}!$. As a result, the time-evolution of the
probability $P(\bm{n},t)$ for the population to be $\bm{n}$ at time $t$ is
given by the Chemical Master Equation [53, 54]

$\displaystyle\partial_{t}P(\bm{n},t)$
$\displaystyle=\sum_{\rho\in\mathcal{R}}w_{+\rho}(\bm{n}-\Delta\bm{\nu}_{\rho})P(\bm{n}-\Delta\bm{\nu}_{\rho},t)$
(9)
$\displaystyle~{}~{}~{}+\sum_{\rho\in\mathcal{R}}w_{-\rho}(\bm{n}+\Delta\bm{\nu}_{\rho})P(\bm{n}+\Delta\bm{\nu}_{\rho},t)$
$\displaystyle~{}~{}~{}-\sum_{\rho\in\mathcal{R}}\\{w_{+\rho}(\bm{n})+w_{-\rho}(\bm{n})\\}P(\bm{n},t)$
$\displaystyle=\hat{\mathcal{L}}P(\bm{n},t).$

Note that conservation laws restrict the set of available states where these
dynamics can take place. Indeed, recall that any basis vector
${\bm{l}}_{\alpha}$ in the left null space of ${\mathsf{S}}$ is by definition
orthogonal to every $\Delta\bm{\nu}_{\rho}$. Thus, chemical reactions, which
change $\bm{n}$ only in increments of $\pm\Delta\bm{\nu}_{\rho}$, can only
move the population through a subspace orthogonal to the left null space of
$\mathsf{S}$. This subspace of accessible states is referred to as a
(nonnegative) stoichiometric compatibility class and is identified by one
particular set of values
$\\{\Lambda_{\alpha}=\bm{l}_{\alpha}\cdot\bm{n}|\alpha=1,2,\cdots,\dim(\ker(\mathsf{S}^{\rm
T}))\\}$ [33]. Note also, that a stoichiometric compatibility class may
consist of several isolated connected components. We can visualize this
microscopic state space as a graph(s) where the vertices are possible
populations of the chemical species ${\bm{n}}$ and the edges represent
allowable chemical reactions that connect accessible configurations (see Fig.
1).

We will assume that $P(\bm{n},t)$ converges to a steady-state distribution
$\pi(\bm{n})$ in the long time limit, given by the solution of
$\hat{\mathcal{L}}\pi(\bm{n})=0$. The total probability of the population
being in an isolated connected component $\Gamma$,
$p_{\Gamma}=\sum_{\bm{n}\in\Gamma}P(\bm{n},0)$, is conserved. This constraint
implies that the steady-state distribution segregates onto each component,
$\pi(\bm{n})=\sum_{\Gamma}p_{\Gamma}\pi_{\Gamma}(\bm{n})$, with conditional
steady-state distribution $\pi_{\Gamma}(\bm{n})$ that we assume to be unique.
With knowledge of this distribution, we can characterize the steady-state
fluctuations in the species numbers ${\bm{n}}$ using the mean and the variance

$\displaystyle\langle n_{i}\rangle$
$\displaystyle=\sum_{\bm{n}}n_{i}\pi(\bm{n})$ (10) $\displaystyle{\rm
Var}\\{n_{i}\\}$ $\displaystyle=\sum_{\bm{n}}(n_{i}-\langle
n_{i}\rangle)^{2}\pi(\bm{n}).$ (11)

Now, while the Chemical Master Equation offers an accurate description of the
steady-state fluctuations, it can quickly become impractical to analyze, both
numerically and analytically, as the system size grows, $\Omega\to\infty$. In
this macroscopic limit, the fluctuations are suppressed and the limiting
behavior of the populations is determined by the evolution of the
concentrations of the chemical species,
$\bm{n}/\Omega\to[\bm{X}]=([X_{1}],\cdots,[X_{N_{S}}])^{\rm T}$, where
$[X_{i}]$ is the number of species $X_{i}$ per unit volume [55, 56]. In this
case, the dynamics are governed by the deterministic Chemical Rate Equation,

$\begin{split}\frac{d}{dt}[\bm{X}]&=\sum_{\rho\in\mathcal{R}}\Delta\bm{\nu}_{\rho}\left(k_{+\rho}[\bm{X}]^{\bm{\nu}^{+}_{\rho}}-k_{-\rho}[\bm{X}]^{\bm{\nu}^{-}_{\rho}}\right)\\\
&={\mathsf{S}}\cdot{\bm{J}}([{\bm{X}}])\end{split}$ (12)

where $\bm{a}^{\bm{b}}$ is a shorthand for the component-wise product
$\prod_{i\in\mathcal{S}}a_{i}^{b_{i}}$ with the convention $0^{0}=1$. The
column vector $\bm{J}([\bm{X}])$ has components
$J_{\rho}([\bm{X}])=k_{+\rho}[\bm{X}]^{\bm{\nu}_{\rho}^{+}}-k_{-\rho}[\bm{X}]^{\bm{\nu}_{\rho}^{-}}$
that represent the net reaction flux through each reaction channel. We note
that the same rate constants appear in the microscopic (9) and the macroscopic
dynamics (12) due to the inclusion of the volume-dependent scaling factor in
the microscopic kinetics (8).

Again, in the absence of intervention the concentrations will relax to a
steady-state value $[\bm{X}]_{\rm ss}$ given as a solution of
$\mathsf{S}\cdot\bm{J}([\bm{X}]_{\rm ss})=\bm{0}$. However, for a single
isolated connected component of the microscopic state space, in sharp contrast
to the unique steady-state distribution of the Chemical Master Equation (9),
the Chemical Rate Equation (12) may have multiple stable fixed points for
certain collections of the rate constants, resulting in multiple potential
steady-state values $[\bm{X}]_{\rm ss}$. Which steady-state is observed as
$t\to\infty$ depends on the initial conditions. This distinction comes from
the ordering of limits: $t\to\infty$ versus $\Omega\to\infty$. Choosing the
order of limits by taking $t\to\infty$ first, we can identify the unique
steady-state concentration that is placed at the dominant peak of the
probability distribution [57, 58], which is one of the stable steady-state
solutions of the deterministic equation (12) in multistable systems [58]. From
now on this most-likely steady-state concentration will be simply denoted as
$[\bm{X}]_{\rm ss}$, and we will only study the response and fluctuations
around this configuration.

The Chemical Rate Equation is a deterministic equation, and therefore does not
offer information about the fluctuations. However, having identified the
steady-state concentration $[{\bm{X}}]_{\rm ss}$, we measure the small
fluctuations about this configuration from the scaled variance

$\displaystyle D_{i}=\lim_{\Omega\to\infty}\frac{1}{\Omega}{\rm
Var}\\{n_{i}\\},$ (13)

which requires some knowledge of the microscopic dynamics. In this article, we
will use two methods to calculate these fluctuations in the macroscopic limit
for different models. When available, we will use the large deviation function
$I([X]_{i})=-\lim_{\Omega\to\infty}\Omega^{-1}\ln\pi(n_{i}=\Omega[X]_{i})$ to
analytically determine the scaled variance, via its curvature at
$[{\bm{X}}]_{\rm ss}$ [57]. Otherwise, we will determine the fluctuations
numerically using the Linear Noise Approximation to the Chemical Master
Equation, which results in a computationally-tractable Langevin-type
stochastic equation for the Gaussian fluctuations about $[{\bm{X}}]_{\rm ss}$
[59, 54]. Details can be found in Appendix A.

###  II.3 Steady-State Response

For a fixed CRN, the steady-state mean number of chemical species
$\langle{\bm{n}}\rangle$ is a function of the values of the rate constants.
Our concern here is how sensitively this average responds to changes to those
rate constants. As it will turn out focusing on (dimensionless) logarithmic
changes will facilitate the mathematical analysis. In this case, the most
generic perturbation is the linear combination

$\frac{\partial\langle{\bm{n}}\rangle}{\partial\lambda}=\sum_{\rho\in{\mathcal{R}}}\lambda_{+\rho}k_{+\rho}\frac{\partial\langle{\bm{n}}\rangle}{\partial
k_{+\rho}}+\lambda_{-\rho}k_{-\rho}\frac{\partial\langle{\bm{n}}\rangle}{\partial
k_{-\rho}},$ (14)

with expansion coefficients $\lambda_{\pm\rho}$. Similarly, in the macroscopic
limit we can address the sensitivity of the steady-state concentration
$\partial_{\lambda}[{\bm{X}}_{\rm ss}]$.

Making predictions for these arbitrary perturbations can be quite challenging.
Therefore inspired by the recent work [30, 31] on the response of
nonequilibrium Markov processes, we will largely focus on two types of
perturbations. The first is perturbing a single rate constant logarithmically:
$\partial/\partial\lambda=k_{\pm\rho}\partial/\partial k_{\pm\rho}$. The
second is a concerted change in the forward and reverse rate constants of a
single reaction channel, which we will compactly denote as

$\frac{\partial}{\partial B_{\rho}}=k_{+\rho}\frac{\partial}{\partial
k_{+\rho}}+k_{-\rho}\frac{\partial}{\partial k_{-\rho}}.$ (15)

An example of such a perturbation is a change in the concentration of a
catalyst, which increases/decreases both rates of a reaction channel
uniformly.

###  II.4 Nonequilibrium Thermodynamics

Driven or nonequilibrium chemical dynamics are characterized by the presence
of nonzero reaction currents, a signal of time-reversal symmetry breaking.
This is only possible when there is an imbalance of reaction rates around
chemical cycles, motivating the introduction of a thermodynamic force that
measures the strength of this asymmetry. Introducing a row vector $\bm{f}$
with components $f_{\rho}=\ln(k_{+\rho}/k_{-\rho})$, we define the
thermodynamic force associated with the chemical cycle ${\bm{r}}$ by [38]

$F_{\bm{r}}=\bm{f}\cdot\bm{r}=\sum_{\rho\in\mathcal{R}}r_{\rho}\ln\left(\frac{k_{+\rho}}{k_{-\rho}}\right).$
(16)

When all the thermodynamic forces are zero, every reaction is balanced by its
reverse, a statistical condition known as detailed balance true of all
equilibrium systems. Turning this around, we can say that a larger force
corresponds to a system driven farther from equilibrium. Note we identify the
thermodynamic forces around chemical cycles (defined via the stoichiometric
matrix), which should be distinguished from the loops in the microscopic
state-space graph. While these microscopic loops are used to identify
thermodynamic forces in stochastic thermodynamics via Schnakenberg’s network
theory [60, 38, 42], in the context of CRNs they are actually repeated copies
of the chemical cycles, and therefore it is sufficient to consider forces just
around the chemical cycles [36, 38, 42].

To each reaction channel, we can associate a maximum thermodynamic force
driving that current,

${\mathcal{F}}=\max_{{\bm{r}}\ni\rho}|F_{\bm{r}}|$ (17)

We will show that this quantity measures the strength of nonequilibrium
driving through the reaction channel, and also constrains the response to
changing the reaction channel’s kinetic rates.

The physical meaning of thermodynamic forces becomes explicit once we take
into account that our CRNs can be open to an environment containing other
chemical species chemostatted to fix concentrations that facilitate the
observed chemical reactions. In this respect, the rate constants $k_{\pm\rho}$
introduced above should be considered pseudo-rate constants, having absorbed
the concentrations of the chemostatted species into their definition. While we
assume that such a thermodynamically-consistent model can be constructed for
the CRNs we study, we will not specify it in order to make the theory as
general as possible.

##  III Birth-death process

We begin our investigation with a class of CRNs with a single chemical species
whose evolution is described by birth-death (or one-step) processes, defined
by the property that every chemical reaction changes the number of chemical
species by one [54]. This class includes CRNs that are both nonlinear and
nonzero deficiency. Despite this complication, there is a known closed-form
expression for the steady-state probability distribution for general reaction
rates, which will allow us to make concrete predictions for the limits of
response in these models.

Consider the CRN composed of a single chemical species $X$ that reacts via the
set of chemical reactions

$z_{\rho}X\xrightleftharpoons[k_{-\rho}]{k_{+\rho}}(z_{\rho}+1)X$ (18)

with $\rho\in\mathcal{R}=\\{1,2,\cdots,N_{R}\\}$ and $z_{\rho}$ are non-
negative integers. In turn, the reaction rates are given by

$\displaystyle w_{+\rho}(n)$
$\displaystyle=\frac{k_{+\rho}}{\Omega^{z_{\rho}-1}}\frac{n!}{(n-z_{\rho})!},$
(19) $\displaystyle w_{-\rho}(n)$
$\displaystyle=\frac{k_{-\rho}}{\Omega^{z_{\rho}}}\frac{n!}{(n-z_{\rho}-1)!}.$

The widely studied Schlögl model fits into this framework with the choices
$N_{R}=2$, $z_{1}=0$, and $z_{2}=2$ [61]. We will use this model to illustrate
the results of this section. A depiction of its microscopic state space and
thermodynamic force can be found in Fig. 2.

Owing to the simple network structure of the microscopic state space,
exemplified in Fig. 2, we can find an exact solution to the steady-state
distribution $\pi(n)$. The key observation is that for a distribution to be
stationary the reaction flux between any pair of neighboring states must be
zero,
$\sum_{\rho\in{\mathcal{R}}}w_{+\rho}(n)\pi(n)=\sum_{\rho\in{\mathcal{R}}}w_{-\rho}(n+1)\pi(n+1)$,
which constrains the ratios of steady-state values all along the chain of
states:

$\frac{\pi(n+1)}{\pi(n)}=\frac{\sum_{\rho\in{\mathcal{R}}}w_{+\rho}(n)}{\sum_{\rho\in{\mathcal{R}}}w_{-\rho}(n+1)}.$
(20)

Iterating this observation using standard methodology (see Ref. [61]) leads to
the steady-state distribution

$\pi(n)=\pi(0)\prod_{i=0}^{n-1}\frac{\sum_{\rho\in\mathcal{R}}w_{+\rho}(i)}{\sum_{\rho\in\mathcal{R}}w_{-\rho}(i+1)},$
(21)

where the probability of the boundary state $n=0$ is fixed by normalization
$\sum_{n=0}^{\infty}\pi(n)=1$.

Figure 2: Microscopic state space of the Schlögl model and its thermodynamics.

Although we have the steady-state distribution, finding a closed-form
expression for the mean number $\langle n\rangle$ cannot be done in general
without specifying the reaction rates. So we analyze the steady-state response
for an arbitrary perturbation directly, before specializing to the uni-
directional and symmetric perturbations. Using the identity
$\sum_{m}\pi(m)\partial_{\lambda}\ln\pi(m)=0$, which follows from probability
conservation, we can rewrite the response $\partial_{\lambda}\langle
n\rangle=\sum_{n}n\partial_{\lambda}\pi(n)$ as

$\frac{\partial\langle
n\rangle}{\partial\lambda}=\sum_{m,m^{\prime}}m\pi(m)\pi(m^{\prime})\frac{\partial}{\partial\lambda}\ln\frac{\pi(m)}{\pi(m^{\prime})}$
(22)

shifting our attention to the sensitivity of ratios of probabilities, which
will allow us to exploit the simple network structure of these models. Indeed,
taking advantage of (20), we have

$\displaystyle\frac{\partial}{\partial\lambda}\ln\frac{\pi(m)}{\pi(m^{\prime})}$
$\displaystyle=\sum_{i=m^{\prime}}^{m-1}\frac{\partial}{\partial\lambda}\ln\frac{\pi(i+1)}{\pi(i)}$
(23)
$\displaystyle=\sum_{i=m^{\prime}}^{m-1}\frac{\partial}{\partial\lambda}\ln\left(\frac{\sum_{\rho\in\mathcal{R}}w_{+\rho}(i)}{\sum_{\rho\in\mathcal{R}}w_{-\rho}(i+1)}\right)$

for $m>m^{\prime}$. The response for two arbitrarily distant states is given
by the sum of the responses for neighboring states.

Let us first consider single-rate perturbations. For two neighboring states,
we have that

$\displaystyle\left|k_{+\rho}\frac{\partial}{\partial
k_{+\rho}}\ln\frac{\pi(n+1)}{\pi(n)}\right|$
$\displaystyle=\frac{w_{+\rho}(n)}{\sum_{\rho^{\prime}\in\mathcal{R}}w_{+\rho^{\prime}}(n)}\leq
1$ (24) $\displaystyle\left|k_{-\rho}\frac{\partial}{\partial
k_{-\rho}}\ln\frac{\pi(n+1)}{\pi(n)}\right|$
$\displaystyle=\frac{w_{-\rho}(n+1)}{\sum_{\rho^{\prime}\in\mathcal{R}}w_{-\rho^{\prime}}(n+1)}\leq
1,$

are bounded by 1. Plugging (24) into (23) and summing over $i$, leads to the
the inequality for non-neighboring states

$\left|k_{\pm\rho}\frac{\partial}{\partial
k_{\pm\rho}}\ln\frac{\pi(m)}{\pi(m^{\prime})}\right|\leq|m-m^{\prime}|.$ (25)

The maximum response for the log-ratio of probabilities at two distant states
is limited by the distance between them, i.e., the number of reactions
required to move between the two configurations. We then obtain a bound on
(22) by first re-expressing the sum as

$\begin{split}&\left|k_{\pm\rho}\frac{\partial\langle n\rangle}{\partial
k_{\pm\rho}}\right|\\\ &\
=\left|\sum_{m>m^{\prime}}(m-m^{\prime})\pi(m)\pi(m^{\prime})k_{\pm\rho}\frac{\partial}{\partial
k_{\pm\rho}}\ln\frac{\pi(m)}{\pi(m^{\prime})}\right|\\\ &\
\leq\sum_{m>m^{\prime}}(m-m^{\prime})^{2}\pi(m)\pi(m^{\prime}),\end{split}$
(26)

having bounded the absolute value using (25) and that the sum is only over
$m>m^{\prime}$. We recognize the result as

$\left|k_{\pm\rho}\frac{\partial\langle n\rangle}{\partial
k_{\pm\rho}}\right|\leq{\rm Var}\\{n\\}.$ (27)

In words, large number fluctuations are necessary for high sensitivity.

Next, we consider the symmetric perturbation, which we anticipate will be
constrained by the thermodynamic force in addition to the number fluctuations.
Again, motivated by (23), we start by analyzing the response of neighboring
states. Differentiating (20), we find

$\displaystyle\frac{\partial}{\partial B_{\rho}}\ln\frac{\pi(n+1)}{\pi(n)}\ $
(28)
$\displaystyle=\frac{w_{+\rho}(n)\bar{w}_{-\rho}(n+1)-w_{-\rho}(n+1)\bar{w}_{+\rho}(n)}{\\{w_{+\rho}(n)+\bar{w}_{+\rho}(n)\\}\\{w_{-\rho}(n+1)+\bar{w}_{-\rho}(n+1)\\}}.$

where for convenience we have introduced the notation
$\bar{w}_{\pm\rho}(n)=\sum_{\rho^{\prime}\in\mathcal{R}\setminus\\{\rho\\}}w_{\pm\rho^{\prime}}(n)$
for the total reaction rate not including the perturbed reaction channel.
Inspired by [30], we factor the numerator as
$(\sqrt{w_{+\rho}\bar{w}_{-\rho}}+\sqrt{w_{-\rho}\bar{w}_{+\rho}})(\sqrt{w_{+\rho}\bar{w}_{-\rho}}-\sqrt{w_{-\rho}\bar{w}_{+\rho}})$
and expand the denominator
$(\sqrt{w_{+\rho}\bar{w}_{-\rho}}+\sqrt{w_{-\rho}\bar{w}_{+\rho}})^{2}+(\sqrt{w_{+\rho}w_{-\rho}}-\sqrt{\bar{w}_{+\rho}\bar{w}_{-\rho}})^{2}$
to place the absolute-value of the response in a form that can readily be
bounded:

$\left|\frac{\partial}{\partial
B_{\rho}}\ln\frac{\pi(n+1)}{\pi(n)}\right|=\frac{1}{1+a_{n}^{2}}\left|\tanh\left(\frac{\chi}{4}\right)\right|$
(29)

with

$a_{n}=\frac{\sqrt{w_{+\rho}(n)w_{-\rho}(n+1)}-\sqrt{\bar{w}_{+\rho}(n)\bar{w}_{-\rho}(n+1)}}{\sqrt{w_{+\rho}(n)\bar{w}_{-\rho}(n+1)}+\sqrt{\bar{w}_{+\rho}(n)w_{-\rho}(n+1)}}.$
(30)

and

$\chi=\ln\left(\frac{\sum_{\rho^{\prime}\in\mathcal{R}\setminus\\{\rho\\}}w_{+\rho}(n)w_{-\rho^{\prime}}(n+1)}{\sum_{\rho^{\prime}\in\mathcal{R}\setminus\\{\rho\\}}w_{+\rho^{\prime}}(n)w_{-\rho}(n+1)}\right).$
(31)

Using $1/(1+a_{n}^{2})\leq 1$ and the monotonicity of the hyperbolic tangent,
we can now bound the response as

$\displaystyle\left|\frac{\partial}{\partial
B_{\rho}}\ln\frac{\pi(n+1)}{\pi(n)}\right|\leq\tanh\left(\frac{|\chi|}{4}\right).$
(32)

The last step is to constrain $\chi$ by the thermodynamic forces driving the
CRN out of equilibrium (17). Importantly, the only forces that will matter are
those around chemical cycles that contain the perturbed reaction channel (that
is $r_{\rho}\neq 0$), which we denote as

$F_{\alpha}=\ln\left(\frac{w_{+\rho}(n)w_{-\alpha}(n+1)}{w_{+\alpha}(n)w_{-\rho}(n+1)}\right)=\ln\left(\frac{k_{+\rho}k_{-\alpha}}{k_{+\alpha}k_{-\rho}}\right).$
(33)

With this definition we can reinterpret (31) as

$|\chi|=\left|\ln\left(\frac{\sum_{\rho^{\prime}\in\mathcal{R}\setminus\\{\rho\\}}w_{+\rho^{\prime}}(n)w_{\rho}(n+1)e^{F_{\rho^{\prime}}}}{\sum_{\rho^{\prime}\in\mathcal{R}\setminus\\{\rho\\}}w_{+\rho^{\prime}}(n)w_{-\rho}(n+1)}\right)\right|\leq{\mathcal{F}},$
(34)

which is bounded by the maximum thermodynamic force
$\mathcal{F}=\max_{\rho^{\prime}}\\{|F_{\rho^{\prime}}|\\}$ through $\rho$.
From which we conclude that

$\displaystyle\left|\frac{\partial}{\partial
B_{\rho}}\ln\frac{\pi(n+1)}{\pi(n)}\right|\leq\tanh\left(\frac{\mathcal{F}}{4}\right).$
(35)

Following the same procedure carried out for single-rate perturbations, we can
now obtain a trade-off between fluctuations, response, and nonequilibrium
driving,

$\left|\frac{\partial\langle n\rangle}{\partial B_{\rho}}\right|\leq{\rm
Var}\\{n\\}\tanh\left(\frac{\mathcal{F}}{4}\right).$ (36)

Note that this bound is tighter than a simple application of the single-rate
bound, $|\partial\langle n\rangle/\partial
B_{\rho}|\leq|k_{+\rho}\partial\langle n\rangle/\partial
k_{+\rho}|+|k_{-\rho}\partial\langle n\rangle/\partial k_{-\rho}|\leq 2{\rm
Var}\\{n\\}$.

Our analytic predictions also have implications for the limits to response in
the macroscopic limit. Dividing (27) and (36) by $\Omega$ and taking the
$\Omega\to\infty$ limit, we can transfer our trade-offs to

$\left|k_{\pm\rho}\frac{\partial[X]_{\rm ss}}{\partial k_{\pm\rho}}\right|\leq
D_{X},$ (37) $\left|\frac{\partial[X]_{\rm ss}}{\partial B_{\rho}}\right|\leq
D_{X}\tanh\left(\frac{\mathcal{F}}{4}\right),$ (38)

where $D_{X}=\lim_{\Omega\to\infty}{\rm Var}\\{n\\}/\Omega$ is the scaled
variance. Remarkably, the response of the _macroscopic_ steady-state
concentration is limited by the _microscopic_ fluctuations, which are
completely absent from the deterministic description. While the information
about the microscopic fluctuations becomes hidden in the macroscopic limit,
its influence remains. A similar observation was made recently in [62], where
the entropy production measured from macroscopic dynamics bounds the
probability of rare fluctuations.

We validate the trade-off (38), by numerically calculating the response and
scaled variance of the Schlögl model. The steady-state concentration of the
Schlögl model satisfies the cubic equation

$k_{+1}-k_{-1}x+k_{+2}x^{2}-k_{-2}x^{3}=0.$ (39)

Depending on the rate constants, the equation has one or three positive
solutions. We will focus on a regime of parameters where there is only a
single fixed point to simplify the analysis. Parameterizing the rate constants
as $k_{+1}=9e^{B_{1}-F}$, $k_{-1}=3e^{B_{1}}$, $k_{+2}=3$, and $k_{-2}=1$, we
find that the equation has a unique solution as long as $F\leq\ln 9$. The
response and number fluctuations are compared in Fig. 3 for $F<\ln 9$ and
$F=\ln 9$. The bound disappears at the critical point where the number
fluctuation diverges [63]. The ratio of the response to the bound demonstrates
that inequality (38) is saturable (insets of Fig. 3).

Figure 3: The response of the steady-state concentration to a symmetric
perturbation in the Schlögl model. The rate constants are
$k_{+1}=9e^{B_{1}-F}$, $k_{-1}=3e^{B_{1}}$, $k_{+2}=3$, and $k_{-2}=1$ with
(a) $F=1$ and (b) $F=\ln 9$. The dashed lines represent $D_{X}$ in (a) and
$D_{X}\tanh(F/4)$ in (b). The insets depict the ratio of the response to the
bound. The scaled variance $D_{X}$ is obtained from the second derivative of
the large deviation function
$I(x)=-\lim_{\Omega\to\infty}\Omega^{-1}\ln\pi(n=\Omega x)$ [57].

##  IV Deficiency Zero

Trade-offs between response and fluctuation can also be shown to exist for
deficiency zero CRNs. We will distinguish two types of CRNs, linear and
nonlinear, based on whether the dynamics in the deterministic Chemical Rate
Equation is linear or nonlinear. For linear CRNs we will be able to derive
simple bounds on the response in terms of number fluctuations, in line with
(2) and (3). For nonlinear CRNs, we will provide numerical studies that
suggest that number fluctuations constrain nonequilibrium response in these
cases as well; however, we cannot prove them. Some analytic progress is
possible though, and other trade-offs between response and fluctuations will
be demonstrated.

Deficiency zero CRNs have a special property that makes their response to
perturbations amenable to theoretical analysis. According to the deficiency
zero theorem [33], for every choice of rate constants, the fixed points of the
Chemical Rate Equation (12) of a (weakly) reversible, deficiency zero CRNs are
such that the net reaction flux between every pair of complexes perfectly
balances, ${\mathsf{C}}\cdot{\bm{J}}([{\bm{X}}]_{\rm ss})={\bm{0}}$, a
property known as complex-balanced (cf. the general condition
${\mathsf{S}}\cdot{\bm{J}}([{\bm{X}}]_{\rm
ss})={\mathsf{Y}}\cdot{\mathsf{C}}\cdot{\bm{J}}([{\bm{X}}]_{\rm
ss})={\bm{0}}$). Furthermore, the complex-balanced steady-state concentration
is unique up to a given stoichiometric compatibility class and stable. There
are a number of implications of this property that will turn out to be
important for the study of response.

The first is that the steady-state distribution of any complex-balanced CRN
has a product form [64, 65]. To be precise, for a state $\bm{n}$ on an
isolated connected component $\Gamma$ of a stoichiometric compatibility class
identified by a set of conservation laws
$\\{\Lambda_{\alpha}=\bm{l}_{\alpha}^{\rm
T}\cdot\bm{n}|\alpha=1,2,\cdots,\dim(\ker(\mathsf{S}^{\rm T}))\\}$, the
conditional steady-state probability is given by

$\pi_{\Gamma}(\bm{n})=\frac{1}{Z_{\Gamma}([\bm{X}]_{\rm
ss})}\frac{(\Omega[\bm{X}]_{\rm ss})^{\bm{n}}}{\bm{n}!},$ (40)

where $Z_{\Gamma}([\bm{X}]_{\rm
ss})=\sum_{\bm{n}\in\Gamma}(\Omega[\bm{X}]_{\rm ss})^{\bm{n}}/\bm{n}!$ is a
normalization factor. Now generically if we were to change the rates of a
complex-balanced CRN there is no guarantee that the resulting dynamics would
be complex-balanced nor would the steady-state distribution maintain its
product form. However, for deficiency zero CRNs, all steady-states are
guaranteed to be complex-balanced, and so a perturbation only has the effect
of shifting the complex-balanced steady-state concentrations while keeping the
product form of (40) intact. As a consequence, the static response of the
conditional average number of any chemical species $\langle
n_{i}\rangle_{\Gamma}=\sum_{\bm{n}\in\Gamma}n_{i}\pi_{\Gamma}(\bm{n})$ can be
determined by differentiating (40),

$\frac{\partial\langle
n_{i}\rangle_{\Gamma}}{\partial\lambda}=\sum_{j\in\mathcal{S}}{\rm
Cov}_{\Gamma}\\{n_{i},n_{j}\\}\frac{\partial\ln[X_{j}]_{\rm
ss}}{\partial\lambda},$ (41)

where ${\rm Cov}_{\Gamma}\\{n_{i},n_{j}\\}=\langle
n_{i}n_{j}\rangle_{\Gamma}-\langle n_{i}\rangle_{\Gamma}\langle
n_{j}\rangle_{\Gamma}$. This equation (41) links the microscopic stochastic
response ($\partial_{\lambda}\langle n_{i}\rangle_{\Gamma}$) with the
macroscopic deterministic response ($\partial_{\lambda}\ln[X_{j}]_{\rm ss}$)
through the fluctuations (${\rm Cov}_{\Gamma}\\{n_{i},n_{j}\\}$).

To make further progress, we need to constrain the macroscopic response, and
again restricting attention to deficiency zero CRNs offers a method. It turns
out that the equation for the complex-balanced steady-state concentration
${\mathsf{C}}\cdot{\bm{J}}([{\bm{X}}]_{\rm ss})={\bm{0}}$ can be cast into a
form that allows for a closed-form graphical solution. To this end, we
introduce the Laplacian matrix $\mathsf{A}$ of the complex-reaction graph with
off-diagonal elements

$A_{lm}=\begin{cases}k_{+\rho}&{\rm
if}~{}\bm{y}_{l}=\bm{\nu}_{\rho}^{-}~{}{\rm
and}~{}\bm{y}_{m}=\bm{\nu}_{\rho}^{+},\\\ k_{-\rho}&{\rm
if}~{}\bm{y}_{l}=\bm{\nu}_{\rho}^{+}~{}{\rm
and}~{}\bm{y}_{m}=\bm{\nu}_{\rho}^{-},\\\ 0&{\rm otherwise},\end{cases}$ (42)

for $l\neq m$ and diagonal elements
$A_{ll}=-\sum_{m\in\mathcal{C}\setminus\\{l\\}}A_{lm}$ [32, 33]. We also
introduce the vector $\bm{\Psi}(\bm{x})=\\{\bm{x}^{\bm{y}_{l}}\\}_{l\in
N_{C}}$ whose components are products of chemical species raised to a power
determined by how often that species appears in a given complex. Despite the
additional complication of introducing these two quantities, it can be shown
that the complex-balanced steady-state concentration now solves the
alternative equation ${\mathsf{C}}\cdot{\bm{J}}([{\bm{X}}]_{\rm
ss})=\mathsf{A}\cdot\bm{\Psi}([\bm{X}]_{\rm ss})={\bm{0}}$. The virtue of this
reformulation is that the right null vectors of the Laplacian matrix can be
given by a graphical solution based on the matrix-tree theorem [66, 67, 30],
which will also allow us to import some of the methodology from [30] to
constrain the response.

###  IV.1 Linear CRNs ($\delta=0$)

Figure 4: Linear CRN with two linkage classes: Each linkage class of a linear
CRN can be regarded as an independent system.

Every linear CRN has deficiency zero. Furthermore, in any linear CRN, each
linkage class can be regarded as an independent system, with chemical species
in separate linkage classes evolving independently from each other (for
example, $A$ and $E$ in Fig. 4). Thus, without loss of generality, for linear
CRNs we need only consider a single linkage class.

For linear CRNs, the results of [30] allow us to directly bound the response
by recognizing that the deterministic $\langle{\bm{n}}\rangle$ satisfies the
same equation $\mathsf{A}\cdot\langle{\bm{n}}\rangle={\bm{0}}$ as the steady-
state distribution of a Markov jump process. For example, in a closed CRN with
a single conservation law, [30] predicts that the single-rate perturbation is
bounded by

$\left|k_{\pm\rho}\frac{\partial\langle n_{i}\rangle}{\partial
k_{\pm\rho}}\right|\leq\langle n_{i}\rangle\left(1-\frac{\langle
n_{i}\rangle}{n_{\rm tot}}\right)$ (43)

Connecting this upper bound to number fluctuations is possible, but requires
information that the steady-state distribution is of product form (40). Below
we take a different approach, building our analysis off the fluctuation-
response equality (41), which allows us to develop in parallel response bounds
for linear and nonlinear deficiency zero networks.

To begin, we need to distinguish whether or not our linear CRN is open to
exchange of chemical species with the surroundings, which would manifest as
one of the complexes being the null vector. If the null vector is a complex,
e.g., as in the right linkage class in Fig. 4, the CRN has no conservation
law. This can be deduced by noting the number of complexes is larger than the
number of species $N_{C}=N_{S}+1$ due to the presence of the null complex.
Furthermore, the number of independent reaction channels is simply one less
than the number of complexes, ${\rm rank}({\mathsf{S}}^{T})=N_{C}-1$, as each
complex is linked to at least one reaction. Then, by the rank-nullity theorem,
there are no conservation laws: ${\rm dim}({\rm
ker}({\mathsf{S}}^{T}))=N_{S}-{\rm rank}({\mathsf{S}}^{T})=N_{S}-(N_{C}-1)=0$.
In this case, the steady-state distribution becomes a multivariate Poisson
distribution [68, 64] over the entire microscopic state space, in which the
numbers of different chemical species are uncorrelated, ${\rm
Cov}\\{n_{i},n_{j}\\}={\rm Var}\\{n_{i}\\}\delta_{ij}$. Then, the fluctuation-
response relation (41) can be written as

$\frac{\partial\langle n_{i}\rangle}{\partial\lambda}={\rm
Var}\\{n_{i}\\}\frac{\partial}{\partial\lambda}\ln[X_{i}]_{\rm ss}.$ (44)

On the other hand, if a linear CRN does not include a null vector as a
complex, then $N_{C}=N_{S}$ and there is a single conservation law, ${\rm
dim}({\rm ker}({\mathsf{S}}^{T}))=N_{S}-{\rm
rank}({\mathsf{S}}^{T})=N_{S}-(N_{C}-1)=1$. The one conservation law is the
total number of chemical species $n_{\rm tot}=\sum_{i\in\mathcal{S}}n_{i}$. In
this case, the steady-state distribution becomes a multinomial distribution
over a stoichiometric compatibility class with a single isolated connected
component, whose mean and covariance are given by $\langle n_{i}\rangle=n_{\rm
tot}[X_{i}]/\sum_{j\in\mathcal{S}}[X_{j}]$ and ${\rm
Cov}\\{n_{i},n_{j}\\}=\langle n_{i}\rangle(\delta_{ij}-\langle
n_{j}\rangle/n_{\rm tot})$ [68, 64]. This allows us to separate the sum in
(41) into two part with $j=i$ and $j\neq i$:

$\displaystyle\sum_{j\in\mathcal{S}\setminus\\{i\\}}{\rm
Cov}\\{n_{i},n_{j}\\}\frac{\partial\ln[X_{j}]_{\rm ss}}{\partial\lambda}$ (45)
$\displaystyle=-{\rm
Var}\\{n_{i}\\}\frac{\partial}{\partial\lambda}\ln\left(\sum_{j\in\mathcal{S}\setminus\\{i\\}}[X_{j}]_{\rm
ss}\right).$

As a result, the fluctuation-response relation (41) becomes

$\frac{\partial\langle n_{i}\rangle}{\partial\lambda}={\rm
Var}\\{n_{i}\\}\frac{\partial}{\partial\lambda}\ln\left(\frac{[X_{i}]_{\rm
ss}}{\sum_{j\in\mathcal{S}\setminus\\{i\\}}[X_{j}]_{\rm ss}}\right).$ (46)

To proceed, we recognize that for linear CRNs there is a one-to-one
correspondence between chemical species and complexes. Thus, we can always
choose a basis for the complexes such that $\Psi_{i}([{\bm{X}}_{\rm
ss}])=[X_{i}]_{\rm ss}$, and the last complex is the null vector
$\Psi_{\bm{\emptyset}}([\bm{X}_{\rm ss}])=([{\bm{X}}]_{\rm ss})^{\bm{0}}=1$
when it is present. In this way, the condition for the complex-balanced
steady-state concentrations becomes $\mathsf{A}\cdot[\bm{X}]_{\rm ss}=\bm{0}$,
where with a slight abuse of notation we write $[X_{N_{S}+1}]=1$ for the null
vector. In this way the steady-state concentrations $[\bm{X}]_{\rm ss}$ have a
graphical representation based on the matrix-tree theorem since $[\bm{X}]_{\rm
ss}$ lies on the right null space of a Laplacian matrix $\mathsf{A}$ [66, 67,
30]. Such a graphical representation has been exploited recently to study the
response of nonequilibrium Markov jump processes [30], whose master equation
has the same structure. The theory developed in [30] is directly applicable to
the present setup. The equation (F1) of [30] can be translated into

$\left|k_{\pm\rho}\frac{\partial}{\partial k_{\pm\rho}}\ln\frac{[X_{i}]_{\rm
ss}}{[X_{j}]_{\rm ss}}\right|\leq 1$ (47)

and

$\left|k_{\pm\rho}\frac{\partial}{\partial
k_{\pm\rho}}\ln\left(\frac{[X_{i}]_{\rm
ss}}{\sum_{j\in\mathcal{S}\setminus\\{i\\}}[X_{i}]_{\rm ss}}\right)\right|\leq
1.$ (48)

Although the counterpart of the former inequality does not appear in [30]
explicitly, it is straightforwardly obtainable by applying the same logic to
prove (F1) of [30]. Similarly, the equations (C12) and (C13) of [30] can be
translated respectively into

$\left|\frac{\partial}{\partial B_{\rho}}\ln\frac{[X_{i}]_{\rm
ss}}{[X_{j}]_{\rm ss}}\right|\leq\tanh\left(\frac{\mathcal{F}}{4}\right)$ (49)

and

$\left|\frac{\partial}{\partial B_{\rho}}\ln\left(\frac{[X_{i}]_{\rm
ss}}{\sum_{j\in\mathcal{S}\setminus\\{i\\}}[X_{j}]_{\rm
ss}}\right)\right|\leq\tanh\left(\frac{\mathcal{F}}{4}\right),$ (50)

Note that only chemical cycles visible in the complex-reaction graph are taken
into account for defining $\mathcal{F}$ since the matrix-tree theorem only
uses information encoded in $\mathsf{A}$, i.e., the connectivity between
complexes. The zero deficiency then ensures all chemical cycles are visible in
the complex-reaction graph.

To complete our analysis, we apply inequalities (47)-(50) to (44) and (46), to
derive trade-offs for linear CRNs

$\left|k_{\pm\rho}\frac{\partial\langle n_{i}\rangle}{\partial
k_{\pm\rho}}\right|\leq{\rm Var}\\{n_{i}\\}$ (51)

and

$\left|\frac{\partial\langle n_{i}\rangle}{\partial B_{\rho}}\right|\leq{\rm
Var}\\{n_{i}\\}\tanh\left(\frac{\mathcal{F}}{4}\right),$ (52)

which hold regardless of the conservation law.

This analysis complements a derivation based directly on the results of [30].

###  IV.2 Nonlinear deficiency zero ($\delta=0$): Theory

The explicit form of the mean and covariance are generally unavailable for
nonlinear CRNs. In addition, different linkage classes are generally not
independent. These facts significantly complicate the analysis, requiring us
to investigate the response of the full steady-state distribution.

While we are not able to prove response is bounded by number fluctuations, we
are able to show it is limited by another measure of fluctuations in chemical
space that we now introduce. Details can be found in Appendix B. To this end,
we define a matrix formed by keeping only the linearly independent columns of
the stoichiometric matrix

$\mathsf{S}^{\prime}=\begin{pmatrix}|&|&&|\\\
\Delta\bm{\nu}_{1}&\Delta\bm{\nu}_{2}&\cdots&\Delta\bm{\nu}_{s}\\\
|&|&&|\end{pmatrix}$ (53)

Since all the columns are linearly independent, $\mathsf{S}^{\prime}$ has a
left inverse matrix $(\mathsf{S}^{\prime})^{-1}$ satisfying
$(\mathsf{S}^{\prime})^{-1}\mathsf{S}^{\prime}=\mathsf{I}$. It is worth
mentioning that ${\mathsf{S}}^{\prime}$ is not unique, in that we could have
chosen a different subset of linearly independent reactions to form
${\mathsf{S}}^{\prime}$.

Applying the matrix tree theorem to the response of the log-ratio of the
steady-state probabilities of two states, we obtain the bounds

$\left|k_{\pm\rho}\frac{\partial\langle n_{i}\rangle_{\Gamma}}{\partial
k_{\pm\rho}}\right|\leq\sum_{\sigma\in\mathcal{R}^{\prime}}\left|\sum_{j\in\mathcal{S}}(\mathsf{S}^{\prime})^{-1}_{\sigma
j}{\rm Cov}_{\Gamma}\\{n_{i},n_{j}\\}\right|$ (54)

and

$\left|\frac{\partial\langle n_{i}\rangle_{\Gamma}}{\partial
B_{\rho}}\right|\leq\sum_{\sigma\in\mathcal{R}^{\prime}}\left|\sum_{j\in\mathcal{S}}(\mathsf{S}^{\prime})^{-1}_{\sigma
j}{\rm
Cov}_{\Gamma}\\{n_{i},n_{j}\\}\right|\tanh\left(\frac{\mathcal{F}}{4}\right).$
(55)

where $\mathcal{R}^{\prime}=\\{1,2,\cdots,s\\}$. These weighted covariances
depend on tangible observable number fluctuations though are coupled through
the network topology. Further insight into the nature of these fluctuations is
offered by considering another representation of the population vectors
${\bm{n}}$, detailed in Appendix B. Any state vector $\bm{n}^{\prime}$ can be
reached from any other state vector $\bm{n}$ through a unique series of the
$s$ independent reactions in the columns of $\mathsf{S}^{\prime}$ (53). If we
denote $\mu_{\sigma}({\bm{n}},{\bm{n}}^{\prime})$ as the number of times
reaction $\sigma$ occurs along the reaction path from ${\bm{n}}$ to
${\bm{n}}^{\prime}$, then the weighted covariances can be re-expressed as

$\sum_{j\in\mathcal{S}}(\mathsf{S}^{\prime})^{-1}_{\sigma j}{\rm
Cov}_{\Gamma}\\{n_{i},n_{j}\\}=\sum_{\bm{n},\bm{n}^{\prime}\in\Gamma}n_{i}\mu_{\sigma}(\bm{n},\bm{n}^{\prime})\pi_{\Gamma}(\bm{n})\pi_{\Gamma}(\bm{n}^{\prime}),$
(56)

in terms of a correlation depending on the distance between the states
${\bm{n}}$ and ${\bm{n}}^{\prime}$ measured in reactions. In this way,
response can be viewed as bounded by number fluctuations or fluctuations in
reaction distance.

Although we were not able to derive the compact forms of trade-offs (2) and
(3), numerical evidence presented in the next subsection suggests they may be
valid for a broad class of nonlinear CRNs as well.

###  IV.3 Nonlinear deficiency zero ($\delta=0$): Numerics

We provide numerical results for several nonlinear, deficiency zero CRNs to
demonstrate the potential validity of the trade-offs in (2) and (3), allowing
for multiple conservation laws and multiple chemical cycles. Without loss of
generality, we consider here only one isolated connected component. This
follows because the total response $\partial\langle
n_{i}\rangle/\partial\lambda$ is a linear combination of the partial responses
$\partial\langle n_{i}\rangle_{\Gamma}/\partial\lambda$. If the trade-offs (2)
and (3) hold for each partial response, then the total response is also
bounded by number fluctuations,

$\left|\frac{\partial\langle
n_{i}\rangle}{\partial\lambda}\right|\leq\sum_{\Gamma}p_{\Gamma}{\rm
Var}_{\Gamma}\\{n_{i}\\}\leq{\rm Var}\\{n_{i}\\}$ (57)

due to the convexity of the variance. We find in each case studied here that
trade-offs (2) and (3) hold with common properties, $\alpha_{1}=\alpha_{2}$
and $g(\mathcal{F})=\tanh(\mathcal{F}/4)$. The values of
$\alpha_{1}=\alpha_{2}$ depend on the chemical species of interest and the
perturbed chemical reaction but not on the size of the stoichiometric
compatibility class. Details of the numerical methods can be found in Appendix
A.

In our numerical analysis, we will distinguish two methodologies for
numerically determining the response and fluctuations based on system size.
When the size of the stoichiometric compatibility class is small enough, say
due to the existence of conservation laws, such that the Chemical Master
Equation (9) is practically solvable, we will use it to directly calculate the
steady-state distribution as well as the steady-state response. Then to
examine the trade-offs (2) and (3), we will plot two quality factors

$\displaystyle{\mathcal{K}}(\pm\rho,X_{i})$ $\displaystyle=\frac{1}{{\rm
Var}\\{n_{i}\\}}\left|k_{\pm\rho}\frac{\partial\langle n_{i}\rangle}{\partial
k_{\pm\rho}}\right|,$ (58) $\displaystyle{\mathcal{B}}(\rho,X_{i})$
$\displaystyle=\frac{1}{{\rm Var}\\{n_{i}\\}}\left|\frac{\partial\langle
n_{i}\rangle}{\partial B_{\rho}}\right|.$ (59)

Other times it will be more feasible to determine the response and
fluctuations in the macroscopic limit, using the Chemical Rate Equation (12)
to determine the steady-state concentration, and the Linear Noise
Approximation to estimate the fluctuations. To distinguish this case, we
introduce modified notation for the quality factors

$\displaystyle{\mathcal{K}}^{\rm M}(\pm\rho,X_{i})$
$\displaystyle=\frac{1}{D_{i}}\left|k_{\pm\rho}\frac{\partial[X_{i}]_{\rm
ss}}{\partial k_{\pm\rho}}\right|,$ (60) $\displaystyle{\mathcal{B}}^{\rm
M}(\rho,X_{i})$ $\displaystyle=\frac{1}{D_{i}}\left|\frac{\partial[X_{i}]_{\rm
ss}}{\partial B_{\rho}}\right|,$ (61)

where the superscript ‘M’ signifies macroscopic. Recalling that the Chemical
Rate Equation can have multiple stable fixed points, we avoid this
complication by only considering the macroscopic limit in situations where
there is a unique stable solution.

####  IV.3.1 Single chemical cycle and single conservation law $(\delta=0)$

Figure 5:  Stoichiometric compatibility class of chemical reaction network
(62) with total number $\Lambda=n_{A}+n_{B}+n_{C}=3$. It consists of four
isolated connected components, including three isolated states $(0,3,0)$,
$(0,0,3)$, and $(1,0,2)$. The shaded area highlights the largest isolated
connected component.

We first consider a deficiency zero, nonlinear CRN with a single chemical
cycle and a single conservation law, with complex-reaction graph

$\schemestart\chemfig{2A}\arrow{<=>[\tiny{\hbox{k_{+2}}}][\tiny{\hbox{k_{-2}}}]}\chemfig{B+C}\arrow{<=>[\tiny{\hbox{k_{-3}}}][\tiny{\hbox{k_{+3}}}]}[125]\chemfig{A+B}\arrow{<=>[\tiny{\hbox{k_{-1}}}][\tiny{\hbox{k_{+1}}}]}[235]\schemestop.$
(62)

This CRN is constructed from a slight modification of the example depicted in
Fig. 1 to make it deficiency zero, but with the same stoichiometric matrix and
conservation law. The total number of chemical species
$\Lambda=n_{A}+n_{B}+n_{C}$ is conserved, and the only chemical cycle is
visible in the complex-reaction graph. For fixed $\Lambda$ the stoichiometric
compatibility class consists of a large isolated connected component and three
isolated states,
$(n_{A},n_{B},n_{C})\in\\{(0,\Lambda,0),~{}(0,0,\Lambda),~{}(1,0,\Lambda-1)\\}$.
The case of $\Lambda=3$ is depicted in Fig. 5. The steady-state probabilities
of the three isolated states are fixed by the initial distribution and do not
change when the dynamics are perturbed. Accordingly, it is sufficient to
consider only the large isolated connected component to examine the trade-offs
since the isolated states do not contribute to response and number
fluctuations.

By choosing a basis where $\\{X_{1}=A,X_{2}=B,X_{3}=C\\}$ and the two
independent chemical reactions $\\{A+B\rightleftharpoons
2A,2A\rightleftharpoons B+C\\}$, we have the reduced stoichiometric matrix

$\mathsf{S}^{\prime}=\begin{pmatrix}1&-2\\\ -1&1\\\ 0&1\end{pmatrix}$ (63)

and its left inverse

$(\mathsf{S}^{\prime})^{-1}=\begin{pmatrix}0&-1&1\\\
-1/3&-1/3&2/3\end{pmatrix}.$ (64)

Plugging the components of $(\mathsf{S}^{\prime})^{-1}$ into (54) and (55) and
using ${\rm Cov}\\{n_{X},\Lambda\\}=0$, we have inequalities

$\displaystyle\left|k_{\pm\rho}\frac{\partial\langle n_{X}\rangle}{\partial
k_{\pm\rho}}\right|$ $\displaystyle\leq|{\rm Cov}\\{n_{X},n_{B}-n_{C}\\}|$
(65) $\displaystyle~{}~{}~{}+|{\rm Cov}\\{n_{X},n_{C}\\}|$

and

$\displaystyle\left|\frac{\partial\langle n_{X}\rangle}{\partial
B_{\rho}}\right|$ $\displaystyle\leq(|{\rm Cov}\\{n_{X},n_{B}-n_{C}\\}|$ (66)
$\displaystyle~{}~{}~{}+|{\rm
Cov}\\{n_{X},n_{C}\\}|)\tanh\left(\frac{\mathcal{F}}{4}\right),$

where $X\in\\{A,B,C\\}$ and $\mathcal{F}$ is the absolute value of the
thermodynamic force associated with the sole chemical cycle (17).

Figure 6: The maximum quality factors among different choices of $\rho$ for
the model in (62) at a fixed
$(X,\Lambda,\\{k_{\pm\rho}\\}_{\rho=1}^{3},\Omega)$. Each open symbol is
obtained from a set of randomly sampled rate constants
$\\{k_{\pm\rho}\\}_{\rho=1}^{3}$ and system size $\Omega$, the logs of which
are uniformly sampled from the range $[-5,5]$, with $\Lambda$ sampled from
$\\{3,4,5,6,7\\}$. The number of open symbols at a fixed $(X,\Lambda)$ in each
panel is $200$. The filled symbols are the maximum values at a fixed
$\mathcal{F}$ and chemical species, which are numerically optimized using
simulated annealing. The perturbed edges used for the optimization are
$A+B\rightharpoonup 2A$ for (a) and $A+B\rightleftharpoons 2A$ for (b), and
$\Lambda$ is fixed by 3. The dashed lines represent $1$ and $2$ in (a),
$\tanh(\mathcal{F}/4)$ and $2\tanh(\mathcal{F}/4)$ in (b).

We numerically calculated the two quality factors ${\mathcal{K}}(\pm\rho,X)$
and ${\mathcal{B}}(\rho,X)$ for each species $X$ and reaction channel $\rho$
for random sets of system parameters
$(\Lambda,\\{k_{\pm\rho}\\}_{\rho=1}^{3},\Omega)$. Then for each random set of
system parameters, we plot in Fig. 6 the maximum quality factors for each
species $X$, maximized over perturbed reaction channels, depicted as open
symbols. The figure suggests the existence of upper bounds on the response
that do not appear to depend on $\Lambda$. To confirm this observation and
determine the shape of the upper bound, we numerically optimized the quality
factors using simulated annealing [69]. The results from the optimization are
drawn as filled symbols in Fig. 6, which align with the trade-offs (2) and (3)
with $g(\mathcal{F})=\tanh(\mathcal{F}/4)$ regardless of $(\rho,X)$. A closer
investigation further reveals that the two prefactors $\alpha_{1}$ and
$\alpha_{2}$ coincide for each pair $(\rho,X)$ (see TABLE 1).

|  chemical reaction ($\rho$)
---|---
chemical species ($X$) |  $\pm 1$ |  $\pm 2$ |  $\pm 3$
$A$ |  (1, 1) |  (1, 1) |  (1, 1)
$B$ |  (2, 2) |  (1, 1) |  (2, 2)
$C$ |  (2, 2) |  (1, 1) |  (2, 2)
Table 1: Numerically estimated values of ($\alpha_{1}$, $\alpha_{2}$) for each
pair $(\rho,X)$ in (62).

To investigate the relationship between our covariance bounds (65) and the
compact trade-offs (2), we compare the tightness of the bounds in Fig. 7 by
plotting the ratios of the response to fluctuations for single rate
perturbations; similar conclusions would be reached by considering the force-
dependent symmetric perturbations bounded in (66) and (3). The numerical
results validate both inequalities as all samples are confined to the unit
square. Moreover we see that one bound is generally not tighter than the other
as samples fall above and below the diagonal, implying that the covariance
bounds (65) and (66) we have derived are not sufficiently tight in this case
to prove the compact trade-offs (2) and (3). Nevertheless, the enveloping
guidelines shown in red based on the inequalities, $(\alpha_{1}/2){\rm
Var}\\{n_{X}\\}\leq|{\rm Cov}\\{n_{X},n_{B}-n_{C}\\}|+|{\rm
Cov}\\{n_{X},n_{C}\\}|\leq 3\alpha_{1}{\rm Var}\\{n_{X}\\}$, imply the bounds
are not completely independent.

Figure 7: Comparison of the tightness of the proven bound (65) and conjectured
bound (2) for the model in (62). The covariance bound in the vertical axis
refers to the right-hand side of (65). Each symbol is obtained from a set of
randomly sampled rate constants $\\{k_{\pm\rho}\\}_{\rho=1}^{3}$ and system
size $\Omega$, the logs of which are uniformly sampled from the range
$[-5,5]$, with $\Lambda=3$. The coefficient $\alpha_{1}$ is used as the
estimated values in TABLE 1. The number of symbols for a fixed pair of
$(X,\pm\rho)$ is $500$. The red guidelines have slopes of $1/3$, $1$, and $2$.

####  IV.3.2 Multiple conservation laws $(\delta=0)$

Figure 8:  The maximum quality factors among different choices of ($\rho$,
$X$) for the model in (67) at a fixed
$(\Lambda_{E},\Lambda_{S},\\{k_{\pm\rho}\\}_{\rho=1}^{4},\Omega)$. Each symbol
is obtained from a set of randomly sampled rate constants
$\\{k_{\pm\rho}\\}_{\rho=1}^{4}$ and system size $\Omega$, the logs of which
are uniformly sampled from $[-5,5]$. The total number of enzymes and
substrates, $(\Lambda_{E},\Lambda_{S})$, take values from
$\\{(1,5),(2,4),(3,3),(4,2),(5,1)\\}$. The number of symbols at a fixed
$(\Lambda_{E},\Lambda_{S})$ in each panel is $200$. The thick lines represent
$1$ in (a), and $\tanh(\mathcal{F}/4)$ in (b).

To see whether and how the number of conservation laws affects the trade-offs,
we consider a CRN with multiple conservation laws, but a single chemical
cycle. Our model is composed of a substrate that can exist in an unmodified
$S$ and modified form $P$. In cells, this modification is often the addition
of a phosphate group, but we will not make that specification here. The
modification and demodification are carried out by a single enzyme $E$ through
two distinct intermediate species $ES_{1}$ and $ES_{2}$, resulting in the
complex-reaction graph

$\schemestart\chemfig{ES_{2}}\arrow{<=>[\tiny{\hbox{k_{-3}}}][\tiny{\hbox{k_{+3}}}]}\chemfig{E+P}\arrow{<=>[\tiny{\hbox{k_{-2}}}][\tiny{\hbox{k_{+2}}}]}[90]\chemfig{ES_{1}}\arrow{<=>[\tiny{\hbox{k_{+1}}}][\tiny{\hbox{k_{-1}}}]}[180]\chemfig{E+S}\arrow{<=>[\tiny{\hbox{k_{-4}}}][\tiny{\hbox{k_{+4}}}]}[270]\schemestop.$
(67)

This CRN was inspired by the Goldbeter-Koshland model for the phosphorylation-
dephosphorylation cycle (PdPC), which is a known mechanism for ultra-
sensitivity. The model was modified to make it reversible and deficiency zero.
Later in Sec. V.2, we will analyze the full (reversible) Goldbeter-Koshland
model.

The modification-demodification cycle in (67) has two independent conservation
laws: the total numbers of enzymes $\Lambda_{E}=n_{E}+n_{ES_{1}}+n_{ES_{2}}$
and substrates $\Lambda_{S}=n_{S}+n_{P}+n_{ES_{1}}+n_{ES_{2}}$. By choosing a
basis where $\\{X_{1}=E,X_{2}=S,X_{3}=P,X_{4}=ES_{1},X_{5}=ES_{2}\\}$ and
three independent reaction channels $\\{E+S\rightleftharpoons
ES_{1},ES_{1}\rightleftharpoons E+P,E+P\rightleftharpoons ES_{2}\\}$, we have
the reduced stoichiometric matrix

$\mathsf{S}^{\prime}=\begin{pmatrix}-1&1&-1\\\ -1&0&0\\\ 0&1&-1\\\ 1&-1&0\\\
0&0&1\end{pmatrix}$ (68)

and its left inverse

$(\mathsf{S}^{\prime})^{-1}=\begin{pmatrix}-1/4&-5/8&3/8&1/8&1/8\\\
0&-1/2&1/2&-1/2&1/2\\\ -1/4&-1/8&-1/8&-3/8&5/8\end{pmatrix}.$ (69)

Plugging the components of $(\mathsf{S}^{\prime})^{-1}$ into (54) and (55) and
using ${\rm Cov}\\{n_{X},\Lambda_{E}\\}={\rm Cov}\\{n_{X},\Lambda_{S}\\}=0$,
we have the inequalities

$\displaystyle\left|k_{\pm\rho}\frac{\partial\langle n_{X}\rangle}{\partial
k_{\pm\rho}}\right|$ $\displaystyle\leq|{\rm Cov}\\{n_{X},n_{S}\\}|+|{\rm
Cov}\\{n_{X},n_{ES_{2}}\\}|$ (70) $\displaystyle~{}~{}~{}+|{\rm
Cov}\\{n_{X},n_{P}+n_{ES_{2}}\\}|$

and

$\displaystyle\left|\frac{\partial\langle n_{X}\rangle}{\partial
B_{\rho}}\right|$ $\displaystyle\leq(|{\rm Cov}\\{n_{X},n_{S}\\}|+|{\rm
Cov}\\{n_{X},n_{ES_{2}}\\}|$ (71) $\displaystyle~{}~{}~{}+|{\rm
Cov}\\{n_{X},n_{P}+n_{ES_{2}}\\}|)\tanh\left(\frac{\mathcal{F}}{4}\right),$

where $X\in\\{E,S,P,ES_{1},ES_{2}\\}$ and $\mathcal{F}$ is the absolute value
of the thermodynamic force associated with the only chemical cycle (17).

We numerically calculated the maximum values ${\mathcal{K}}(\pm\rho,X)$ and
${\mathcal{B}}(\rho,X)$ among different pairs of ($\rho$, $X$) at fixed system
parameters $(\Lambda_{E},\Lambda_{S},\\{k_{\pm\rho}\\}_{\rho=1}^{4},\Omega)$.
Rate constants $\\{k_{\pm\rho}\\}_{\rho=1}^{4}$ and system size $\Omega$ are
randomly chosen. The total numbers of enzymes and substrates are selected from
$(\Lambda_{E},\Lambda_{S})\in\\{(1,5),~{}(2,4),~{}(3,3),~{}(4,2),~{}(5,1)\\}$.
The results are shown in Fig. 8, which demonstrate that the maximum quality
factors are bounded from above by ${\mathcal{K}}\leq 1$ and
${\mathcal{B}}\leq\tanh(\mathcal{F}/4)$. A detailed investigation further
shows that the model (67) obeys trade-offs, $|k_{\pm}\partial\langle
n_{X}\rangle/\partial k_{\pm\rho}|\leq 1$ and $|\partial\langle
n_{X}\rangle/\partial B_{\rho}|\leq\tanh(\mathcal{F}/4)$ no matter the choice
of pair $(\rho,X)$.

In Fig. 9, we numerically compare the covariance bound (70) with the compact
trade-offs (2) for a single rate perturbation by plotting the ratio of
response to the fluctuations. The validity of both bonds is apparent as all
samples fall inside the unit square. Here, however, all samples fall below the
diagonal implying that covariance bound (70) is tighter then the compact
trade-off (2) in contrast to what was observed for the model in (62). Together
the results in Figs. 7 and 9 suggest that the relative tightness of the two
types of bounds depends on the model.

Figure 9: Comparison of the tightness of the proven bound (70) and conjectured
bound (2) for the model in (67). The covariance bound in the vertical axis
refers to the right-hand side of (70). Each symbol is obtained from a set of
randomly sampled rate constants $\\{k_{\pm\rho}\\}_{\rho=1}^{4}$ and system
size $\Omega$, the logs of which are uniformly sampled from the range
$[-5,5]$, with $(\Lambda_{E},\Lambda_{S})=(3,3)$. The coefficient $\alpha_{1}$
is set to $1$ for all pairs of $(X,\pm\rho)$. The number of symbols for a
fixed pair of $(X,\pm\rho)$ is $200$. The guidelines have slopes of $1/3$ and
$1$.

####  IV.3.3 Multiple chemical cycles $(\delta=0)$

Figure 10: The maximum quality factors among different choices of
$\bm{(}\rho,X(\neq W)\bm{)}$ for the model in (72) at a fixed
$(S_{T},\\{k_{\pm\rho}\\}_{\rho=1}^{7})$ (circles). The maximum quality
factors for $X=W$ are calculated separately (squares). Each symbol is obtained
from a set of randomly sampled rate constants $\\{k_{\pm\rho}\\}_{\rho=1}^{7}$
and the total concentration of substrates $S_{T}$, the logs of which are
uniformly sampled within the range $[-5,5]$. The number of each symbol in each
panel is $10^{3}$. The thick lines represent $1/2$ and $1$ in (a), and
$(1/2)\tanh(\mathcal{F}/4)$ and $\tanh(\mathcal{F}/4)$ in (b).

We further alter the modification-demodification cycle to include the effect
of multiple chemical cycles:

$\schemestart\chemfig{ES_{2}}\arrow{<=>[\tiny{\hbox{k_{-3}}}][\tiny{\hbox{k_{+3}}}]}\chemfig{E+P}\arrow{<=>[\tiny{\hbox{k_{-2}}}][\tiny{\hbox{k_{+2}}}]}[90]\chemfig{ES_{1}}\arrow{<=>[\tiny{\hbox{k_{+1}}}][\tiny{\hbox{k_{-1}}}]}[180]\chemfig{E+S}\arrow{<=>[\tiny{\hbox{k_{-4}}}][\tiny{\hbox{k_{+4}}}]}[270]\schemestop\hskip
28.45274pt\schemestart\chemfig{E}\arrow{<=>[\tiny{\hbox{k_{+6}}}][\tiny{\hbox{k_{-6}}}]}\chemfig{2W}\arrow{<=>[\tiny{\hbox{k_{-7}}}][\tiny{\hbox{k_{+7}}}]}[120]\chemfig{\emptyset}\arrow{<=>[\tiny{\hbox{k_{-5}}}][\tiny{\hbox{k_{+5}}}]}[240]\schemestop.$
(72)

The enzyme $E$ can now also be reversibly converted into a waste form $W$, and
both the enzyme and the waste species are in contact with chemical reservoirs.
Two chemical cycles are visible in the complex-reaction graph. We denote the
associated thermodynamic forces in the left and right linkage classes as
$F_{l}$ and $F_{r}$. The maximum thermodynamic force $\mathcal{F}$ (17) is
equal to $|F_{l}|$ ($|F_{r}|$) if the perturbed reaction is in the left
(right, resp.) linkage class. The total number of enzymes is no longer
conserved, while the conservation law for substrates,
$\Lambda_{S}=n_{S}+n_{P}+n_{ES_{1}}+n_{ES_{2}}$, remains. Since the numbers of
$E$ and $W$ are not limited, the stoichiometric compatibility class is
infinitely large even if $\Lambda_{S}$ is small. At the same time, the zero
deficiency of the CRN ensures that the steady-state (positive) solution of the
Chemical Rate Equations is unique [33]. Instead of solving the Chemical Master
Equation with small $\Lambda_{S}$, we examine the validity of the macroscopic
version of trade-offs by considering quality factors ${\mathcal{K}}^{\rm
M}(\pm\rho,X)$ and ${\mathcal{B}}^{\rm M}(\rho,X)$. For the macroscopic
concentrations, we denote the conserved total concentration of substrates by
$S_{T}=[S]+[P]+[ES_{1}]+[ES_{2}]$.

We numerically calculated the maximum values of ${\mathcal{K}}^{\rm
M}(\pm\rho,X)$ and ${\mathcal{B}}^{\rm M}(\rho,X)$ for different pairs
$(\rho,X)$ for randomly-sampled combinations of system parameters
$(S_{T},\\{k_{\pm\rho}\\}_{\rho=1}^{7})$. The numerical results are presented
in Fig. 10 in which the cases $X=W$ and $X\neq W$ are treated separately. The
maximum quality factor of the wasted form of enzyme $W$ is bounded from above
by 1/2 while those of other chemical species are bounded by 1. Figure 10 and
TABLE 2 suggest the model in (72) obeys the macroscopic version of trade-offs
${\mathcal{K}}^{\rm M}\leq\alpha_{1}$ and ${\mathcal{B}}^{\rm
M}\leq\alpha_{2}\tanh(\mathcal{F}/4)$ with $\alpha_{1}=\alpha_{2}=1/2$ for
$X=W$ and $\alpha_{1}=\alpha_{2}=1$ for $X\neq W$.

The suppressed response of $W$ is a nonlinear effect. Since the steady-state
distribution is a product form but no conservation law constrains $n_{W}$, the
fluctuations of $n_{W}$ are decoupled from other chemical species and thus
${\rm Cov}_{\Gamma}\\{n_{W},n_{j}\\}={\rm
Var}_{\Gamma}\\{n_{W}\\}\delta_{Wj}$. As a result, the fluctuation-response
relation (41) for $W$ simplifies to

$\frac{\partial\langle n_{W}\rangle_{\Gamma}}{\partial\lambda}=\frac{1}{2}{\rm
Var}_{\Gamma}\\{n_{W}\\}\frac{\partial}{\partial\lambda}\ln[W].$ (73)

which in turn leads to trade-offs with $\alpha_{1}=\alpha_{2}=1/2$. Since a
compatibility class is composed of a single connected component in this model,
applying the matrix tree theorem, we can actually prove trade-offs for $E$ and
$W$, whose numbers are not conserved for any $\Lambda_{S}$.

The above statements can be generalized. For a deficiency zero CRN, if the
number of a chemical species $A$ is not conserved and a complex $zA$ with a
positive integer $z$ is included in the CRN, the fluctuation-response relation
(41) for $A$ can be written as $\partial_{\lambda}\langle
n_{A}\rangle_{\Gamma}=(1/z){\rm
Var}_{\Gamma}\\{n_{A}\\}\partial_{\lambda}\ln[A]^{z}$. By noticing
$\Psi_{zA}([\bm{X}]_{\rm ss})=[A]^{z}$ and $\Psi_{\emptyset}([\bm{X}]_{\rm
ss})=1$ and applying the matrix-tree theorem, we can derive trade-offs
$|k_{\pm\rho}\partial\langle n_{A}\rangle_{\Gamma}/\partial
k_{\pm\rho}|\leq{\rm Var}_{\Gamma}\\{n_{A}\\}/z$ and $|\partial\langle
n_{A}\rangle_{\Gamma}/\partial B_{\rho}|\leq({\rm
Var}_{\Gamma}\\{n_{A}\\}/z)\tanh(\mathcal{F}/4)$, provided the complexes
$\emptyset$ and $zA$ are in the same linkage class.

|  chemical reaction $(\rho)$
---|---
$X$ | $\pm 1$ |  $\pm 2$ |  $\pm 3$ |  $\pm 4$ |  $\pm 5$ |  $\pm 6$ |  $\pm 7$
$E$ | (0, 0) |  (0, 0) |  (0, 0) |  (0, 0) |  (1, 1) |  (1, 1) |  (1, 1)
$S$ | (1, 1) |  (1, 1) |  (1, 1) |  (1, 1) |  (1, 1) |  (1, 1) |  (1, 1)
$P$ | (1, 1) |  (1, 1) |  (1, 1) |  (1, 1) |  (1, 1) |  (1, 1) |  (1, 1)
$ES_{1}$ | (1, 1) |  (1, 1) |  (1, 1) |  (1, 1) |  (1, 1) |  (1, 1) |  (1, 1)
$ES_{2}$ | (1, 1) |  (1, 1) |  (1, 1) |  (1, 1) |  (1, 1) |  (1, 1) |  (1, 1)
$W$ | (0, 0) |  (0, 0) |  (0, 0) |  (0, 0) |  ($\frac{1}{2}$, $\frac{1}{2}$) |  ($\frac{1}{2}$, $\frac{1}{2}$) |  ($\frac{1}{2}$, $\frac{1}{2}$)
Table 2: Numerically estimated values of ($\alpha_{1}$, $\alpha_{2}$) to each
given pair $(\rho,X)$ for the model in (72).

##  V Deficiency Nonzero Chemical Reaction Networks

Analytical approaches are further limited for deficiency nonzero CRNs since
the general form of the steady-state distribution is unknown. So, we
exclusively rely on numerical methods here. The validity of trade-offs (2) and
(3) are examined by calculating the quality factors defined in the previous
section for two models with different features.

###  V.1 Multiple chemical cycles $(\delta\neq 0)$

The Schögl model illustrated in Sec. III is a deficiency nonzero model with a
single chemical cycle, for which the trade-offs between fluctuation and
response are proven. In order to see the effect of multiple thermodynamic
forces, we consider a deficiency nonzero model with two independent chemical
cycles, whose complex-reaction graph is given by

$\schemestart\chemfig{X}\arrow{<=>[\tiny{\hbox{k_{+2}}}][\tiny{\hbox{k_{-2}}}]}\chemfig{2X}\arrow{<=>[\tiny{\hbox{k_{-3}}}][\tiny{\hbox{k_{+3}}}]}[120]\chemfig{\emptyset}\arrow{<=>[\tiny{\hbox{k_{-1}}}][\tiny{\hbox{k_{+1}}}]}[240]\schemestop.$
(74)

Apart from the cycle visible in the complex-reaction graph with thermodynamic
force $F_{1}=\ln\bm{(}(k_{+1}k_{+2}k_{+3})/(k_{-1}k_{-2}k_{-3})\bm{)}$, an
independent cycle with thermodynamic force
$F_{2}=\ln\bm{(}(k_{-1}k_{+2})/(k_{+1}k_{-2})\bm{)}$ is further identified by
analysis of the stoichiometric matrix. Along with $F_{1}$ and $F_{2}$, two
more cycles characterized by thermodynamics forces
$F_{1}-F_{2}=\ln\bm{(}(k_{+1}^{2}k_{+3})/(k_{-1}^{2}k_{-3})\bm{)}$ and
$F_{1}+F_{2}=\ln\bm{(}(k_{+2}^{2}k_{+3})/(k_{-2}^{2}k_{-3})\bm{)}$ constitute
the entire set of cycles that reside in the microscopic state space. The
maximum thermodynamic force $\mathcal{F}$ depends on the perturbed edge (17).
For example, if the chemical reaction $\emptyset\rightleftharpoons X$ is
perturbed, $\mathcal{F}=\max\\{|F_{1}|,|F_{2}|,|F_{1}-F_{2}|\\}$.

The number of chemical species $X$ is not conserved with no upper bound. Thus
the only stoichiometric compatibility class is infinitely large. We examine
the validity of the macroscopic version of the trade-offs (2) and (3). The
steady-state concentration $[X]_{\rm ss}$ is the unique positive solution of
the quadratic equation

$(k_{+1}+2k_{-3})+(-k_{-1}+k_{+2})x-(k_{-2}+2k_{+3})x^{2}=0.$ (75)

We numerically calculate the quality factors ${\mathcal{K}}^{\rm
M}(\pm\rho,X)$ and $\mathcal{B}^{\rm M}(\rho,X)$ for a fixed set of rate
constants $\\{k_{\pm\rho}\\}_{\rho=1}^{3}$. The logs of the rate constants are
randomly sampled from a uniform distribution over the range $[-10,10]$. From
these $10^{6}$ data points, the maximum response is singled out from each
interval $\mathcal{F}\in[x,x+0.5]$ with $x\in\\{0,0.5,1,\cdots,19.5\\}$ and
plotted in Fig. 11. Numerical results are consistent with the macroscopic
version of the trade-offs. The prefactors and the function $g(\mathcal{F})$
vary for different choices of the perturbed chemical reaction. For example,
when the chemical reaction $2X\rightleftharpoons\emptyset$ is perturbed, the
quality factors do not exceed $2/3$ and $g(\mathcal{F})$ differs from
$\tanh(\mathcal{F}/4)$. Unlike in the deficiency zero CRNs without a
conservation law (see Sec. IV.3.3), the effect of nonlinearity is not to
simply reduce the response function by an integer factor.

Figure 11: Quality factors for the model in (74). Raw data are calculated from
$10^{6}$ sets of randomly sampled rate constants, the logs of which are
uniformly sampled over the range $[-10,10]$. The symbols represent the maximum
values of quality factors selected from the sampled data within range
$\mathcal{F}\in[x,x+0.5]$ where $x\in\\{0,0.5,1,\cdots,19.5\\}$. The thick
lines represent $\\{1/2,2/3,1\\}$ in (a) and
$\\{(2/3)\tanh(\mathcal{F}/4),\tanh(\mathcal{F}/4)\\}$ in (b).

###  V.2 Multiple conservation laws $(\delta\neq 0)$

Our inspiration here is the Goldbeter-Koshland model for the PdPC. In its
original formulation, a substrate can exist in an unmodified $S$ form and a
modified or phosphorylated form $P$. Phosphorylation is catalyzed by a kinase
enzyme $E_{1}$, whereas a phosphatase enzyme $E_{2}$ catalyzes
dephosphorylation. By introducing two intermediate species $E_{1}S$ and
$E_{2}P$, Goldbeter and Koshland considered two chains of irreversible
chemical reactions [1]

$\displaystyle E_{1}+S\rightleftharpoons E_{1}S\rightharpoonup E_{1}+P,$ (76)
$\displaystyle E_{2}+P\rightleftharpoons E_{2}P\rightharpoonup E_{2}+S$

as a mechanism to explain switch-like behavior observed in some biological
scenarios. The high sensitivity, called zero-order ultrasensitivity [70], is
observed in a regime where the total numbers of both enzymes are much lower
than that of the substrates and the chemical cycle is highly irreversible [2].
We will examine the emergence of this high sensitivity from the point of view
of the trade-offs.

The irreversible chemical reactions of the model imply infinite entropy
production and are thermodynamically implausible. In line with [2], we
consider a reversible version with complex-reaction graph

$\displaystyle
E_{1}+S\xrightleftharpoons[k_{-1}]{k_{+1}}E_{1}S\xrightleftharpoons[k_{-2}]{k_{+2}}E_{1}+P,$
(77) $\displaystyle
E_{2}+P\xrightleftharpoons[k_{-3}]{k_{+3}}E_{2}P\xrightleftharpoons[k_{-4}]{k_{+4}}E_{2}+S.$

The deficiency of the model is one. The stoichiometric compatibility class is
specified by three conservation laws: the total numbers of the respective
enzymes, $\Lambda_{E_{1}}=n_{E_{1}}+n_{E_{1}S}$ and
$\Lambda_{E_{2}}=n_{E_{2}}+n_{E_{2}P}$, as well as the total number of
substrates, $\Lambda_{S}=n_{S}+n_{P}+n_{E_{1}S}+n_{E_{2}P}$. The thermodynamic
force associated with the only chemical cycle invisible in the complex-
reaction graph is $F=\ln\bm{(}\prod_{i=1}^{4}(k_{+i}/k_{-i})\bm{)}$, which is
the chemical potential difference between ATP and ADP+Pi in the unit of
$k_{B}T$.

In order to examine the validity of the trade-offs (2) and (3), we numerically
calculate the quality factors $\mathcal{K}(\pm\rho,X)$ and
$\mathcal{B}(\rho,X)$ at $10^{6}$ fixed values of
$(\\{k_{\pm\rho}\\}_{\rho=1}^{4},\Omega)$. The logs of the rate constants
$\\{k_{\pm\rho}\\}_{\rho=1}^{4}$ and system volume $\Omega$ are randomly
sampled from a uniform distribution over the range $[-5,5]$, from which each
$\mathcal{K}(\pm\rho,X)$ and $\mathcal{B}(\rho,X)$ are determined for all
pairs $(\rho,X)$. For each value of the system parameters, we then select the
maximum of $\mathcal{K}(\pm\rho,X)$ and $\mathcal{B}(\rho,X)$ over all pairs
$(\rho,X)$ to obtain one sample of the max-quality factors. Plotted in Figure
12 is then the maximum quality factor within each interval
$\mathcal{F}\in[x,x+0.5]$ for $x\in\\{0,0.5,1,\cdots,19.5\\}$. The response
depends on the total numbers of enzymes and substrates,
$(\Lambda_{E_{1}},\Lambda_{E_{2}},\Lambda_{S})\in\\{(1,1,3),~{}(1,1,5),~{}(3,3,3),~{}(1,1,7)\\}$,
which also determine the respective size of the stoichiometric compatibility
classes, $\\{12,20,20,28\\}$. Two features distinct from deficiency zero CRNs
can be observed in Fig. 12.

We first compare the cases
$(\Lambda_{E_{1}},\Lambda_{E_{2}},\Lambda_{S})=(1,1,5)$ and $(3,3,3)$ (orange
and green), which have the same numbers of microscopic states in the
stoichiometric compatibility class, 20. Although the numbers of microscopic
states are the same, the structural differences in the microscopic state space
lead to a different quantitative relationship between response and
fluctuations. While the maximum values of the quality factors are slightly
greater than 1 when the number of enzymes is comparable with that of the
substrates, they appear to grow to a much larger value with thermodynamic
force when the enzymes are scarce. This high response in Fig. 12 observed when
the enzymes are scarce is reminiscent of zero-order ultrasensitivity; however,
in the biochemical literature ultrasensitivity is defined by the response to a
change in the total numbers of enzymes $\Lambda_{E_{1}}$ and/or
$\Lambda_{E_{2}}$ [70], which is different from what we consider here. With
that caveat in mind, a recent theoretical study has shown that in the limit
that substrate concentration goes to infinity (infinitely scarce enzymes), the
log-sensitivity to changes in the total enzyme concentration can actually grow
as fast as exponential with thermodynamic driving [5]. This prediction appears
consistent with our observation that our quality factors grow over a wider
range of thermodynamic forces when there are less enzymes.

Next, we compare the cases
$(\Lambda_{E_{1}},\Lambda_{E_{2}},\Lambda_{S})=(1,1,3),~{}(1,1,5)$, and
$(1,1,7)$. On the one hand, the quality factors appear to obey the trade-offs
(2) and (3) for a fixed $\Lambda_{S}$. On the other hand, as the total number
of substrates increases, the maximum values of the quality factors increase as
well. This suggests that the prefactors $\alpha_{1,2}$ in the trade-offs
depend on system details, such as the ratio between the numbers of substrates
and enzymes, which deform the structure of the stoichiometric compatibility
class.

Figure 12: The maximum quality factors among different choices of ($\rho$,
$X$) for the model in (77) with
$(\Lambda_{E_{1}},\Lambda_{E_{2}},\Lambda_{S})\in\\{(1,1,3),~{}(1,1,5),~{}(3,3,3),~{}(1,1,7)\\}$
at a fixed $(\\{k_{\pm\rho}\\}_{\rho=1}^{4},\Omega)$. Raw data are calculated
from $10^{6}$ sets of randomly sampled rate constants
$\\{k_{\pm\rho}\\}_{\rho=1}^{4}$ and system size $\Omega$, the logs of which
are uniformly sampled within the range $[-5,5]$. The symbols represent the
maximum values of quality factors selected out of raw data within range
$\mathcal{F}\in[x,x+0.5]$ where $x\in\\{0,0.5,1,\cdots,19.5\\}$ The thick
lines represent $1$ in (a), and $\tanh(\mathcal{F}/4)$ in (b).

##  VI Concluding Remarks

In this work, we investigated fluctuation-response trade-offs for CRNs in
nonequilibrium steady-states. We proved the trade-offs, (2) and (3), for
linear CRNs (Sec. IV) and a class of nonlinear CRNs with a single chemical
species (Sec. III). For deficiency zero nonlinear CRNs, the response is still
bounded by fluctuations and thermodynamic forces via the inequalities (54) and
(55), though we have not been able to connect them to the compact trade-offs
(2) and (3). Nevertheless, numerical calculations for nonlinear models suggest
that the trade-offs (2) and (3) may continue to hold. To support this
observation, we have analyzed various models with differing aspects, such as
whether there are conservation laws, the number of independent thermodynamic
forces, the number of chemical species, and the deficiency of the CRN.
Regardless of different aspects, deficiency zero CRNs appear to share the
following properties: (i) $\alpha_{1}=\alpha_{2}$, (ii)
$g(\mathcal{F})=\tanh(\mathcal{F}/4)$, and (iii) the prefactors do not depend
on the size of the stoichiometric compatibility class. In contrast, numerical
results for deficiency-nonzero CRNs show that the prefactors depend on system
details that affect the structure of the stoichiometric compatibility class
(Sec. V).

Trade-offs between number fluctuations and response suggest limits to the
accuracy and sensitivity of biochemical systems. Thus, our results pave the
way for identifying fundamental limitations to nonequilibrium biochemical
processes. Moving forward new analytic techniques may be required, especially
for nonlinear CRNs where the general form of the steady-state distribution is
unknown. Such advances would allow us to discern how characteristics of CRNs,
such as conservation laws or deficiency, restrict the topological structure of
the microscopic state space and thus affect the response.

## Acknowledgements

H.-M.C. was supported in part by a KIAS Individual Grant (PG089401) at Korea
Institute for Advanced Study. J.M.H. acknowledges support from the National
Science Foundation under Grant No. 2142466.

##  Appendix A Numerical Methods

When the size of the isolated connected component of a stoichiometric
compatibility class is sufficiently small, we first numerically solve the
Chemical Master Equation (9) for the steady-state distribution,
$\hat{\mathcal{L}}\pi_{\Gamma}(\bm{n})=0$ with normalization condition
$\sum_{\bm{n}\in\Gamma}\pi_{\Gamma}({\bm{n}})=1$, which is a system of linear
equations. To determine the steady-state response
$\partial_{\lambda}\pi_{\Gamma}({\bm{n}})$, we differentiate the Chemical
Master Equation with respect to the perturbation parameter $\lambda\in\\{\ln
k_{+\rho},\ln k_{-\rho},B_{\rho}\\}$, which leads to a system of inhomogeneous
linear equations for the response

$\hat{\mathcal{L}}\frac{\partial\pi_{\Gamma}(\bm{n})}{\partial\lambda}+\frac{\partial\hat{\mathcal{L}}}{\partial\lambda}\pi_{\Gamma}(\bm{n})=0,$
(78)

which are solved numerically together with the condition
$\sum_{\bm{n}\in\Gamma}\partial_{\lambda}\pi_{\Gamma}({\bm{n}})=0$ that is
required by normalization of the steady-state distribution.

When the system size is large, numerical solutions to the Chemical Master
Equation (9) are infeasible. In these scenarios, we first numerically solve
the Chemical Rate Equation for the steady-state concentration $[{\bm{X}}_{\rm
ss}]$. To determine the steady-state response
$\partial_{\lambda}[{\bm{X}}_{\rm ss}]$, we develop a nonlinear equation for
it by differentiating the Chemical Rate Equation

$\displaystyle\sum_{\rho\in\mathcal{R}}\Delta\bm{\nu}_{\rho}\left(\frac{\partial
k_{+\rho}}{\partial\lambda}[\bm{X}]_{\rm
ss}^{\bm{\nu}^{+}_{\rho}}-\frac{\partial
k_{-\rho}}{\partial\lambda}[\bm{X}]_{\rm ss}^{\bm{\nu}^{-}_{\rho}}\right)$
(79)
$\displaystyle+\sum_{\rho\in\mathcal{R}}\Delta\bm{\nu}_{\rho}\left(k_{+\rho}\frac{\partial[\bm{X}]_{\rm
ss}^{\bm{\nu}^{+}_{\rho}}}{\partial\lambda}-k_{-\rho}\frac{\partial[\bm{X}]_{\rm
ss}^{\bm{\nu}^{-}_{\rho}}}{\partial\lambda}\right)=\bm{0},$

which allows us to numerically determine the response. However, the Chemical
Rate Equation does not encode any information about the fluctuations.

To ascertain the scaled variance $D_{i}$, we use Van Kampen’s system size
expansion, also known as the Linear Noise Approximation [54], to approximate
the Chemical Master Equation as a diffusion process, which we review here. In
this approximation, the numbers of chemical species are assumed to be split as
$\bm{n}=\Omega[\bm{X}]+\Omega^{1/2}\delta\bm{n}$ in terms of the deterministic
solution $[\bm{X}]$ of the Chemical Rate Equation (12) and Gaussian
fluctuations $\delta\bm{n}$. Such a decomposition of $\bm{n}$ allows us to
approximate the mean scaled chemical reaction rates as

$\displaystyle\frac{\langle w_{\pm\rho}(\bm{n})\rangle}{\Omega}$
$\displaystyle=\frac{\langle
w_{\pm\rho}(\Omega[\bm{X}]+\Omega^{1/2}\delta\bm{n})\rangle}{\Omega}$ (80)
$\displaystyle\approx\frac{w_{\pm\rho}(\Omega[\bm{X}])}{\Omega}+\left.\frac{\langle\delta\bm{n}\rangle\cdot\nabla_{\bm{n}}w_{\pm\rho}(\bm{n})}{\Omega^{1/2}}\right|_{\bm{n}=\Omega[\bm{X}]}$
$\displaystyle\approx
k_{\pm\rho}[\bm{X}]^{\bm{\nu}_{\rho}^{\pm}}+\Omega^{-1/2}\langle\delta\bm{n}\rangle\cdot\bm{\nabla}_{[\bm{X}]}k_{\pm\rho}[\bm{X}]^{\bm{\nu}^{\pm}_{\rho}}.$

As a result, we can approximate the time-evolution of the mean and covariance
of $\delta\bm{n}$ as

$\displaystyle\frac{d\langle\delta\bm{n}_{t}\rangle}{dt}$
$\displaystyle=\Omega^{1/2}\left(\frac{1}{\Omega}\frac{d\langle\bm{n}_{t}\rangle}{dt}-\frac{d[\bm{X}]_{t}}{dt}\right)$
(81)
$\displaystyle=\Omega^{1/2}\sum_{\rho\in\mathcal{R}}\Delta\bm{\nu}_{\rho}\left(\frac{\langle
w_{+\rho}(\bm{n}_{t})\rangle-\langle
w_{-\rho}(\bm{n}_{t})\rangle}{\Omega}\right)$
$\displaystyle~{}~{}~{}-\Omega^{1/2}\sum_{\rho\in\mathcal{R}}\Delta\bm{\nu}_{\rho}\left(k_{+\rho}[\bm{X}]^{\bm{\nu}_{\rho}^{+}}-k_{-\rho}[\bm{X}]^{\bm{\nu}_{\rho}^{-}}\right)$
$\displaystyle\approx\langle\delta\bm{n}_{t}\rangle\cdot\bm{\nabla}_{[\bm{X}_{t}]}\sum_{\rho\in\mathcal{R}}\Delta\bm{\nu}_{\rho}\left(k_{+\rho}[\bm{X}]_{t}^{\bm{\nu}^{+}_{\rho}}-k_{-\rho}[\bm{X}]_{t}^{\bm{\nu}^{-}_{\rho}}\right)$
$\displaystyle\equiv\mathsf{L}\cdot\langle\delta\bm{n}_{t}\rangle$

and

$\displaystyle\frac{d{\rm Cov}_{t}\\{\delta n_{i},\delta
n_{j}\\}}{dt}=\frac{1}{\Omega}\frac{d{\rm Cov}_{t}\\{n_{i},n_{j}\\}}{dt}$ (82)
$\displaystyle=\frac{1}{\Omega}\sum_{\rho\in\mathcal{R}}S_{i\rho}S_{j\rho}\left\\{w_{+\rho}(\bm{n}_{t})+w_{-\rho}(\bm{n}_{t})\right\\}$
$\displaystyle\approx\sum_{\rho\in\mathcal{R}}S_{i\rho}S_{j\rho}\left(k_{+\rho}[\bm{X}]_{t}^{\bm{\nu}^{+}_{\rho}}+k_{-\rho}[\bm{X}]_{t}^{\bm{\nu}^{-}_{\rho}}\right)\equiv
N_{ij}.$

The last equalities define the relaxation matrix $\mathsf{L}$ and noise matrix
$\mathsf{N}$. Together these equations imply that the Gaussian dynamics of
$\delta\bm{n}$ can be captured with the linear Langevin equation

$\frac{d(\delta\bm{n}_{t})}{dt}=\mathsf{L}\cdot\delta\bm{n}_{t}+\bm{\zeta}_{t},$
(83)

where $\bm{\zeta}_{t}$ is a multivariate zero-mean Gaussian white noise with
covariance $\mathsf{N}$. As a result, the steady-state covariance $V_{ij}={\rm
Cov}\\{n_{i},n_{j}\\}$ of a linear Langevin equation can be obtained as the
solution of the algebraic equation [71]

$\mathsf{L}\cdot\mathsf{V}+\mathsf{V}\cdot\mathsf{L}^{\rm T}+\mathsf{N}=0,$
(84)

which we solve numerically. The scaled variance is then obtained as
$D_{i}={\rm Var}\\{n_{i}\\}/\Omega={\rm Var}\\{\delta n_{i}\\}$.

##  Appendix B Trade-offs for General Deficiency Zero CRNs

The structure of the equation $\mathsf{A}\cdot\bm{\Psi}([\bm{X}]_{\rm
ss})=\bm{0}$ that determines $[\bm{X}]_{\rm ss}$ allows us to find
restrictions on the response functions for the deterministic concentrations.
The vector $\bm{\Psi}([\bm{X}]_{\rm ss})$ lies in the right null space of the
Laplacian matrix $\mathsf{A}$. One useful basis of the right null space of the
Laplacian matrix is provided by the matrix-tree theorem (MTT) [66, 67, 30].

To exploit this connection, a few graph-theoretical notions are needed. A
complex-reaction graph may consist of multiple linkage classes, which we
denote as $g_{p}$ with $p\in\\{1,2,\cdots,\ell\\}$. The MTT is then built out
of rooted spanning trees of each $g_{p}$, which are connected subgraphs that
contain no cycles, include all vertices, and have each edge oriented towards a
specific vertex, $\bm{y}_{l}$. A weight is then assigned to each tree by
noting that each directed edge is endowed with a rate constant from the
corresponding chemical reaction: the weight of the subgraph is then product of
the associated rate constants. We denote the sum of the weights of every
spanning tree rooted at vertex $\bm{y}_{l}$ by $T_{l}$. The MTT then states
that the vectors $\bm{T}^{(p)}$ with components

$T^{(p)}_{l}=\begin{cases}T_{l}&{\rm if}~{}\bm{y}_{l}\in g_{p},\\\ 0&{\rm
otherwise},\end{cases}$ (85)

lie in the right null space of the Laplacian matrix, i.e.,
$\mathsf{A}\cdot\bm{T}^{(p)}=\bm{0}$ for all $p\in\\{1,2,\cdots,\ell\\}$ [66,
67, 30]. Since the $\bm{T}^{(p)}$ with different $p$ are linearly independent
of each other, the collection of $\bm{T}^{(p)}$ for
$p\in\\{1,2,\cdots,\ell\\}$ can serve as a basis of the right null space of
$\mathsf{A}$. It follows that $\Psi_{l}([\bm{X}]_{\rm
ss})/\Psi_{m}([\bm{X}]_{\rm ss})=T_{l}/T_{m}$ if $\bm{y}_{l}$ and $\bm{y}_{m}$
are in the same linkage class. The ratios for two vertices belonging to
different linkage classes are determined by constraints imposed by
conservation laws.

Now, due to the product form of the steady-state distribution (40), the static
response of the log-ratio of steady-state probabilities to a perturbation is
given by

$\frac{\partial}{\partial\lambda}\ln\frac{\pi_{\Gamma}(\bm{n})}{\pi_{\Gamma}(\bm{n}^{\prime})}=\sum_{i\in\mathcal{S}}(n_{i}-n_{i}^{\prime})\frac{\partial}{\partial\lambda}\ln[X_{i}]_{\rm
ss}.$ (86)

In order to take advantage of the graph-theoretic tools, (86) needs to be
rearranged to a form involving the $\Psi_{i}([\bm{X}]_{\rm ss})$, which can be
replaced with spanning trees. To this end, we change the basis from number,
and instead identify each microscopic state by the number of (independent)
reactions required to arrive at it. We accomplish this by noting the matrix
$\mathsf{S}^{\prime}$, defined in (53), converts the difference of state
vectors as
$\bm{n}-\bm{n}^{\prime}=\mathsf{S}^{\prime}\cdot\bm{\mu}(\bm{n},\bm{n}^{\prime})$
into $\bm{\mu}(\bm{n},\bm{n}^{\prime})$ vector whose components are the number
of each reaction required to move from ${\bm{n}}^{\prime}$ to ${\bm{n}}$. Note
that the decomposition of independent and dependent reactions used to
construct $\mathsf{S}^{\prime}$ is not unique. Substituting this change of
basis into (86), coupled with the identity
$\sum_{j\in\mathcal{S}}S_{j\rho}^{\prime}\ln[X_{j}]=\ln[\bm{X}]_{\rm
ss}^{\Delta\bm{\nu}_{\rho}}=\ln([\bm{X}]_{\rm
ss}^{\bm{\nu}_{\rho}^{-}}/[\bm{X}]_{\rm ss}^{\bm{\nu}_{\rho}^{+}})$, leads to

$\displaystyle\frac{\partial}{\partial\lambda}\ln\frac{\pi_{\Gamma}(\bm{n})}{\pi_{\Gamma}(\bm{n}^{\prime})}$
$\displaystyle=\sum_{\rho\in\mathcal{R}^{\prime}}\mu_{\rho}(\bm{n},\bm{n}^{\prime})\frac{\partial}{\partial\lambda}\ln\frac{[\bm{X}]_{\rm
ss}^{\bm{\nu}_{\rho}^{-}}}{[\bm{X}]_{\rm ss}^{\bm{\nu}_{\rho}^{+}}}$ (87)
$\displaystyle=\sum_{\rho\in\mathcal{R}^{\prime}}\mu_{\rho}(\bm{n},\bm{n}^{\prime})\frac{\partial}{\partial\lambda}\ln\frac{\Psi_{{\bm{\nu}}_{\rho}^{-}}([{\bm{X}}_{\rm
ss}])}{\Psi_{{\bm{\nu}}_{\rho}^{+}}([{\bm{X}}_{\rm ss}])}$
$\displaystyle=\sum_{\rho\in\mathcal{R}^{\prime}}\mu_{\rho}(\bm{n},\bm{n}^{\prime})\frac{\partial}{\partial\lambda}\ln\frac{T_{{\bm{\nu}}_{\rho}^{-}}}{T_{{\bm{\nu}}_{\rho}^{+}}}$

where sum runs over independent reactions,
$\mathcal{R}^{\prime}=\\{1,2,\cdots,s\\}$. Note the complexes
$\bm{\nu}_{\rho}^{-}$ and $\bm{\nu}_{\rho}^{+}$ are always in the same linkage
class since they are connected via the $\rho$-th chemical reaction.

As a result, the analysis reduces to bounding the sensitivity of ratios of
spanning tree weights. The bounds to the specific perturbations
$\partial/\partial\lambda=k_{\pm\rho}\partial/\partial k_{\pm\rho}$ and
$\partial/\partial\lambda=\partial/\partial B_{\rho}$, now follow directly
from the methods developed in [30], with the results

$\displaystyle\left|k_{\pm\rho}\frac{\partial}{\partial
k_{\pm\rho}}\ln\frac{T_{\bm{\nu}_{\sigma^{-}}}}{T_{\bm{\nu}_{\sigma^{+}}}}\right|\leq
1,$ (88) $\displaystyle\left|\frac{\partial}{\partial
B_{\rho}}\ln\frac{T_{\bm{\nu}_{\sigma^{-}}}}{T_{\bm{\nu}_{\sigma^{+}}}}\right|\leq\tanh\left(\frac{\mathcal{F}}{4}\right).$
(89)

Applying the inequalities (88) and (89) to (87), we obtain the bounds,

$\displaystyle\left|k_{\pm\rho}\frac{\partial}{\partial
k_{\pm\rho}}\ln\frac{\pi_{\Gamma}(\bm{n})}{\pi_{\Gamma}(\bm{n}^{\prime})}\right|\leq\sum_{\rho\in\mathcal{R}^{\prime}}\left|\mu_{\rho}(\bm{n},\bm{n}^{\prime})\right|,$
(90) $\displaystyle\left|\frac{\partial}{\partial
B_{\rho}}\ln\frac{\pi_{\Gamma}(\bm{n})}{\pi_{\Gamma}(\bm{n}^{\prime})}\right|\leq\sum_{\rho\in\mathcal{R}^{\prime}}\left|\mu_{\rho}(\bm{n},\bm{n}^{\prime})\right|\tanh\left(\frac{\mathcal{F}}{4}\right).$
(91)

The bound (90) is reminiscent of (25) since
$\mu_{\rho}(\bm{n},\bm{n}^{\prime})$ is a measure of the distance between two
states $\bm{n}$ and $\bm{n}^{\prime}$. If the dynamics are a birth-death
process, then the chemical reaction changes the number of chemical species by
one, and the bound (90) coincides with (25). Nevertheless, the derivations are
independent. The class of CRNs considered in Sec. III may have a nonzero
deficiency (e.g. the Schlögl model has a deficiency of one), whereas (90)
applies only to deficiency zero CRNs.

Lastly, we derive trade-offs between the response of the mean number of a
chemical species and fluctuations. The response of the conditional mean number
$\langle
n_{i}\rangle_{\Gamma}=\sum_{\bm{n}\in\Gamma}n_{i}\pi_{\Gamma}(\bm{n})$ on an
isolated connected component $\Gamma$ is given by multiplying
$n_{i}\pi_{\Gamma}(\bm{n})\pi_{\Gamma}(\bm{n}^{\prime})$ by both sides of (87)
and summing over $\bm{n},\bm{n}^{\prime}\in\Gamma$. Using the identity
$\sum_{\bm{n}^{\prime}}\pi_{\Gamma}(\bm{n}^{\prime})\partial\ln\pi_{\Gamma}(\bm{n}^{\prime})/\partial\lambda=0$,
we have

$\frac{\partial\langle
n_{i}\rangle_{\Gamma}}{\partial\lambda}=\sum_{\rho\in\mathcal{R}^{\prime}}\langle
n_{i}\mu_{\rho}(\bm{n},\bm{n}^{\prime})\rangle_{\Gamma}\frac{\partial}{\partial\lambda}\ln\frac{T_{{\bm{\nu}}_{\rho}^{-}}}{T_{{\bm{\nu}}_{\rho}^{+}}}$
(92)

with

$\langle
n_{i}\mu_{\rho}(\bm{n},\bm{n}^{\prime})\rangle_{\Gamma}=\sum_{\bm{n},\bm{n}^{\prime}\in\Gamma}n_{i}\mu_{\rho}(\bm{n},\bm{n}^{\prime})\pi_{\Gamma}(\bm{n})\pi_{\Gamma}(\bm{n}^{\prime}).$
(93)

The correlation function $\langle
n_{i}\mu_{\rho}(\bm{n},\bm{n}^{\prime})\rangle_{\Gamma}$ measures the
fluctuations in the numbers of chemical species through the definition of
$\mu_{\rho}(\bm{n},\bm{n}^{\prime})$. Since all the columns are linearly
independent, the matrix $\mathsf{S}^{\prime}$ has a left inverse matrix
$(\mathsf{S}^{\prime})^{-1}$ satisfying
$(\mathsf{S}^{\prime})^{-1}\mathsf{S}^{\prime}=\mathsf{I}$. The correlation
function $\langle n_{i}\mu_{\rho}(\bm{n},\bm{n}^{\prime})\rangle_{\Gamma}$ can
be written as

$\langle
n_{i}\mu_{\rho}(\bm{n},\bm{n}^{\prime})\rangle_{\Gamma}=\sum_{j\in\mathcal{S}}(\mathsf{S}^{\prime})^{-1}_{\rho
j}{\rm Cov}_{\Gamma}\\{n_{i},n_{j}\\}.$ (94)

Note that the fluctuation-response relation in (41) can be restored by
noticing

$\sum_{\rho\in\mathcal{R}^{\prime}}(\mathsf{S})_{\rho
j}^{-1}\ln\frac{[\bm{X}]^{\bm{\nu}_{\rho}^{-}}_{\rm
ss}}{[\bm{X}]^{\bm{\nu}_{\rho}^{+}}_{\rm ss}}=\ln[X_{j}]_{\rm ss}.$ (95)

Applying the inequalities (88) and (89) to (92) together with (94), we obtain
trade-offs between the response and fluctuations,

$\left|k_{\pm\rho}\frac{\partial\langle n_{i}\rangle_{\Gamma}}{\partial
k_{\pm\rho}}\right|\leq\sum_{\sigma\in\mathcal{R}^{\prime}}\left|\sum_{j\in\mathcal{S}}(\mathsf{S}^{\prime})^{-1}_{\sigma
j}{\rm Cov}_{\Gamma}\\{n_{i},n_{j}\\}\right|$ (96)

and

$\left|\frac{\partial\langle n_{i}\rangle_{\Gamma}}{\partial
B_{\rho}}\right|\leq\sum_{\sigma\in\mathcal{R}^{\prime}}\left|\sum_{j\in\mathcal{S}}(\mathsf{S}^{\prime})^{-1}_{\sigma
j}{\rm
Cov}_{\Gamma}\\{n_{i},n_{j}\\}\right|\tanh\left(\frac{\mathcal{F}}{4}\right).$
(97)

## References

  * Goldbeter and Koshland [1981] A. Goldbeter and D. E. Koshland, An amplified sensitivity arising from covalent modification in biological systems, Proc. Natl. Acad. Sci. USA 78, 6840 (1981).
  * Qian [2003] H. Qian, Thermodynamic and kinetic analysis of sensitivity amplification in biological signal transduction, Biophys. Chem. 105, 585 (2003).
  * Ge and Qian [2008] H. Ge and M. Qian, Sensitivity amplification in the phosphorylation-dephosphorylation cycle: Nonequilibrium steady states, chemical master equation, and temporal cooperativity, J. Chem. Phys. 129, 015104 (2008).
  * Mallory  _et al._ [2019] J. D. Mallory, A. B. Kolomeisky, and O. A. Igoshin, Trade-offs between error, speed, noise, and energy dissipation in biological processes with proofreading, J. Phys. Chem. B 123, 4718 (2019).
  * [5] J. A. Owen, P. Talla, J. W. Biddle, and J. Gunawardena, Thermodynamic bounds on ultrasenitivity in covalent switching, arXiv:2211.10472.
  * Estrada  _et al._ [2016] J. Estrada, F. Wong, A. DePace, and J. Gunawardena, Information integration and energy expenditure in gene regulation, Cell 166, 234 (2016).
  * Tran  _et al._ [2018] H. Tran, J. Desponds, C. A. P. Romero, M. Coppey, C. Fradin, N. Dostatni, and A. M. Walczak, Precision in a rush: Trade-offs between reproducibility and steepnees of the hunchmback expression pattern, PLoS Comput. Biol. 14, e1006513 (2018).
  * Falasco  _et al._ [2019] G. Falasco, T. Cossetto, E. Penocchio, and M. Esposito, Negative differential response in chemical reactions, New J. Phys. 21, 073005 (2019).
  * Nguyen  _et al._ [2018] B. Nguyen, U. Seifert, and A. C. Barato, Phase transition in thermodynamically consistent biochemical oscillators, J. Chem. Phys. 149, 045101 (2018).
  * Callen and Welton [1951] H. B. Callen and T. A. Welton, Irreversibility and generalized noise, Phys. Rev. 83, 34 (1951).
  * Kubo [1966] R. Kubo, The fluctuation-dissipation theorem, Rep. Prog. Phys. 29, 255 (1966).
  * Forster [1975] D. Forster,  _Hydrodynamic fluctuations, broken symmetry, and correlation functions_ (W. A. Benjamin Inc., Massachusetts, 1975).
  * Marconi  _et al._ [2008] U. M. B. Marconi, A. Puglisi, L. Rondoni, and A. Vulpiani, Fluctuation–dissipation: response theory in statistical physics, Phys. Rep. 461, 111 (2008).
  * Agarwal [1972] G. S. Agarwal, Fluctuation-dissipation theorems for systems in non-thermal equilibrium and applications, Z. Phys. A 252, 25 (1972).
  * Harada and Sasa [2005] T. Harada and S. Sasa, Equality connecting energy dissipation with a violation of the fluctuation-response relation, Physical Review Letters 95, 130602 (2005).
  * Speck and Seifert [2006] T. Speck and U. Seifert, Restoring a fluctuation-dissipation theorem in a nonequilibrium steady state, Europhys. Lett. 74, 391 (2006).
  * Baiesi  _et al._ [2009] M. Baiesi, C. Maes, and B. Wynants, Fluctuations and response of nonequilibrium states, Phys. Rev. Lett. 103, 010602 (2009).
  * Prost  _et al._ [2009] J. Prost, J.-F. Joanny, and J. M. Parrondo, Generalized fluctuation-dissipation theorem for steady-state systems, Phys. Rev. Lett. 103, 090601 (2009).
  * Chetrite and Gawedzki [2009] R. Chetrite and K. Gawedzki, Eulerian and lagrangian pictures of non-equilibrium diffusions, J. Stat. Phys. 137, 890 (2009).
  * Seifert and Speck [2010] U. Seifert and T. Speck, Fluctuation-dissipation theorem in nonequilibrium steady states, Europhys. Lett. 89, 10007 (2010).
  * Baiesi and Maes [2013] M. Baiesi and C. Maes, An update on the nonequilibrium linear response, New J. Phys. 15, 013004 (2013).
  * Wang  _et al._ [2016] S.-W. Wang, K. Kawaguchi, S. I. Sasa, and L. H. Tang, Entropy produciton of nanosystems with time scale seperation, Phys. Rev. Lett. 117, 070601 (2016).
  * Lubensky [2010a] D. K. Lubensky, Equilibriumlike behavior in chemical reaction networks far from equilibrium, Phys. Rev. E 81, 060102(R) (2010a).
  * Caprini [2021] L. Caprini, Generalized fluctuation-dissipation relations holding in non-equilibrium dynamics, J. Stat. Mech. , 063202 (2021).
  * Graham [1977] R. Graham, Covariant formulation of non-equilibrium statistical thermodynamics, Z. Phys. B Condens. matter 26, 397 (1977).
  * Asheichyk  _et al._ [2019] K. Asheichyk, A. P. Solon, C. M. Rohwer, and M. Krüger, Response of active brownian particles to shear flow, J. Chem. Phys. 150, 144111 (2019).
  * Dal Cengio  _et al._ [2019] S. Dal Cengio, D. Levis, and I. Pagonabarraga, Linear response theory and green-kubo relations for active matter, Phys. Rev. Lett. 123, 238003 (2019).
  * Chun  _et al._ [2021] H.-M. Chun, Q. Gao, and J. M. Horowitz, Nonequilibrium green-kubo relations for hydrodynamic transport from an equilibrium-like fluctuation-response equality, Phys. Rev. Res. 3, 043172 (2021).
  * Dechant and Sasa [2020] A. Dechant and S.-i. Sasa, Fluctuation–response inequality out of equilibrium, Proc. Natl. Acad. Sci. U.S.A. 117, 6430 (2020).
  * Owen  _et al._ [2020] J. A. Owen, T. R. Gingrich, and J. M. Horowitz, Universal thermodynamic bounds on nonequilibrium response with biochemical applications, Phys. Rev. X 10, 011066 (2020).
  * Gao  _et al._ [2022] Q. Gao, H.-M. Chun, and J. M. Horowitz, Thermodynamic constraints on the nonequilibrium response of one-dimensional diffusions, Phys. Rev. E 105, L012102 (2022).
  * Gunawardena [2003] J. Gunawardena, Chemical reaction network theory for in-silico biologists, Notes available for download at http://vcp.med.harvard.edu/papers/crnt.pdf (2003).
  * Feinberg [2019] M. Feinberg,  _Foundations of chemical reaction network theory_ (Springer, 2019).
  * Nicolis and Prigogine [1977] G. Nicolis and I. Prigogine,  _Self-Organization in Nonequilibrium Systems: From Dissipative Structures to Order through Fluctuations_ (Wiley, 1977).
  * De Groot and Mazur [2011] S. R. De Groot and P. Mazur,  _Non-Equilibrium Thermodynamics_ (Dover Publications, 2011).
  * Polettini and Esposito [2014] M. Polettini and M. Esposito, Irreversible thermodynamics of open chemical networks. i. emergent cycles and broken conservation laws, J. Chem. Phys. 141, 024117 (2014).
  * Ge and Qian [2016] H. Ge and H. Qian, Mesoscopic kinetic basis of macroscopic chemical thermodynamics: A mathematical theory, Phys. Rev. E 94, 052150 (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).
  * Avanzini  _et al._ [2020] F. Avanzini, G. Falasco, and M. Esposito, Thermodynamics of non-elementary chemical reaction networks, New J. Phys. 22, 093040 (2020).
  * Avanzini  _et al._ [2021] F. Avanzini, E. Penocchio, G. Falasco, and M. Esposito, Nonequilibrium thermodynamics of non-ideal chemical reaction networks, J. Chem. Phys. 154, 094114 (2021).
  * Schmiedl and Seifert [2007] T. Schmiedl and U. Seifert, Stochastic thermodynamics of chemical reaction networks, J. Chem. Phys. 126, 044101 (2007).
  * Rao and Esposito [2018] R. Rao and M. Esposito, Conservation laws and work fluctuation relations in chemical reaction networks, J. Chem. Phys. 149, 245101 (2018).
  * Gaspard [2004] P. Gaspard, Fluctuation theorem for nonequilibrium reactions, J. Chem. Phys. 120, 8898 (2004).
  * Andrieux and Gaspard [2004] D. Andrieux and P. Gaspard, Fluctuation theorem and onsager reciprocity relations, J. Chem. Phys. 121, 6167 (2004).
  * Andrieux and Gaspard [2007] D. Andrieux and P. Gaspard, Fluctuation theorem for currents and schnakenberg network theory, J. Stat. Phys. 127, 107 (2007).
  * Yoshimura and Ito [2021a] K. Yoshimura and S. Ito, Information geometric inequalities of chemical thermodynamics, Phys. Rev. Res. 3, 013175 (2021a).
  * Yoshimura and Ito [2021b] K. Yoshimura and S. Ito, Thermodynamic uncertainty relation and thermodynamic speed limit in deterministic chemical reaction networks, Phys. Rev. Lett. 127, 160601 (2021b).
  * Craciun  _et al._ [2006] G. Craciun, T. Tang, and M. Feinberg, Understanding bistability in complex enzyme-driven reaction networks, Proc. Natl. Acad. Sci. U.S.A. 103, 8697 (2006).
  * Tyson  _et al._ [2008] J. J. Tyson, R. Albert, A. Goldbeter, P. Ruoff, and S. J., Biological switches and clocks, J. R. Soc. Interface 5, S1 (2008).
  * Polettini  _et al._ [2015] M. Polettini, A. Wachtel, and M. Esposito, Dissipation in noisy chemical networks: The role of deficiency, J. Chem. Phys. 143, 184103 (2015).
  * Cornish-Bowden [2013] A. Cornish-Bowden,  _Fundamentals of enzyme kinetics_ (John Wiley & Sons, 2013).
  * Wachtel  _et al._ [2018] A. Wachtel, R. Rao, and M. Esposito, Thermodynamically consistent coarse graining of biocatalysts beyond michaelis–menten, New J. Phys. 20, 042002 (2018).
  * Gillespie [1992] D. T. Gillespie, A rigorous derivation of the chemical master equation, Phys. A: Stat. Mech. Appl. 188, 404 (1992).
  * Van Kampen [2007] N. G. Van Kampen,  _Stochastic processes in physics and chemistry_ , 3rd ed. (North Holland, 2007).
  * Kurtz [1972] T. G. Kurtz, The relationship between stochastic and deterministic models for chemical reactions, J. Chem. Phys. 57, 2976 (1972).
  * Gillespie [2009] D. T. Gillespie, Deterministic limit of stochastic chemical kinetics, The Journal of Physical Chemistry B 113, 1640 (2009).
  * Vellela and Qian [2009] M. Vellela and H. Qian, Stochastic dynamics and non-equilibrium thermodynamics of a bistable chemical system: the schlögl model revisited, J. R. Soc. Interface 6, 925 (2009).
  * Ge and Qian [2009] H. Ge and H. Qian, Thermodynamic limit of a nonequilibrium steady state: Maxwell-type construction for a bistable biochemical system, Phys. Rev. Lett. 103, 148103 (2009).
  * Gillespie [2000] D. T. Gillespie, The chemical langevin equation, J. Chem. Phys. 113, 297 (2000).
  * Schnakenberg [1976] J. Schnakenberg, Network theory of microscopic and macroscopic behavior of master equation systems, Reviews of Modern physics 48, 571 (1976).
  * Schlögl [1972] F. Schlögl, Chemical reaction models for non-equilibrium phase transitions, Zeitschrift für physik 253, 147 (1972).
  * Freitas and Esposito [2022] N. Freitas and M. Esposito, Emergent second law for non-equilibrium steady states, Nat. Commun. 13, 5084 (2022).
  * Nicolis and Turner [1977] G. Nicolis and J. Turner, Stochastic analysis of a nonequilibrium phase transition: Some exact results, Phys. A: Stat. Mech. Appl. 89, 326 (1977).
  * Anderson  _et al._ [2010] D. F. Anderson, G. Craciun, and T. G. Kurtz, Product-form stationary distributions for deficiency zero chemical reaction networks, Bull. Math. Biol. 72, 1947 (2010).
  * Lubensky [2010b] D. K. Lubensky, Equilibriumlike behavior in chemical reaction networks far from equilibrium, Phys. Rev. E 81, 060102 (2010b).
  * Hill [1966] T. L. Hill, Studies in irreversible thermodynamics iv. diagrammatic representation of steady state fluxes for unimolecular systems, J. Theor. Biol. 10, 442 (1966).
  * Gunawardena [2012] J. Gunawardena, A linear framework for time-scale separation in nonlinear biochemical systems, PloS one 7, e36321 (2012).
  * Heuett and Qian [2006] W. J. Heuett and H. Qian, Grand canonical markov model: a stochastic theory for open nonequilibrium biochemical networks, J. Chem. Phys. 124, 044110 (2006).
  * Kirkpatrick  _et al._ [1983] S. Kirkpatrick, C. D. Gelatt Jr, and M. P. Vecchi, Optimization by simulated annealing, science 220, 671 (1983).
  * Ferrell Jr and Ha [2014] J. E. Ferrell Jr and S. H. Ha, Ultrasensitivity part i: Michaelian responses and zero-order ultrasensitivity, Trends Biochem. Sci. 39, 496 (2014).
  * Gardiner [2009] C. Gardiner,  _Stochastic Methods: A Handbook for the Natural and Social Sciences_ , 4th ed. (Springer, New York, 2009).