Reaction Network Analysis
of Metabolic Insulin Signaling

Patrick Vincent N. Lubenia Systems and Computational Biology Research Unit, Center for Natural Sciences and Environmental Research, 2401 Taft Avenue, Manila, 0922, Metro Manila, Philippines Eduardo R. Mendoza Systems and Computational Biology Research Unit, Center for Natural Sciences and Environmental Research, 2401 Taft Avenue, Manila, 0922, Metro Manila, Philippines Department of Mathematics and Statistics, De La Salle University, 2401 Taft Avenue, Manila, 0922, Metro Manila, Philippines Max Planck Institute of Biochemistry, Am Klopferspitz 18, 82152, Martinsried near Munich, Germany Angelyn R. Lao Systems and Computational Biology Research Unit, Center for Natural Sciences and Environmental Research, 2401 Taft Avenue, Manila, 0922, Metro Manila, Philippines Department of Mathematics and Statistics, De La Salle University, 2401 Taft Avenue, Manila, 0922, Metro Manila, Philippines

Abstract

Absolute concentration robustness (ACR) and concordance are novel concepts in the theory of robustness and stability within Chemical Reaction Network Theory. In this paper, we have extended Shinar and Feinberg’s reaction network analysis approach to the insulin signaling system based on recent advances in decomposing reaction networks. We have shown that the network with 20 species, 35 complexes, and 35 reactions is concordant, implying at most one positive equilibrium in each of its stoichiometric compatibility class. We have obtained the system’s finest independent decomposition consisting of 10 subnetworks, a coarsening of which reveals three subnetworks which are not only functionally but also structurally important. Utilizing the network’s deficiency-oriented coarsening, we have developed a method to determine positive equilibria for the entire network. Our analysis has also shown that the system has ACR in 8 species all coming from a deficiency zero subnetwork. Interestingly, we have shown that, for a set of rate constants, the insulin-regulated glucose transporter GLUT4 (important in glucose energy metabolism), has stable ACR.

Keywords: independent decomposition, metabolic insulin signaling, positive equilibria, reaction network, subnetwork

1 Introduction

Between 2010 and 2015, Shinar and Feinberg published a series of papers based on the novel concepts of absolute concentration robustness (ACR) and concordance, which one may view as the beginnings of a theory of robustness within Chemical Reaction Network Theory (CRNT) [17, 24, 25, 26, 27]. ACR characterizes the invariance of the concentrations of a species at all positive equilibria of a kinetic system and, from experimental observations in Escherichia coli subsystems, the authors extracted sufficient (mathematical) conditions for the property. They related ACR to bifunctionality of enzymes and viewed the condition as a “structural source” [24] or “design principle” [25] of robustness. Concordance, on the other hand, is seen to indicate “architectures that, by their very nature, enforce duller, more restrictive behavior despite what might be great intricacy in the interplay of many species, even independent of values that kinetic parameters might take” [26]. Shinar and Feinberg provided small models of biological systems such as the osmotic pressure response system EnvZ-OmpR of Escherichia coli for ACR and calcium dynamics of olfactory cilia for concordance. The goal of this paper is to explore the extension of their reaction network analysis approach to larger models of biochemical systems based on recent advances in decomposing reaction networks [12, 13, 14] using the example of a widely used model of the insulin signaling system [23].

The insulin signaling system is an important metabolic system that, upon binding of insulin to its receptor at the cell surface, initiates the uptake of glucose into the cell. The analysis of this process in muscle cells, hepatocytes, cells of adipose tissue, and (most recently) neurons is crucial for understanding the underlying mechanisms of insulin resistance. This reduced ability of cells to use available insulin for energy metabolism is viewed as a common factor in diseases such as obesity, type 2 diabetes, metabolic syndrome, and cancer. More recently, brain insulin resistance has received increased attention from neuroscientists in connection with mild cognitive impairment and Alzheimer’s Disease (AD) [2, 22]. Our studies of the metabolic aspects of AD, in fact, provided the motivation for this analysis of a mathematical model of the insulin signaling system [28].

The complexity of the insulin signaling system, both in terms of the number of molecular components involved as well as the intricate combination of positive and negative feedback loops, clearly warrants the application of mathematical modeling and computational tools. A natural starting point for our studies of such models is the seminal work of Sedaghat et al [23]. This widely-cited work has been applied in various contexts and the authors conveniently provided WinPP source files for the system. We utilized the Hars-Tóth criterion presented in [4] to derive a chemical reaction network underlying the system of ordinary differential equations (ODEs) of the Sedaghat et al model, leading to its realization as a mass action system with 20 species, 35 complexes, and 35 reactions.

The first main result of our analysis of the Sedaghat et al system is that the underlying network is concordant. Concordance is a property abstracted from the continuous flow stirred tank reactor model widely used in chemical engineering and its occurrence in complex biological systems is not straightforward. To date, only two smaller systems—the previously-mentioned calcium signaling and the Wnt signaling [8]—have been shown to have the property. Concordance has numerous structural and kinetic consequences including monostationarity for the insulin signaling system (see Section 2.3). A detailed discussion of concordance properties can be found in Chapter 10 of Feinberg’s recent book [8].

The remaining main results derive from the systematic use of decomposition theory.

Decomposition theory was initiated by Feinberg in his 1987 review [7] where he also introduced the important concept of an independent decomposition, a decomposition wherein the stoichiometric subspace of the whole network is the direct sum of the stoichiometric subspaces of the subnetworks. He highlighted its importance by showing that in an independent decomposition the set of positive equilibria of the whole network is the intersection of the equilibria sets of the subnetworks. Hence, if a species has ACR in a subnetwork of an independent decomposition, it also has ACR in the whole network since the latter’s equilibria set is contained in that of the subnetwork. Recent results of Hernandez and De la Cruz [14] provided a criterion and procedure for determining the existence of a nontrivial independent decomposition (the trivial independent decomposition is the network itself).

The second main result is the existence of nontrivial independent decompositions of the network. This property allows insights into structural relationships between positive equilibria of the whole network and its subnetworks. Since the algorithm from Hernandez and De la Cruz [14] determines the finest independent decomposition and independence is invariant under decomposition coarsening, all independent decompositions can be generated from the result. We developed a Matlab code for the said algorithm and applied it to the insulin signaling system.

Third, it is shown that the networks of three functional modules—the insulin receptor binding and recycling subsystem, the postreceptor signaling subsystem, and the GLUT4 translocation subsystem—discussed by Sedaghat et al also form an independent decomposition. These subsystems, hence, not only are functionally but also structurally significant, i.e., positive equilibria of the whole system come from equilibria of these systems. Further details can be found in Section 3.1.

Fourth is the existence of a large weakly reversible, deficiency zero subnetwork constituting a well-understood part of the system and which is also the source of all 8 ACR species.

For the ACR analysis, we implemented the algorithm of Fontanil et al [12] in Matlab and discovered Shinar-Feinberg reaction pairs in appropriate low-deficiency subnetworks of coarsenings of the finest independent decomposition. We found that the system has ACR in 8 (out of 20) species. This restricts the variability of the positive equilibria and suggests that this “structural source of robustness” may be an important factor in the system’s overall robustness, a property which, according to a previous study by Dexter et al [5], is not common in biochemical systems. Overall, however, there is still a high variation in equilibria composition due to the lack of ACR species in the deficiency 7 subnetwork of the network’s deficiency-oriented decomposition. To our knowledge, this is the first large kinetic system for which an ACR assessment has been documented.

The last main result is the discovery of ACR of the essential glucose transporter GLUT4 which, coupled with adequate glucose supply, enables reliable cellular energy production. Further details are discussed in Section 4.3.

The paper is organized as follows: The construction of the metabolic insulin signaling reaction network, as well as its basic properties are described in Section 2. Section 3 discusses the finest nontrivial independent decomposition of the insulin signaling system and its deficiency-oriented coarsening. In Section 4, a brief review of relevant ACR results is followed by ACR analyses based on decompositions of the network. A short summary and the future direction of our research conclude the paper in Section 5. Basic concepts and notations related to CRNT, as well as some detailed computations relevant to the paper, are provided in the Appendix.

2 Model Realization as a Mass Action System

In this section, we present the ODE system of Sedaghat et al [23] and its mass action system realization. We then analyze various properties of the underlying chemical reaction network (CRN), including its positive dependency, $t$ -minimality, nonconservativity, and concordance.

2.1 Insulin Signaling Model

We consider the system of ODEs in Sedaghat et al [23] which has only the state variables $X_{2}$ to $X_{21}$ . We streamlined the notation to clearly identify the 35 reactions involved in the metabolic insulin signaling network (see Appendix B for the description and units of the variables):

	$\displaystyle\dot{X}_{2}=k_{2}X_{3}+k_{6}X_{5}-k_{1}X_{2}+k_{8}X_{6}-k_{7}X_{2}$
	$\displaystyle\dot{X}_{3}=k_{1}X_{2}-k_{2}X_{3}-k_{5}X_{3}$
	$\displaystyle\dot{X}_{4}=k_{3}X_{5}-k_{4}X_{4}+k_{10}X_{7}-k_{9}X_{4}$
	$\displaystyle\dot{X}_{5}=k_{5}X_{3}+k_{4}X_{4}-k_{3}X_{5}-k_{6}X_{5}+k_{12}X_{8}-k_{11}X_{5}$
	$\displaystyle\dot{X}_{6}=k_{13}-k_{14}X_{6}+k_{15}X_{7}+k_{16}X_{8}+k_{7}X_{2}-k_{8}X_{6}$
	$\displaystyle\dot{X}_{7}=k_{9}X_{4}-k_{10}X_{7}-k_{15}X_{7}$
	$\displaystyle\dot{X}_{8}=k_{11}X_{5}-k_{12}X_{8}-k_{16}X_{8}$
	$\displaystyle\dot{X}_{9}=k_{19}X_{10}-k_{17}X_{9}X_{4}-k_{18}X_{9}X_{5}$
	$\displaystyle\dot{X}_{10}=k_{17}X_{9}X_{4}+k_{18}X_{9}X_{5}+k_{21}X_{12}-k_{19}X_{10}-k_{20}X_{11}X_{10}$
	$\displaystyle\dot{X}_{11}=k_{21}X_{12}-k_{20}X_{10}X_{11}$
	$\displaystyle\dot{X}_{12}=k_{20}X_{10}X_{11}-k_{21}X_{12}$
	$\displaystyle\dot{X}_{13}=k_{22}X_{12}X_{14}+k_{24}X_{15}-k_{23}X_{13}-k_{25}X_{13}$
	$\displaystyle\dot{X}_{14}=k_{23}X_{13}-k_{22}X_{12}X_{14}$
	$\displaystyle\dot{X}_{15}=k_{25}X_{13}-k_{24}X_{15}$
	$\displaystyle\dot{X}_{16}=k_{27}X_{17}-k_{26}X_{13}X_{16}$
	$\displaystyle\dot{X}_{17}=k_{26}X_{13}X_{16}-k_{27}X_{17}$
	$\displaystyle\dot{X}_{18}=k_{29}X_{19}-k_{28}X_{13}X_{18}$
	$\displaystyle\dot{X}_{19}=k_{28}X_{13}X_{18}-k_{29}X_{19}$
	$\displaystyle\dot{X}_{20}=k_{31}X_{21}-k_{30}X_{20}-k_{32}X_{17}X_{20}-k_{33}X_{19}X_{20}+k_{34}-k_{35}X_{20}$
	$\displaystyle\dot{X}_{21}=k_{30}X_{20}+k_{32}X_{17}X_{20}+k_{33}X_{19}X_{20}-k_{31}X_{21}$

Sedaghat et al used a biochemical map to derive the ODE system and an incomplete extract of the map is shown in their paper. For better visual orientation of the reader, we combined information from the paper’s text and the ODE system to reconstruct the complete biochemical map shown in Figure 1.

Refer to caption — Figure 1: Biochemical map of the insulin signaling system where $X_{2},\ldots,X_{21}$ are the species of the network (see Appendix B for the description and units of the variables); $k_{1},\ldots,k_{35}$ are the rate constants of the reactions; and solid lines represent mass transfer reactions while broken lines represent regulatory reactions

2.2 Application of the Hars-Tóth Criterion for Mass Action System Realization

Theorem 4 of Chellaboina et al [4], called the Hars-Tóth criterion, guarantees a mass action system realization such that the ODE system with species $X_{1},\ldots,X_{m}$ and $r$ reactions is given by $\dot{X}=f(X)=(B-A)^{\top}(k\circ X^{A})$ where $\circ$ represents componentwise multiplication, $A=[a_{ij}]$ , $B=[b_{ij}]$ , $k=[k_{1},\ldots,k_{r}]^{\top}$ , $X=[X_{1},\ldots,X_{m}]^{\top}$ , and $X^{A}$ is the element of $\mathbb{R}^{m}$ with $i$ th component $X_{1}^{a_{i1}}\cdots X_{m}^{a_{im}}$ for $i=1,\ldots,r$ and $j=1,\ldots,m$ .

In the insulin signaling model $\dot{X}=f(X)$ in Section 2.1, it is easy to check that the function $f_{j}(X_{2},\ldots,X_{j-1},0,X_{j+1},\ldots,X_{21})$ is a multivariate polynomial with nonnegative coefficients for each $j=2,\ldots,21$ . For instance, for $j=2$ and $j=3$ , we have

	$\displaystyle f_{2}(0,X_{3},X_{4}\ldots,X_{21})=k_{2}X_{3}+k_{6}X_{5}+k_{8}X_{6}$
	$\displaystyle f_{3}(X_{2},0,X_{4},\ldots,X_{21})=k_{1}X_{2}.$

