# Comparative Analysis of Kinetic
Realizations of 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

Several studies have developed dynamical models to understand the underlying
mechanisms of insulin signaling, a signaling cascade that leads to the
translocation of glucose, the human body’s main source of energy. Fortunately,
reaction network analysis allows us to extract properties of dynamical systems
without depending on their model parameter values. This study focuses on the
comparison of insulin signaling in healthy state (INSMS or INSulin Metabolic
Signaling) and in type 2 diabetes (INRES or INsulin RESistance) using reaction
network analysis. The analysis uses network decomposition to identify the
different subsystems involved in insulin signaling (e.g., insulin receptor
binding and recycling, GLUT4 translocation, and ERK signaling pathway, among
others). Furthermore, results show that INSMS and INRES are similar with
respect to some network, structo-kinetic, and kinetic properties. Their
differences, however, provide insights into what happens when insulin
resistance occurs. First, the variation in the number of species involved in
INSMS and INRES suggests that when irregularities occur in the insulin
signaling pathway, other complexes (and, hence, other processes) get involved,
characterizing insulin resistance. Second, the loss of concordance exhibited
by INRES suggests less restrictive interplay between the species involved in
insulin signaling, leading to unusual activities in the signaling cascade.
Lastly, GLUT4 losing its absolute concentration robustness in INRES may
signify that the transporter has lost its reliability in shuttling glucose to
the cell, inhibiting efficient cellular energy production. This study also
suggests possible applications of the equilibria parametrization and network
decomposition, resulting from the analysis, to potentially establish absolute
concentration robustness in a species.

Keywords: insulin resistance, insulin signaling, kinetic realization, reaction
network, subnetwork

##  1 Introduction

Insulin signaling plays a crucial role in the human body’s energy metabolism.
Insulin reception triggers other processes that lead to the translocation of
glucose, the body’s main source of energy [36]. However, irregularities that
decrease insulin reception can lead to insulin resistance in cells.
Researchers have studied the role of insulin resistance in the development of
many metabolic disorders such as type 2 diabetes, and of other phenomena in
the human body including oxidative stress, inflammation, insulin receptor
mutation, and mitochondrial dysfunction [22, 24, 31]. To this day, the
underlying cause of insulin resistance has not been fully understood since the
phenomenon is highly complicated and the underlying mechanisms in healthy and
insulin-resistant cells have not yet been identified by researchers dealing
with the difficult problem.

Mathematical models have helped the scientific community to better understand
insulin signaling in healthy cells [15, 25, 26, 35] and diseased cells [2, 3,
23]. However, the complexity of metabolic insulin signaling have also
naturally led to complicated dynamical models. The results generated from
these models depend on simulations that rely on numerous parameter values that
need to be identified from literature or derived from experimental data.
Chemical Reaction Network Theory (CRNT) can address this challenge inherent in
many dynamical systems. Hence, analysis using CRNT is a useful approach in
studying biological systems.

CRNT translates a system of differential equations into a chemical reaction
network (CRN) representation or kinetic realization of the system’s network of
reactions. A CRN is composed of nonempty finite sets of species, complexes,
and reactions. A complex is one which appears as a reactant or a product in a
reaction. Upon identification of a kinetics (i.e., an assignment of rate
functions to the reactions), one can then talk about the chemical kinetic
system.

Constructing a particular kinetic realization of a system and identifying its
kinetics allow the extraction of formal properties of the system, independent
from parameter values of the model. Working with kinetic realizations also
allows us to compare different instances of the same system, e.g., insulin
signaling in healthy and insulin-resistant cells. One can also decompose the
CRN representation into subnetworks to determine subgroups in the system, a
task that is difficult to do with differential equations and biochemical maps.
Several papers have taken advantage of the usefulness of CRNT in dealing with
complex networks [1, 11, 13, 18, 19, 20, 27, 28, 32, 33, 34].

In this study, we use two systems of ordinary differential equations (ODEs): a
model of insulin signaling in healthy cells by Sedaghat et al. [26] and a
model of insulin signaling in type 2 diabetes by Nyman et al. [23]. Lubenia et
al. [19] have constructed a kinetic realization of the Sedaghat et al. mass
action model and have already performed a reaction network analysis of the
system. Their key findings are as follows: $(i)$ the underlying network is
concordant; $(ii)$ the network of three functional modules discussed by
Sedaghat and colleagues form an independent decomposition; and $(iii)$ the
network has 8 species with absolute concentration robustness (ACR), one of
which is the essential glucose transporter GLUT4 which, coupled with adequate
glucose supply, enables reliable cellular energy production. We denote the
kinetic realization of the insulin signaling in healthy cell by INSMS (INSulin
Metabolic Signaling).

A natural next step is to check what happens to insulin signaling when there
is insulin resistance. To this end, we use the Nyman et al. model of insulin
signaling in type 2 diabetes. We then perform a comparative analysis of the
kinetic realizations of the two states of insulin signaling, an analysis
which, to our knowledge, is the first of its kind regarding metabolic insulin
signaling. Similar to the Sedaghat et al. model, the Nyman et al. model is
also a mass action system, making our reaction network analysis parameter-
free. We denote the kinetic realization of the insulin signaling in type 2
diabetes by INRES (INsulin RESistance).

In our comparative analysis, we utilize decomposition theory to identify the
subsystems involved in insulin signaling. Decomposing the networks into
subnetworks allows us to identify in INRES the three functional modules that
also appear in INSMS: insulin receptor binding and recycling, postreceptor
signaling, and GLUT4 translocation. Furthermore, the decomposition shows
further subsystems that are present in INRES but not in INSMS: the S6 and S6K
formation from mTORC, ERK signaling pathway, and nuclear transcription. More
importantly, decomposition theory helps identify the subnetworks of INRES
which correspond to the “single mechanism” that explains insulin resistance in
type 2 diabetes according to Nyman et al.: an attenuated positive feedback
from mTORC to IRS1.

We also look at the essential properties of INSMS and INRES. The properties
considered are essential in the sense that they do not depend on the
parameters, i.e., rate constants, used in the mass action models they were
derived from. A network (or structural) property is one that can be specified
in terms of the network components alone (e.g., concordance). Properties which
may change under different dynamic equivalences are called structo-kinetic
properties (e.g., monostationarity) while those which are invariant under any
dynamic equivalence are called (purely) kinetic properties (e.g., ACR). INSMS
and INRES are similar with respect to some network, structo-kinetic, and
kinetic properties. But more importantly, this study confirms our expectation
that kinetic system differences reveal notable biological differences of the
same system under differing states, i.e., healthy and insulin-resistant.

First, the (molecular) species involved in insulin signaling (and consequently
the functional modules they constitute) differ strongly. This variation in the
number of species involved in INSMS and INRES suggests that when
irregularities occur in the insulin signaling pathway, other complexes (and,
hence, other processes) get involved, characterizing insulin resistance. This
and other properties point to a generally higher complexity of signal
processing in the insulin-resistant case.

Second, INRES is discordant, i.e., the concordance of INSMS is lost. According
to Shinar and Feinberg [30], concordant networks 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”. Thus, the loss of concordance
exhibited by INRES suggests less restrictive interplay between the species
involved in insulin signaling, leading to unusual activities in the signaling
cascade.

Finally, there are no ACR species in INRES, marking the loss of ACR in 6 of
its common species with INSMS, including the critical GLUT4 which is
responsible for transporting glucose into the glycolytic system. The loss of
ACR in GLUT4 suggests that, in an insulin-resistant cell, the transporter has
lost its reliability in shuttling glucose to the cell, inhibiting efficient
cellular energy production. We discuss in Section 5.5.2 possible applications
of a network’s equilibria parametrization and its finest independent
decomposition to potentially establish ACR in a species.

The similarities and differences in the features of the two networks suggest
that insulin signaling in healthy and insulin-resistant cells follow different
pathways leading to the glucose transporter GLUT4 and its translocation of
glucose.

