Robustness in power law kinetic systems with reactant-determined interactions

TL;DR

This paper extends the concept of absolute concentration robustness (ACR) from mass-action systems to more general power-law kinetic systems with reactant-determined interactions, providing new insights into system robustness.

Contribution

It generalizes ACR conditions to power-law kinetic systems with reactant-determined interactions, broadening the applicability of robustness analysis beyond mass-action models.

Findings

Power-law kinetic systems with reactant-determined interactions can exhibit ACR.

A power-law approximation of the global carbon cycle model shows ACR properties.

The extended framework helps identify robustness components in complex chemical networks.

Abstract

Robustness against the presence of environmental disruptions can be observed in many systems of chemical reaction network. However, identifying the underlying components of a system that give rise to robustness is often elusive. The influential work of Shinar and Feinberg established simple yet subtle network-based conditions for absolute concentration robustness (ACR), a phenomena in which a species in a mass-action system has the same concentration for any steady state the network may admit. In this contribution, we extend this result to embrace kinetic systems more general than mass-action systems, namely, power-law kinetic systems with reactant-determined interactions (denoted by "PL-RDK"). In PL-RDK, the kinetic order vectors (which we call "interactions") of reactions with the same reactant complex are identical. As illustration, we considered a scenario in the pre-industrial state of global carbon cycle. A power-law approximation of the dynamical system of this scenario is found to be dynamically equivalent to an ACR-possessing PL-RDK system.

Tables

Table 1. Table A.1: Power law approximation of the process rates.

Carbon Flux	Function	Power law approx.
$A_2\to A_1$	$K_1=r_{tc}P\left(t\right)A_1\left[1-\frac{A_1}{k}\right]$	$k_1{A}_1^{p_1}{A}_2^{q_1}$
$A_1\to A_2$	$K_2=r_{tc}R\left(t\right)A_1\left[1-\frac{A_1}{k}\right]+\alpha A_1$	$k_2{A}_1^{p_2}{A}_2^{q_2}$
$A_2\to A_3$	$K_3=a_m{A}_2$	$a_m{A}_2$
$A_3\to A_2$	$K_4=a_m\beta A_3$	$a_m\beta A_3$

Equations

$$S:=\text{span }\{y^{\prime }-y\in {\mathbb{R}}^{\mathcal{S}}|y\to y^{\prime }\in \mathcal{R}\}.$$

$$A_{k}x:=\sum _{y\to y^{\prime }\in \mathcal{R}}{k}_{y\to y^{\prime }}{x}_y({\omega }_{y^{\prime }}-{\omega }_y),$$

$$\dim (\text{Ker }Y{A}_k)\le \delta +t.$$

$$K=[K_1,K_2,...,K_r{]}^{\top }:{\Omega }_K\to {\mathbb{R}}_{\ge 0}^{\mathcal{R}}.$$

$$\text{For each reaction }{R}_j:y_j\to y_j^{\prime }\in \mathcal{R},K_j(c)>0\text{ if and only if }\text{supp }{y}_j\subset \text{supp }c.$$

$$f(c)=NK(c)=\sum _{y_j\to y_j^{\prime }\in \mathcal{R}}{K}_j(c)(y_j^{\prime }-y_j).$$

$$K_i(x)=k_i{x}^{F_{i,\cdot }}\forall i=1,\dots ,r$$

$$K_j(c)=k_j{x}^{Y_{.,j}}\text{for all }{R}_j:y_j\to y_j^{\prime }\in \mathcal{R}$$

$$\frac{dx}{dt}=Y\circ A_k\circ {\Psi }_k,$$

$$X_1$$

$$\left.\begin{array}[]{rl}\dot{X}_{1}&=k_{1}X_{2}^{0.8}-k_{2}X_{1}^{0.5}X_{2}^{0.8}\\ \dot{X}_{2}&=-k_{1}X_{2}^{0.8}+k_{2}X_{1}^{0.5}X_{2}^{0.8}\end{array}\right.$$

$$X_1={\left(\frac{k_1}{k_2}\right)}^2\text{and}{X}_2=\Gamma -{\left(\frac{k_1}{k_2}\right)}^2.$$

$$\left(T_{\cdot ,y}-T_{\cdot ,y^{\prime }}\right)\cdot \log \left(\frac{c^{\ast \ast }}{c^{\ast }}\right)=0.$$

$$\widetilde{Y}=\left[\begin{array}[]{c|c}\text{ }T\text{ }&\text{ }0\text{ }\\ \end{array}\right],$$

$$\sum _{y\to y^{\prime }\in \mathcal{R}}{k}_{y\to y^{\prime }}(c^{\ast }{)}^{\tilde{y}}(y^{\prime }-y)=0.$$

$${\kappa }_{y\to y^{\prime }}:=k_{y\to y^{\prime }}(c^{\ast }{)}^{\tilde{y}}.$$