Similar computations for $j=4,\ldots,21$ yield the same conclusion. Therefore, by the Hars-Tóth criterion, there is a reaction network of the form $AX\xrightarrow{k}BX$ such that $f(X)=(B-A)^{\top}(k\circ X^{A})$ where $A$ and $B$ have nonnegative integer entries (see Appendix C for the detailed computation). The CRN corresponding to the insulin signaling model is:

	$\displaystyle R_{1}:X_{2}\rightarrow X_{3}$
	$\displaystyle R_{2}:X_{3}\rightarrow X_{2}$
	$\displaystyle R_{3}:X_{5}\rightarrow X_{4}$
	$\displaystyle R_{4}:X_{4}\rightarrow X_{5}$
	$\displaystyle R_{5}:X_{3}\rightarrow X_{5}$
	$\displaystyle R_{6}:X_{5}\rightarrow X_{2}$
	$\displaystyle R_{7}:X_{2}\rightarrow X_{6}$
	$\displaystyle R_{8}:X_{6}\rightarrow X_{2}$
	$\displaystyle R_{9}:X_{4}\rightarrow X_{7}$
	$\displaystyle R_{10}:X_{7}\rightarrow X_{4}$
	$\displaystyle R_{11}:X_{5}\rightarrow X_{8}$
	$\displaystyle R_{12}:X_{8}\rightarrow X_{5}$
	$\displaystyle R_{13}:0\rightarrow X_{6}$
	$\displaystyle R_{14}:X_{6}\rightarrow 0$
	$\displaystyle R_{15}:X_{7}\rightarrow X_{6}$
	$\displaystyle R_{16}:X_{8}\rightarrow X_{6}$
	$\displaystyle R_{17}:X_{9}+X_{4}\rightarrow X_{10}+X_{4}$
	$\displaystyle R_{18}:X_{9}+X_{5}\rightarrow X_{10}+X_{5}$
	$\displaystyle R_{19}:X_{10}\rightarrow X_{9}$
	$\displaystyle R_{20}:X_{10}+X_{11}\rightarrow X_{12}$
	$\displaystyle R_{21}:X_{12}\rightarrow X_{10}+X_{11}$
	$\displaystyle R_{22}:X_{14}+X_{12}\rightarrow X_{13}+X_{12}$
	$\displaystyle R_{23}:X_{13}\rightarrow X_{14}$
	$\displaystyle R_{24}:X_{15}\rightarrow X_{13}$
	$\displaystyle R_{25}:X_{13}\rightarrow X_{15}$
	$\displaystyle R_{26}:X_{16}+X_{13}\rightarrow X_{17}+X_{13}$
	$\displaystyle R_{27}:X_{17}\rightarrow X_{16}$
	$\displaystyle R_{28}:X_{18}+X_{13}\rightarrow X_{19}+X_{13}$
	$\displaystyle R_{29}:X_{19}\rightarrow X_{18}$
	$\displaystyle R_{30}:X_{20}\rightarrow X_{21}$
	$\displaystyle R_{31}:X_{21}\rightarrow X_{20}$
	$\displaystyle R_{32}:X_{20}+X_{17}\rightarrow X_{21}+X_{17}$
	$\displaystyle R_{33}:X_{20}+X_{19}\rightarrow X_{21}+X_{19}$
	$\displaystyle R_{34}:0\rightarrow X_{20}$
	$\displaystyle R_{35}:X_{20}\rightarrow 0.$

2.3 CRNT Analysis of the Mass Action System: Key Basic Properties

From this point onward, we denote the CRN constructed in the previous section for the metabolic insulin signaling network as $\mathscr{N}=(\mathscr{S},\mathscr{C},\mathscr{R})$ with mass action kinetics $K$ , stoichiometric subspace $S$ , set of species $\mathscr{S}=\{X_{2},\ldots,X_{21}\}$ , set of complexes $\mathscr{C}$ , and set of reactions $\mathscr{R}=\{R_{1},\ldots,R_{35}\}$ .

Table 1 presents the network numbers of $\mathscr{N}$ . These numbers provide a cursory analysis of the structural and dynamical properties of the CRN. Properties inferred by Table 1 include non-weak reversibility (since $s\ell=24\neq 13=\ell$ ), $t$ -minimality (since $t=13=\ell$ ), branching (since $r=35>24=n_{r}$ ), and non-point terminality and non-cycle terminality (since $t\neq n-n_{r}=11$ and $n-n_{r}\neq 0$ , respectively).

Table 1: Network numbers of

\mathscr{N}

Network numbers
Species	$m$	20
Complexes	$n$	35
Reactant complexes	$n_{r}$	24
Reversible reactions	$r_{\text{rev}}$	10
Irreversible reactions	$r_{\text{irrev}}$	15
Reactions	$r$	35
Linkage classes	$\ell$	13
Strong linkage classes	$s\ell$	24
Terminal strong linkage classes	$t$	13
Rank	$s$	15
Reactant rank	$q$	20
Deficiency	$\delta$	7
Reactant deficiency	$\delta_{p}$	4

Recall that deficiency measures the linear dependence of the network’s reactions, i.e., the higher the deficiency, the higher the linear dependence. The insulin signaling system is of deficiency $\delta=7$ , indicative of high complexity of the network.

2.3.1 Positive Dependency

CRNToolbox [10] results show that the metabolic insulin signaling network is positive dependent, implying that a set of rate constants exists for which the mass action system can admit a positive equilibrium.

Proposition 1 shows the key implication of positive dependency in a general form for rate constant-interaction map decomposable (RID) kinetic systems (a kinetics $K$ is RID if it assigns to each $i$ th reaction a function $K_{i}(x)=k_{i}I_{i}(X)$ with rate constant $k_{i}>0$ and interaction map $I_{i}(X)\in\mathbb{R}^{\mathscr{R}}$ for $X\in\mathbb{R}_{\geq 0}^{\mathscr{S}}$ ). More information about RID kinetic systems can be found in Nazareno et al [21] where they were introduced. The proposition with the implication of positive dependency was shown for mass action systems by Feinberg [7].

Proposition 1.

For any RID kinetics $K$ on a positive dependent network $\mathscr{N}$ , there are rate constants such that the set of positive equilibria $E_{+}(\mathscr{N},K)\neq\varnothing$ .

Proof.

Since $\mathscr{N}$ is positive dependent, then for each reaction $i:C_{i}\rightarrow C_{i}^{\prime}$ there is a positive number $\alpha_{i}$ such that

\sum_{i}\alpha_{i}(C_{i}^{\prime}-C_{i})=0.

For a positive vector $X^{*}$ , $I_{i}(X^{*})>0$ by definition of RID kinetics. Set $\displaystyle k_{i}^{*}=\frac{\alpha_{i}}{I_{i}(X^{*})}$ . Then

	$\displaystyle f(X^{*})$	$\displaystyle=\sum_{i}K_{i}(X^{*})(C_{i}^{\prime}-C_{i})$
		$\displaystyle=\sum_{i}k_{i}^{}I_{i}(X^{})(C_{i}^{\prime}-C_{i})$
		$\displaystyle=\sum_{i}\alpha_{i}(C_{i}^{\prime}-C_{i})$
		$\displaystyle=0$

i.e., $X^{*}\in E_{+}(\mathscr{N},K)$ . ∎

Remark 1.

Proposition 1 is a generalization of Lemma 3.5.3 of Feinberg [8].

2.3.2 $t$ -minimality

Table 1 shows that there is only one terminal strong linkage class in each linkage class since $t=13=\ell$ , i.e., the network is $t$ -minimal. Feinberg and Horn [9] showed that this implies the coincidence of the kinetic and stoichiometric subspaces.

The noncoincidence of the kinetic and stoichiometric subspaces is closely related to the degeneracy of equilibria. In [11], the authors showed that if the two subspaces differ, then all equilibria are degenerate. In his book, Feinberg describes anomalies that can occur if the two subspaces do not coincide (Section 3.A.1 of [8]).

Thus, the $t$ -minimality of the metabolic insulin signaling network is a necessary condition for the existence of nondegenerate equilibria as detailed in Remark 4.11 of [11].

2.3.3 Nonconservativity

Recall that a reaction network $(\mathscr{S},\mathscr{C},\mathscr{R})$ is conservative whenever the orthogonal complement $S^{\perp}$ of its stoichiometric subspace $S$ contains a strictly positive member of $\mathbb{R}^{\mathscr{S}}$ , i.e., $S^{\perp}\cap\mathbb{R}_{>0}^{\mathscr{S}}\neq\varnothing$ . Otherwise, the network is called nonconservative.

The metabolic insulin signaling network has $m=20$ species and rank $s=15$ . Thus, $S^{\perp}$ has dimension 5. A basis for $S^{\perp}$ is the set

	$\displaystyle\{$	$\displaystyle[0,0,0,0,0,0,0,1,1,-1,0,0,0,0,0,0,0,0,0,0]^{\top},$
		$\displaystyle[0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0]^{\top},$
		$\displaystyle[0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0]^{\top},$
		$\displaystyle[0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0]^{\top},$
		$\displaystyle[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0]^{\top}\}$

which implies that $S^{\perp}\cap\mathbb{R}_{>0}^{\mathscr{S}}=\varnothing$ . Therefore, the metabolic insulin signaling network is nonconservative. The positive implications of this property for ACR will be discussed in Section 4.

2.3.4 Concordance

The most important specific result of the CRNT analysis is the network’s concordance.

Concordance is closely related to weakly monotonic kinetics as shown by Proposition 4.8 of [26]. It shows that, for any weakly monotonic kinetic system, injectivity (and, therefore, the absence of distinct stoichiometrically compatible equilibria, at least one of which is positive) can be precluded merely by establishing concordance of the underlying reaction network. The converse of the said proposition is not true, i.e., there can be a weakly monotonic kinetic system (in fact, a mass action system) that is injective even when its underlying reaction network is not concordant.

Theorem 4.11 of [26] shows that the class of concordant networks is precisely the class of networks that are injective for every assignment of a weakly monotonic kinetics.

In summary, concordance implies injectivity (of the mass action kinetics) and hence monostationarity, i.e., there is at most one positive equilibrium in each stoichiometric compatibility class. We ran the Concordance Test in the CRNToolbox which showed that the metabolic insulin signaling network is concordant. Running the Mass Action Injectivity Test and the Higher Deficiency Test confirms that the kinetics of the insulin signaling system is injective and that the system is monostationary, respectively.

3 Independent Decompositions of the Insulin Signaling System

In this section, we present the finest nontrivial independent decomposition of the metabolic insulin signaling network. We also present a deficiency-oriented coarsening of the decomposition which will be useful in the analysis of the network’s ACR.

3.1 The Finest Nontrivial Independent Decomposition

Hernandez and De la Cruz [14] provide an algorithm for constructing an independent decomposition. Their procedure utilizes a coordinate graph that represents the network using the reaction vectors of the CRN. A nontrivial independent decomposition is generated if the coordinate graph is not connected and, in this case, each connected component of the graph constitutes a partition of the set of reaction vectors of the CRN. Following their procedure, we developed a Matlab code to run the said algorithm and got the independent decomposition $\mathscr{N}=\{R_{1},\ldots,R_{35}\}$ consisting of 10 subnetworks:

	$\displaystyle\mathscr{N}_{1}=\{R_{1},\ldots,R_{12},R_{15},R_{16}\}$
	$\displaystyle\mathscr{N}_{2}=\{R_{13},R_{14}\}$
	$\displaystyle\mathscr{N}_{3}=\{R_{17},R_{18},R_{19}\}$
	$\displaystyle\mathscr{N}_{4}=\{R_{20},R_{21}\}$
	$\displaystyle\mathscr{N}_{5}=\{R_{22},R_{23}\}$
	$\displaystyle\mathscr{N}_{6}=\{R_{24},R_{25}\}$
	$\displaystyle\mathscr{N}_{7}=\{R_{26},R_{27}\}$
	$\displaystyle\mathscr{N}_{8}=\{R_{28},R_{29}\}$
	$\displaystyle\mathscr{N}_{9}=\{R_{30},\ldots,R_{33}\}$
	$\displaystyle\mathscr{N}_{10}=\{R_{34},R_{35}\}.$

Furthermore, the independent decomposition obtained is precisely the finest independent decomposition of $\mathscr{N}$ [16]. The network numbers of the subnetworks are presented in Table 2.

Table 2: Network numbers of the subnetworks

\mathscr{N}_{i}

of the finest independent decomposition of the metabolic insulin signaling network

\mathscr{N}