This paper is organized as follows: in Section 2, we summarize the results of
the reaction network analysis of insulin signaling in healthy cells by Lubenia
et al. [19]. In Section 3, we construct the kinetic realization of the model
of insulin signaling in type 2 diabetes by Nyman et al. [23]. Section 4
details the comparative analysis of network properties of INSMS and INRES
while Section 5 considers the networks’ structo-kinetic and kinetic
properties. Finally, we conclude the paper with a summary in Section 6. An
Appendix with background on CRNT and other pertinent details of the models
considered is provided for the reader.

##  2 Reaction Network Analysis of Insulin Signaling in
Healthy Cells

This section revisits the results of the reaction network analysis of Lubenia
et al. [19] of a model of insulin signaling in healthy cells. They constructed
a chemical reaction network (CRN) with mass action kinetics based on the
insulin signaling model by Sedaghat et al. [26]. We denote this kinetic
realization by INSMS (INSulin Metabolic Signaling). Tables 1 and 2 give the
network numbers and an overview of the properties of INSMS. The reader may
refer to Appendix A for a review of some basic concepts in Chemical Reaction
Network Theory (CRNT).

Table 1: Network numbers of INSMS The high values of species, (reactant)
complexes, reactions, (reactant) rank, and (reactant) deficiency quantify the
high complexity of insulin signaling in healthy cells; the network numbers
indicate that INSMS is a closed network with high reactant diversity; the
network is also branching, $t$-minimal, not weakly reversible, and its
terminal classes contain points and cycles
Characteristic | Value
---|---
Species | 20
Complexes | 35
Reactant complexes | 24
Reactions | 35
Linkage classes | 13
Strong linkage classes | 24
Terminal strong linkage classes | 13
Rank | 15
Reactant rank | 20
Deficiency | 7
Reactant deficiency | 4
Table 2: Overview of properties of INSMS Some of the network properties listed
are described in Table 1; INSMS has a positive equilibrium for some sets of
rate constants (since it is positive dependent), its reactions have a high
degree of linear dependence among each other (due to high deficiency), and is
nonconservative and concordant; the key structo-kinetic properties of INSMS
include its monostationarity (it cannot have multiple positive equilibria for
a set of rate constants) and injectivity; its kinetic properties include the
nondegeneracy of its equilibria and the existence of 8 ACR species in the
system
Property Class | INSMS
---|---
Network | Closed | Not (weakly) reversible
$t$-minimal | Branching
Terminal classes include points and cycles
Positive dependent | Deficiency 7
Nonconservative | Concordant
Structo-Kinetic | Monostationary | Injective
Kinetic | Nondegenerate equilibria | 8 ACR species

INSMS is clearly a complex network as it involves 20 species, 35 complexes,
and 35 reactions. It is a closed network with high reactant diversity. Viewed
as a digraph, INSMS is branching, $t$-minimal, not (weakly) reversible, and
its terminal classes contain points and cycles. Based on its stoichiometry,
the network is positive dependent, of deficiency 7, nonconservative, and
concordant. The authors also observed that INSMS has a nontrivial finest
independent decomposition composed of 10 subnetworks. The biological
significance of these and further properties are discussed in the succeeding
chapters.

One of the key properties implied by the concordance of INSMS is its
injectivity. The weak momotonicity of the mass action system implies that the
network is also monostationary. Furthermore, all positive equilibria of INSMS
are nondegenerate. And one of the most important results of Lubenia and
colleagues was the discovery that 8 out of 20 species, including the critical
GLUT4, in the insulin signaling in healthy cells has absolute concentration
robustness (ACR).

##  3 Kinetic Realization of the Nyman et al. Model

In this section, we discuss the model of insulin signaling in type 2 diabetes
by Nyman et al. [23]. After an overview of their biological findings, we
construct a CRN representation of the dynamical system with mass action
kinetics and discuss some of the basic properties of this kinetic realization.

###  3.1 Summary of Novel Biological Insights from the Model

Using data based on human adipocytes from healthy and diabetic individuals,
Nyman and colleagues identified three hallmarks of insulin resistance in type
2 diabetes:

  1. 1.

Diminished concentration of insulin receptors;

  2. 2.

Low concentration of the glucose transporter GLUT4; and

  3. 3.

Impaired feedback from mammalian target of rapamycin complex 1 (mTORC1) to
insulin receptor substrate 1 (IRS1).

######  Remark 1.

Hallmarks 1 and 2 are reflected in the loss of concentration robustness of the
ACR species in INSMS (discussed in Section 5.5).

The model the authors used is a modification of an existing system of ordinary
differential equations (ODEs) representing insulin signaling in type 2
diabetes constructed by Brännmark et al. [3]. Nyman et al. added the
extracellular signal-regulated kinase (ERK) signaling pathway, an important
signaling branch in type 2 diabetes. After comparing observed data and
generating simulations, the authors determined that, among the three main
hallmarks of insulin resistance, the reduced positive feedback from mTORC1 to
IRS1 could explain insulin resistance in all parts of the insulin signaling
pathway.

###  3.2 The Reaction Network of the Nyman et al. Model

The Nyman et al. model is composed of 32 ODEs (see Appendix B for the system
of ODEs and the description of the variables). All 44 reactions in their model
are modeled using mass action kinetics.

For easy comparison, the species occurring in both Sedaghat et al. and Nyman
et al. models are denoted by the same variable. Variables for those occurring
solely in Nyman et al.’s model continue their numbering from those in the
Sedaghat et al. model.

For better visual orientation of the reader, we reconstruct the complete
biochemical map of the Nyman et al. mass action system in Figure 1. The
supplementary materials of Brännmark et al. [3] and Nyman et al. [23] contain
details of the construction of the ODEs.

Figure 1: Biochemical map of insulin signaling in type 2 diabetes $X_{2}$,
$X_{3}$, $X_{4}$, $X_{6}$, $X_{7}$, $X_{9}$, $X_{10}$, $X_{20},\ldots,X_{44}$
are the species of the network (see Appendix B for the description of the
variables); $k_{1},\ldots,k_{44}$ are the rate constants of the reactions; and
solid lines represent mass transfer reactions while broken lines represent
regulatory reactions; identifying subnetworks in insulin signaling in type 2
diabetes is difficult to do with biochemical maps but can be easily done using
kinetic realizations

Applying the the Hars-Tóth criterion for mass action system realization, as
used in Lubenia et al. [19] and detailed in Chellaboina et al. [4], a kinetic
realization of the Nyman et al. model is as follows:

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

This realization coincides with the biochemical map in Figure 1. Moreover, we
denote this kinetic realization by INRES (INsulin RESistance) (also
$\mathscr{N}=(\mathscr{S},\mathscr{C},\mathscr{R})$ with mass action kinetics
$K$, set $\mathscr{S}$ of 32 species, set $\mathscr{C}$ of 70 complexes, and
set $\mathscr{R}$ of 44 reactions). Tables 3 and 4 present the network numbers
and an overview of the properties of INRES.

Table 3: Network numbers of INRES The high values of species, (reactant)
complexes, reactions, (reactant) rank, and (reactant) deficiency quantify the
high complexity of insulin signaling in insulin-resistant cells; the network
numbers indicate that INRES is a closed network ($s<m$) with high reactant
diversity ($n_{r}>s$); the network is also branching ($r>n_{r}$), $t$-minimal
($t=\ell$), not weakly reversible ($s\ell\neq\ell$), and its terminal classes
contain points ($t\neq n-n_{r}$) and cycles ($n\neq n_{r}$)
Characteristic | Notation | Value
---|---|---
Species | $m$ | 32
Complexes | $n$ | 70
Reactant complexes | $n_{r}$ | 41
Reactions | $r$ | 44
Linkage classes | $\ell$ | 31
Strong linkage classes | $s\ell$ | 65
Terminal strong linkage classes | $t$ | 31
Rank | $s$ | 20
Reactant rank | $q$ | 32
Deficiency | $\delta$ | 19
Reactant deficiency | $\delta_{p}$ | 9
Table 4: Overview of properties of INRES Some of the network properties listed
are described in Table 3; INRES has a positive equilibrium for some sets of
rate constants (since it is positive dependent), its reactions have a high
degree of linear dependence among each other (due to high deficiency), and is
conservative and discordant; the key structo-kinetic properties of INRES
include its monostationarity (it cannot have multiple positive equilibria for
a set of rate constants) and non-injectivity; its kinetic properties include
the nondegeneracy of its equilibria, and its lack of ACR species
Property Class | INRES
---|---
Network | Closed | Not (weakly) reversible
$t$-minimal | Branching
Terminal classes include points and cycles
Positive dependent | Deficiency 19
Conservative | Discordant
Structo-Kinetic | Monostationary | Non-injective
Kinetic | Nondegenerate equilibria | No ACR species