$$\sum _{y\to y^{\prime }\in \mathcal{R}}{\kappa }_{y\to y^{\prime }}(y^{\prime }-y)=0.$$

$$\sum _{y\to y^{\prime }\in \mathcal{R}}{k}_{y\to y^{\prime }}(c^{\ast \ast }{)}^{\tilde{y}}(y^{\prime }-y)=0.$$

$$\mu :=\log c^{\ast \ast }-\log c^{\ast }.$$

$$\sum _{y\to y^{\prime }\in \mathcal{R}}{\kappa }_{y\to y^{\prime }}{e}^{\tilde{y}\cdot \mu }(y^{\prime }-y)=0.$$

$$1^{\mathcal{C}}=\sum _{y\in \mathcal{C}}{\omega }_y.$$

$$Y{A}_{\kappa }{1}^{\mathcal{C}}=0,\text{ and }Y{A}_{\kappa }\left(\sum _{y\in \mathcal{C}}{e}^{\tilde{y}\cdot \mu }{\omega }_y\right)=0.$$

$$1^{\mathcal{C}}\in \text{Ker }Y{A}_{\kappa },\text{ and}$$

$$\sum _{y\in \mathcal{C}}{e}^{\tilde{y}\cdot \mu }{\omega }_y\in \text{Ker }Y{A}_{\kappa }.$$

$$\dim (\text{Ker }Y{A}_{\kappa })\le 1+t$$

$$1^{\mathcal{C}}\notin \text{Span }\{b^1,b^2,\dots ,b^t\}.$$

$$\sum _{y\in \mathcal{C}}{e}^{\tilde{y}\cdot \mu }{\omega }_y={\lambda }_0{1}^{\mathcal{C}}+\sum _{i=1}^t{\lambda }_i{b}^i.$$

$$\tilde{y}\cdot \mu =\widetilde{y^{\prime }}\cdot \mu .$$

$$T_{\cdot ,y}\cdot \mu =T_{\cdot ,y^{\prime }}\cdot \mu ,$$

$$T_{\cdot ,y}-T_{\cdot ,y^{\prime }}=a{X}_i$$

Full text

Keywords. Absolute concentration robustness , Chemical reaction network, Power law kinetics, Reactant-determined interactions, Carbon cycle model

1 Introduction

Robustness may be generally defined [15, 26] as a system-level dynamical property that allows a system to sustain its functions despite changes in internal and external conditions. This feature, in fact, is fundamental and ubiquitous in many biological processes, including cellular networks and entire organisms [2, 15, 26]. One type of robust behavior is “concentration robustness,” wherein some quantity involving the concentrations of the different species in a network is fixed at equilibrium [10]. In a well-cited paper published in Science, Shinar and Feinberg [26] introduced absolute concentration robustness (ACR), a condition in which the concentration of a species in a network attains the same value in every positive steady state set by parameters and does not depend on initial conditions.

Shinar and Feinberg presented sufficient structure-based conditions for a chemical reaction network (CRN) to display ACR on a particular species through a structural index called the deficiency. This non-negative parameter has been the center of many powerful results in Chemical Reaction Network Theory (CRNT), a theoretical body of work that associates the structure of a CRN to the dynamical behaviour of the system [11, 12]. CRNT employs mathematical methods from graph theory, linear algebra, group theory and the theory of ordinary differential equations. In CRNT, chemical reaction networks are viewed as digraphs whose vertices (called complexes) are mapped to non-negative vectors representing compositions of chemical species and whose arcs represent chemical reactions between them. The Shinar-Feinberg Theorem on ACR holds for systems whose evolution are modelled by ordinary differential equations with mass-action kinetics (MAK), and is stated as follows:

Consider a mass-action system that admits a positive steady state and suppose that the deficiency of the underlying reaction network is one. If there are two nonterminal nodes in the network that differ only in species $S$, then the system has absolute concentration robustness in $S$.

Here, we show that this result extends to systems endowed with power law kinetics (PLK), which generalize mass-action kinetics [7, 14]. Several experiments have shown that the kinetic order of a reaction with respect to a given reactant is a function of the geometry within which the reaction occurs [16, 17, 18, 22, 25]. In the case of reactions occurring within a three-dimensional homogenous space (as in mass-action systems), the kinetic order is the same as the number of molecules entering into the reaction. However, for systems characterized by molecular overcrowding (e.g., when other molecules deny the reactants from the supposedly allowable space, and to stickiness, when the reactants are found along the surfaces of the reaction vessel) the kinetic orders for the reactions can exhibit non-integer values [24] found in power law formalism [23, 28, 29]. For instance, in intracellular environments, which are highly structured and characterized by molecular crowding, reactions in vivo are likely to take place on membranes or channels and as such, reactions follow fractal-like kinetics [6, 8, 19, 25]. The presence of power law kinetics in reaction systems thus motivated CRN-based studies on PLK systems ([9, 13, 20, 27] among others), some of which are extensions or modifications of existing results on MAK systems.