Network numbers		$\mathscr{N}_{1}$	$\mathscr{N}_{2}$	$\mathscr{N}_{3}$	$\mathscr{N}_{4}$	$\mathscr{N}_{5}$	$\mathscr{N}_{6}$	$\mathscr{N}_{7}$	$\mathscr{N}_{8}$	$\mathscr{N}_{9}$	$\mathscr{N}_{10}$
Species	$m$	7	1	4	3	3	2	3	3	4	1
Complexes	$n$	7	2	6	2	4	2	4	4	6	2
Reactant complexes	$n_{r}$	7	2	3	2	2	2	2	2	4	2
Reversible reactions	$r_{\text{rev}}$	5	1	0	1	0	1	0	0	1	1
Irreversible reactions	$r_{\text{irrev}}$	4	0	3	0	2	0	2	2	2	0
Reactions	$r$	14	2	3	2	2	2	2	2	4	2
Linkage classes	$\ell$	1	1	3	1	2	1	2	2	3	1
Strong linkage classes	$s\ell$	1	1	6	1	4	1	4	4	5	1
Terminal strong linkage classes	$t$	1	1	3	1	2	1	2	2	3	1
Rank	$s$	6	1	1	1	1	1	1	1	1	1
Reactant rank	$q$	7	1	3	2	2	2	2	2	4	1
Deficiency	$\delta$	0	0	2	0	1	0	1	1	2	0
Reactant deficiency	$\delta_{p}$	0	1	0	0	0	0	0	0	0	1

Consider the coarsening $\mathscr{N}=\mathscr{N}_{1}^{*}\cup\mathscr{N}_{2}^{*}\cup\mathscr{N}_{3}^{*}$ where $\mathscr{N}_{1}^{*}=\mathscr{N}_{1}\cup\mathscr{N}_{2}$ , $\mathscr{N}_{2}^{*}=\mathscr{N}_{3}\cup\ldots\cup\mathscr{N}_{8}$ , and $\mathscr{N}_{3}^{*}=\mathscr{N}_{9}\cup\mathscr{N}_{10}$ . This decomposition shows the three subsystems considered by Sedaghat et al in the construction of their system: the insulin receptor binding and recycling subsystem ( $\mathscr{N}_{1}^{*}$ ), the postreceptor signaling subsystem ( $\mathscr{N}_{2}^{*}$ ), and the GLUT4 translocation subsystem ( $\mathscr{N}_{3}^{*}$ ). It is significant to note that CRNT analysis using independent decomposition reveals that these subnetworks of the insulin signaling system are not just functionally but also structurally important, i.e., positive equilibria of the whole system come from equilibria of these systems.

Recall that a decomposition is said to be bi-independent if it is both independent and incidence independent.

Proposition 2.

The decomposition $\mathscr{N}=\mathscr{N}_{1}\cup\ldots\cup\mathscr{N}_{10}$ is bi-independent.

Proof.

Table 2 shows that $\delta=\delta_{1}+\ldots+\delta_{10}$ where $\delta_{i}$ is the deficiency of subnetwork $\mathscr{N}_{i}$ . Hence, by Proposition 8 of [6], the decomposition is bi-independent. ∎

Note that since Proposition 7 of [6] implies that every coarsening is also incidence independent, then every coarsening is bi-independent as well.

3.2 Deficiency-Oriented Coarsening

We next consider a deficiency-oriented coarsening of the independent decomposition. The deficiency-oriented coarsening is $\mathscr{N}=\mathscr{N}_{A}\cup\mathscr{N}_{B}$ where $\mathscr{N}_{A}=\mathscr{N}_{1}\cup\mathscr{N}_{2}\cup\mathscr{N}_{4}\cup\mathscr{N}_{6}\cup\mathscr{N}_{10}$ and $\mathscr{N}_{B}=\mathscr{N}_{3}\cup\mathscr{N}_{5}\cup\mathscr{N}_{7}\cup\mathscr{N}_{8}\cup\mathscr{N}_{9}$ . Let $\mathscr{S}_{i}$ , $\mathscr{C}_{i}$ , $K_{i}$ , and $S_{i}$ be the set of species, set of complexes, kinetics, and stoichiometric subspace, respectively, of $\mathscr{N}_{i}$ for $i=A,B$ .

In $\mathscr{N}_{A}$ , we consider together all subnetworks of the finest independent decomposition which have a deficiency of 0. On the other hand, we put together all nonzero-deficiency subnetworks in $\mathscr{N}_{B}$ . The network numbers of $\mathscr{N}_{A}$ and $\mathscr{N}_{B}$ are presented in Table 3. Note that the set of common species of the subnetworks is $\mathscr{S}_{A}\cap\mathscr{S}_{B}=\{X_{4},X_{5},X_{10},X_{12},X_{13},X_{20}\}$ while its set of common complexes $\mathscr{C}_{A}\cap\mathscr{C}_{B}=\{X_{13},X_{20}\}$ .

Table 3: Network numbers of the subnetworks of the deficiency-oriented coarsening of the metabolic insulin signaling network

\mathscr{N}=\mathscr{N}_{A}\cup\mathscr{N}_{B}

Network numbers		$\mathscr{N}_{A}$	$\mathscr{N}_{B}$
Species	$m$	13	13
Complexes	$n$	13	24
Reactant complexes	$n_{r}$	13	13
Reversible reactions	$r_{\text{rev}}$	9	1
Irreversible reactions	$r_{\text{irrev}}$	4	11
Reactions	$r$	22	13
Linkage classes	$\ell$	3	12
Strong linkage classes	$s\ell$	3	23
Terminal strong linkage classes	$t$	3	12
Rank	$s$	10	5
Reactant rank	$q$	12	11
Deficiency	$\delta$	0	7
Reactant deficiency	$\delta_{p}$	1	2

The two subnetworks contrast in further properties besides deficiency: $\mathscr{N}_{A}$ is weakly reversible and nonconservative, while $\mathscr{N}_{B}$ is not weakly reversible and conservative. However, they are both $t$ -minimal and concordant (hence monostationary). The Deficiency Zero Theorem for mass action systems, in fact, guarantees that $\mathscr{N}_{A}$ has a unique complex balanced equilibrium in every stoichiometric compatibility class. Because $\mathscr{N}_{B}$ is not weakly reversible, it has no complex balanced equilibria and there may be stoichiometric compatibility classes without an equilibrium.

Lemma 3 paves the way for the usefulness of the deficiency-oriented decomposition.

Lemma 3.

Let $\mathscr{N}=(\mathscr{S},\mathscr{C},\mathscr{R})$ be a CRN. If $\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}_{2}$ is an independent decomposition of $\mathscr{N}$ with $\mathscr{S}_{i}$ , $S_{i}$ , and $K_{i}$ the set of species, stoichiometric subspace, and kinetics, respectively, of $\mathscr{N}_{i}$ for $i=1,2$ , then:

$(i)$

If $\mathscr{S}_{1}\cap\mathscr{S}_{2}=\varnothing$ , then for $i=1,2$ , the projection maps $p_{i}:\mathbb{R}^{\mathscr{S}}\rightarrow\mathbb{R}^{\mathscr{S}_{i}}$ induce an isomorphism of $\mathbb{R}^{\mathscr{S}}/S\rightarrow\mathbb{R}^{\mathscr{S}_{1}}/S_{1}\times\mathbb{R}^{\mathscr{S}_{2}}/S_{2}$ where $\mathscr{S}_{1}\cup\mathscr{S}_{2}=\mathscr{S}$ ; and
$(ii)$

If $\mathscr{S}_{1}\cap\mathscr{S}_{2}\neq\varnothing$ , let $p_{\mathscr{CS}}:\mathbb{R}^{\mathscr{S}}\rightarrow\mathbb{R}^{\mathscr{CS}}$ be the projection to $\mathscr{CS}:=\mathscr{S}_{1}\cap\mathscr{S}_{2}$ . Then $x\in E_{+}(\mathscr{N}_{i},K_{i})$ is in $E_{+}(\mathscr{N},K)$ if and only if $p_{\mathscr{CS}}(x)\in p_{\mathscr{CS}}(E_{+}(\mathscr{N}_{j},K_{j}))$ , $j\neq i$ .

Proof.

$(i)$ Let $p(x+S)=(p_{1}(x)+S_{1},p_{2}(x)+S_{2})$ . It is well-defined because we have the following: $x-x^{\prime}\in S\Rightarrow p_{i}(x)-p_{i}(x^{\prime})=p_{i}(x-x^{\prime})\in S_{i}$ . Since $\mathscr{S}_{1}\cap\mathscr{S}_{2}=\varnothing$ , $x=p_{1}(x)+p_{2}(x)\in S$ if the projections are in $S_{1}$ and $S_{2}$ , respectively, showing the injectivity. Furthermore, the number of species $m_{1}+m_{2}=m$ (since $\mathscr{S}_{1}\cap\mathscr{S}_{2}=\varnothing$ ) and the rank $s_{1}+s_{2}=s$ (due to direct sum). Hence, the domain and codomain have the same dimensions, showing the isomorphism.

$(ii)$ This follows directly from Feinberg’s Decomposition Theorem that if a decomposition is independent, then $E_{+}(\mathscr{N},K)=E_{+}(\mathscr{N}_{1},K_{1})\cap E_{+}(\mathscr{N}_{2},K_{2})$ . ∎

Remark 2.

1.

Statement $(i)$ of Lemma 3 can be generalized to the case of $\mathscr{S}_{1}\cap\mathscr{S}_{2}\neq\varnothing$ but the isomorphism obtained is not directly relevant to our equilibria analysis. Hence, we have relegated it to Appendix D.
2.

Although statement $(ii)$ of Lemma 3 is an easy consequence of Feinberg’s Decomposition Theorem, it turns out to be useful in cases where at least one of the subnetworks has a well-understood set of positive equilibria, as in the case of the metabolic insulin signaling network. This is shown in Proposition 4 and Example 3.2.

Proposition 4.

Let $(\mathscr{N},K)$ be the metabolic insulin signaling network and $\mathscr{N}=\mathscr{N}_{A}\cup\mathscr{N}_{B}$ its deficiency-oriented decomposition. Then:

$(i)$

The map $\epsilon:E_{+}(\mathscr{N},K)\rightarrow\mathbb{R}^{\mathscr{S}}/S$ given by $\epsilon(x):=x+S$ is injective.
$(ii)$

The map $\epsilon_{A}:E_{+}(\mathscr{N}_{A},K_{A})\rightarrow\mathbb{R}^{\mathscr{S}}/S$ given by $\epsilon_{A}(x):=x+S$ is surjective.
$(iii)$

$Az+S\in\text{Im}(\epsilon)$ if and only if there is $x\in\epsilon_{A}^{-1}(z+S)$ such that $p_{\mathscr{CS}}(x)\in p_{\mathscr{CS}}(E_{+}(\mathscr{N}_{B},K_{B}))$ .

Proof.

$(i)$ is equivalent to the statement that $(\mathscr{N},K)$ is monostationary. $(ii)$ follows from the fact that the map is the composition of $\tilde{\epsilon}_{A}:E_{+}(\mathscr{N},K)\rightarrow\mathbb{R}^{\mathscr{S}}/S_{A}$ (surjective due to the Deficiency Zero Theorem) and the surjective map $x+S_{A}\rightarrow x+S$ . Finally, $(iii)$ is the application of statement $(ii)$ of Lemma 3 to the metabolic insulin signaling network. ∎

Example 1. The following example illustrates the use of statement $(iii)$ of Proposition 4.

First, we express the positive equilibria of $\mathscr{N}_{A}$ in terms of the rate constants:

	$\displaystyle X_{2}=\frac{BCDF-ACEG-BCDF+ACDH}{A^{2}EG+ABDF-A^{2}DH}$
	$\displaystyle X_{3}=\frac{k_{1}}{k_{2}+k_{5}}X_{2}$
	$\displaystyle X_{4}=\frac{CDEF}{ADEG+BD^{2}F-AD^{2}H}$
	$\displaystyle X_{5}=\frac{CDF}{ADH-BDF-AEG}$
	$\displaystyle X_{6}=\frac{k_{13}}{k_{14}}$
	$\displaystyle X_{7}=\frac{k_{9}}{k_{10}+k_{15}}X_{4}$
	$\displaystyle X_{8}=\frac{k_{11}}{k_{12}+k_{16}}X_{5}$
	$\displaystyle X_{12}=\frac{k_{20}}{k_{21}}X_{10}X_{11}$
	$\displaystyle X_{15}=\frac{k_{25}}{k_{24}}X_{13}$
	$\displaystyle X_{20}=\frac{k_{34}}{k_{35}}$

	$\displaystyle A=\frac{k_{1}k_{2}}{k_{2}+k_{5}}-k_{1}-k_{7}$
	$\displaystyle B=k_{6}$
	$\displaystyle C=\frac{k_{8}k_{13}}{k_{14}}$
	$\displaystyle D=\frac{k_{9}k_{10}}{k_{10}+k_{15}}-k_{4}-k_{9}$
	$\displaystyle E=k_{3}$
	$\displaystyle F=\frac{k_{1}k_{5}}{k_{2}+k_{5}}$
	$\displaystyle G=k_{4}$
	$\displaystyle H=\frac{k_{11}k_{12}}{k_{12}+k_{16}}-k_{3}-k_{6}-k_{11}.$

Observe that in subnetwork $\mathscr{N}_{A}$ , the equilibria of $X_{10}$ , $X_{11}$ , and $X_{13}$ are independent of the rate constants. Note also that there are restrictions on the rate constant values since the $X_{i}$ ’s have to be positive. Using the rate constants from [23], we obtain a positive equilibrium for $\mathscr{N}_{A}$ :

	$\displaystyle X_{2}=0.3700$
	$\displaystyle X_{3}=8.8788\times 10^{-4}$
	$\displaystyle X_{4}=3.2844$
	$\displaystyle X_{5}=10.9491$
	$\displaystyle X_{6}=10.0$
	$\displaystyle X_{7}=0.0150$
	$\displaystyle X_{8}=0.0499$
	$\displaystyle X_{10}=7.0929\times 10^{-27}$
	$\displaystyle X_{11}=0.00105$
	$\displaystyle X_{12}=5.2614\times 10^{-19}$
	$\displaystyle X_{13}=0.3100$
	$\displaystyle X_{15}=0.2900$
	$\displaystyle X_{20}=96.$