The complexity of the network is evident in the network’s high number of
species, complexes, reactant complexes, and reactions, together with its high
deficiency ($\delta=19$). INRES is also a closed network (since $s=20<32=m$)
with high reactant diversity (since $n_{r}=41>20=s$).

Viewed as a digraph, the network is not weakly reversible (since $s\ell=65\neq
31=\ell$) but branching (since $r=44>41=n_{r}$). Other properties inferred by
the numbers include its $t$-minimality (since $t=31=\ell$), and non-point
terminality and non-cycle terminality (since $t\neq n-n_{r}=29$ and
$n-n_{r}\neq 0$, respectively).

The software CRNToolbox [8] can be used to analyze mass action systems, as is
the case for INRES. Using the application, we find that INRES is positive
dependent, meaning there is a set of rate constants wherein the network has a
positive equilibrium. The software’s Basic Report also indicates that the
network is conservative. The Mass Action Injectivity Report shows that the
system is not injective, implying that it is also discordant (this is
confirmed by CRNToolbox’s Concordance Report). Finally, the Higher Deficiency
Report concludes that INRES is monostationary, i.e., for a given set of rate
constants, the network cannot admit multiple equilibria (moreover, an
equilibrium of INRES cannot be degenerate).

###  3.3 The Finest Independent Decomposition

We now present the finest independent decomposition of
$\mathscr{N}=\\{R_{1},\ldots,R_{44}\\}$ (where $R_{1},\ldots,R_{44}$ are the
reactions of INRES) which will be helpful later in our analyses. Applying the
algorithm of Hernandez and De la Cruz [12] to determine the finest independent
decomposition of a network, we find that this decomposition has 12
subnetworks:

$\displaystyle\mathscr{N}_{1}=\\{R_{1},\ldots,R_{7}\\}$
$\displaystyle\mathscr{N}_{2}=\\{R_{8},\ldots,R_{14}\\}$
$\displaystyle\mathscr{N}_{3}=\\{R_{15},R_{16}\\}$
$\displaystyle\mathscr{N}_{4}=\\{R_{17},\ldots,R_{22}\\}$
$\displaystyle\mathscr{N}_{5}=\\{R_{23},R_{24},R_{25}\\}$
$\displaystyle\mathscr{N}_{6}=\\{R_{26},R_{27}\\}$
$\displaystyle\mathscr{N}_{7}=\\{R_{28},R_{29},R_{30}\\}$
$\displaystyle\mathscr{N}_{8}=\\{R_{31},R_{32}\\}$
$\displaystyle\mathscr{N}_{9}=\\{R_{33},R_{34}\\}$
$\displaystyle\mathscr{N}_{10}=\\{R_{35},R_{36},R_{37}\\}$
$\displaystyle\mathscr{N}_{11}=\\{R_{38},\ldots,R_{42}\\}$
$\displaystyle\mathscr{N}_{12}=\\{R_{43},R_{44}\\}.$

Table 5 presents the network numbers of the subnetworks of INRES. The analysis
and biological implications of the decomposition are discussed in Section 4.2.

Table 5: Network numbers of the subnetworks of INRES The finest independent decomposition of INRES has 12 subnetworks representing different processes involved in insulin signaling in type 2 diabetes Characteristics | $\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}$ | $\mathscr{N}_{11}$ | $\mathscr{N}_{12}$
---|---|---|---|---|---|---|---|---|---|---|---|---
Species | 6 | 6 | 3 | 7 | 4 | 3 | 4 | 3 | 3 | 4 | 5 | 3
Complexes | 7 | 8 | 4 | 10 | 6 | 4 | 6 | 4 | 4 | 6 | 7 | 4
Reactant complexes | 5 | 6 | 2 | 6 | 3 | 2 | 3 | 2 | 2 | 3 | 5 | 2
Reactions | 7 | 7 | 2 | 6 | 3 | 2 | 3 | 2 | 2 | 3 | 5 | 2
Linkage classes | 2 | 3 | 2 | 4 | 3 | 2 | 3 | 2 | 2 | 3 | 3 | 2
Strong linkage classes | 5 | 7 | 4 | 10 | 6 | 4 | 6 | 4 | 4 | 6 | 5 | 4
Terminal strong linkage classes | 2 | 3 | 2 | 4 | 3 | 2 | 3 | 2 | 2 | 3 | 3 | 2
Rank | 4 | 3 | 1 | 3 | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 1
Reactant rank | 5 | 6 | 2 | 6 | 3 | 2 | 3 | 2 | 2 | 3 | 5 | 2
Deficiency | 1 | 2 | 1 | 3 | 2 | 1 | 2 | 1 | 1 | 2 | 2 | 1
Reactant deficiency | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0

##  4 Comparative Analysis of Network Properties of
INSMS and INRES

This section compares the network properties of INSMS and INRES. A network (or
structural) property is one that can be specified in terms of the network
components alone, i.e., there is no need for specifying a kinetics.

We highlight three differing network properties:

  1. 1.

The sets of species (consequently, functional modules) of the networks are
different;

  2. 2.

There is a much higher level of complexity in INRES; and

  3. 3.

The concordance of INSMS is lost, i.e., INRES is discordant.

We conclude the section with a discussion of further differences between the
network properties of the two kinetic realizations.

###  4.1 Common Network Properties

Comparing Tables 2 and 4, we observe that INSMS and INRES are both branching,
$t$-minimal, positive dependent, not (weakly) reversible, non-point terminal,
and non-cycle terminal. Based on Sections 2 and 3.2, both are also closed
networks with high reactant diversity.

Lemma 3.5.4 of Feinberg [7] and Proposition 1 of Lubenia et al. [19] show that
a positive dependent network has a set of rate constants such that the
corresponding mass action system has a positive equilibrium. This result is
important when we check INRES for ACR species in Section 5.5.

###  4.2 Differences in Species Sets and Functional Modules

Among the 20 species in INSMS, only nine appear in INRES as well. Seven of
these are involved in the initial signaling steps while the other two are
significant in the process of glucose transfer to the glycolytic system. There
are 23 unique species in INRES (i.e., those not in INSMS). The presence of
these unique species suggests that other complexes get involved when
irregularities occur in the insulin signaling pathway, characterizing insulin
resistance. Moreover, this difference in species sets clearly shows that the
functional modules active in the intermediate steps of the two networks are
entirely different.

Lubenia et al.’s reaction network analysis of INSMS revealed the three
functional modules used by Sedaghat et al. in the construction of their model.
A coarsening of the finest independent decomposition of INRES also shows these
functional modules. A coarsening of a decomposition has some of the original
subnetworks combined to form fewer subnetworks. We utilize this concept here
to identify the various processes involved in insulin signaling in type 2
diabetes.

Consider the coarsening described in Table 6. $\mathscr{N}_{1}^{*}$,
$\mathscr{N}_{2}^{*}$, and $\mathscr{N}_{3}^{*}$ correspond to the subsystems
used to construct the Sedaghat et al. model (insulin receptor binding and
recycling, postreceptor signaling, and GLUT4 translocation, respectively),
albeit with different actors. The last subnetwork $\mathscr{N}_{4}^{*}$
appears in the Nyman et al. model only; it is a signaling branch (composed of
S6 and S6K formation from mTORC, ERK signaling pathway, and nuclear
transcription) that the authors deemed significant in insulin signaling in
type 2 diabetes.