This contribution specifically shows that the result of Shinar and Feinberg on ACR applies to a class of PLK system called* power law kinetic systems with reactant-determined interactions* (denoted by “PL-RDK”). PL-RDK systems are kinetic systems with power law rate functions whose kinetic order vectors are identical for reactions with the same reactant complex. Since the kinetic orders of the mass-action rate functions are precisely the stoichiometric coefficients of the reactant complex, one can see that MAK is a special case of PL-RDK.

As an application, we employ the theorem to a power law approximation of the ODE system corresponding to a specific scenario in the pre-industrial carbon cycle model developed by Anderies et al. [1]. Particulary, for the pre-industrial scenario where there are anthropogenic causes that reduce the capacity of terrestrial carbon pool to store carbon, the power law approximation leads to an ACR-possessing PL-RDK system.

The rest of the paper is organized as follows: Section 2 assembles preliminary concepts in Chemical Reaction Network Theory required in stating and proving the results. Section 3 discusses the extension of the Shinar-Feinberg Theorem on ACR for PL-RDK systems. Section 4 applies the main result obtained from the previous section to a carbon cycle model. In Section 5, we summarize our results and outline some research perspectives.

2 Fundamentals of Chemical Reaction Networks and Kinetic Systems

We recall some fundamental notions about chemical reaction networks (CRNs) and chemical kinetic systems (CKS) assembled in [5, 27]. Some concepts introduced by Feinberg in [11, 12] are also reviewed.

Notation. We denote the real numbers by $\mathbb{R}$, the non-negative real numbers by $\mathbb{R}_{\geq 0}$ and the positive real numbers by $\mathbb{R}_{>0}$. Objects in the reaction systems are viewed as members of vector spaces. Suppose $\mathscr{I}$ is a finite index set. By $\mathbb{R}^{\mathscr{I}}$, we mean the usual vector space of real-valued functions with domain $\mathscr{I}$. For $x\in\mathbb{R}^{\mathscr{I}}$, the $i^{\text{th}}$ coordinate of $x$ is denoted by $x_{i}$, where $i\in\mathscr{I}$. The sets $\mathbb{R}_{\geq 0}^{\mathscr{I}}$ and $\mathbb{R}_{>0}^{\mathscr{I}}$ are called the non-negative and positive orthants of $\mathbb{R}^{\mathscr{I}}$, respectively. Addition, subtraction, and scalar multiplication in $\mathbb{R}^{\mathscr{I}}$are defined in the usual way. If $x\in\mathbb{R}_{>0}^{\mathscr{I}}$ and $y\in\mathbb{R}^{\mathscr{I}}$, we define $x^{y}\in\mathbb{R}_{>0}$ by $x^{y}=\prod_{i\in\mathscr{I}}x_{i}^{y_{i}}.$ The vector $\log x\in\mathbb{R}^{\mathscr{I}}$,where $x\in\mathbb{R}_{>0}^{\mathscr{I}}$, is given by $(\log x)_{i}=\log x_{i},\text{ for all }i\in\mathscr{I}.$ If $x,y\in\mathbb{R}^{\mathscr{I}}$, the standard scalar product $x\cdot y\in\mathbb{R}$ is defined by $x\cdot y=\sum_{i\in\mathscr{I}}x_{i}y_{i}.$ By the support of $x\in\mathbb{R}^{\mathscr{I}}$, denoted by $\text{supp }x$, we mean the subset of $\mathscr{I}$ assigned with non-zero values by $x$. That is, $\text{supp }x:=\{i\in\mathscr{I}|x_{i}\neq 0\}.$

Definition 1**.**

A chemical reaction network (CRN) $\mathscr{N}$ is a triple $(\mathscr{S},\mathscr{C},\mathscr{R})$ of three finite sets:

a set $\mathscr{S}=\{X_{1},X_{2},\dots,X_{m}\}$ of species; 2. 2.

a set $\mathscr{C}\subset\mathbb{R}^{\mathscr{S}}_{\geq 0}$ of complexes; 3. 3.

a set $\mathscr{R}=\{R_{1},R_{2},\dots,R_{r}\}\subset\mathscr{C}\times\mathscr{C}$ of reactions such that $(y,y)\notin\mathscr{R}$ for any $y\in\mathscr{C}$, and for each $y\in\mathscr{C}$, there exists $y^{\prime}\in\mathscr{C}$ such that either $(y,y^{\prime})\in\mathscr{R}$ or $(y^{\prime},y)\in\mathscr{R}$.

We denote the number of species with $m$, the number of complexes with $n$ and the number of reactions with $r$