Next, we substitute the values of the species common to $\mathscr{N}_{A}$ and $\mathscr{N}_{B}$ into the ODEs for $\mathscr{N}_{B}$ and look for a positive solution (this implements statement $(iii)$ of Proposition 4):

	$\displaystyle\dot{X}_{9}=k_{19}X_{10}-k_{17}X_{9}X_{4}-k_{18}X_{9}X_{5}$
	$\displaystyle\dot{X}_{10}=k_{17}X_{9}X_{4}+k_{18}X_{9}X_{5}-k_{19}X_{10}$
	$\displaystyle\dot{X}_{13}=k_{22}X_{12}X_{14}-k_{23}X_{13}$
	$\displaystyle\dot{X}_{14}=k_{23}X_{13}-k_{22}X_{12}X_{14}$
	$\displaystyle\dot{X}_{16}=k_{27}X_{17}-k_{26}X_{13}X_{16}$
	$\displaystyle\dot{X}_{17}=k_{26}X_{13}X_{16}-k_{27}X_{17}$
	$\displaystyle\dot{X}_{18}=k_{29}X_{19}-k_{28}X_{13}X_{18}$
	$\displaystyle\dot{X}_{19}=k_{28}X_{13}X_{18}-k_{29}X_{19}$
	$\displaystyle\dot{X}_{20}=k_{31}X_{21}-k_{30}X_{20}-k_{32}X_{17}X_{20}-k_{33}X_{19}X_{20}$
	$\displaystyle\dot{X}_{21}=k_{30}X_{20}+k_{32}X_{17}X_{20}+k_{33}X_{19}X_{20}-k_{31}X_{21}$

Solving $\dot{X}=0$ , we get the following equilibrium for the other species in $\mathscr{N}_{B}$ :

	$\displaystyle X_{9}=1.5\times 10^{-40}$
	$\displaystyle X_{14}=99.3$
	$\displaystyle X_{17}=(5.1269\times 10^{-9})X_{16}$
	$\displaystyle X_{19}=(5.1269\times 10^{-9})X_{18}$
	$\displaystyle X_{21}=(6.7676\times 10^{-9})X_{16}+(2.7070\times 10^{-8})X_{18}+4.0$

Note that the equilibrium values of $X_{17}$ , $X_{19}$ , and $X_{21}$ depend on any values of $X_{16}$ and $X_{18}$ . If we (randomly) choose $X_{16}=99.9$ and $X_{18}=99.9$ , we get the following:

	$\displaystyle X_{17}=5.1218\times 10^{-7}$
	$\displaystyle X_{19}=5.1218\times 10^{-7}$
	$\displaystyle X_{21}=4.0000.$

Remark 3.

Example 3.2 provides an alternative method for solving for positive equilibria using smaller subnetworks instead of the (more complex) entire network.

4 ACR of Species in Insulin Signaling

ACR denotes the invariance of the concentrations of a species at all positive equilibria of a kinetic system. Shinar and Feinberg introduced the concept in 2010 [24] and, from experimental observations in Escherichia coli subsystems, extracted sufficient (mathematical) conditions for the property. In this section, after a brief review of relevant results on ACR, we present an analysis of the property in the insulin signaling system.

4.1 Absolute Concentration Robustness

A pair of reactant complexes $\{C,C^{\prime}\}$ is a Shinar-Feinberg pair (SF-pair) in species $X$ if their kinetic order vectors (i.e., the corresponding row in the kinetic order matrix) differ only in $X$ . In mass action systems, the stoichiometric coefficients of the complexes play the role of the kinetic order vectors.

An SF-pair is called nonterminal if both complexes are not in terminal strong linkage classes. The pair is said to be linked if both complexes are in a linkage class.

Proposition 5.7 of [18] presents a framework for ACR using SF-pairs. A special case of condition $(ii)$ of the said proposition is any weakly reversible power law system with reactant-determined kinetics, zero deficiency, and a linked SF-pair in $X$ . A proof of ACR in $X$ in this case was already provided in Theorem 6 in the Appendix of [13]. We refer the reader to [18] for the detailed discussion of positive equilibria log-parametrized systems and a proof of the general framework. We apply Proposition 5.7 of [18] only in the special case of the deficiency zero subnetwork of the metabolic insulin signaling network.

We will often call subnetworks with the property $(i)$ or $(ii)$ of Proposition 5.7 of [18] “(low deficiency) ACR building blocks”. A computational approach to the framework was developed in [12]. For the ACR analysis, we implemented this algorithm using Matlab.

Another theorem that we will apply to the ACR analysis of the metabolic insulin signaling network is Theorem 5.5 of [20] regarding the stable ACR for one-dimensional networks which we will refer to as the “Meshkat et al criterion”. As defined by Meshkat et al [20], a kinetic system has stable ACR in a species $X$ (for a set of rate constants) if all the steady states of the kinetic system (for the set of rate constants) are stable.

Remark 4.

The sufficient conditions for ACR in a species in the works of Lao et al [18] and Meshkat et al [20] establish the property for all rate constants for which the system has a positive equilibrium. Meshkat et al have proposed the convention of “vacuous ACR” for the case in which no positive equilibrium exists: in such a scenario, ACR in all species is assumed. This enables the more convenient terminology that the above sufficient conditions provide ACR “for all rate constants” and we adopt this terminology in the following analysis.

4.2 ACR Analysis of Rank 1 Subnetworks of the Finest Independent Decomposition

The finest independent decomposition using the Hernandez-De la Cruz algorithm consists of one subnetwork of rank 6 ( $\mathscr{N}_{1}$ ) and 9 subnetworks of rank 1 ( $\mathscr{N}_{2},\ldots,\mathscr{N}_{10}$ ) (see Table 2). It is natural to first try to apply the ACR criterion of Meshkat et al to the latter set of subnetworks.

Condition $2(b)$ of the Meshkat et al criterion says that all reactions, taken pairwise, must be SF-pairs in species $X$ . Table 4 shows the result of the verification.

Table 4: Verification of condition

2(b)

of the Meshkat et al criterion

Subnetwork	Reactions	Non-SF-pair, e.g.	$2(b)$ satisfied for species $X_{i}$
$\mathscr{N}_{2}$	$R_{13}$ , $R_{14}$	none	Yes for $X_{6}$
$\mathscr{N}_{3}$	$R_{17}$ , $R_{18}$ , $R_{19}$	$R_{17}$ , $R_{18}$	No
$\mathscr{N}_{4}$	$R_{20}$ , $R_{21}$	$R_{20}$ , $R_{21}$	No
$\mathscr{N}_{5}$	$R_{22}$ , $R_{23}$	$R_{22}$ , $R_{23}$	No
$\mathscr{N}_{6}$	$R_{24}$ , $R_{25}$	$R_{24}$ , $R_{25}$	No
$\mathscr{N}_{7}$	$R_{26}$ , $R_{27}$	$R_{26}$ , $R_{27}$	No
$\mathscr{N}_{8}$	$R_{28}$ , $R_{29}$	$R_{28}$ , $R_{29}$	No
$\mathscr{N}_{9}$	$R_{30}$ , $R_{31}$ , $R_{32}$ , $R_{33}$	$R_{30}$ , $R_{31}$	No
$\mathscr{N}_{10}$	$R_{34}$ , $R_{35}$	none	Yes for $X_{20}$

Subnetworks $\mathscr{N}_{2}$ and $\mathscr{N}_{10}$ satisfy condition $2(a)$ of the Meshkat et al criterion as well since they each consist of a reversible pair and are identical to the required embedded network. Hence, $X_{6}$ and $X_{20}$ have ACR for all rate constants in the whole network since the decomposition is independent (see Proposition 4.4 of [18]).

Remark 5.

1.

The Meshkat et al criterion ensures stability of the equilibria only within the subnetworks and not for the whole network so we cannot infer any claim about the property for $X_{6}$ and $X_{20}$ .
2.

In a similar vein, non-ACR for a species in the subnetwork does not necessarily imply non-ACR in the whole network since the set of positive equilibria of the latter is generally smaller than that of the former.

4.3 ACR Analysis of Subnetworks in the Deficiency-Oriented Independent Decomposition

Implementing the algorithm of Fontanil et al [12] and using Proposition 5.7 of [18], we discover Shinar-Feinberg reaction pairs in appropriate low-deficiency subnetworks of coarsenings of the finest independent decomposition in Section 3.1. We find that the system has ACR in 8 species (from the zero-deficiency building blocks): $X_{2}$ , $X_{3}$ , $X_{4}$ , $X_{5}$ , $X_{6}$ , $X_{7}$ , $X_{8}$ , $X_{20}$ . The other species are not identified as having ACR because the reactant complexes of their associated SF-pairs are not nonterminal in the subnetwork generated by the building block. Table 5 provides an overview.

It is particularly significant that the system’s output to glucose energy metabolism ( $X_{20}=\text{intracellular GLUT4}$ ) has ACR. GLUT4 is a key transporter of glucose into neurons. Under healthy conditions (i.e., no insulin resistance), the insulin signaling system works such that GLUT4, which determines the amount of glucose transported into a cell, is held constant. GLUT4, coupled with adequate glucose supply, enables reliable cellular energy production. Keeping the value of GLUT4 is very important for glucose energy metabolism, showing real robustness in the system. Energy processing of neurons works properly due to this robustness.

Table 5: Species with ACR

Species	SF-pair	ACR building block	Deficiency	Comments
$X_{2}$	$R_{1}$ , $R_{13}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{1}$ , $R_{34}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{7}$ , $R_{13}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{7}$ , $R_{34}$	$\mathscr{N}_{A}$	0	linked SF-pair
$X_{3}$	$R_{2}$ , $R_{13}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{2}$ , $R_{34}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{5}$ , $R_{13}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{5}$ , $R_{35}$	$\mathscr{N}_{A}$	0	linked SF-pair
$X_{4}$	$R_{4}$ , $R_{13}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{4}$ , $R_{34}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{9}$ , $R_{13}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{9}$ , $R_{34}$	$\mathscr{N}_{A}$	0	linked SF-pair
$X_{5}$	$R_{3}$ , $R_{13}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{3}$ , $R_{34}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{6}$ , $R_{13}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{6}$ , $R_{34}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{11}$ , $R_{13}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{11}$ , $R_{34}$	$\mathscr{N}_{A}$	0	linked SF-pair
$X_{6}$	$R_{8}$ , $R_{13}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{8}$ , $R_{34}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{13}$ , $R_{14}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{14}$ , $R_{34}$	$\mathscr{N}_{A}$	0	linked SF-pair
$X_{7}$	$R_{10}$ , $R_{13}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{10}$ , $R_{34}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{13}$ , $R_{15}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{15}$ , $R_{34}$	$\mathscr{N}_{A}$	0	linked SF-pair
$X_{8}$	$R_{12}$ , $R_{13}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{12}$ , $R_{34}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{13}$ , $R_{16}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{16}$ , $R_{34}$	$\mathscr{N}_{A}$	0	linked SF-pair
$X_{13}$	$R_{13}$ , $R_{23}$	$\mathscr{N}_{2}\cup\mathscr{N}_{5}$	1
	$R_{23}$ , $R_{34}$	$\mathscr{N}_{5}\cup\mathscr{N}_{10}$	1
$X_{14}$	$R_{21}$ , $R_{22}$	$\mathscr{N}_{4}\cup\mathscr{N}_{5}$	1
$X_{16}$	$R_{25}$ , $R_{26}$	$\mathscr{N}_{6}\cup\mathscr{N}_{7}$	1
$X_{17}$	$R_{13}$ , $R_{27}$	$\mathscr{N}_{2}\cup\mathscr{N}_{7}$	1
	$R_{27}$ , $R_{34}$	$\mathscr{N}_{7}\cup\mathscr{N}_{10}$	1
$X_{18}$	$R_{25}$ , $R_{28}$	$\mathscr{N}_{6}\cup\mathscr{N}_{8}$	1
$X_{19}$	$R_{13}$ , $R_{29}$	$\mathscr{N}_{2}\cup\mathscr{N}_{8}$	1
	$R_{29}$ , $R_{34}$	$\mathscr{N}_{8}\cup\mathscr{N}_{10}$	1
$X_{20}$	$R_{13}$ , $R_{35}$	$\mathscr{N}_{A}$	0	linked SF-pair
	$R_{34}$ , $R_{35}$	$\mathscr{N}_{A}$	0	linked SF-pair

Remark 6.

The ACR in 8 species of the network is inferred from ACR in the deficiency zero subnetwork $\mathscr{N}_{A}$ . A necessary condition for this occurrence is the nonconservativity of $\mathscr{N}$ . If $\mathscr{N}$ were conservative, i.e., there was a positive vector in $S^{\perp}$ , then $S^{\perp}=(S_{A}+S_{B})^{\perp}=S_{A}^{\perp}\cap S_{B}^{\perp}$ . Hence, $\mathscr{N}_{A}$ would be a conservative deficiency zero mass action network and have no ACR in any species (Theorem 9.7.1 of [8]).

4.4 ACR Analysis of the Metabolic Insulin Signaling Network for a Set of Rate Constants