Table 6: Coarsening of the finest independent decomposition of INRES in
relation to the Sedaghat et al. model $\mathscr{N}_{1}^{*}$,
$\mathscr{N}_{2}^{*}$, and $\mathscr{N}_{3}^{*}$ correspond to the insulin
receptor binding and recycling subsystem, postreceptor signaling subsystem,
and GLUT4 translocation subsystem, respectively; $\mathscr{N}_{4}^{*}$ is a
subsystem found in the Nyman et al. model representing a signaling branch
(composed of S6 and S6K formation from mTORC, ERK signaling pathway, and
nuclear transcription) that the authors deemed significant in insulin
signaling in type 2 diabetes
$\mathscr{N}=\mathscr{N}_{1}^{*}\cup\mathscr{N}_{2}^{*}\cup\mathscr{N}_{3}^{*}\cup\mathscr{N}_{4}^{*}$
---
$\mathscr{N}_{1}^{*}$ | $\mathscr{N}_{1}$
$\mathscr{N}_{2}^{*}$ | $\mathscr{N}_{2}\cup\ldots\cup\mathscr{N}_{6}$
$\mathscr{N}_{3}^{*}$ | $\mathscr{N}_{7}\cup\mathscr{N}_{8}$
$\mathscr{N}_{4}^{*}$ | $\mathscr{N}_{9}\cup\ldots\cup\mathscr{N}_{12}$

Table 7 presents another coarsening of the 12 subnetworks of INRES which
reveals the 9 subsystems that Nyman and colleagues considered in the
construction of their model. Subsystems $\mathscr{N}_{1}^{{}^{\prime}}$ to
$\mathscr{N}_{9}^{{}^{\prime}}$ correspond to the insulin receptor signaling,
IRS1 dynamics, negative feedback to insulin receptors, PKB dynamics, mTORC
dynamics, GLUT4 translocation, S6 and S6K formation from mTORC, ERK signaling
pathway, and nuclear transcription, respectively.

Table 7: Coarsening of the finest independent decomposition of INRES in
relation to the Nyman et al. model $\mathscr{N}_{1}^{{}^{\prime}}$ to
$\mathscr{N}_{9}^{{}^{\prime}}$ correspond to the insulin receptor signaling,
IRS1 dynamics, negative feedback to insulin receptors, PKB dynamics, mTORC
dynamics, GLUT4 translocation, S6 and S6K formation from mTORC, ERK signaling
pathway, and nuclear transcription, respectively
$\mathscr{N}=\mathscr{N}_{1}^{{}^{\prime}}\cup\mathscr{N}_{2}^{{}^{\prime}}\cup\mathscr{N}_{3}^{{}^{\prime}}\cup\mathscr{N}_{4}^{{}^{\prime}}\cup\mathscr{N}_{5}^{{}^{\prime}}\cup\mathscr{N}_{6}^{{}^{\prime}}\cup\mathscr{N}_{7}^{{}^{\prime}}\cup\mathscr{N}_{8}^{{}^{\prime}}\cup\mathscr{N}_{9}^{{}^{\prime}}$
---
$\mathscr{N}_{1}^{{}^{\prime}}$ | $\mathscr{N}_{1}$
$\mathscr{N}_{2}^{{}^{\prime}}$ | $\mathscr{N}_{2}$
$\mathscr{N}_{3}^{{}^{\prime}}$ | $\mathscr{N}_{3}$
$\mathscr{N}_{4}^{{}^{\prime}}$ | $\mathscr{N}_{4}$
$\mathscr{N}_{5}^{{}^{\prime}}$ | $\mathscr{N}_{5}\cup\mathscr{N}_{6}$
$\mathscr{N}_{6}^{{}^{\prime}}$ | $\mathscr{N}_{7}\cup\mathscr{N}_{8}$
$\mathscr{N}_{7}^{{}^{\prime}}$ | $\mathscr{N}_{9}\cup\mathscr{N}_{10}$
$\mathscr{N}_{8}^{{}^{\prime}}$ | $\mathscr{N}_{11}$
$\mathscr{N}_{9}^{{}^{\prime}}$ | $\mathscr{N}_{12}$

In the finest independent decomposition of INRES, subnetwork $\mathscr{N}_{2}$
is significant because it contains reaction $R_{11}$ which exhibits the mTORC
feedback to IRS1. Affecting this feedback are reactions $R_{23}$ and $R_{24}$
in subnetwork $\mathscr{N}_{5}$: these are the reactions for the activation of
mTORC. These reactions together form the “single mechanism” that explains
insulin resistance in type 2 diabetes according to Nyman et al.: an attenuated
positive feedback from mTORC to IRS1.

These analyses show how decomposition theory aids in identifying the
functional modules, and crucial processes and reactions involved in insulin
signaling in type 2 diabetes.

Aside from the number of species and functional modules, the number of
reactions and the value of the deficiency of INRES also indicate that it is
more complex than INSMS. Although INRES has only nine more reactions (44 vs
35), its deficiency is nearly triple that of INSMS (19 vs 7), signifying a
much higher level of (linear) dependence among the reactions. On the other
hand, the reactant deficiency more than doubles “only” from 4 to 9. This
appears to be mainly due to the (maximally) high reactant rank of INRES (32 vs
20 for INSMS). The length of the finest independent decomposition of INRES (12
vs 10 for INSMS) gives further testament to its higher complexity compared to
INSMS. Although this difference is relatively moderate, this is consistent
with the rank relation (20 for INRES vs 15 for INSMS).

###  4.3 Concordance vs Discordance

The most significant difference between INSMS and INRES in terms of network
properties is the concordance of the former and the discordance of the latter.
Discordance marks the loss of a network characteristic which Shinar and
Feinberg [30] consider 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”. The discordance of INRES suggests that the
interactions between the species in the system is less restrictive, leading to
irregularities in the insulin signaling process. The root causes of
discordance will be discussed in Section 5.2.

######  Remark 2.

While large classes of biochemical systems have been shown to be discordant,
only a handful of systems have been published as having the concordance
property. For example, Fariñas et al. [5] have shown that the independent
kinetic realization of any S-system in two or more variables is discordant. To
our knowledge, only three concordant biochemical networks have been published
to date: a kinetic realization of the Wnt signaling pathway in Chapter 10.5.4
of Feinberg [7], of the insulin signaling in healthy cells by Lubenia et al.
[19], and of the Schmitz carbon cycle model by Fortun and Mendoza [11].

###  4.4 Further Differences in Network Properties

Table 8 provides an overview of further differences between INSMS and INRES.

Table 8: Further differences between INSMS and INRES The nonconservativity of
INSMS, and its deficiency zero subnetworks and rank 1 subnetworks are
significant in the existence of ACR species in the network; INRES, on the
other hand, is conservative and, as shown in Section 5.5, has no ACR species
INSMS | INRES | Comment
---|---|---
Nonconservative | Conservative |
5 (of 10) subnetworks | 0 (of 12) subnetworks | Large subnetwork ($s=10$) with
with $\delta=0$ | with $\delta=0$ |  $\delta=0$ in INSMS
9 (of 10) subnetworks | 8 (of 12) subnetworks |
with $s=1$ | with $s=1$ |
0 (of 10) subnetworks | 12 (of 12) subnetworks |
with $\delta_{\rho}=0$ | with $\delta_{\rho}=0$ |

The analysis of Lubenia et al. showed that INSMS is nonconservative while
CRNToolbox results reveal that INRES is conservative. The nonconservativity of
INSMS was a necessary condition for the ACR in 8 species in INSMS to be
inferred from ACR in its deficiency zero subnetwork. Meanwhile, INRES has no
deficiency zero subnetwork and, as will be shown in Section 5.5, has no ACR
species as well.