A CRN can be viewed as a digraph $(\mathscr{C},\mathscr{R})$ with vertex-labelling. In particular, it is a digraph where each vertex $y\in\mathscr{C}$ has positive degree and stoichiometry, i.e., there is a finite set $\mathscr{S}$ of species such that $\mathscr{C}$ is a subset of $\mathbb{R}^{\mathscr{S}}_{\geq 0}$. The vertices are the complexes whose coordinates are in $\mathbb{R}^{\mathscr{S}}_{\geq 0}$, which are the stoichiometric coefficients. The arcs are precisely the reactions.

We use the convention that an element $R_{j}=(y_{j},y_{j}^{\prime})\in\mathscr{R}$ is denoted by $R_{j}:y_{j}\rightarrow y_{j}^{\prime}$. In this reaction, we say that $y_{j}$ is the reactant complex and $y^{\prime}_{j}$ is the product complex. Connected components of a CRN are called linkage classes, strongly connected components are called strong linkage classes, and strongly connected components without outgoing arcs are called terminal strong linkage classes. We denote the number of linkage classes with $\ell$, that of the strong linkage classes with $s\ell$, and that of terminal strong linkage classes with $t$. A complex is called terminal if it belongs to a terminal strong linkage class; otherwise, the complex is called nonterminal.

With each reaction $y\rightarrow y^{\prime}$, we associate a reaction vector obtained by subtracting the reactant complex $y$ from the product complex $y^{\prime}$. The stoichiometric subspace $S$ of a CRN is the linear subspace of $\mathbb{R}^{\mathscr{S}}$ defined by

[TABLE]

The rank of the CRN, $s$, is defined as $s=\dim S$.

Many features of CRNs can be examined by working in terms of finite dimensional spaces $\mathbb{R}^{\mathscr{S}}$ (species space), $\mathbb{R}^{\mathscr{C}}$ (complex space), and $\mathbb{R}^{\mathscr{R}}$ (reaction space). Suppose the set $\{\omega_{i}\in\mathbb{R}^{\mathscr{I}}\mid i\in\mathscr{I}\}$ forms the standard basis for $\mathbb{R}^{\mathscr{I}}$ where $\mathscr{I}=\mathscr{S,C}$ or $\mathscr{R}$. We recall four maps relevant in the study of CRNs: map of complexes, incidence map, stoichiometric map and Laplacian map.

Definition 2**.**

Let $\mathscr{N}=(\mathscr{S,C,R})$ be a CRN.

The map of complexes $\displaystyle{Y:\mathbb{R}^{\mathscr{C}}\rightarrow\mathbb{R}^{\mathscr{S}}}$ maps the basis vector $\omega_{y}$ to the complex $y\in\mathscr{C}$. 2. 2.

The incidence map $\displaystyle{I_{a}:\mathbb{R}^{\mathscr{R}}\rightarrow\mathbb{R}^{\mathscr{C}}}$ is the linear map defined by mapping for each reaction $\displaystyle{R_{j}:y_{j}\rightarrow y_{j}^{\prime}\in\mathscr{R}}$, the basis vector $\omega_{j}$ to the vector $\omega_{y_{j}^{\prime}}-\omega_{y_{j}}\in\mathscr{C}$. 3. 3.

The stoichiometric map $\displaystyle{N:\mathbb{R}^{\mathscr{R}}\rightarrow\mathbb{R}^{\mathscr{S}}}$ is defined as $N=Y\circ I_{a}$. 4. 4.

For each $k\in\mathbb{R}^{\mathscr{R}}_{>0}$ , the linear transformation $A_{k}:\mathbb{R}^{\mathscr{C}}\rightarrow\mathbb{R}^{\mathscr{C}}$ called Laplacian map is the mapping defined by

[TABLE]

where $x_{y}$ refers to the $y^{\text{th}}$ component of $x\in\mathbb{R}^{\mathscr{C}}$ relative to the standard basis.

The following result, named as the Structure Theorem of the Laplacian Kernel (STLK) by Arceo et al. in [5], is crucial in deriving important results in CRNT [11, 12].

Proposition 1** (Structure Theorem of the Laplacian Kernel (STLK), Prop. 4.1 [11]).**

Let $\mathscr{N}=\mathscr{(S,C,R)}$ be a CRN with terminal strong linkage classes $\mathscr{C}^{1},\mathscr{C}^{2},\dots,\mathscr{C}^{t}$. Let $k\in\mathbb{R}^{\mathscr{R}}_{>0}$ and $A_{k}$ its associated Laplacian. Then $\text{Ker }A_{k}$ has a basis $b^{1},b^{2},\dots,b^{t}$ such that $\text{supp }b^{i}=\mathscr{C}^{i}$ for all $i=1,2,\dots,t$.