Table 5 shows that there are SF-pairs in deficiency one subnetworks in further independent coarsenings of the metabolic insulin signaling network. The failure of the three sufficient conditions for ACR for all rate constants (Meshkat et al criterion) led us to suspect that the remaining species (i.e., $X_{13}$ , $X_{14}$ , $X_{16}$ , $X_{17}$ , $X_{18}$ , and $X_{19}$ ) do not have this property. The test is done using the set of rate constants available from [23].

Example 3.2 in Section 3.2 presents a positive equilibrium of the metabolic insulin signaling network. It is evident from the example that positive equilibria for the metabolic insulin signaling network have the same value for $X_{2},\ldots,X_{8}$ and $X_{20}$ . For the other 12 species, their equilibrium value is dependent on the following variables: $X_{10}$ , $X_{11}$ , $X_{13}$ , $X_{16}$ , and $X_{18}$ . This suggests that there are infinitely many possible positive equilibria.

Our computations verify that, for the set of rate constants used above, species $X_{2},\ldots,X_{8}$ and $X_{20}$ have stable ACR while the remaining species do not have ACR (since the 5 species identified as “independent variables” can take various values and the remaining species will vary according to the variation of the variables).

These computations indeed confirm that species $X_{13}$ , $X_{14}$ , $X_{16}$ , $X_{17}$ , $X_{18}$ , and $X_{19}$ do not have ACR for all rate constants.

Remark 7.

The number of species exhibiting ACR (which we call “ACR species” for short) is an inverse measure of the variation in the equilibria composition: the more ACR species there are, the less the variation. An extreme case is ACR in all species, which is equivalent to the system having a unique equilibrium (in the entire species space). With 8 ACR species among 20, the metabolic insulin signaling network has a fairly high variation in equilibria composition. The deficiency-oriented decomposition reveals the source of this variation. The deficiency zero subnetwork $\mathscr{N}_{A}$ contains all the ACR species. Furthermore, since the insulin signaling system is a log-parametrized system, it follows from the results of Lao et al [18] that the number of ACR species (among the total 13 species in the subnetwork) is bounded by the subnetwork’s rank ( $s_{A}=10$ ). This shows that the high variation in equilibria composition is caused primarily by the lack of ACR species among the 13 species in the deficiency 7 subnetwork $\mathscr{N}_{B}$ .

5 Summary and Outlook

The insulin signaling system is an important metabolic system. In this study, we derived a CRN of the Sedaghat et al insulin signaling model with 20 species, 35 complexes, and 35 reactions. We have shown that it is a nonconservative, non-weakly reversible, and high deficiency ( $\delta=7$ ) system. The positive dependence of the reaction network ensures the existence of rate constants under which the mass action system has positive equilibria (Proposition 1). Additionally, the network’s $t$ -minimality implies that the kinetic and stoichiometric subspaces coincide, which is necessary for the existence of nondegenerate equilibria. Moreover, the network is concordant, which implies that the system’s species formation rate function, when restricted to any stoichiometric compatibility class, is injective. It follows, then, that the kinetic system is monostationary, i.e., there is at most one positive equilibrium in each stoichiometric compatibility class.

We obtained the finest independent decomposition of the metabolic insulin signaling network $\mathscr{N}=\{R_{1},\ldots,R_{35}\}$ consisting of 10 subnetworks which we have shown to be bi-independent (Proposition 2). CRNT analysis using a coarsening of the decomposition revealed three subnetworks of the metabolic insulin signaling network which are not only functionally but also structurally important. Upon considering a deficiency-oriented coarsening of the finest decomposition, we have shown how a binary decomposition can be viewed in relation to a network’s set of positive equilibria (Lemma 3). We have also developed a method of determining positive equilibria of the metabolic insulin signaling network using its deficiency-oriented coarsening (Proposition 4). This provides an alternative method for solving positive equilibria analytically using smaller subnetworks instead of the (more complex) entire network.

For the ACR analysis, we have shown that subnetworks $\mathscr{N}_{2}$ and $\mathscr{N}_{10}$ satisfy the Meshkat et al criterion, but species $X_{6}$ and $X_{20}$ (which come from the said subnetworks) have ACR for all rate constants only in their respective subnetwork. This implies that the stability of the equilibria is only within the subnetworks and not for the whole network. Upon implementing the algorithm of Fontanil et al [12], however, we found that the system had ACR in 8 species. This restricts the variability of the positive equilibria and also suggests that this “structural source of robustness” may be an important factor in the system’s overall robustness. However, overall there is still a high variation in equilibria composition due to the lack of ACR species in the deficiency 7 subnetwork of the network’s deficiency-oriented coarsening. For the rate constants used in the study of Sedaghat et al [23], we have verified that species $X_{2},\ldots,X_{8}$ , and $X_{20}$ indeed have stable ACR. Interestingly, $X_{20}$ , i.e., the insulin-regulated glucose transporter GLUT4, plays an important role in glucose energy metabolism.

Our analysis of the Sedaghat et al model is the first part of a two-step research effort on metabolic insulin signaling. In their paper, Sedaghat et al constructed the model from data of healthy cells. As a second step, we have started a reaction network analysis of a model of metabolic insulin signaling by Brännmark et al [3] based on cell data from type 2 diabetes patients, i.e., cells with insulin resistance. Preliminary results already indicate significant differences to our main results on the Sedaghat et al model [19].

References

[1] Arceo C, Jose E, Lao A, Mendoza E (2017) Reaction networks and kinetics of biochemical systems. Math Biosci 283:13–29. https://doi.org/10.1016/j.mbs.2016.10.004
[2] Arnold S, Arvanitakis Z, Macauley-Rambach S, Koenig A, Wang H, Ahima R, Craft S, Gandy S, Buettner C, Stoeckel L, Holtzman D, Nathan D (2018) Brain insulin resistance in type 2 diabetes and Alzheimer disease: concepts and conundrums. Nat Rev Neurol 14(3):168–181. https://doi.org/10.1038/nrneurol.2017.185
[3] Brännmark C, Nyman E, Fagerholm S, Bergenholm L, Ekstrand E, Cedersund G, Strålfors P (2013) Insulin signaling in type 2 diabetes: experimental and modeling analyses reveal mechanisms of insulin resistance in human adipocytes. J Biol Chem 288(14):9867–9880.
https://doi.org/10.1074/jbc.m112.432062
[4] Chellaboina V, Bhat S, Haddad W, Bernstein D (2009) Modeling and analysis of mass-action kinetics. IEEE Control Syst 29(4):60–78.
https://doi.org/10.1109/MCS.2009.932926
[5] Dexter J, Dasgupta T, Gunawardena, J (2015) Invariants reveal multiple forms of robustness in bifunctional enzyme systems. Integr Biol 7(8):883–894.
https://doi.org/10.1039/c5ib00009b
[6] Fariñas H, Mendoza E, Lao A (2021) Chemical reaction network decompositions and realizations of S-systems. Philipp Sci Lett 14(1):147–157
[7] Feinberg M (1987) Chemical reaction network structure and the stability of complex isothermal reactors: I. The deficiency zero and deficiency one theorems. Chem Eng Sci 42(10):2229–2268. https://doi.org/10.1016/0009-2509(87)80099-4
[8] Feinberg M (2019) Foundations of chemical reaction network theory. Springer, Switzerland,
https://doi.org/10.1007/978-3-030-03858-8
[9] Feinberg M, Horn J (1977) Chemical mechanism structure and the coincidence of the stoichiometric and kinetic subspaces. Arch Ration Mech Anal 66(1):83–97.
https://doi.org/10.1007/bf00250853
[10] Feinberg M, Ellison P, Ji H, Knight D (2018) The Chemical Reaction Network Toolbox Version 2.35. https://doi.org/10.5281/zenodo.5149266
[11] Feliu E, Wiuf C (2012) Preclusion of switch behavior in networks with mass-action kinetics. Appl Math Comput 219(14):1449–1467. https://doi.org/10.1016/j.amc.2012.07.048
[12] Fontanil L, Mendoza E, Fortun N (2021) A computational approach to concentration robustness in power law kinetic systems of Shinar-Feinberg type. MATCH Commun Math Comput Chem 86(3):489–516
[13] Fortun N, Mendoza E (2021) Absolute concentration robustness in power law kinetic systems. MATCH Commun Math Comput Chem 85(3):669–691
[14] Hernandez B, De la Cruz R (2021) Independent decompositions of chemical reaction networks. Bull Math Biol 83(76):1–23. https://doi.org/10.1007/s11538-021-00906-3
[15] Hernandez B, Mendoza E (2021) Positive equilibria of Hill-type kinetic systems. J Math Chem 59:840–870. https://doi.org/10.1007/s10910-021-01230-w
[16] Hernandez B, Amistas D, De la Cruz R, Fontanil L, de los Reyes V A, Mendoza E (2022) Independent, incidence independent and weakly reversible decompositions of chemical reaction networks. MATCH Commun Math Comput Chem 87(2):367–396.
https://doi.org/10.46793/match.87-2.367H
[17] Knight D, Shinar G, Feinberg M (2015) Sharper graph-theoretical conditions for the stabilization of complex reaction networks. Math Biosci 262:10–27.
https://doi.org/10.1016/j.mbs.2015.01.002
[18] Lao A, Lubenia P, Magpantay D, Mendoza E (2022) Concentration robustness in LP kinetic systems. MATCH Commun Math Comput Chem 88(1):29–66.
https://doi.org/10.46793/match.88-1.029L
[19] Lubenia P, Mendoza E, Lao A (2022) Comparison of kinetic realizations of insulin signaling systems (in preparation)
[20] Meshkat N, Shiu A, Torres A (2021) Absolute concentration robustness in networks with low-dimensional stoichiometric subspace. Vietnam J Math
https://doi.org/10.1007/s10013-021-00524-5
[21] Nazareno A, Eclarin R, Mendoza E, Lao A (2019) Linear conjugacy of chemical kinetic systems. Math Biosci Eng 16(6):8322–8355. https://doi.org/10.3934/mbe.2019421
[22] Neth B, Craft S (2017) Insulin resistance and Alzheimer’s disease: Bioenergetic linkages. Front Aging Neurosci 9(345):1–20. https://doi.org/10.3389/fnagi.2017.00345
[23] Sedaghat A, Sherman A, Quon M (2002) A mathematical model of metabolic insulin signaling pathways. Am J Physiol Endocrinol Metab 283(5):E1084–E1101.
https://doi.org/10.1152/ajpendo.00571.2001
[24] Shinar G, Feinberg M (2010) Structural sources of robustness in biochemical reaction networks. Science 327(5971):1389–1391. https://doi.org/10.1126/science.1183372
[25] Shinar G, Feinberg M (2011) Design principles for robust biochemical reaction networks: What works, what cannot work, and what might almost work. Math Biosci 231(1):39–48.
https://doi.org/10.1016/j.mbs.2011.02.012
[26] Shinar G, Feinberg M (2012) Concordant chemical reaction networks. Math Biosci 240(2):92–113. https://doi.org/10.1016/j.mbs.2012.05.004
[27] Shinar G, Feinberg M (2012) Concordant chemical reaction networks and the species-reaction graph. Math Biosci 241(1):1–23. https://doi.org/10.1016/j.mbs.2012.08.002
[28] Villasis O, Tan D, Mendoza E, Lao A (2022) Systems biology approach of understanding insulin resistance: Linkage between type 2 diabetes & Alzheimer’s disease.
bioRxiv 2022.04.13.463890

Appendix A Notations and Definition of Terms

In this section, we lay the foundation of Chemical Reaction Network Theory (CRNT) by discussing the definition of terms used in the paper. After discussing the fundamentals of chemical reaction networks and kinetic systems, we review important terminologies related to decomposition theory.

A.1 Chemical Reaction Networks

A chemical reaction network (CRN) $\mathscr{N}$ is a triple $(\mathscr{S},\mathscr{C},\mathscr{R})$ of nonempty finite sets $\mathscr{S}$ , $\mathscr{C}$ , and $\mathscr{R}$ of $m$ species, $n$ complexes, and $r$ reactions, respectively, where $\mathscr{C}\subseteq\mathbb{R}_{\geq 0}^{\mathscr{S}}$ and $\mathscr{R}\subset\mathscr{C}\times\mathscr{C}$ satisfying the following:

( $i$ )

$(C,C)\notin\mathscr{R}$ for any $C\in\mathscr{C}$ ; and
( $ii$ )

For each $C_{i}\in\mathscr{C}$ , there exists $C_{j}\in\mathscr{C}$ such that $(C_{i},C_{j})\in\mathscr{R}$ or $(C_{j},C_{i})\in\mathscr{R}$ .

Given a set $I$ , $\mathbb{R}^{I}$ refers to the vector space of real-valued functions with domain $I$ . If $x$ is a vector in $\mathbb{R}^{I}$ , we use the symbol $x_{i}$ to denote the number that $x$ assigns to $i\in I$ . In the context of a CRN, $\mathbb{R}^{\mathscr{S}}$ , $\mathbb{R}^{\mathscr{C}}$ , and $\mathbb{R}^{\mathscr{R}}$ are referred to as the species space, complex space, and reaction space, respectively.