Lubenia et al. also observered that INSMS has nine (out of 10) subnetworks
whose rank is 1. In INRES, eight (out of 12) subnetworks have rank 1. The
Meshkat et al. criterion [21] was applied to the INSMS rank 1 subnetworks to
determine two ACR species. We tried applying the said criterion to rank 1
subnetworks of INRES to determine ACR species but the effort yielded no
results.

Finally, another difference between INSMS and INRES is terms of the number of
subnetworks with zero reactant deficiency: INSMS has none while every
subnetwork of INRES has a reactant deficiency of zero.

##  5 Comparative Analysis of Structo-kinetic and
Kinetic Properties of INSMS and INRES

In this section, we discuss the structo-kinetic and kinetic properties of
INSMS and INRES, and explore their relationships to insulin signaling in
healthy and insulin-resistant cells. Structo-kinetic properties are properties
that depend on both network and kinetic properties while (purely) kinetic
properties depend on the kinetics alone.

We also highlight in this section the most striking difference in kinetic
properties of the two networks: INSMS has 8 ACR species while INRES has none.

###  5.1 Common Structo-kinetic Properties

The first common structo-kinetic property of INSMS and INRES is the absence of
complex balanced equilibria since both networks are not weakly reversible.

Secondly, the $t$-minimality of both networks implies the coincidence of their
respective kinetic and stoichiometric subspaces [9].

A notable third structo-kinetic commonality of INSMS and INRES is their
monostationarity for any mass action kinetics: the network cannot admit
multiple positive equilibria for a given a set of rate constants. For INSMS,
this property derives from a general property of weakly monotonic kinetics on
a concordant network. For INRES, monostationarity was concluded by the Higher
Deficiency Report of CRNToolbox and is considered an infrequent phenomenon in
view of the underlying network’s discordance.

###  5.2 Injective vs Non-injective

The reaction network analysis of INSMS showed that it is injective since the
network is concordant. On the other hand, the non-injectivity of INRES is
shown by the Mass Action Injectivity Report of the CRNToolbox.

Based on Theorem 4.11 of Shinar and Feinberg [29], the class of concordant
networks is precisely the class of networks that are injective for every
assignment of a weakly monotonic kinetics. Since every mass action kinetics is
weakly monotonic, the root cause of the discordance of INRES is the
combination of two factors: the weak monotonicity and non-injectivity of its
kinetics.

###  5.3 Monostationarity in All Weakly Monotonic Systems vs
Multistationarity for Some Weakly Monotonic Systems

Sections 5.1 and 5.2 show that while the injectivity of INSMS, induced by a
weakly monotonic kinetics on a concordant network, is lost on the discordant
INRES, monostationarity is maintained.

Since INSMS is a concordant network, then it has injectivity in all weakly
monotonic kinetic systems derived from it (Theorem 4.11 of Shinar and Feinberg
[29]). Moreover, together with Remark 4.4 in Shinar and Feinberg [29], we
conclude that INSMS is monostationary in all weakly monotonic monotonic
systems.

Despite its monostationarity for any mass action kinetics, it follows from
Theorem 10.5.10 of Feinberg [7] that for some weakly monotonic kinetics, INRES
will have the capacity for multiple positive equilibria in a stoichiometric
class, i.e., the network can display multistationarity. The root of this is,
hence, the network property of discordance in combination with the weak
monotonicity of its kinetics.

###  5.4 Common Kinetic Properties

We observed only a single common kinetic property of INSMS and INRES: the non-
degeneracy of all their positive equilibria. This property can be derived from
the $t$-minimality of both networks: $t$-minimality is a necessary condition
for the existence of nondegenerate equilibria, as detailed in Remark 4.11 of
Feinberg and Horn [9].

######  Remark 3.

  1. 1.

This kinetic property also implies the common structo-kinetic property of
coincidence of the kinetic and stoichiometric subsets.

  2. 2.

The nondegeneracy of all equilibria of INSMS was a consequence of concordance
of the underlying network: if a concordant network has weakly monotonic
kinetics, then all equilibria are degenerate. On the other hand, for INRES,
the nondegeneracy of its equilibria came from CRNToolbox results.

###  5.5 Presence of ACR in 8 of 20 Species vs. Absence of ACR in
32 Species

ACR is the invariance of the concentrations of a species at all positive
equilibria of a kinetic system [27, 28]. In Section 4.1, we showed that INSMS
and INRES have a positive equilibrium for a set of rate constants. Thus, we
can check for ACR species in both networks. Lubenia and colleagues already
showed that INSMS has 8 species (out of 20) with ACR and it is noteworthy that
this set of species includes GLUT4, the glucose transporter essential for
feeding the glycolytic process with glucose. ACR of a species suggests that
the species is reliable in maintaining the efficiency of the process it
belongs to.

Applying the algorithm of Fontanil et al. [10] and the Meshkat et al.
criterion [21], we find that, remarkably, none of the 32 species in INRES have
ACR. Significantly, GLUT4 lost its concentration robustness in the insulin-
resistant case, suggesting that it lost its reliability in shuttling glucose
to the cell, a crucial process in sustaining efficient energy production.

######  Remark 4.

Aside from GLUT4, five of the other species from INSMS, which lost their ACR
in INRES, are pooled insulin receptors. Hence, the first two key biological
findings of Nyman et al. on lower concentration of receptors and GLUT4 (see
Section 3.1) are also reflected in the loss of ACR.

####  5.5.1 Equilibria Parametrization Method Applied to the Sedaghat et al.
Model

Hernandez et al. [13] developed a method to determine the parameterized
equilibrium of a mass action system. They did this using the fact that the set
of positive equilibria of the whole network is the intersection of the
equilibria sets of the subnetworks [6], together with network translation
methods developed by Johnston and colleagues [16, 17].

The analysis of INSMS showed that the number of ACR species
$m_{\mathrm{ACR}}\geq 8$ for all rate constants and $m_{\mathrm{ACR}}=8$ for
some. Subsequently, Hernandez et al. [13] sharpened this to
$m_{\mathrm{ACR}}=8$ for all rate constants using their method applied to the
Sedaghat et al. model.

####  5.5.2 Equilibria Parametrization Method Applied to the Nyman et al.
Model

To verify that INRES indeed has no ACR species, we tried using the Hernandez
et al. method to derive the parametrized equilibrium of the Nyman et al.
model. However, since their method did not deliver full results, we used
network decomposition to compute manually the equilibria relationships and
obtained the following solution where $x_{i}$ denotes the concentration of
species $X_{i}$:

$\displaystyle x_{2}=\frac{k_{3}}{k_{1}}x_{3}$ $\displaystyle
x_{4}=\frac{k_{3}(k_{2}+k_{1})}{k_{1}(k_{4}+k_{6})}x_{3}$ $\displaystyle
x_{6}=\frac{k_{4}k_{3}(k_{2}+k_{1})}{k_{7}k_{1}(k_{4}+k_{6})}x_{3}$
$\displaystyle
x_{9}=\left(\frac{k_{10}}{k_{8}}+\frac{k_{13}k_{11}}{k_{8}(k_{12}+k_{13})}x_{31}\right)\frac{x_{10}}{x_{7}}$
$\displaystyle
x_{20}=\frac{k_{32}}{\left(\frac{k_{31}k_{19}k_{17}}{k_{30}(k_{18}+k_{19}x_{33})}x_{10}x_{26}x_{33}x_{34}\right)\left[\left(1+\frac{k_{20}k_{11}}{k_{22}(k_{12}+k_{13})}x_{10}x_{31}\right)\frac{k_{28}}{k_{21}}+\frac{k_{29}}{k_{22}}\right]}x_{21}$
$\displaystyle x_{22}=\frac{k_{11}}{k_{12}+k_{13}}x_{10}x_{31}$ $\displaystyle
x_{23}=\frac{k_{9}k_{10}}{k_{14}k_{8}}\frac{x_{10}}{x_{7}}+\left(\frac{k_{9}}{k_{8}}\frac{1}{x_{7}}+1\right)\frac{k_{13}k_{11}}{k_{14}(k_{12}+k_{13})}x_{10}x_{31}$
$\displaystyle
x_{24}=\frac{k_{16}k_{4}k_{3}(k_{2}+k_{1})}{k_{15}k_{5}k_{1}(k_{4}+k_{6})}\frac{x_{3}}{x_{7}x_{10}}$
$\displaystyle
x_{25}=\frac{k_{4}k_{3}(k_{2}+k_{1})}{k_{5}k_{1}(k_{4}+k_{6})}\frac{x_{3}}{x_{7}}$
$\displaystyle x_{27}=\frac{k_{17}}{k_{18}+k_{19}x_{33}}x_{10}x_{26}$
$\displaystyle
x_{28}=\frac{k_{19}k_{17}}{k_{22}(k_{18}+k_{19}x_{33})}x_{10}x_{26}x_{33}$
$\displaystyle
x_{29}=\left(1+\frac{k_{20}k_{11}}{k_{22}(k_{12}+k_{13})}x_{10}x_{31}\right)\frac{k_{19}k_{17}}{k_{21}(k_{18}+k_{19}x_{33})}x_{10}x_{26}x_{33}$
$\displaystyle
x_{30}=\frac{k_{25}}{\left(\frac{k_{17}}{k_{18}+k_{19}x_{33}}x_{10}x_{26}\right)\left[\left(1+\frac{k_{20}k_{11}}{k_{22}(k_{12}+k_{13})}x_{10}x_{31}\right)\frac{k_{23}k_{19}}{k_{21}}x_{33}+k_{24}\right]}x_{31}$
$\displaystyle x_{32}=\frac{k_{27}}{k_{26}}\frac{x_{33}}{x_{7}}$
$\displaystyle
x_{35}=\left(\frac{k_{19}k_{17}}{k_{30}(k_{18}+k_{19}x_{33})}x_{10}x_{26}x_{33}x_{34}\right)\left[\left(1+\frac{k_{20}k_{11}}{k_{22}(k_{12}+k_{13})}x_{10}x_{31}\right)\frac{k_{28}}{k_{21}}+\frac{k_{29}}{k_{22}}\right]$
$\displaystyle x_{37}=\frac{k_{33}}{k_{34}}x_{31}x_{36}$ $\displaystyle
x_{39}=\frac{k_{35}k_{33}}{k_{37}k_{34}}x_{31}x_{36}x_{38}+\frac{k_{36}}{k_{37}}x_{38}x_{41}$
$\displaystyle
x_{40}=\frac{k_{41}}{k_{38}x_{7}+\frac{k_{39}k_{11}}{k_{12}+k_{13}}x_{10}x_{31}+k_{40}}x_{41}$
$\displaystyle x_{42}=\frac{k_{41}}{k_{42}}x_{41}$ $\displaystyle
x_{44}=\frac{k_{43}}{k_{44}}x_{41}x_{43}$

Observe that the parametrized equilibrium of the species above depend on other
species $x_{3}$, $x_{7}$, $x_{10}$, $x_{21}$, $x_{26}$, $x_{31}$, $x_{33}$,
$x_{34}$, $x_{36}$, $x_{38}$, $x_{41}$, and $x_{43}$, implying that there are
infinitely many possible equilibrium solutions for the system. We, therefore,
conclude that INRES has no ACR species.

The parametrized equilibria can potentially be used in collaboration with
biologists. For example, let us focus on the equilibrium parametrization of
the concentration of the intracellular glucose transporter GLUT4 ($x_{20}$).
Measurements of the amount of GLUT4 in experimental data were used by Sedaghat
et al. [26] and Nyman et al. [23] in the construction of their models of
insulin signaling in healthy and insulin-resistant cells. The ACR of GLUT4 in
the healthy cell model suggests that keeping the amount of GLUT4
(approximately) constant could be valuable against insulin resistance as
glucose can be efficiently transported to the cell. The equilibria
parametrization shows that the concentration of GLUT4 is dependent on the
concentrations of other species ($x_{10}$, $x_{21}$, $x_{26}$, $x_{31}$,
$x_{33}$, and $x_{34}$). One can possibly work with biologists to determine if
the concentration of these species can be modified (e.g., using medical
interventions) to ensure that the amount of GLUT4 remains approximately
constant or at least bounded, resulting in continued efficient glucose
translocation into the cell.

Another possible application of our analysis in relation to ACR is the usage
of a network’s finest independent decomposition. The idea is to start by
looking at the subnetwork/s where a species of interest is present and finding
a way to establish ACR of the species in the small subnetwork. One can then
try to extend it to bigger subnetworks and, eventually, to the entire network.
Working with biologists is, of course, essential in this synthetic biology
approach to ensure that results are biologically sound.

##  6 Summary and Conclusion

To study and understand the underlying mechanism of insulin signaling, we have
constructed, analyzed, and compared INSMS and INRES. Section 2 summarizes the
reaction network analysis of Lubenia et al. on a kinetic realization of a
model of metabolic insulin signaling in healthy cells (INSMS). In Section 3,
we constructed a CRN of a model of insulin signaling in an insulin-resistant
cell (INRES).

We showed how to use decomposition theory to identify the subsystems involved
in insulin signaling. Subnetworks of INRES revealed the functional modules
that are also present in INSMS: insulin receptor binding and recycling,
postreceptor signaling, and GLUT4 translocation. Furthermore, additional
subsystems are present in INRES but not in INSMS: the S6 and S6K formation
from mTORC, ERK signaling pathway, and nuclear transcription. More notably,
decomposition theory allowed us to identify the subnetworks of INRES which
correspond to the “single mechanism” that explains insulin resistance in type
2 diabetes according to Nyman et al.: an attenuated positive feedback from
mTORC to IRS1.

We then compared the different network, structo-kinetic, and kinetic
properties of the two reaction networks. We have summarized in a table in
Appendix C all the properties of INSMS and INRES discussed in this paper.
INSMS and INRES are similar with respect to some properties. But more
importantly, we highlighted that kinetic system differences revealed notable
biological differences of the same system under differing states, i.e.,
healthy and insulin-resistant.

First, the set of species involved in insulin signaling (and consequently the
functional modules they constitute) are very different. This variation in the
number of species involved in INSMS and INRES suggests that when
irregularities occur in the insulin signaling pathway, other complexes (and,
hence, other processes) get involved, characterizing insulin resistance.

Second, the concordance of INSMS is lost in INRES. This suggests less
restrictive interplay between the species involved in insulin signaling,
leading to unusual activities in the signaling cascade.

Finally, INSMS has 8 ACR species while INRES has none. This marked the loss of
ACR in 6 of its common species with INSMS, including the critical GLUT4 which
is responsible for transporting glucose into the glycolytic system. The loss
of ACR in GLUT4 suggests that, in an insulin-resistant cell, the transporter
has lost its reliability in shuttling glucose to the cell, inhibiting
efficient cellular energy production. We also talked about possible ways of
using information gained from a network’s equilibria parametrization and its
finest independent decomposition to potentially establish ACR in a species. We
emphasized that consultation with biologist is essential to ensure proper
biological interpretation of results.

These findings, from our view, have shown the usefulness of doing a kinetic
realization analysis, beyond what we already know from the analysis of the
dynamical systems of Sedaghat et al. and Nyman et al.

## 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] Braatz, E, Coleman, R (2015) A mathematical model of insulin resistance in Parkinson’s disease. Comput Biol Chem 56:84–97.  https://doi.org/10.1016/j.compbiolchem.2015.04.003
  * [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] Fariñas H, Mendoza E, Lao A (2021) Chemical reaction network decompositions and realizations of S-systems. Philipp Sci Lett 14(1):147–157
  * [6] 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
  * [7] Feinberg M (2019) Foundations of chemical reaction network theory. Springer, Switzerland,