A non-negative integer, called the deficiency, can be associated to each CRN. The deficiency of a CRN, denoted by $\delta$, is the integer defined by $\delta=n-\ell-s$. This index has been the center of many studies in CRNT due to its relevance in the dynamic behaviour of the system. In [11], Feinberg provided a geometric interpretation of deficiency: $\delta=\dim(\text{Ker }Y\cap\text{Im }I_{a})$. From this fact and the STLK, the following result follows.

Corollary 1** (Cor. 4.12 [11]).**

Let $\mathscr{N}=(\mathscr{S,C,R})$ be a CRN with deficiency $\delta$ and $t$ terminal strong linkage classes. Then for each $k\in\mathbb{R}^{\mathscr{R}}_{>0}$,

[TABLE]

By kinetics of a CRN, we mean the assignment of a rate function to each reaction in the CRN. It is defined formally as follows.

Definition 3**.**

A kinetics of a CRN $\mathscr{N}=(\mathscr{S},\mathscr{C},\mathscr{R})$ is an assignment of a rate function $\displaystyle{K_{j}:\Omega_{K}\to\mathbb{R}_{\geq 0}}$ to each reaction $R_{j}\in\mathscr{R}$, where $\Omega_{K}$ is a set such that $\mathbb{R}^{\mathscr{S}}_{>0}\subseteq\Omega_{K}\subseteq{\mathbb{R}}^{\mathscr{S}}_{\geq 0}$. A kinetics for a network $\mathscr{N}$ is denoted by

[TABLE]

The pair $(\mathscr{N},K)$ is called the chemical kinetic system (CKS).

The above definition is adopted from [30]. It is expressed in a more general context than what one typically finds in CRNT literature. For power law kinetic systems, one sets $\Omega_{K}=\mathbb{R}^{\mathscr{S}}_{>0}$. Here, we focus on the kind of kinetics relevant to our context:

Definition 4**.**

A chemical kinetics is a kinetics $K$ satisfying the positivity condition:

[TABLE]

Once a kinetics is associated with a CRN, we can determine the rate at which the concentration of each species evolves at composition $c\in\mathbb{R}^{\mathscr{S}}_{>0}$.

Definition 5**.**

The species formation rate function of a chemical kinetic system is the vector field

[TABLE]

The equation $dc/dt=f(c)$ is the ODE or dynamical system of the CKS. A positive equilibrium or steady state $c^{*}$ is an element of $\mathbb{R}^{\mathscr{S}}_{>0}$ for which $f(c^{*})=0$. The set of positive equilibria of a chemical kinetic system is denoted by $E_{+}(\mathscr{N},K)$.

Power law kinetics is defined by an $r\times m$ matrix $F=[F_{ij}]$, called the kinetic order matrix, and vector $k\in\mathbb{R}^{\mathscr{R}}$, called the rate vector.

Definition 6**.**

A kinetics $K:\mathbb{R}^{\mathscr{S}}_{>0}\rightarrow\mathbb{R}^{\mathscr{R}}$ is a power law kinetics (PLK) if

[TABLE]

with $k_{i}\in\mathbb{R}_{>0}$ and $F_{ij}\in\mathbb{R}$. A PLK system has reactant-determined kinetics (of type PL-RDK) if for any two reactions $R_{i}$, $R_{j}\in\mathscr{R}$ with identical reactant complexes, the corresponding rows of kinetic orders in $F$ are identical, i.e., $F_{ik}=F_{jk}$ for $k=1,...,m$.

An example of PL-RDK is the well-known mass-action kinetics (MAK), where the kinetic order matrix is the transpose of the matrix representation of the map of complexes $Y$ [11]. That is, a kinetics is a MAK if

[TABLE]

where $k_{j}\in\mathbb{R}_{>0}$, called rate constants. Note that $Y_{.,j}$ pertains to the stoichiometric coefficients of a reactant complex $y_{j}\in\mathscr{C}$.

Remark 1**.**

In [5], Arceo et al. discussed several sets of kinetics of a network and drew a “kinetic landscape”. They identified two main sets: the complex factorizable kinetics and its complement, the non-complex factorizable kinetics. Complex factorizable kinetics generalize the key structural property of MAK – that is, the species formation rate function decomposes as

[TABLE]

where $Y$ is the map of complexes, $A_{k}$ is the Laplacian map, and $\Psi_{k}:\mathbb{R}^{\mathscr{S}}_{\geq 0}\rightarrow\mathbb{R}^{\mathscr{C}}_{\geq 0}$ such that $I_{a}\circ K(x)=A_{k}\circ\Psi_{k}(x)$ for all $x\in\mathbb{R}^{\mathscr{S}}_{\geq 0}$. In the set of power law kinetics, PL-RDK is the subset of complex-factorizable kinetics.