In a CRN, we denote the species as $X_{1},\ldots,X_{m}$ . This way, $X_{i}$ can be identified with the vector in $\mathbb{R}^{m}$ with 1 in the $i$ th coordinate and zero elsewhere. We denote the reactions as $R_{1},\ldots,R_{r}$ . We denote the complexes as $C_{1},\ldots,C_{n}$ where the manner in which the complexes are numbered play no essential role. A complex $C_{i}\in\mathscr{C}$ is given as $\displaystyle C_{i}=\sum_{j=1}^{m}c_{ij}X_{j}$ or as the vector $(c_{i1},\ldots,c_{im})\in\mathbb{R}^{m}$ . The coefficient $c_{ij}$ is called the stoichiometric coefficient of species $X_{j}$ in complex $C_{i}$ . Stoichiometric coefficients are all nonnegative numbers. We define the zero complex as the zero vector in $\mathbb{R}^{m}$ . A reaction $0\rightarrow X$ is called an inflow reaction while a reaction $X\rightarrow 0$ is called an outflow reaction. The ordered pair $(C_{i},C_{j})$ corresponds to the familiar notation $C_{i}\rightarrow C_{j}$ which indicates the reaction where complex $C_{i}$ reacts to complex $C_{j}$ . We call $C_{i}$ the reactant complex and $C_{j}$ the product complex. We denote the number of reactant complexes as $n_{r}$ . A reaction $C_{i}\rightarrow C_{j}$ is called reversible if it is accompanied by its reverse reaction $C_{j}\rightarrow C_{i}$ . Otherwise, it is called irreversible.

Example A.1. Consider the reaction

2X_{1}+X_{2}\rightarrow 2X_{3}.

$X_{1}$ , $X_{2}$ , and $X_{3}$ are the species, the reactant complex is $2X_{1}+X_{2}$ , and $2X_{3}$ is the product complex. The stoichiometric coefficients are 2, 1, and 2 for $X_{1}$ , $X_{2}$ , and $X_{3}$ , respectively.

Let $\mathscr{N}=(\mathscr{S},\mathscr{C},\mathscr{R})$ be a CRN. For each reaction $C_{i}\rightarrow C_{j}\in\mathscr{R}$ , we associate the reaction vector $C_{j}-C_{i}\in\mathbb{R}^{m}$ . The linear subspace of $\mathbb{R}^{m}$ spanned by the reaction vectors is called the stoichiometric subspace of $\mathscr{N}$ , defined as $S=\text{span}\{C_{j}-C_{i}\in\mathbb{R}^{m}\mid C_{i}\rightarrow C_{j}\in\mathscr{R}\}$ . The rank of $\mathscr{N}$ is given by $s=\text{dim}(S)$ , i.e., the rank of the network is the rank of its set of reaction vectors. In this paper, we sometimes use the notation $\mathscr{N}=\{R_{1},\ldots,R_{r}\}$ where we loosely use the notation $R_{i}$ to refer to either reaction $i$ or its corresponding reaction vectors.

Two vectors $x^{*},x^{**}\in\mathbb{R}^{m}$ are said to be stoichiometrically compatible if $x^{*}-x^{**}$ is an element of the stoichiometric subspace $S$ . Stoichiometric compatibility is an equivalence relation that induces a partition of $\mathbb{R}_{\geq 0}^{\mathscr{S}}$ or $\mathbb{R}_{>0}^{\mathscr{S}}$ into equivalence classes called the stoichiometric compatibility classes or positive stoichiometric compatibility classes, respectively, of the network. In particular, the stoichiometric compatibility class containing $x\in\mathbb{R}_{\geq 0}^{\mathscr{S}}$ is the set $(x+S)\cap\mathbb{R}_{\geq 0}^{\mathscr{S}}$ where $x+S$ is the left coset of $S$ containing $x$ . Similarly, the positive stoichiometric compatibility class containing $x\in\mathbb{R}_{>0}^{\mathscr{S}}$ is the set $(x+S)\cap\mathbb{R}_{>0}^{\mathscr{S}}$ .

The molecularity matrix $Y$ is an $m\times n$ matrix whose entry $Y_{ij}$ is the stoichiometric coefficient of species $X_{i}$ in complex $C_{j}$ . The incidence matrix $I_{a}$ is an $n\times r$ matrix whose entry $(I_{a})_{ij}$ is defined as follows:

(I_{a})_{ij}=\left\{\begin{array}[]{rl}-1&\text{if $C_{i}$ is the reactant complex of reaction $R_{j}$}\\ 1&\text{if $C_{i}$ is the product complex of reaction $R_{j}$}\\ 0&\text{otherwise}\\ \end{array}\right..

The stoichiometric matrix $N$ is the $m\times r$ matrix given by $N=YI_{a}$ . The columns of $N$ are the reaction vectors of the system. From the definition of stoichiometric subspace, we can see that $S$ is the image of $N$ , written as $S=\text{Im}(N)$ . Observe that $s=\text{dim}(S)=\text{dim}(\text{Im}(N))=\text{rank}(N)$ .

Example A.2. Consider the following CRN:

	$\displaystyle R_{1}:2X_{1}\rightarrow X_{3}$
	$\displaystyle R_{2}:X_{2}+X_{3}\rightarrow X_{3}$
	$\displaystyle R_{3}:X_{3}\rightarrow X_{2}+X_{3}$
	$\displaystyle R_{4}:3X_{4}\rightarrow X_{2}+X_{3}$
	$\displaystyle R_{5}:2X_{1}\rightarrow 3X_{4}.$

The set of species and complexes are $\mathscr{S}=\{X_{1},X_{2},X_{3},X_{4}\}$ and $\mathscr{C}=\{2X_{1},X_{2}+X_{3},X_{3},3X_{4}\}$ , respectively. Thus, there are $m=4$ species, $n=4$ complexes, $n_{r}=4$ reactant complexes, and $r=5$ reactions. The network’s molecularity matrix, incidence matrix, and stoichiometric matrix are as follows:

Y=\blockarray{ccccc}2X_{1}&X_{2}+X_{3}X_{3}3X_{4}\\ \block{[cccc]l}2000{\bigstrut[t]}X_{1}\\ 0100X_{2}\\ 0110X_{3}\\ 0003{\bigstrut[b]}X_{4}\\

I_{a}=\blockarray{rrrrrr}R_{1}&R_{2}R_{3}R_{4}R_{5}{\color[rgb]{1,1,1}1}\\ \block{[rrrrr]l}{\color[rgb]{1,1,1}1}-1000-1{\color[rgb]{1,1,1}1}{\bigstrut[t]}2X_{1}\\ 0-1110{\color[rgb]{1,1,1}1}X_{2}+X_{3}\\ 11-100{\color[rgb]{1,1,1}1}X_{3}\\ 000-11{\color[rgb]{1,1,1}1}{\bigstrut[b]}3X_{4}\\

N=YI_{a}=\left[\begin{array}[]{rrrrr}-2&0&0&0&-2\\ 0&-1&{\color[rgb]{1,1,1}-}1&1&0\\ 1&0&0&1&0\\ 0&0&0&-3&3\\ \end{array}\right].

The network has rank $s=\text{rank}(N)=3$ .

CRNs can be viewed as directed graphs where the complexes are represented by vertices and the reactions by edges. The linkage classes of a CRN are the subnetworks of its reaction graph where for any complexes $C_{i}$ and $C_{j}$ of the subnetwork, there is a path between them. The number of linkage classes is denoted by $\ell$ . The linkage class is said to be a strong linkage class if there is a directed path from $C_{i}$ to $C_{j}$ , and vice versa, for any complexes $C_{i}$ and $C_{j}$ of the subnetwork. The number of strong linkage classes is denoted by $s\ell$ . Moreover, terminal strong linkage classes, the number of which is denoted as $t$ , are the maximal strongly connected subnetworks where there are no edges (reactions) from a complex in the subgraph to a complex outside the subnetwork. Complexes belonging to terminal strong linkage classes are called terminal; otherwise, they are called nonterminal.

In Example A.1, the number of linkage classes is $\ell=1$ : $\{2X_{1},X_{3},X_{2}+X_{3},3X_{4}\}$ ; the number of strong linkage classes is $s\ell=3$ : $\{X_{3},X_{2}+X_{3}\},\{2X_{1}\},\{3X_{4}\}$ ; and the number of terminal strong linkage classes is $t=1$ : $\{X_{3},X_{2}+X_{3}\}$ . $X_{3}$ and $X_{2}+X_{3}$ are terminal complexes while $2X_{1}$ and $3X_{4}$ are nonterminal complexes.

A CRN is called weakly reversible if $s\ell=\ell$ , $t$ -minimal if $t=\ell$ , point terminal if $t=n-n_{r}$ , and cycle terminal if $n-n_{r}=0$ . The deficiency of a CRN is given by $\delta=n-\ell-s$ .

For a CRN $\mathscr{N}$ , the linear subspace of $\mathbb{R}^{m}$ generated by the reactant complexes is called the reactant subspace of $\mathscr{N}$ , defined as $R=\text{span}\{C_{i}\in\mathbb{R}^{m}\mid C_{i}\rightarrow C_{j}\in\mathscr{R}\}$ . The reactant rank of $\mathscr{N}$ is given by $q=\text{dim}(R)$ , i.e., the reactant rank of the network is the rank of its set of complexes. The reactant deficiency of $\mathscr{N}$ is given by $\delta_{p}=n_{r}-q$ .

To make sense of the reactant subspace $R$ , write the incidence matrix as $I_{a}=I_{a}^{+}-I_{a}^{-}$ where $I_{a}^{+}$ consists only of the 0’s and 1’s in $I_{a}$ while $I_{a}^{-}$ contains only the 0’s and absolute values of the $-1$ ’s. We form the reactant matrix $N^{-}$ (size $m\times r$ ) given by $N^{-}=YI_{a}^{-}$ . The columns of $N^{-}$ contains the reactant complexes of the system. From the definition of reactant subspace, we can see that $R$ is the image of $N^{-}$ , written as $R=\text{Im}(N^{-})$ . Observe that $q=\text{dim}(R)=\text{dim}(\text{Im}(N^{-}))=\text{rank}(N^{-})$ .

The incidence matrix of the network in Example A.1 can be written as

	$\displaystyle I_{a}$	$\displaystyle=I_{a}^{+}-I_{a}^{-}$
	$\displaystyle\left[\begin{array}[]{rrrrr}-1&0&0&0&-1\\ 0&-1&1&1&0\\ 1&1&-1&0&0\\ 0&0&0&-1&1\\ \end{array}\right]$	$\displaystyle=\left[\begin{array}[]{rrrrr}0&0&0&0&0\\ 0&0&1&1&0\\ 1&1&0&0&0\\ 0&0&0&0&1\\ \end{array}\right]-\left[\begin{array}[]{rrrrr}1&0&0&0&1\\ 0&1&0&0&0\\ 0&0&1&0&0\\ 0&0&0&1&0\\ \end{array}\right]$

allowing us to form the reactant matrix

N^{-}=YI_{a}^{-}=\left[\begin{array}[]{rrrrr}2&0&0&0&2\\ 0&1&0&0&0\\ 0&1&1&0&0\\ 0&0&0&3&0\\ \end{array}\right].

The network has reactant rank $q=\text{rank}(N^{-})=4$ and reactant deficiency $\delta_{p}=n_{r}-q=4-4=0$ .

A.2 Chemical Kinetic Systems

A kinetics $K$ for a CRN $\mathscr{N}=(\mathscr{S},\mathscr{C},\mathscr{R})$ is an assignment to each reaction $C_{i}\rightarrow C_{j}\in\mathscr{R}$ of a rate function $K_{C_{i}\rightarrow C_{j}}:\mathbb{R}_{\geq 0}^{\mathscr{S}}\rightarrow\mathbb{R}_{\geq 0}$ such that

K_{C_{i}\rightarrow C_{j}}(X)>0\text{ if and only if }\text{supp}(C_{i})\subset\text{supp}(X).

The system $(\mathscr{N},K)$ is called a chemical kinetic system (CKS). The support of complex $C_{i}\in\mathscr{C}$ is $\text{supp}(C_{i})=\{X_{j}\in\mathscr{S}\mid c_{ij}\neq 0\}$ , i.e, it is the set of all species that have nonzero stoichiometric coefficients in complex $C_{i}$ .

A kinetics gives rise to two closely related objects: the species formation rate function and the associated ordinary differential equation system.

The species formation rate function (SFRF) of a CKS is given by

f(X)=\sum_{C_{i}\rightarrow C_{j}}K_{C_{i}\rightarrow C_{j}}(X)(C_{j}-C_{i})

where $X$ is the vector of species in $\mathscr{S}$ and $K_{C_{i}\rightarrow C_{j}}$ is the rate function assigned to reaction $C_{i}\rightarrow C_{j}\in\mathscr{R}$ . The SFRF is simply the summation of the reaction vectors for the network, each multiplied by the corresponding rate function. The kinetic subspace $\mathcal{K}$ for a CKS is the linear subspace of $\mathbb{R}^{\mathscr{S}}$ defined by $\mathcal{K}=\text{span}\{\text{Im}(f)\}.$ Note that the SFRF can be written as $f(X)=NK(X)$ where $K$ the vector of rate functions. The equation $\dot{X}=f(X)$ is the ordinary differential equation (ODE) system or dynamical system of the CKS.

The ODE system of the CRN in Example A.1 can be written as

\dot{X}=\left[\begin{array}[]{c}\dot{X}_{1}\\ \dot{X}_{2}\\ \dot{X}_{3}\\ \dot{X}_{4}\\ \end{array}\right]=\left[\begin{array}[]{rrrrr}-2&0&0&0&-2\\ 0&-1&{\color[rgb]{1,1,1}-}1&1&0\\ 1&0&0&1&0\\ 0&0&0&-3&3\\ \end{array}\right]\left[\begin{array}[]{l}k_{1}X_{1}^{f_{11}}\\ k_{2}X_{2}^{f_{22}}X_{3}^{f_{23}}\\ k_{3}X_{3}^{f_{33}}\\ k_{4}X_{4}^{f_{44}}\\ k_{5}X_{1}^{f_{51}}\\ \end{array}\right]=NK(X).

A zero of the SFRF is called an equilibrium or a steady state of the system. If $f$ is differentiable, an equilibrium $X^{*}$ is called degenerate if $\text{Ker}(J_{X^{*}}(f))\cap S\neq\{0\}$ where $J_{X^{*}}(f)$ is the Jacobian of $f$ evaluated at $X^{*}$ and Ker is the kernel function; otherwise, the equilibrium is said to be nondegenerate.

A vector $X\in\mathbb{R}_{>0}^{\mathscr{S}}$ is called complex balanced if $K(X)$ is contained in $\text{Ker}(I_{a})$ where $I_{a}$ is the incidence matrix. Furthermore, if $X$ is a positive equilibrium, then we call it a complex balanced equilibrium. A CKS is called complex balanced if it has a complex balanced equilibrium.

The reaction vectors of a CRN are positively dependent if, for each reaction $C_{i}\rightarrow C_{j}\in\mathscr{R}$ , there exists a positive number $\alpha_{C_{i}\rightarrow C_{j}}$ such that

\sum_{C_{i}\rightarrow C_{j}}\alpha_{C_{i}\rightarrow C_{j}}(C_{j}-C_{i})=0.

A CRN with positively dependent reaction vectors is said to be positive dependent. Shinar and Feinberg [26] showed that a CKS can admit a positive equilibrium only if its reaction vectors are positively dependent. The set of positive equilibria of a CKS is given by

E_{+}(\mathscr{N},K)=\{X\in\mathbb{R}_{>0}^{\mathscr{S}}\mid f(X)=0\}.

A CRN is said to admit multiple (positive) equilibria if there exist positive rate constants such that the ODE system admits more than one stoichiometrically compatible equilibria. Analogously, the set of complex balanced equilibria of a CKS $(\mathscr{N},K)$ is given by

Z_{+}(\mathscr{N},K)=\{X\in\mathbb{R}_{>0}^{\mathscr{S}}\mid I_{a}K(X)=0\}\subseteq E_{+}(\mathscr{N},K).

Let $F$ be an $r\times m$ matrix of real numbers. Define $X^{F}$ by $\displaystyle(X^{F})_{i}=\prod_{j=1}^{m}X_{j}^{f_{ij}}$ for $i=1,\ldots,r$ . A power law kinetics (PLK) assigns to each $i$ th reaction a function

K_{i}(X)=k_{i}(X^{F})_{i}

with rate constant $k_{i}>0$ and kinetic order $f_{ij}\in\mathbb{R}$ . The vector $k\in\mathbb{R}^{r}$ is called the rate vector and the matrix $F$ is called the kinetic order matrix. We refer to a CRN with PLK as a power law system. The PLK becomes the well-known mass action kinetics (MAK) if the kinetic order matrix consists of stoichiometric coefficients of the reactants. We refer to a CRN with MAK as a mass action system.

In the ODE system of Example A.1, we assumed PLK so that the kinetic order matrix is

F=\left[\begin{array}[]{cccc}f_{11}&0&0&0\\ 0&f_{22}&f_{23}&0\\ 0&0&f_{33}&0\\ 0&0&0&f_{44}\\ f_{51}&0&0&0\\ \end{array}\right]

where $f_{ij}\in\mathbb{R}$ . If we assume MAK, the kinetic order matrix is

F=\left[\begin{array}[]{cccc}2&0&0&0\\ 0&1&1&0\\ 0&0&1&0\\ 0&0&0&3\\ 2&0&0&0\\ \end{array}\right].

The reactions $R_{i},R_{j}\in\mathscr{R}$ are called branching reactions if they have the same reactant complex. One way to check if we have identified all branching reactions is through the formula

r-n_{r}=\sum_{C_{i}}(|R_{C_{i}}|-1)

where $C_{i}$ is the reactant complex in $C_{i}\rightarrow C_{j}\in\mathscr{R}$ , $R_{C_{i}}$ is the set of branching reactions of $C_{i}$ , and $|R_{C_{i}}|$ is the cardinality of $R_{C_{i}}$ . $r-n_{r}=0$ if and only if all reactant complexes are nonbranching. A CRN is called branching if $r>n_{r}$ .

We can classify a power law system based on the kinetic orders assigned to its branching reactions. A power law system has reactant-determined kinetics (of type PL-RDK) if, for any two branching reactions $R_{i},R_{j}\in\mathscr{R}$ , their corresponding rows of kinetic orders in $F$ are identical, i.e., $f_{ik}=f_{jk}$ for $k=1,\ldots,m$ . Otherwise, a power law system has non-reactant-determined kinetics (of type PL-NDK). From Proposition 12 in Arceo et al [1], a nonbranching power law system is PL-RDK.

Note that the stoichiometric subspace $S$ is just the set of all linear combinations of the reaction vectors, i.e., the set of all vectors in $\mathbb{R}^{\mathscr{S}}$ can be written in the form

\sum_{C_{i}\rightarrow C_{j}}\alpha_{C_{i}\rightarrow C_{j}}(C_{j}-C_{i}).

Let $L:\mathbb{R}^{\mathscr{R}}\rightarrow S$ be the linear map defined by

L(\alpha)=\sum_{C_{i}\rightarrow C_{j}}\alpha_{C_{i}\rightarrow C_{j}}(C_{j}-C_{i}).

$\text{Ker}(L)$ is the set of all vectors $\alpha\in\mathbb{R}^{\mathscr{R}}$ such that $L(\alpha)=0$ .

We say that a CRN is concordant if there do not exist an $\alpha\in\text{Ker}(L)$ and a nonzero $\sigma\in S$ having the following properties:

( $i$ )

For each $C_{i}\rightarrow C_{j}\in\mathscr{R}$ such that $\alpha_{C_{i}\rightarrow C_{j}}\neq 0$ , $\text{supp}(C_{i})$ contains a species $X$ for which $\text{sgn}(\sigma_{X})=\text{sgn}(\alpha_{C_{i}\rightarrow C_{j}})$ where $\sigma_{X}$ denotes the term in $\sigma$ involving the species $X$ and $\text{sgn}(\cdot)$ is the signum function.
( $ii$ )

For each $C_{i}\rightarrow C_{j}\in\mathscr{R}$ such that $\alpha_{C_{i}\rightarrow C_{j}}=0$ , either $\sigma_{X}=0$ for all $X\in\text{supp}(C_{i})$ , or else $\text{supp}(C_{i})$ contains species $X$ and $X^{\prime}$ for which $\text{sgn}(\sigma_{X})=-\text{sgn}(\sigma_{X^{\prime}})$ , but not zero.

A network that is not concordant is discordant.

A CKS is injective if, for each pair of distinct stoichiometrically compatible vectors $X^{*},X^{**}\in\mathbb{R}_{\geq 0}^{\mathscr{S}}$ , at least one of which is positive,

\sum_{C_{i}\rightarrow C_{j}}K_{C_{i}\rightarrow C_{j}}(X^{**})(C_{j}-C_{i})\neq\sum_{C_{i}\rightarrow C_{j}}K_{C_{i}\rightarrow C_{j}}(X^{*})(C_{j}-C_{i}).

Clearly, an injective kinetic system cannot admit two distinct stoichiometrically compatible equilibria, at least one of which is positive.

A kinetics for a CRN is weakly monotonic if, for each pair of vectors $X^{*},X^{**}\in\mathbb{R}_{\geq 0}^{\mathscr{S}}$ , the following implications hold for each reaction $C_{i}\rightarrow C_{j}\in\mathscr{R}$ such that $\text{supp}(C_{i})\subset\text{supp}(X^{*})$ and $\text{supp}(C_{i})\subset\text{supp}(X^{**})$ :

( $i$ )

$K_{C_{i}\rightarrow C_{j}}(X^{**})>K_{C_{i}\rightarrow C_{j}}(X^{*})$ implies that there is a species $X_{k}\in\text{supp}(C_{i})$ with $X_{k}^{**}>X_{k}^{*}$ .
( $ii$ )

$K_{C_{i}\rightarrow C_{j}}(X^{**})=K_{C_{i}\rightarrow C_{j}}(X^{*})$ implies that $X_{k}^{**}=X_{k}^{*}$ for all $X_{k}\in\text{supp}(C_{i})$ or else there are species $X_{k},X_{k}^{\prime}\in\text{supp}(C_{i})$ with $X_{k}^{**}>X_{k}^{*}$ and $(X_{k}^{\prime})^{**}<(X_{k}^{\prime})^{*}$ .

We say that a CKS is weakly monotonic when its kinetics is weakly monotonic.

Example A.3. Every MAK is weakly monotonic.

A.3 Decomposition Theory

A covering of a CRN is a collection of subsets $\{\mathscr{R}_{1},\ldots,\mathscr{R}_{k}\}$ whose union is $\mathscr{R}$ . A covering is called a decomposition of $\mathscr{N}$ if the sets $\mathscr{R}_{i}$ form a partition of $\mathscr{R}$ . $\mathscr{R}_{i}$ defines a subnetwork $\mathscr{N}_{i}$ of $\mathscr{N}$ where $\mathscr{N}_{i}=(\mathscr{S}_{i},\mathscr{C}_{i},\mathscr{R}_{i})$ such that $\mathscr{C}_{i}$ consists of all complexes occurring in $\mathscr{R}_{i}$ and $\mathscr{S}_{i}$ has all the species occurring in $\mathscr{C}_{i}$ . In this paper, we will denote a decomposition as a union of the subnetworks: $\mathscr{N}=\mathscr{N}_{1}\cup\ldots\cup\mathscr{N}_{k}$ . We refer to a “decomposition” with a single “subnetwork” $\mathscr{N}=\mathscr{N}_{1}$ as the trivial decomposition. Furthermore, when a network has been decomposed into subnetworks, we can refer to the said network as the parent network.

The most widely used decomposition of a reaction network is the set of linkage classes. Linkage classes have the special property that they not only partition the set of reactions but also the set of complexes.

A decomposition $\mathscr{N}=\mathscr{N}_{1}\cup\ldots\cup\mathscr{N}_{k}$ is independent if the parent network’s stoichiometric subspace $S$ is the direct sum of the subnetworks’ stoichiometric subspaces $S_{i}$ . Equivalently, the sum is direct when the rank of the parent network is equal to the sum of the ranks of the individual subnetworks, i.e.,

s=\sum_{i=1}^{k}s_{i}\text{ where }s_{i}=\text{dim}(S_{i}).

A network decomposition $\mathscr{N}=\mathscr{N}_{1}\cup\ldots\cup\mathscr{N}_{k}$ is a refinement of $\mathscr{N}=\mathscr{N}_{1}^{\prime}\cup\ldots\cup\mathscr{N}_{k^{\prime}}^{\prime}$ (and the latter a coarsening of the former) if it is induced by a refinement $\{\mathscr{R}_{1},\ldots,\mathscr{R}_{k}\}$ of $\{\mathscr{R}_{1}^{\prime},\ldots,\mathscr{R}_{k^{\prime}}^{\prime}\}$ .

Example A.4. If $\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}_{2}\cup\mathscr{N}_{3}$ and $\mathscr{N}^{\prime}=\mathscr{N}_{2}\cup\mathscr{N}_{3}$ , then

•

$\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}_{2}\cup\mathscr{N}_{3}$ is a refinement of $\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}^{\prime}$ ; and
•

$\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}^{\prime}$ is a coarsening of $\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}_{2}\cup\mathscr{N}_{3}$ .

A decomposition is said to be incidence independent if its incidence matrix is the direct sum of the incidence matrices of the subnetworks. Consider the network decomposition $\mathscr{N}=\mathscr{N}_{1}\cup\ldots\cup\mathscr{N}_{k}$ . If $\mathscr{N}$ ( $\mathscr{N}_{i}$ ) has $n$ ( $n_{i}$ ) complexes and $\ell$ ( $\ell_{i}$ ) linkage classes, then incidence independence is equivalent to

n-\ell=\sum_{i=1}^{k}(n_{i}-\ell_{i}).

Despite its early founding by Feinberg in 1987, decomposition theory has received little attention in the Chemical Reaction Network Theory (CRNT) community. Recently, however, its usefulness expanded beyond results on equilibria existence and parametrization, e.g., in the relatively new field of concentration robustness. Interesting results for large and high deficiency systems have been derived for more general kinetic systems such as power law systems [12, 13] and Hill-type systems [15].

Appendix B Variable Descriptions and Units