https://doi.org/10.1007/978-3-030-03858-8

  * [8] Feinberg M, Ellison P, Ji H, Knight D (2018) The Chemical Reaction Network Toolbox Version 2.35.  https://doi.org/10.5281/zenodo.5149266
  * [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] 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
  * [11] Fortun N, Mendoza E (2023) Comparative analysis of carbon cycle models via kinetic representations. J Math Chem 61:896–932.  https://doi.org/10.1007/s10910-022-01442-8
  * [12] 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
  * [13] Hernandez B, Lubenia P, Johnston M, Kim J (2023) A framework for deriving analytic steady states of biochemical reaction networks. PLoS Comput Biol 19(4):e1011039.  https://doi.org/10.1371/journal.pcbi.1011039
  * [14] 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
  * [15] Huang C, Wu M, Du J, Liu D, Chan C (2014) Systematic modeling for the insulin signaling network mediated by IRS1 and IRS2. J Theor Biol 355:40–52.  https://doi.org/10.1016/j.jtbi.2014.03.030
  * [16] Johnston M, Burton E (2019) Computing weakly reversible deficiency zero network translations using elementary flux modes. Bull Math Biol 81:1613–1644.  https://doi.org/10.1007/s11538-019-00579-z
  * [17] Johnston M, Müller S, Pantea C (2019) A deficiency-based approach to parameterizing positive equilibria of biochemical reaction systems. Bull Math Biol 81:1143–1172.
https://doi.org/10.1007/s11538-018-00562-0

  * [18] Karp R, Pérez-Millán M, Dasgupta T, Dickenstein A, Gunawardena J (2012) Complex-linear invariants of biochemical networks. J Theor Biol 311:130–138.
https://doi.org/10.1016/j.jtbi.2012.07.004

  * [19] Lubenia P, Mendoza E, Lao A (2022) Reaction network analysis of metabolic insulin signaling. Bull Math Biol 84(129):1–22.  https://doi.org/10.1007/s11538-022-01087-3
  * [20] MacLean A, Rosen Z, Byrne H, Harrington H (2015) Parameter-free methods distinguish Wnt pathway models and guide design of experiments. Proc Natl Acad Sci USA 112(9):2652–2657.  https://doi.org/10.1073/pnas.1416655112
  * [21] Meshkat N, Shiu A, Torres A (2022) Absolute concentration robustness in networks with low-dimensional stoichiometric subspace. Vietnam J Math
https://doi.org/10.1007/s10013-021-00524-5

  * [22] Nguyen T, Ta Q, Nguyen T, Le T, Vo V (2020) Role of insulin resistance in the Alzheimer’s disease progression. Neurochem Res 45:1481–1491.  https://doi.org/10.1007/s11064-020-03031-0
  * [23] Nyman E, Rajan M, Fagerholm S, Brännmark C, Cedersund G, Strålfors P (2014) A single mechanism can explain network-wide insulin resistance in adipocytes from obese patients with type 2 diabetes. J Biol Chem 289(48):33215–33230.  https://doi.org/10.1074/jbc.M114.608927
  * [24] Ormazabal V, Nair S, Elfeky O, Aguayo C, Salomon C, Zuñiga F (2018) Association between insulin resistance and the development of cardiovascular disease. Cardiovasc Diabetol 17:1–14.  https://doi.org/10.1186/s12933-018-0762-4
  * [25] Quon, Michael J. and Campfield, L. Arthur (1991) A mathematical model and computer simulation study of insulin receptor regulation. J Theor Biol 150(1):59–72.  https://doi.org/10.1016/S0022-5193(05)80475-8
  * [26] 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

  * [27] 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
  * [28] 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

  * [29] 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
  * [30] Shinar G, Feinberg M (2013) 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
  * [31] Tanase D, Gosav E, Costea C, Ciocoiu M, Lacatusu C, Maranduca M, Ouatu A, Floria M (2020) The intricate relationship between type 2 diabetes mellitus (T2DM), insulin resistance (IR), and nonalcoholic fatty liver disease (NAFLD). J Diabetes Res 2020:1–16.  https://doi.org/10.1155/2020/3920196
  * [32] Thomson M, Gunawardena J (2009) The rational parameterisation theorem for multisite post-translational modification systems. J Theor Biol 261(4):626–636.  https://doi.org/10.1016/j.jtbi.2009.09.003
  * [33] Thomson M, Gunawardena J (2009) Unlimited multistability in multisite phosphorylation systems. Nature 460, 274–277.  https://doi.org/10.1038/nature08102
  * [34] Villar J, Lubenia P, Mendoza E, Arceo C (2019) Structural stability analysis of models of dopamine synthesis and D1 receptor trafficking in RPT cells using CRNT. Philipp J Sci 148(3):523–533
  * [35] Wanant S, Quon M (2000) Insulin receptor binding kinetics: Modeling and simulation studies. J Theor Biol 205(3):355–364.  https://doi.org/10.1006/jtbi.2000.2069
  * [36] Yaribeygi H, Farrokhi F, Butler A, Sahebkar A (2019) Insulin resistance: Review of the underlying molecular mechanisms. J Cell Physiol 234(6):8152–8161.  https://doi.org/10.1002/jcp.27603

##  Appendix A Notations and Definition of Terms

In this section, we lay the foundation of Chemical Reaction Network Theory 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.

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}_{\geq 0}^{m}$ (the subscript $\geq 0$
means we consider only the nonnegative real numbers). 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}$. 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}$.

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}^{m}$ or $\mathbb{R}_{>0}^{m}$ 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}^{m}$ is the set $(x+S)\cap\mathbb{R}_{\geq 0}^{m}$
where $x+S$ is the left coset of $S$ containing $x$. Similarly, the positive
stoichiometric compatibility class containing $x\in\mathbb{R}_{>0}^{m}$ is the
set $(x+S)\cap\mathbb{R}_{>0}^{m}$.

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}^{m}\rightarrow\mathbb{R}_{\geq 0}$. The system $(\mathscr{N},K)$ is called
a chemical kinetic system (CKS).

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 concentrations 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}^{m}$ 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}^{m}$ 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 (i.e., all coordinates are
positive), 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 [29] 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}^{m}\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}^{m}\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}$.

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}^{m}$ 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}^{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}^{r}$ such that
$L(\alpha)=0$.

We define 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}$. Alternatively, we can also define the support of a vector
$x\in\mathbb{R}^{n}$ as the set of indices of the coordinates of $x$ from the
index set $\\{1,\dots,n\\}$ for which $x_{i}\neq 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:

  1. ($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.

  2. ($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}^{m}$, 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}^{m}$, 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^{**})$:

  1. ($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}^{*}$.

  2. ($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.

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}$.

##  Appendix B The Nyman et al. Model

This section presents the system of ODEs of the insulin signaling in type 2
diabetes model of Nyman et al. [23]. It also lists in detail the definition of
all the variables involved.

The system of ODEs of the Nyman et al. model are as follows:

$\displaystyle\dot{x}_{2}=-k_{1}x_{2}-k_{2}x_{2}+k_{6}x_{4}+k_{7}x_{6}$
$\displaystyle\dot{x}_{4}=k_{2}x_{2}+k_{3}x_{3}-k_{4}x_{4}-k_{6}x_{4}$
$\displaystyle\dot{x}_{3}=k_{1}x_{2}-k_{3}x_{3}$
$\displaystyle\dot{x}_{7}=k_{4}x_{4}-k_{5}x_{7}x_{25}$
$\displaystyle\dot{x}_{6}=k_{5}x_{7}x_{25}-k_{7}x_{6}$
$\displaystyle\dot{x}_{9}=-k_{8}x_{7}x_{9}-k_{9}x_{9}+k_{10}x_{10}+k_{14}x_{23}$
$\displaystyle\dot{x}_{10}=k_{8}x_{7}x_{9}-k_{10}x_{10}-k_{11}x_{10}x_{31}+k_{12}x_{22}$
$\displaystyle\dot{x}_{22}=k_{11}x_{10}x_{31}-k_{12}x_{22}-k_{13}x_{22}$
$\displaystyle\dot{x}_{23}=k_{9}x_{9}+k_{13}x_{22}-k_{14}x_{23}$
$\displaystyle\dot{x}_{24}=-k_{15}x_{10}x_{24}+k_{16}x_{25}$
$\displaystyle\dot{x}_{25}=k_{15}x_{10}x_{24}-k_{16}x_{25}$
$\displaystyle\dot{x}_{26}=-k_{17}x_{10}x_{26}+k_{18}x_{27}+k_{22}x_{28}$
$\displaystyle\dot{x}_{27}=k_{17}x_{10}x_{26}-k_{18}x_{27}-k_{19}x_{27}x_{33}$
$\displaystyle\dot{x}_{28}=-k_{20}x_{22}x_{28}+k_{21}x_{29}-k_{22}x_{28}$
$\displaystyle\dot{x}_{29}=k_{19}x_{27}x_{33}+k_{20}x_{22}x_{28}-k_{21}x_{29}$
$\displaystyle\dot{x}_{30}=-k_{23}x_{29}x_{30}-k_{24}x_{27}x_{30}+k_{25}x_{31}$
$\displaystyle\dot{x}_{31}=k_{23}x_{29}x_{30}+k_{24}x_{27}x_{30}-k_{25}x_{31}$
$\displaystyle\dot{x}_{32}=-k_{26}x_{7}x_{32}+k_{27}x_{33}$
$\displaystyle\dot{x}_{33}=k_{26}x_{7}x_{32}-k_{27}x_{33}$
$\displaystyle\dot{x}_{34}=-k_{28}x_{29}x_{34}-k_{29}x_{28}x_{34}+k_{30}x_{35}$
$\displaystyle\dot{x}_{35}=k_{28}x_{29}x_{34}+k_{29}x_{28}x_{34}-k_{30}x_{35}$
$\displaystyle\dot{x}_{20}=-k_{31}x_{20}x_{35}+k_{32}x_{21}$
$\displaystyle\dot{x}_{21}=k_{31}x_{20}x_{35}-k_{32}x_{21}$
$\displaystyle\dot{x}_{36}=-k_{33}x_{31}x_{36}+k_{34}x_{37}$
$\displaystyle\dot{x}_{37}=k_{33}x_{31}x_{36}-k_{34}x_{37}$
$\displaystyle\dot{x}_{38}=-k_{35}x_{37}x_{38}-k_{36}x_{38}x_{41}+k_{37}x_{39}$
$\displaystyle\dot{x}_{39}=k_{35}x_{37}x_{38}+k_{36}x_{38}x_{41}-k_{37}x_{39}$
$\displaystyle\dot{x}_{40}=-k_{38}x_{7}x_{40}-k_{39}x_{22}x_{40}-k_{40}x_{40}+k_{42}x_{42}$
$\displaystyle\dot{x}_{41}=k_{38}x_{7}x_{40}+k_{39}x_{22}x_{40}+k_{40}x_{40}-k_{41}x_{41}$
$\displaystyle\dot{x}_{42}=k_{41}x_{41}-k_{42}x_{42}$
$\displaystyle\dot{x}_{43}=-k_{43}x_{41}x_{43}+k_{44}x_{44}$
$\displaystyle\dot{x}_{44}=k_{43}x_{41}x_{43}-k_{44}x_{44}$

where the following are concentrations:

$\displaystyle x_{2}=\text{Inactive insulin receptors}$ $\displaystyle
x_{3}=\text{Insulin-bound receptors}$ $\displaystyle x_{4}=\text{Tyrosine-
phosphorylated receptors}$ $\displaystyle x_{6}=\text{Internalized
dephosphorylated receptors}$ $\displaystyle x_{7}=\text{Tyrosine-
phosphorylated and internalized receptors}$ $\displaystyle
x_{9}=\text{Inactive IRS-1}$ $\displaystyle x_{10}=\text{Tyrosine-
phosphorylated IRS-1}$ $\displaystyle x_{20}=\text{Intracellular GLUT4}$
$\displaystyle x_{21}=\text{Cell surface GLUT4}$ $\displaystyle
x_{22}=\text{Combined tyrosine/serine 307-phosphorylated IRS-1}$
$\displaystyle x_{23}=\text{Serine 307-phosphorylated IRS-1}$ $\displaystyle
x_{24}=\text{Inactive negative feedback}$ $\displaystyle x_{25}=\text{Active
negative feedback}$ $\displaystyle x_{26}=\text{Inactive PKB}$ $\displaystyle
x_{27}=\text{Threonine 308-phosphorylated PKB}$ $\displaystyle
x_{28}=\text{Serine 473-phosphorylated PKB}$ $\displaystyle
x_{29}=\text{Combined threonine 308/serine 473-phosphorylated PKB}$
$\displaystyle x_{30}=\text{mTORC1}$ $\displaystyle x_{31}=\text{mTORC1
involved in phosphorylation of IRS-1 at serine 307}$ $\displaystyle
x_{32}=\text{mTORC2}$ $\displaystyle x_{33}=\text{mTORC2 involved in
phosphorylation of PKB at threonine 473}$ $\displaystyle x_{34}=\text{AS160}$
$\displaystyle x_{35}=\text{AS160 phosphorylated at threonine 642}$
$\displaystyle x_{36}=\text{S6K}$ $\displaystyle x_{37}=\text{Activated S6K
phosphorylated at threonine 389}$ $\displaystyle x_{38}=\text{S6}$
$\displaystyle x_{39}=\text{Activated S6 phosphorylated at serine 235 and
serine 236}$ $\displaystyle x_{40}=\text{ERK}$ $\displaystyle x_{41}=\text{ERK
phosphorylated at threonine 202 and tyrosine 204}$ $\displaystyle
x_{42}=\text{ERK sequestered in an inactive pool}$ $\displaystyle
x_{43}=\text{Elk1}$ $\displaystyle x_{44}=\text{Elk1 phosphorylated at serine
383}.$

Note that all the species concentrations are measured in %. The total
concentration of each of the following sets of species is 100% at any given
time: $\\{x_{2},x_{3},x_{4},x_{6},x_{7}\\}$,
$\\{x_{9},x_{10},x_{22},x_{23}\\}$, $\\{x_{20},x_{21}\\}$,
$\\{x_{24},x_{25}\\}$, $\\{x_{26},\ldots,x_{29}\\}$, $\\{x_{30},x_{31}\\}$,
$\\{x_{32},x_{33}\\}$, $\\{x_{34},x_{35}\\}$, $\\{x_{36},x_{37}\\}$,
$\\{x_{38},x_{39}\\}$, $\\{x_{40},x_{41},x_{42}\\}$, and
$\\{x_{43},x_{44}\\}$.

We observe that the system has mass action kinetics since the exponent of
every species concentration in the terms of the ODEs correspond to their
stoichiometry in their respective reaction.

##  Appendix C Summary of the Properties of
INSMS and INRES

We provide in this section a table showing the different properties discussed
in this paper.

Table C.1 summarizes the properties of INSMS and INRES.

Table C.1: Summary of the properties of INSMS and INRES
Property Class | INSMS | INRES
---|---|---
Network | Branching
Closed
High reactant diversity
Non-cycle terminal
Non-point terminal
Not (weakly) reversible
Positive dependent
$t$-minimal
10 subnetworks | 12 subnetworks
Nonconservative | Conservative
Concordant | Discordant
5 subnetworks with $\delta=0$ | 0 subnetworks with $\delta=0$
9 subnetworks with $s=1$ | 8 subnetworks with $s=1$
0 subnetworks with $\delta_{\rho}=0$ | 12 subnetworks with $\delta_{\rho}=0$
Structo-Kinetic | No complex balanced equilibria
Kinetic and stoichiometric subspaces coincide
Monostationary
Injective | Non-injective
Monostationary in all weakly | Multistationary for some weakly
monotonic systems |  monotonic systems
Kinetic | Nondegenerate equilibria
8 ACR species | No ACR species