We recall the definition of the $m\times n$ matrix $\widetilde{Y}$ from the work of Müller and Regensburger [20, 21]: For each reactant complex, the associated column of $\widetilde{Y}$ is the transpose of the kinetic order matrix row of the complex’s reaction, otherwise (i.e., for non-reactant complexes), the column is 0. We form the $\bm{T}$-matrix of a PL-RDK system by truncating away the columns of the non-reactant complexes in $\widetilde{Y}$, obtaining an $m\times n_{r}$ matrix, where $n_{r}$ denotes the number of reactant complexes [27].

3 Absolute Concentration Robustness in PL-RDK Systems

To illustrate absolute concentration robustness, we consider the following toy model:

[TABLE]

The map depicts a biochemical system involving transfer of material from two pools: $X_{2}$ to $X_{1}$ and $X_{1}$ to $X_{2}$, but with $X_{2}$ regulating the second process. Suppose the system evolves according to the following set of ODEs:

[TABLE]

The positive equilibrium of the system is attained when

[TABLE]

where $\Gamma$ is the conserved amount of total material. These equations indicate that whenever $\Gamma>(k_{1}/k_{2})^{2}$, a positive steady state exists. Furthermore, since $X_{1}$ has the same value in any steady state, the system exhibits ACR in $X_{1}$.

We define absolute concentration robustness in PL-RDK systems as follows:

Definition 7**.**

A PL-RDK system $(\mathscr{N},K)$ has absolute concentration robustness(ACR) in species $X_{i}\in\mathscr{S}$ if there exists $c^{*}\in E_{+}(\mathscr{N},K)$ and for every other $c^{**}\in E_{+}(\mathscr{N},K)$, we have $c^{**}_{i}=c^{*}_{i}$.

The following proposition adapts Theorem S3.15 found in supplementary online material of the paper of Shinar and Feinberg [26] to deal with PL-RDK systems.

Proposition 2**.**

Let $\mathscr{N}=(\mathscr{S,C,R})$ be a deficiency-one CRN. Suppose that $(\mathscr{N},K)$ is a PL-RDK system which admits a positive equilibrium $c^{*}$. If $y,y^{\prime}\in\mathscr{C}$ are nonterminal complexes, then each positive equilibrium $c^{**}$ of the system satisfies the equation

[TABLE]

We largely reproduce the proof of Shinar and Feinberg in the said supplementary material of their paper. Since in their proof, the sums are often taken over all complexes, we use the notation of Müller and Regensburger in [20, 21]:

[TABLE]

adjoining $n-n_{r}$ zero columns for the non-reactant complexes, where $n_{r}$ denotes the number of reactant complexes. Furthermore, we write $\widetilde{y}$ for $\widetilde{Y}_{\cdot,y}$.

Proof.

Assume that $c^{*}$ is a positive steady state of the PL-RDK system $(\mathscr{N},K)$. That is,

[TABLE]

For each $y\rightarrow y^{\prime}\in\mathscr{R}$, define the positive number $\kappa_{y\rightarrow y^{\prime}}$ by

[TABLE]

Thus, we obtain

[TABLE]

Suppose that $c^{**}$ is also a positive equilibrium of the system. Hence,

[TABLE]

Define

[TABLE]

With $\kappa\in\mathbb{R}^{\mathscr{R}}_{>0}$ given by Equation (3.6) and $\mu$ given by Equation (3.9), it follows from Equation (3.8) that

[TABLE]

Let $\bm{1}^{\mathscr{C}}\in\mathbb{R}^{\mathscr{C}}$ such that

[TABLE]

Observe that Equations (3.7) and (3.10) can be respectively written as

[TABLE]

Equivalently,

[TABLE]

[TABLE]

Therefore, $c^{*}$ and $c^{**}$ are positive equilibria of the PL-RDK system $(\mathscr{N},K)$ if and only if Equations (3.11) and (3.12) hold. From Corollary 1, we have

[TABLE]

for the CRN under consideration. Let $\{b^{1},b^{2},\dots,b^{t}\}\subset\mathbb{R}^{\mathscr{C}}_{\geq 0}$ be a basis for $\text{Ker }A_{\kappa}$ as in Proposition 1 (STLK). Since $\text{Ker }A_{\kappa}\subseteq\text{Ker }YA_{\kappa}$, this basis of $\text{Ker }A_{\kappa}$ can be extended to form a basis of $\text{Ker }YA_{\kappa}$. Recall from Equation (3.11) that $\bm{1}^{\mathscr{C}}$ is in $\text{Ker }YA_{\kappa}$. We assert that the set $\{\bm{1}^{\mathscr{C}},b^{1},b^{2},\dots,b^{t}\}$ is a basis for $\text{Ker }YA_{\kappa}$ (and hence, equality holds in Equation (3.13)). This follows if

[TABLE]