The following are the variables used in the metabolic insulin signaling network together with their descriptions and unit of measurement:

	$\displaystyle X_{2}=\text{Unbound surface insulin receptors (in molar)}$
	$\displaystyle X_{3}=\text{Unphosphorylated once-bound surface receptors (in molar)}$
	$\displaystyle X_{4}=\text{Phosphorylated twice-bound surface receptors (in molar)}$
	$\displaystyle X_{5}=\text{Phosphorylated once-bound surface receptors (in molar)}$
	$\displaystyle X_{6}=\text{Unbound unphosphorylated intracellular receptors (in molar)}$
	$\displaystyle X_{7}=\text{Phosphorylated twice-bound intracellular receptors (in molar)}$
	$\displaystyle X_{8}=\text{Phosphorylated once-bound intracellular receptors (in molar)}$
	$\displaystyle X_{9}=\text{Unphosphorylated IRS-1 (in molar)}$
	$\displaystyle X_{10}=\text{Tyrosine-phosphorylated IRS-1 (in molar)}$
	$\displaystyle X_{11}=\text{Unactivated PI 3-kinase (in molar)}$
	$\displaystyle X_{12}=\text{Tyrosine-phosphorylated IRS-1/activated PI 3-kinase complex (in molar)}$
	$\displaystyle X_{13}=\text{PI(3,4,5)P${}_{3}$ out of the total lipid population (in \%)}$
	$\displaystyle X_{14}=\text{PI(4,5)P${}_{2}$ out of the total lipid population (in \%)}$
	$\displaystyle X_{15}=\text{PI(3,4)P${}_{2}$ out of the total lipid population (in \%)}$
	$\displaystyle X_{16}=\text{Unactivated Akt (in \%)}$
	$\displaystyle X_{17}=\text{Activated Akt (in \%)}$
	$\displaystyle X_{18}=\text{Unactivated PKC-$\zeta$ (in \%)}$
	$\displaystyle X_{19}=\text{Activated PKC-$\zeta$ (in \%)}$
	$\displaystyle X_{20}=\text{Intracellular GLUT4 (in \%)}$
	$\displaystyle X_{21}=\text{Cell surface GLUT4 (in \%)}.$

Appendix C Determining the CRN of the Metabolic Insulin Signaling Network

In Section 2.2, it is determined by the Hars-Tóth criterion that there is a reaction network of the form $AX\xrightarrow{k}BX$ such that $f(X)=(B-A)^{\top}(k\circ X^{A})$ where $A$ and $B$ have nonnegative integer entries. In this section, we show the detailed computation on how to get matrices $A$ and $B$ .

To determine matrix $A$ , we write

k\circ X^{A}=\left[\begin{array}[]{l}k_{1}X_{2}\\ k_{2}X_{3}\\ k_{3}X_{5}\\ {\color[rgb]{1,1,1}k_{n}l}\vdots\\ k_{34}\\ k_{35}X_{20}\end{array}\right]

which gives us

A=\blockarray{ccccccl}X_{2}&X_{3}X_{4}\ldots X_{20}X_{21}\\ \block{[cccccc]l}100\ldots 00{\bigstrut[t]}R_{1}\\ 010\ldots 00R_{2}\\ 000\ldots 00R_{3}\\ \vdots\vdots\vdots\cdots\vdots\vdots{\color[rgb]{1,1,1}1}\vdots\\ 000\ldots 00R_{34}\\ 000\ldots 10{\bigstrut[b]}R_{35}\\

To get $B$ , observe first that we can write the ODE system as

\dot{X}=f(X)=\left[\begin{array}[]{rrrrrr}-1&1&0&\ldots&0&0\\ 1&-1&0&\ldots&0&0\\ 0&0&{\color[rgb]{1,1,1}-}1&\ldots&{\color[rgb]{1,1,1}-}0&0\\ \vdots&\vdots&\vdots&\cdots&\vdots&\vdots\\ 0&0&0&\ldots&1&-1\\ 0&0&0&\ldots&0&0\end{array}\right]\left[\begin{array}[]{l}k_{1}X_{2}\\ k_{2}X_{3}\\ k_{3}X_{5}\\ {\color[rgb]{1,1,1}k_{n}l}\vdots\\ k_{34}\\ k_{35}X_{20}\end{array}\right]

providing us with

(B-A)^{T}=\left[\begin{array}[]{rrrrrr}-1&1&0&\ldots&0&0\\ 1&-1&0&\ldots&0&0\\ 0&0&{\color[rgb]{1,1,1}-}1&\ldots&{\color[rgb]{1,1,1}-}0&0\\ \vdots&\vdots&\vdots&\cdots&\vdots&\vdots\\ 0&0&0&\ldots&1&-1\\ 0&0&0&\ldots&0&0\end{array}\right]

allowing us to easily get $B$ :

B=\blockarray{ccccccl}X_{2}&X_{3}X_{4}\ldots X_{20}X_{21}\\ \block{[cccccc]l}010\ldots 00{\bigstrut[t]}R_{1}\\ 100\ldots 00R_{2}\\ 001\ldots 00R_{3}\\ \vdots\vdots\vdots\cdots\vdots\vdots{\color[rgb]{1,1,1}1}\vdots\\ 000\ldots 10R_{34}\\ 000\ldots 00{\bigstrut[b]}R_{35}\\ .

Matrices $A$ and $B$ are used to determine the CRN corresponding to the ODE system $\dot{X}=f(X)$ in Section 2.1.

Appendix D Generalization of Lemma 3

In this section, we present two propositions related to Lemma 3 in Section 3.2. The first can be viewed as a corollary of statement $(i)$ of Lemma 3. The rest of the section leads to a generalization of this corollary to the case of a nonempty intersection of the species sets.

Proposition 5.

Let $\mathscr{N}=(\mathscr{S},\mathscr{C},\mathscr{R})$ be a CRN with kinetics $K$ and subsets $\mathscr{S}_{1}$ and $\mathscr{S}_{2}$ such that $\mathscr{S}_{1}\cup\mathscr{S}_{2}=\mathscr{S}$ . If $\mathscr{S}_{1}\cap\mathscr{S}_{2}=\varnothing$ , then for any kinetics, if each subsystem has an equilibrium in every stoichiometric class, then so does $(\mathscr{N},K)$ . If the equilibria are unique positive equilibria, the same holds for $(\mathscr{N},K)$ .

Proof.

By Feinberg’s Decomposition Theorem, any equilibrium of the whole system is also an equilibrium in each subnetwork. Hence, such an equilibrium in a stoichiometric compatibility class corresponds to its projection pair in the subnetwork stoichiometric compatibility classes. ∎

The next proposition describes the relationship between the stoichiometric classes of a CRN and those of a binary decomposition. This lays the foundation for the generalization of statement $(i)$ of Lemma 3 to the nonempty intersection case.

Proposition 6.

Let $\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}_{2}$ be a decomposition and $S$ , $S_{1}$ , and $S_{2}$ the corresponding stoichiometric subspaces. Let $d:\mathbb{R}^{\mathscr{S}}\rightarrow\mathbb{R}^{\mathscr{S}}/S_{1}\times\mathbb{R}^{\mathscr{S}}/S_{2}$ be the linear map $d(x)=(x+S_{1},x+S_{2})$ and $\Delta:\mathbb{R}^{\mathscr{S}}/S_{1}\times\mathbb{R}^{\mathscr{S}}/S_{2}\rightarrow\mathbb{R}^{\mathscr{S}}/S$ the linear map $\Delta(x_{1}+S_{1},x_{2}+S_{2})=(x_{1}-x_{2})+S$ . Then

$(i)$

$\Delta$ is surjective.
$(ii)$

$Ker(\Delta)=\{(x_{1}+S_{1},x_{2}+S_{2})\mid(x_{1}+S_{1})\cap(x_{2}+S_{2})\neq\varnothing\}$
$(iii)$

$Im(d)=Ker(\Delta)$

If, in addition, the decomposition is independent, then

$(iv)$

$d:\mathbb{R}^{\mathscr{S}}\rightarrow Ker(\Delta)$ is an isomorphism.
$(v)$

$|(x_{1}+S_{1})\cap(x_{2}+S_{2})|\leq 1$

Proof.

It is easily verified that the maps are well-defined.

$(i)$ For any $z+S\in\mathbb{R}^{\mathscr{S}}/S$ , write $z=z^{\prime}+z^{\prime\prime}$ where $z^{\prime}\in S$ and $z^{\prime\prime}\in S^{\perp}$ . Let $z^{\prime}=z^{\prime}_{1}+z^{\prime}_{2}$ where $z^{\prime}_{i}=p_{i}(z^{\prime})$ (the projection map $p_{i}$ is defined in statement $(i)$ of Lemma 3). Then $\Delta(\frac{z^{\prime\prime}}{2}+z^{\prime}_{2}+S_{1},-\frac{z^{\prime\prime}}{2}-z^{\prime}_{1}+S_{2})=z^{\prime\prime}+z^{\prime}_{2}+z^{\prime}_{1}+S$ .

$(ii)$ By definition, we have $\text{Ker}(\Delta)=\{(x_{1}+S_{1},x_{2}+S_{2})\mid x_{1}-x_{2}\in S\}$ . This implies that $x_{1}-x_{2}=s_{1}+s_{2}$ where $s_{i}\in S_{i}$ . Therefore, $x=x_{1}-s_{1}=x_{2}+s_{2}\in(x_{1}+S_{1})\cap(x_{2}+S_{2})$ . Conversely, $x=x_{1}+s_{1}=x_{2}+s_{2}$ implies that $x_{1}-x_{2}=(-s_{1})+s_{2}\in S$ .

$(iii)$ Clearly, $\text{Im}(d)\subseteq\text{Ker}(\Delta)$ . On the other hand, for any $(x_{1}+S_{1},x_{2}+S_{2})\in\text{Ker}(\Delta)$ and any $x\in(x_{1}+S_{1})\cap(x_{2}+S_{2})$ , $x+S_{1}=x_{1}+S_{1}$ and $x+S_{2}=x_{2}+S_{2}$ . Hence, $(x_{1}+S_{1},x_{2}+S_{2})=d(x)$ .

$(iv)$ $d(x)=(S_{1},S_{2})\Leftrightarrow x\in S_{1}$ and $x\in S_{2}\Rightarrow x=0$ since the sum is direct.

$(v)$ Suppose the intersection is nonempty with common elements $x$ and $x^{\prime}$ . Then $x-x^{\prime}\in S_{1}$ and $x-x^{\prime}\in S_{2}$ implying that $x-x^{\prime}=0$ since the sum is direct. ∎

Hence, for any binary independent decomposition, we have

\mathbb{R}^{\mathscr{S}}/S_{1}\times\mathbb{R}^{\mathscr{S}}/S_{2}=\text{Ker}(\Delta)\times\mathbb{R}^{\mathscr{S}}/S.

Since $\mathbb{R}^{\mathscr{S}}/S_{i}=\mathbb{R}^{\mathscr{S}}/S_{i}\times\mathbb{R}^{m-m_{i}}$ where $m$ and $m_{i}$ are the numbers of species in $\mathscr{S}$ and $\mathscr{S}_{i}$ , respectively, we have

	$\displaystyle\mathbb{R}^{\mathscr{S}}/S_{1}\times\mathbb{R}^{\mathscr{S}}/S_{2}$	$\displaystyle=\mathbb{R}^{\mathscr{S}}/S_{1}\times\mathbb{R}^{\mathscr{S}}/S_{2}\times\mathbb{R}^{m-m_{1}}\times\mathbb{R}^{m-m_{2}}$
		$\displaystyle=\mathbb{R}^{\mathscr{S}}/S_{1}\times\mathbb{R}^{\mathscr{S}}/S_{2}\times\mathbb{R}^{m-c}$

where $2m-m_{1}-m_{2}=m-c$ ( $c$ is the number of common species). Since $d$ is an isomorphism between $\mathbb{R}^{m}$ and $\text{Ker}(\Delta)$ , we have

\mathbb{R}^{\mathscr{S}}/S_{1}\times\mathbb{R}^{\mathscr{S}}/S_{2}=\mathbb{R}^{c}\times\mathbb{R}^{\mathscr{S}}/S.

This reduces to statement $(i)$ of Lemma 3 when $c=0$ .

Reaction Network Analysis of Metabolic Insulin Signaling

Abstract

1 Introduction

2 Model Realization as a Mass Action System

2.1 Insulin Signaling Model

2.2 Application of the Hars-Tóth Criterion for Mass Action System Realization

2.3 CRNT Analysis of the Mass Action System: Key Basic Properties

2.3.1 Positive Dependency

Proposition 1.

Proof.

Remark 1.

2.3.2 tt-minimality

2.3.3 Nonconservativity

2.3.4 Concordance

3 Independent Decompositions of the Insulin Signaling System

3.1 The Finest Nontrivial Independent Decomposition

Proposition 2.

Proof.

3.2 Deficiency-Oriented Coarsening

Lemma 3.

Proof.

Remark 2.

Proposition 4.

Proof.

Remark 3.

4 ACR of Species in Insulin Signaling

4.1 Absolute Concentration Robustness

Remark 4.

4.2 ACR Analysis of Rank 1 Subnetworks of the Finest Independent Decomposition

Remark 5.

4.3 ACR Analysis of Subnetworks in the Deficiency-Oriented Independent Decomposition

Remark 6.

4.4 ACR Analysis of the Metabolic Insulin Signaling Network for a Set of Rate Constants

Remark 7.

5 Summary and Outlook

References

Appendix A Notations and Definition of Terms

A.1 Chemical Reaction Networks

A.2 Chemical Kinetic Systems

A.3 Decomposition Theory

Appendix B Variable Descriptions and Units

Appendix C Determining the CRN of the Metabolic Insulin Signaling Network

Appendix D Generalization of Lemma 3

Proposition 5.

Proof.

Proposition 6.

Proof.

Reaction Network Analysis
of Metabolic Insulin Signaling

2.3.2 $t$ -minimality