From Proposition 1, every element of $\text{Ker }A_{\kappa}$ must have its support contained entirely in the set of terminal complexes. However, the support of $\bm{1}^{\mathscr{C}}$ consists of all complexes. By assumption, there are nonterminal complexes and hence, $\bm{1}^{\mathscr{C}}$ cannot lie in $\text{Ker }A_{\kappa}$ (i.e., Equation (3.14) holds).

From Equation (3.12), there exist scalars $\lambda_{0},\lambda_{1},\dots,\lambda_{t}$ such that

[TABLE]

Observe that each vector $b^{i}$, $i=0,1,\dots,t$, has its support entirely on terminal complexes. This fact, along with Equation (3.15), implies that for each pair of nonterminal complexes $y\in\mathscr{C}$ and $y^{\prime}\in\mathscr{C}$, we have

[TABLE]

Since $y$ and $y^{\prime}$ are nonterminal, they are reactant complexes. Hence, Equation (3.16) may be written as

[TABLE]

which is equivalent to Equation (3.4) in Theorem 2. ∎

The extension of the Shinar-Feinberg Theorem on ACR to PL-RDK systems is stated as follows.

Theorem 1**.**

Let $\mathscr{N}=(\mathscr{S,C,R})$ be a deficiency-one CRN and suppose that $(\mathscr{N},K)$ is a PL-RDK system which admits a positive equilibrium. If $y,y^{\prime}\in\mathscr{C}$ are nonterminal complexes whose kinetic order vectors differ only in species $X_{i}$, then the system has ACR in $X_{i}$.

Proof.

Suppose $c^{*}$ and $c^{**}$ are positive equilibria of the PL-RDK system $(\mathscr{N},K)$. Observe that since $y,y^{\prime}\in\mathscr{C}$ are nonterminal complexes whose kinetic order vectors differ only in species $X_{i}$, we have

[TABLE]

for some nonzero $a$. Thus Equation (3.4) reduces to

[TABLE]

It follows that

[TABLE]

That is, the system has ACR in species $X_{i}$. ∎

The ODE system in Equation (3.2) can be translated into a dynamically equivalent CRN with associated kinetic order matrix by employing the notion of total CRN representation of Generalized Mass Action (GMA) systems, proposed by Arceo et al. [5]. GMA system is a canonical framework used in Biochemical Systems Theory (BST) wherein every mass transfer rate is approximated separately with a power law term, and these terms are added together, with a plus sign for incoming fluxes and a minus sign for outgoing fluxes [28, 29]. For BST-related concepts, the reader may refer to the BST tutorial in the Appendix of Arceo et al. [3].

The total CRN representation of a GMA system allows for the CRN-based analysis of the dynamical system. Viewed as a GMA system, the set of ODEs in (3.2) has the following total CRN representation:

[TABLE]

with associated kinetic order matrix $F$ given by

[TABLE]

The CRN in (3.18) is a deficiency-one network with nonterminal complexes $X_{1}+X_{2}$ and $X_{2}$ whose kinetic order rows differ only in $X_{1}$. The previous theorem indicates ACR in $X_{1}$, which agrees with the computation in (3.3).

The following simple proposition provides some examples for the ACR theorem for PL-RDK systems. As preparation, we recall some notions from [4, 27] which are used in the result. A PL-RDK is said to be reactant set linear independent (of type PL-RLK) if the columns of $T$ are linearly independent. We also recall the reactant matrix $Y_{\text{res}}$, which is obtained from the matrix representation of $Y$ by removing the columns corresponding to non-reactant complexes. Its image $\text{Im }Y_{\text{res}}$ is called the reactant subspace $R$, whose dimension $q$ is called the reactant rank of the CRN. The reactant deficiency $\delta_{\rho}$ is the difference between the number of reactant complexes $n_{r}$ and the reactant rank $q$.

Proposition 3**.**

Let $(\mathscr{S,C,R})$ be a deficiency one reaction network, which with PL-RDK, admits a positive equilibrium. Suppose the network has zero reactant deficiency, two nonterminal complexes $y,y^{\prime}\in\mathscr{C}$ differing only in $X_{j}$ and the map

[TABLE]

is given by

[TABLE]

Then the system is PL-RLK and has ACR in $X$.

Proof.

Since $\widehat{y}$ is an isomorphism, $T=\widehat{y}\circ Y_{\text{res}}$ is also an isomorphism. This implies that the system is PL-RLK. The kinetic order vector difference of $y$ and $y^{\prime}$ is $(0,\dots,ka_{j},\dots,0)$ for some nonzero real $k$ so that Theorem 1’s condition is fulfilled. ∎

4 Application to a Carbon Cycle Model

The pre-industrial carbon cycle model of Anderies et al. [1] is a simple mass balance which involves three interacting carbon pools: land, atmosphere and ocean. Pictorially, the system can be depicted using a biochemical map comprised of nodes that represent carbon pools, solid arrows that indicate transfer of carbon, and dashed arrows that indicate if a pool affects or modulates a process. Figure 4.1 presents the biochemical map of the model of interest. \

In our previous work [13], we reviewed the model’s design and underlying assumptions and described the parameters and ODEs present in the pre-industrial state of the carbon cycle model. We also approximated all rate processes by products of power law functions in order to obtain a GMA system approximation of the original system. The resulting ODEs of the approximation is given in (4.1):

[TABLE]

We also obtained in [13], using total CRN representation of [5], the following deficiency-one CRN representation for the model:

[TABLE]

Its associated kinetic order matrix is the transpose of the following $T$-matrix:

[TABLE]

In the Appendix, it is shown that there is a scenario in the pre-industrial state leading to a GMA system approximation such that the kinetic order vectors of the nonterminal vertices $A_{1}+2A_{2}$ and $A_{1}+A_{2}$ differ only in $A_{2}$; that is, $p_{1}-p_{2}=0$ and $q_{1}-q_{2}\neq 0$. In particular, this occurs when the human terrestrial carbon off-take term (which accounts for human activities that reduce the capacity of terrestrial pool to capture carbon such deforestation and land-use change) vanishes. Assuming the existence of a steady state, Theorem 1 indicates that the system has ACR in $A_{2}$. In fact, when $p_{1}=p_{2}$, steady state computation of the system in (4.1) yields the following equilibria set for the system:

[TABLE]

where $A_{0}=$ total conserved carbon at pre-industrial state.

5 Conclusion and Outlook

In conclusion, we summarize our results and outline some perspectives for further research.

We modified the Shinar-Feinberg Theorem on ACR for mass-action systems to include PL-RDK systems, a kinetic system more general than mass-action systems. 2. 2.

The theorem is applied to a power law approximation of Anderies et al.’s Earth’s carbon cycle in its pre-industrial state. The analysis reveals that there is a scenario in the pre-industrial state which yields a power law approximation where there is ACR in the atmospheric carbon pool. Specifically, the power law approximation leads to an ACR-possessing PL-RDK system when the human off-take coefficient, which accounts for the which accounts for human activities that reduce the capacity of terrestrial pool to sequester carbon, vanishes. 3. 3.

The investigation of other forms of “concentration robustness” identified by Dexter et al. [10] for PL-RDK systems offers a further interesting research perspective. 4. 4.

The extension of the stochastic analysis of CRNs with ACR of Anderson et al. [2] for PL-RDK systems is another promising area for further investigation.

Acknowledgements

NTF acknowledges the support of the Department of Science and Technology-Science Education Institute (DOST-SEI), Philippines through the ASTHRDP Scholarship grant and Career Incentive Program (CIP). ARL and LFR held research fellowships from De La Salle University and would like to acknowledge the support of De La Salle University’s Research Coordination Office.

Appendix A Pre-industrial Carbon Cycle Model of Anderies et al.

The complete set of ODEs for the pre-industrial state is given by

[TABLE]

where

[TABLE]

For the description of the parameters, the reader is referred to [1] and the Appendix of [13]. The parameter values are identical to the values used in [13] but with $\alpha=0$. This particular parameter is assigned as the human terrestrial carbon off-take rate. It is associated to human activities such as clearing, burning or farming, which reduce the capacity of land to capture carbon.

A power law approximation of the ODE system at an operating point is obtained to generate a Generalized Mass Action (GMA) System [28, 29]. Mathematically, GMA system approximation is equivalent to Taylor approximation up to the linear term in logarithmic space. The function $V(X_{1},X_{2},\dots,X_{m})$ can be approximated by $\displaystyle{V=\alpha X_{1}^{p_{1}}X_{2}^{p_{2}}\cdots X_{m}^{p_{m}}}$ at an operating point where

[TABLE]

Table LABEL:tab:power_law_approx presents the four carbon fluxes present in the pre-industrial state of the Anderies et al. model, and their corresponding rate functions. Furthermore, the last column lists their respective target power law approximation. The last two functions, $a_{m}A_{2}$ and $a_{m}\beta A_{3}$, are already in the desired format and are thus, kept as is. To compute for the kinetic orders (and rate constants), we apply (A.2). By taking the parameter values used in [13] but with $\alpha=0$, and assuming the initial values to be $A_{1}=2850/4500$, $A_{2}=750/4500$ and $A_{3}=900/4500$ (as in[1]), the ODE system in (A.1) reaches the following steady state: $A_{1}=0.7,\quad A_{2}=0.15\text{ and }A_{3}=0.15.$

The algebraic calculations are implemented in Mathematica as shown in Figure A.1. When $\alpha=0$ (i.e., the human off-take term vanishes),

[TABLE]

For the power law approximation, we choose values close to the equilibrium point as operating point: $A_{1}=0.69$, $A_{2}=0.155$ and $A_{3}=0.155$. Consequently, we obtain

[TABLE]