definitionDefinition[section] remark[definition]Remark proofProof potProof of
Theorem LABEL:thm

[1]

1]organization=Centro de Bioinformática Médica, Universidade Federal de São
Paulo, addressline=Edifício de Pesquisas 2, city=São Paulo,
postcode=04039-032, state=SP, country=Brazil

[cor1]Corresponding author

2]organization=Department of Mathematics, The Ohio State University,
addressline=231 W 18th Ave, city=Columbus, postcode=43210, state=OH,
country=USA

3]organization=Mathematics Institute, University of Warwick,
addressline=Zeeman Building, city=Coventry, citysep=, postcode=CV4 7AL,
country=UK

# Homeostasis in Input-Output Networks:
Structure, Classification and Applications

Fernando Antoneli<EMAIL_ADDRESS>[  Martin Golubitsky
<EMAIL_ADDRESS>[  Jiaxin Jin<EMAIL_ADDRESS>Ian Stewart
<EMAIL_ADDRESS>[

###### Abstract

Homeostasis is concerned with regulatory mechanisms, present in biological
systems, where some specific variable is kept close to a set value as some
external disturbance affects the system. Many biological systems, from gene
networks to signaling pathways to whole tissue/organism physiology, exhibit
homeostatic mechanisms. In all these cases there are homeostatic regions where
the variable is relatively to insensitive external stimulus, flanked by
regions where it is sensitive. Mathematically, the notion of homeostasis can
be formalized in terms of an input-output function that maps the parameter
representing the external disturbance to the output variable that must be kept
within a fairly narrow range. This observation inspired the introduction of
the notion of infinitesimal homeostasis, namely, the derivative of the input-
output function is zero at an isolated point. This point of view allows for
the application of methods from singularity theory to characterize
infinitesimal homeostasis points (i.e. critical points of the input-output
function). In this paper we review the infinitesimal approach to the study of
homeostasis in input-output networks. An input-output network is a network
with two distinguished nodes ‘input’ and ‘output’, and the dynamics of the
network determines the corresponding input-output function of the system. This
class of dynamical systems provides an appropriate framework to study
homeostasis and several important biological systems can be formulated in this
context. Moreover, this approach, coupled to graph-theoretic ideas from
combinatorial matrix theory, provides a systematic way for classifying
different types of homeostasis (homeostatic mechanisms) in input-output
networks, in terms of the network topology. In turn, this leads to new
mathematical concepts, such as, homeostasis subnetworks, homeostasis patterns,
homeostasis mode interaction. We illustrate the usefulness of this theory with
several biological examples: biochemical networks, chemical reaction networks
(CRN), gene regulatory networks (GRN), Intracellular metal ion regulation and
so on.

###### keywords:

Infinitesimal Homeostasis Input-Output Networks Perfect Adaptation

##  1 Introduction

The idea of homeostasis has its roots in the work of the French physiologist
Claude Bernard [17] who observed a kind of regulation in the ‘milieu
intérieur’ (internal environment) of human organs such as the liver and
pancreas. The term ’homeostasis’, coined by the American physiologist Walter
Cannon in 1926 [24, 25], derives from the Greek language and refers to a
capacity of maintaining _similar stasis_. In this context, ‘homeostasis’
refers to certain regulatory mechanisms by which a feature is maintained at a
steady state despite disturbances caused by changes in the environment. At the
scale of a whole organism, homeostasis manifests itself in many forms: body
temperature, blood sugar level, concentration of ions in body fluids with
changes in the external environment.

More recently, the study of homeostasis in biological systems through
mathematical models has gained prominence. For example, the extensive work of
Nijhout, Reed, Best and collaborators [107, 18, 106, 105, 104] consider
biochemical networks associated with metabolic signaling pathways. Further
examples include regulation of cell number and size [87], control of sleep
[135], and expression level regulation in housekeeping genes [7].

In order to mathematically model ‘homeostasis’ is necessary to have a
framework where a precise definition can be given. The simplest mathematical
setup, very often used to model biological phenomena, is the theory of
ordinary differential equations. In this setting ‘homeostasis’ can be
interpreted in two mathematically distinct ways. One boils down to the
existence of a ‘globally stable equilibrium’. Here, changes in the environment
are considered to be perturbations of _initial conditions_ [128]. A stronger
(and, in our view, more appropriate) usage works with a parametrized family of
differential equations, with a corresponding family of stable equilibria. Now
‘homeostasis’ means that some quantity defined through this equilibrium
changes by a relatively small amount when the parameter varies by a much
larger amount. Homeostasis does not imply that the whole system remains
invariant with change in external variables. In fact, changes in external
variables can cause changes in certain internal variables, while other
internal variables remain almost unchanged.

A precise definition of homeostasis can be give as follows. Consider system of
ODEs depending on an _input parameter_ ${\mathcal{I}}$, which varies over a
range of external stimuli. Suppose there is a family of equilibrium points
$X({\mathcal{I}})$ and an observable $\phi$, such that the _input-output
function_ $z({\mathcal{I}})=\phi(X({\mathcal{I}}))$ is well-defined on the
range of ${\mathcal{I}}$. In this situation, we say that the system exhibits
_homeostasis_ if, under variation of the input parameter ${\mathcal{I}}$, the
input-output function $z({\mathcal{I}})$ remains approximately constant over
the interval of external stimuli (see section 2). Golubitsky and Stewart [54]
observe that homeostasis on some neighborhood of a specific value
${\mathcal{I}}_{0}$ follows from infinitesimal homeostasis, where
$z^{\prime}({\mathcal{I}}_{0})=0$ and ′ indicates differentiation with respect
to ${\mathcal{I}}$. This observation is essentially the well-known idea that
the value of a function changes most slowly near a stationary (or critical)
point.

Infinitesimal homeostasis is a _sufficient_ condition for homeostasis over
some interval of parameters, but it is not a necessary condition. A function
can vary slowly without having a stationary point. In applications, the
quantity that experiences homeostasis, represented by the observable $\phi$
can be a function of several internal variables, such as a sum of
concentrations, or the period of an oscillation (see subsection 4.4).

There are other variations on the formulation of homeostasis. Control-
theoretic models of homeostasis often require _perfect homeostasis_ , also
known as _perfect adaptation_ , in which the input-output function is exactly
constant over the parameter range [80, 81, 47]. The requirement of perfect
homeostasis is very strong, demanding that the model equations have a somewhat
restricted functional form [10].

Another important idea of [54] is to take advantage of fact that infinitesimal
homeostasis is suitable for analysis using methods from singularity theory
(see subsection 2.2). A motivational example for the use of singularity theory
to analyze homeostasis is the regulation of the output ‘body temperature’ in
an opossum, when the input ‘environmental temperature’ varies [101]. In [57,
Figs. 1–2] the graphs of body temperature against environmental temperature
${\mathcal{I}}$ are approximately linear, with nonzero slope, when
${\mathcal{I}}$ is either small or large, while in between is a broad flat
region, where homeostasis occurs. This general shape is called a ‘chair’ by
Nijhout and Reed [106] (see also [103, 107]) and suggest that there is a chair
in the body temperature data of opossums [101] given by such a piecewise
linear function [103, Fig. 1]. Golubitsky and Stewart [54] take a singularity-
theoretic point of view and suggest that chairs are better described locally
by a homogeneous cubic function (that is, like
$z({\mathcal{I}})\approx{\mathcal{I}}^{3}$) rather than by the previous
piecewise linear description. Based on this suggestion, [54] goes on to show
that the singularities of input-output functions are described by _Elementary
Catastrophe Theory_ [126, 127, 142].

The last contribution of [54] is that homeostasis can be thought of as a
network concept. This is motivated by the abundance of networks in biology,
specially connected to ODE modeling, see e.g. [103, 115, 46]. In fact, almost
all ODE models in biology come together with a ‘wiring diagram’, or a network,
describing the interactions among the elements in the model. The recently
published book [56] gives an exposition of a formal framework for studying
networks of coupled ODEs, that have been developed the authors and
collaborators in the past decade. More precisely, a network of coupled ODEs is
determined by directed graph whose nodes and edges are classified into types.
Nodes, or cells, represent the variables of a component ODE. Edges, or arrows,
between nodes represent couplings from the tail node to the head node. Nodes
of the same type have the same phase spaces (up to a canonical
identification); edges of the same type represent identical couplings.

In the setting of network coupled ODEs there is a preferable class of
observables, namely the output of the node variables (coordinate functions).
Network systems are distinguished from large systems by the ability to keep
track of the output from each node individually. Now homeostasis can be
naturally defined as the fact that the output of a node variable (the ‘output
node’) is held approximately constant as other variables (other nodes) vary
(perhaps wildly) under variation of an input parameter that affects another
node (the ‘input node’). Placing homeostasis in the general context of network
dynamics leads naturally to the methods reviewed here.

A special kind of network of coupled ODE, with two distinguished nodes (one
input node and one output node) is called an _input-output network_. It was
introduced in [58] in their study of homeostasis on $3$-node networks. Wang et
al. [134] extended the notion of _input-output network_ to arbitrary large
networks and developed a combinatorial theory for the classification of
‘homeostasis types’ in such networks. A _homeostasis type_ is essentially a
‘combinatorial mechanism’ that causes homeostasis in the full network and is
represented by specific subnetworks, called _homeostasis subnetworks_. The
homeostasis subnetworks can be subdivided into to classes called _structural_
and _appendage_.

The motivation for the term _structural homeostasis_ comes from [115], where
the authors identify the feedforward loop as one of the homeostastic motifs in
$3$-node biochemical networks, that is, a generalized ‘feedforward loop’ The
intuition behind the term _appendage homeostasis_ is that homeostasis is
generated by a cycle of regulatory nodes, that is, a generalized ‘feedback
loop’. The structural and appendage classes are abstract generalizations of
the usual ‘feedforward’ and ‘feedback’ mechanisms. A striking outcome the
approach of [134] is that they do not specify any homeostasis generating
mechanisms at the outset and they find _a posteriori_ that there are
essentially only the two types of homeostasis generating mechanisms:
_generalized feedback_ and _generalized feedforward_.

In Duncan et al. [40] the combinatorial formalism of [134] id further refined
to allow one to determine, for each homeostasis subnetwork, which nodes are
(and are not) going to be simultaneously homeostatic (besides the designated
output node) when that particular homeostasis subnetwrok is the ‘homeostasis
trigger’, i.e. the one that causes homeostasis in the network. These
collections of simultaneous homeostatic nodes are called _homeostasis
patterns_. The general picture that emerges from all these results is the
following: the homeostasis subnetworks and the homeostasis patterns of network
correspond uniquely to each other.

##  2 Structure

In this section we give the basic definition of homeostasis in a parametrized
system of ordinary differential equations. We start with a very general
definition that will be specialized to more ‘concrete’ situations throughout
the paper.

###  2.1 A Dynamical Formalism for Homeostasis

Golubitsky and Stewart [54, 55] proposed a mathematical framework for the
study of homeostasis based on dynamical systems theory (see [57] and the book
[56]).

In this framework one considers a system of differential equations

$\dot{X}=F(X,{\mathcal{I}})$ (2.1)

with state variable $X=(x_{1},\ldots,x_{n})\in\mathbb{R}^{n}$ and input
parameters
${\mathcal{I}}=({\mathcal{I}}_{1},\ldots,{\mathcal{I}}_{N})\in\mathbb{R}^{N}$
representing the external input to the system.

Suppose that $(X^{*},{\mathcal{I}}^{*})$ is a linearly stable equilibrium of
(2.1). By the implicit function theorem, there is a function
$\tilde{X}({\mathcal{I}})$ defined in a neighborhood of ${\mathcal{I}}^{*}$
such that $\tilde{X}({\mathcal{I}}^{*})=X^{*}$ and
$F(\tilde{X}({\mathcal{I}}),{\mathcal{I}})\equiv 0$.

A smooth function $\phi:\mathbb{R}^{n}\to\mathbb{R}$ is called an
_observable_. Define the _input-output function_
$z:\mathbb{R}^{N}\to\mathbb{R}$ associated to $\phi$ and $\tilde{X}$ as
$z({\mathcal{I}})=\phi(\tilde{X}({\mathcal{I}}))$. The input-output function
allows one to formulate several definitions that capture the notion of
homeostasis (see [89, 4, 125, 54, 55]).

######  Definition 2.1.

Let $z({\mathcal{I}})$ be an input-output function. We say that
$z({\mathcal{I}})$ exhibits

  1. [(a)]

  2. 1.

_Perfect Homeostasis_ on an open set $\Omega\subseteq\operatorname{dom}(z)$ if

$\nabla z({\mathcal{I}})=0\qquad\text{for all}\;{\mathcal{I}}\in\Omega$ (2.2)

That is, $z$ is constant on $\Omega$.

  3. 2.

_Near-perfect Homeostasis_ relative to a _set point_
${\mathcal{I}}_{0}\in\Omega\subseteq\operatorname{dom}(z)$ if, for fixed
$\delta$,

$\|z({\mathcal{I}})-z({\mathcal{I}}_{0})\|\leqslant\delta\qquad\text{for
all}\;{\mathcal{I}}\in\Omega$ (2.3)

That is, $z$ stays within the range $z({\mathcal{I}}_{0})\pm\delta$ over
$\Omega$.

  4. 3.

_Infinitesimal Homeostasis_ at the point
${\mathcal{I}}_{0}\in\operatorname{dom}(z)$ if

$\nabla z({\mathcal{I}}_{0})=0$ (2.4)

That is, ${\mathcal{I}}_{0}$ is a _critical point_ of $z$.  $\Diamond$

It is clear that perfect homeostasis implies near-perfect homeostasis, but the
converse does not hold. Inspired by Nijhout, Reed, Best et al. [107, 18, 103],
Golubitsky and Stewart [54, 55] introduced the notion of infinitesimal
homeostasis which is intermediate between perfect and near-perfect
homeostasis. It is obvious that perfect homeostasis implies infinitesimal
homeostasis. On the other hand, it follows from Taylor’s theorem that
infinitesimal homeostasis implies near-perfect homeostasis in a neighborhood
of ${\mathcal{I}}_{0}$ (see [56] for details). It is easy to see that the
converse to both implications is not generally valid (see Reed et al. [115]).

The geometric interpretation of infinitesimal homeostasis is that, if
${\mathcal{I}}_{0}$ is a critical point of the input-output function, then
$z({\mathcal{I}})$ differs from $z({\mathcal{I}}_{0})$ in a manner that
depends quadratically (or to higher order) on
$\|{\mathcal{I}}-{\mathcal{I}}_{0}\|$. This makes the graph of
$z({\mathcal{I}})$ flatter than any growth rate with a nonzero linear term.

Definition 2.1 is very general in the sense that it allows for an arbitrary
number of input parameters and and arbitrary (smooth) function of state
variables as the output. In the engineering literature, a system of inputs and
outputs can be described as one of four types: SISO (single input, single
output), SIMO (single input, multiple output), MISO (multiple input, single
output), or MIMO (multiple input, multiple output). By this terminology,
Definition 2.1 describes MIMO systems. However, as in bifurcation theory, the
most important situation is when a single scalar input parameter is considered
at a time. Moreover, as we will define latter there is a natural class of
dynamical systems, called network dynamical systems or coupled cell systems,
that comes with a distinguished set of observables [56]. By adjusting the
notion of homeostasis with a single input parameter to work in the class of
network dynamical systems we obtain a very general description of SISO
systems, called input-output networks. Remarkably, it is in this class of
systems that we can obtain the most complete theory of homeostasis (see
Subsection 2.3).

Nevertheless, we can still say something in the most general setting. We can
use implicit differentiation and the cofactor formula for the inverse of a
matrix to write a formula for $\nabla z$. Let us denote the Jacobian matrix of
$F$ at the equilibrium $(X(\mathcal{I}),\mathcal{I})$ (we will omit the ~ over
$X$ from now on) by

$J=DF_{(X({\mathcal{I}}),{\mathcal{I}})}$

The partial derivatives of $F=(f_{1},\ldots,f_{n})$ with respect to $x_{i}$
are denoted by $f_{j,x_{i}}$ and the partial derivatives with respect to
${\mathcal{I}}_{M}$ are denoted by $f_{j,{\mathcal{I}}_{M}}$. Let
$F_{{\mathcal{I}}_{m}}$ denote the vector of partial derivatives of $F$ with
respect to ${\mathcal{I}}_{M}$, that is,
$F_{{\mathcal{I}}_{M}}=(f_{1,{\mathcal{I}}_{M}},\ldots,f_{n,{\mathcal{I}}_{M}})$.

######  Lemma 2.2.

Let $z(\mathcal{I})$ be an input-output function associated to a system of
differential equations (2.1). The gradient of $\nabla z({\mathcal{I}})$ with
respect to the multiple input ${\mathcal{I}}$ is given by

$\nabla z=\frac{-1}{\det
J}\bigg{(}\langle\nabla\phi\,,\operatorname{adj}(J)F_{{\mathcal{I}}_{M}}\rangle\bigg{)}_{M=1}^{N}$
(2.5)

where $\operatorname{adj}(J)$ is the _adjugate matrix_ or _classical adjoint_
of $J$ (defined as the transpose of the cofactor matrix of $J$). In
particular, ${\mathcal{I}}_{0}$ is an infinitesimal homeostasis point of $z$
if and only if

$\langle\nabla\phi\,,\operatorname{adj}(J)F_{{\mathcal{I}}_{M}}\rangle\big{|}_{(X({\mathcal{I}}_{0}),{\mathcal{I}}_{0})}=0\;,\quad\forall\;M=1,\ldots,N$
(2.6)

######  Proof 2.3.

The gradient of $z$ is

$\nabla z=\bigg{(}\frac{\partial
z}{\partial{\mathcal{I}}_{1}},\ldots,\frac{\partial
z}{\partial{\mathcal{I}}_{N}}\bigg{)}$ (2.7)

The partial derivative of $z(\mathcal{I})=\phi(X(\mathcal{I}))$ with respect
to ${\mathcal{I}}_{M}$ (for each $M=1,\ldots,N$) is

$\frac{\partial
z}{\partial{\mathcal{I}}_{M}}=\frac{\partial}{\partial{\mathcal{I}}_{j}}\phi(X({\mathcal{I}}))=\langle\nabla\phi\,,X_{{\mathcal{I}}_{M}}\rangle$
(2.8)

where $X_{{\mathcal{I}}_{M}}$ is the partial derivative of $X$ with respect to
${\mathcal{I}}_{M}$. Implicit differentiation of the equation
$F(X(\mathcal{I}),\mathcal{I})=0$, with respect to ${\mathcal{I}}_{j}$, yields
the linear system

$J\,X_{{\mathcal{I}}_{M}}=-F_{{\mathcal{I}}_{M}}$ (2.9)

From (2.9) and the fact that $X({\mathcal{I}})$ is assumed to be a linearly
stable equilibrium, for all ${\mathcal{I}}\in\Omega$, and hence $\det(J)\neq
0$ over $\Omega$, it follows that

$X_{{\mathcal{I}}_{M}}=-\frac{1}{\det
J}\operatorname{adj}(J)F_{{\mathcal{I}}_{M}}$ (2.10)

Therefore, substituting (2.10) into (2.8) and then into (2.7) we obtain (2.5).
∎

Lemma 2.2 allows us to draw some general remarks about the function $\nabla
z$. The first observation is that the term
$\operatorname{adj}(J)F_{{\mathcal{I}}_{M}}$ is a multivariate polynomial in
the partial derivatives $f_{j,x_{i}}$ and $f_{j,{\mathcal{I}}_{M}}$, linear in
$f_{j,{\mathcal{I}}_{M}}$. Hence, the defining equations for homeostasis (2.6)
are polynomial functions in the partial derivatives $f_{j,x_{i}}$,
$f_{j,{\mathcal{I}}_{M}}$ and $\phi_{x_{i}}$, linear in both
$f_{j,{\mathcal{I}}_{M}}$ and $\phi_{x_{i}}$. This will become manifest in the
several examples discussed in the paper.

The second observation is that there are two ‘trivial’ ways in which
homeostasis can occur: (i) by the existence of a critical point of $\phi$,
that is, $\nabla\phi=0$ at $(X({\mathcal{I}}_{0}),{\mathcal{I}}_{0})$, and
(ii) the simultaneous vanishing of all $F_{{\mathcal{I}}_{M}}$ at
$(X({\mathcal{I}}_{0}),{\mathcal{I}}_{0})$. The first case depends on the
choice of the observable $\phi$ and it indeed occurs in some applications, for
instance [107] (see Example 2.4 below). The second case is a highly non-
generic situation which essentially says that the system is homeostatic “at
the input source”. Henceforth, unless stated otherwise, we always assume the
following genericity condition: the derivative of the vector field $F$ with
respect to the input parameters generically satisfies

$F_{{\mathcal{I}}_{M}}\neq 0\;,\quad\forall\;M=1,\ldots,N$ (2.11)

The third observation is that the formula (2.5) is valid as long as the
equilibrium $X({\mathcal{I}})$ is linearly stable, since $\det(J)\neq 0$ at
$(X({\mathcal{I}}),{\mathcal{I}})$. Moreover, $\det(J)=0$ at some
$(X({\mathcal{I}}_{0}),{\mathcal{I}}_{0}))$ if and only if ${\mathcal{I}}_{0}$
is a steady-state bifurcation point. Therefore, the input-output function is
non-differentiable at the steady-state bifurcation points of
$X({\mathcal{I}})$. Although it may be continuously extended beyond those
points and become differentiable again.

The final observation is concerned with the fact that, very often, model
equations depend on several parameters (besides the input parameters). Because
of this dependence the equilibrium family $\tilde{X}({\mathcal{I}})$ and
input-output function may also depend on some of these parameters
(generically, they depend on all of them), but they are suppressed from
$z({\mathcal{I}})$ most the time.

More precisely, if we write the vector field in (2.1) as
$F(X,{\mathcal{I}},\alpha)$, where $\alpha=(\alpha_{1},\ldots,\alpha_{\ell})$
is the vector of all scalar parameters that appear in the definition of $F$
then the implicit function theorem applied to a point
$(X^{*},{\mathcal{I}}^{*},\alpha^{*})$ gives the function
$\tilde{X}({\mathcal{I}},\alpha)$ defined in a neighborhood of
$({\mathcal{I}}^{*},\alpha^{*})$ such that
$\tilde{X}({\mathcal{I}}^{*},\alpha^{*})=X^{*}$ and
$F(\tilde{X}({\mathcal{I}},\alpha),{\mathcal{I}},\alpha)\equiv 0$. Therefore,
the input-output function $z$ can be written as $z({\mathcal{I}},\alpha)$ and
infinitesimal homeostasis may depend on a fixed $\alpha_{0}$. That is,
infinitesimal homeostasis occurs at ${\mathcal{I}}_{0}$ only when
$\alpha=\alpha_{0}$.

It is here that notion of infinitesimal homeostasis places the study of
homeostasis in the context of singularity theory. For instance, the dependence
on the parameters $\alpha$ for the occurrence of infinitesimal homeostasis is
related to higher degeneracy conditions, in addition to (2.4), which leads to
distinct ‘forms’ of infinitesimal homeostasis.

For example, let us consider the single parameter case
${\mathcal{I}}\in\mathbb{R}$. We say that $z$ exhibits _simple homeostasis_ if

$z^{\prime}({\mathcal{I}}_{0})=0\quad\text{and}\quad
z^{\prime\prime}({\mathcal{I}}_{0})\neq 0$ (2.12)

In this case, ${\mathcal{I}}_{0}$ is called a Morse singularity (or a
nondegenerate critical point) and it is persistent under any small
perturbation of $\alpha$ (see subsection 2.2). In fact, a Morse singularity is
persistent under any small perturbation of $z$ in the space of smooth
functions (with the appropriate topology) [51]. We say that $z$ exhibits
_chair homeostasis_ if

$z^{\prime}({\mathcal{I}}_{0})=0\quad\text{and}\quad
z^{\prime\prime}({\mathcal{I}}_{0})=0\quad\text{and}\quad
z^{\prime\prime\prime}({\mathcal{I}}_{0})\neq 0$ (2.13)

In this case, there is at least one parameter $\alpha_{j}$ (or a combination
of several components of $\alpha$) that should be fixed together with
${\mathcal{I}}_{0}$ to have infinitesimal homeostasis. However, as we will see
below, singularity theory implies that small perturbations of $z$ (that is,
variation of the suppressed parameters) change the ‘almost flat region’ only
slightly. Following Nijhout et al. [103] we define a _plateau_ as a region of
$\Omega$ over which $z$ is approximately constant. In other words, the
infinitesimal homeostasis points are not persistent, but the family of
functions $z({\mathcal{I}},\alpha)$ have a plateau for all
$({\mathcal{I}},\alpha)$ in a neighborhood of $({\mathcal{I}}_{0},\alpha_{0})$
(see Definition 2.10 in subsection 2.2).

The observation that an unstable structure ‘determines’ (unfolds) all the
neighboring stable behaviors goes back to René Thom’s notion of a singularity
as an _organizing center_ [127]. Even though the singularity itself is
unstable under perturbation, it ‘organizes’ all its small perturbations into
the universal unfolding. Unlike the singularity, the universal unfolding is
structurally stable and thus contains all the possible persistent behaviors,
as well as the number of parameters to (qualitatively) model all those
behaviors.

This is important in applications because it may be very difficult to
calculate exactly the values of $({\mathcal{I}}_{0},\alpha_{0})$ where
infinitesimal homeostasis occurs and then one must resort to numerical
methods. It is possible to compute parts of the graph of input-output
functions and find the plateau using numerical methods for continuation of
equilibrium points, such as Auto from XPPAut [44] (for example, see Figure 7
in subsection 2.4 and Figure 25 in subsection 4.3).

Finally, we should mention that there is another notion of ‘robustness’ often
referred in the control theoretic literature, see [4, 125, 10, 9] for example.
This notion is closely related to other concepts such as _structural
identifiability_ [98, 99, 59, 60, 20] and _dynamical compensation_ [77, 124,
132].

####  2.1.1 Examples of Homeostasis

Here we present some examples of input-output functions associated to
biological models.

######  Example 2.4 (Folate Cycle).

Nijhout et al. [107] developed a model for the folate cycle based on standard
biochemical kinetics and used the model to provide new insights into several
different mechanisms of folate homeostasis. Their model is described by 6 ODEs
for the concentrations of the following substrates [107, Eqs. 4–9]: (1)
tetrahydrofolate (THF), (2) 5-methyltetrahydrofolate (5mTHF), (3)
dihydrofolate (DHF), (4) 5,10-methylenetetrahydrofolate (5,10-CH2-THF), (5)
5,10-methenyltetrahydro-folate (5,10-CH=THF) and (6) 10-formyltetrahydrofolate
(10f-THF). The coupling structure of the model is represented by a diagram
containing the 6 substrates and additional 12 enzymes [107, Fig. 1]. In the
model, folate enters and leaves the cell as 5mTHF as indicated by the presence
of the input parameter
$\mathcal{I}=\mathcal{F}_{\mathrm{in}}-\mathcal{F}_{\mathrm{out}}$ in the
corresponding differential equation Here, $\mathcal{F}_{\mathrm{in}}$ and
$\mathcal{F}_{\mathrm{out}}$ are the rates at which 5mTHF enters and leaves
the cell, respectively. The authors consider observables of the form
$V_{\mathrm{ENZYME}}=\phi(S,F)$ representing the velocity of a reaction
associated to ENZYME, called ‘fluxes’ [107, Tab. III], in terms of substrates.
The smooth functions $\phi$ are given by certain expressions involving
Michaelis-Menten functions [107, Eqs. 1–3]. The plots of the fluxes as
functions of the input parameter were numerically computed and homeostasis is
shown by the fact that velocities do not decline substantially until total
folate is close to $0$ [107, Figs. 4–6]. It can be checked that the functions
$\phi$ have a critical point completely determined by the kinetic parameters
entering its definition.  $\Diamond$

######  Example 2.5 (Kidney Flow).

Sgouralis and Layton [119] show by using a mathematical model how well the
effect of the myogenic mechanism affects the afferent arteriole flow rate
(output) in response to pressure variation (input). The myogenic response in
the smooth muscle enables nephrons to regulate the flow in afferent arterioles
for a wide range of pressure values The solid curves in [119, Fig. 4] (see
also [103, Fig. 6]) shows the percentage change at steady state in their model
as the pressure in the afferent arteriole is varied from $50\,mm\mathrm{Hg}$
to $230\,mm\mathrm{Hg}$. The flow rate is remarkably stable between
$90\,mm\mathrm{Hg}$ and $190\,mm\mathrm{Hg}$. Below $90\,mm\mathrm{Hg}$ the
flow drops quickly and above $190\,mm\mathrm{Hg}$ the flow begins to increase
substantially giving the same chair-shaped curve that we have seen in the
Introduction. The dashed curves shows what the flow rate would be if the
myogenic response were turned off and the walls of the arteriole responded
passively.  $\Diamond$

######  Example 2.6 (Glycolysis Pathway).

Mulukutla et al. [102] use a mathematical model based on reported mechanisms
for the allosteric regulations of the enzymes in the flux of glycolysis. They
show that glycolysis exhibits multiple steady state behavior segregating
glucose metabolism into high flux and low flux states. Here, the input
parameter is the glucose consumption and the output variable is the specific
lactate production. In [102, Fig. 3] the authors show a phenomenon of
bistability of homeostasis in cultured HeLa cells (see also [41, Fig. 4]). The
data suggests the existence of two homeostasis points on the lower branch: one
which is apparent in the figure, and another which we would expect to see if
it were extended further. The similar glucose consumption rates of both types
of cells in very high and very low glucose environments indicate two switches
on the border of the plateaus. The homeostasis points and the bistable
behavior are suggestive of the behavior depicted [102, Fig. 2] obtained by
simulation of the model equations. Duncan and Golubitsky [41, Fig. 3] propose
a mechanism combining homeostasis and bifurcation that qualitatively
reproduces the bistability observed in [102].  $\Diamond$

######  Example 2.7 (Extracellular Dopamine).

Best et al. [18, 103] propose a biological network to study homeostasis of
extracellular dopamine (eDA) in response to variation in the activities of the
enzyme tyrosine hydroxylase (TH) and the dopamine transporters (DAT). Hence,
in this model we have two input parameters (eDA,TH) and one output variable
DAT. The authors derive a differential equation model for the biochemical
network [18, Fig. 1]. In [103] the authors fix reasonable values for all
parameters in the model with the exception of the concentrations of TH and
DAT. In [103, Fig. 8] the authors show the equilibrium value of eDA as a
function of TH and DAT in their model. The white dots indicate the predicted
eDA values for the observationally determined values of TH and DAT in the
wild-type genotype (large white disk) and the polymorphisms observed in human
populations (small white disks). These polymorphisms raise or lower the
activity of TH and lower the activity of the DAT as indicated in [103, Tab.
1]. This result is scientifically important because almost all of these
genetic polymorphisms (white disks) lie on the plateau and that indicates
homeostasis of eDA. Presumably, this is because polymorphisms that would
result in large changes in eDA are deleterious and would, therefore, be
removed from the population. Thus, although the polymorphisms have large
effects on the activities of TH and DAT, they have only a very small effect on
the phenotypic variable eDA, and thus can be considered to be cryptic genetic
variation. Note that the plateau contains a line from left to right at about
eDA $=0.9$. In this respect, the surface graph in [103, Fig. 8] appears to
resemble that of a non-singular perturbed hyperbolic umbilic (see Table 1) and
Figure 2. See also the level contours of the hyperbolic umbilic in Figure 2.
This figure shows that the hyperbolic umbilic is the only low-codimension
singularity that contains a single line in its zero set, see [55] for more
details.  $\Diamond$

###  2.2 Singularity Theory of Input-Output Functions: Elementary Catastrophe
Theory

As discussed in the Introduction, Nijhout et al. [103] observe that
homeostasis appears in many applications through the notion of a chair.
Golubitsky and Stewart [55] observed that a chair can be thought of as a
singularity of a scalar input-output function, one where $z({\mathcal{I}})$
‘looks like’ a homogeneous cubic, e.g.
$z({\mathcal{I}})\approx{\mathcal{I}}^{3}$. More precisely, the mathematics of
singularity theory [112, 50] replaces ‘looks like’ by ‘up to a change of
coordinates.’

######  Definition 2.8.

A smooth function $z:\mathbb{R}^{N}\to\mathbb{R}$ is singular at a point
${\mathcal{I}}_{0}\in\mathbb{R}^{N}$ if

$\nabla z({\mathcal{I}}_{0})=0$ (2.14)

A point ${\mathcal{I}}_{0}$ satisfying (2.14) is called a singularity, or
critical point, of $z$. A point ${\mathcal{I}}_{0}$ is a nondegenerate
singularity of $z$ if, in addition to (2.14), its Hessian matrix is
invertible. Otherwise, ${\mathcal{I}}_{0}$ is a called a degenerate
singularity.  $\Diamond$

######  Definition 2.9.

Two smooth functions $z,w:\mathbb{R}^{N}\to\mathbb{R}$ are right equivalent on
a neighborhood of ${\mathcal{I}}_{0}\in\mathbb{R}^{N}$ if

$w({\mathcal{I}})=z(\Lambda({\mathcal{I}}))+K$

where $\Lambda:\mathbb{R}^{N}\to\mathbb{R}^{N}$ is an _invertible change of
coordinates_ , or simply a _diffeomorphism_ , on a neighborhood of
${\mathcal{I}}_{0}$ and $K\in\mathbb{R}^{N}$ is a constant.  $\Diamond$

The simplest classification theorem of singularity theory states that
$z:\mathbb{R}^{N}\to\mathbb{R}$ is right equivalent to
$w({\mathcal{I}}_{1},\ldots,{\mathcal{I}}_{N})=\sum_{i=1}^{N}\pm{\mathcal{I}}_{i}^{2}$
on a neighborhood of a singularity if and only if the singularity is
nondegenerate. This is the content of the classical Morse Lemma.
Geometrically, this says that the graph of $z$ ‘looks like’ a quadratic
function near a nondegenerate singularity. Now, the the main goal of
singularity theory is the study of degenerate singularities.

The transformations of the input-output map $z({\mathcal{I}})$ given in
Definition 2.9 are just the standard change of coordinates in elementary
catastrophe theory [50, 112, 142]. In [54, 55] it is shown that a special
class of transformations on the vector field induces the class of right
equivalences on the input-output functions. Although their proof is formulated
for a special class of observables (the projections onto the coordinate
variables) the use of right equivalence is completely justified. The proof
that the action of the diffeomorphism $\Lambda$ on the vector field induces
the correct action on the input-output function works in the general setting
and the addition of the constant $K$ is justified because we are looking for
plateaus on which $z$ is approximately constant but not specifically what that
constant is.

Now we can therefore use standard results from elementary catastrophe theory
to find normal forms and universal unfoldings of $z$, as we now explain.
Informally, the codimension of a singularity is the number of conditions on
derivatives that determine it. This is also the minimum number of extra
variables required to specify all small perturbations of the singularity, up
to changes of coordinates. These perturbations can be organized into a family
of maps called the universal unfolding, which has that number of extra
variables.

######  Definition 2.10.

A smooth function $W({\mathcal{I}},a)$, $a\in\mathbb{R}^{\ell}$, is an
unfolding of $w({\mathcal{I}})$ if $W({\mathcal{I}},0)=w({\mathcal{I}})$.
Here, $a\in\mathbb{R}^{\ell}$ is called unfolding parameter. An unfolding $W$
of $w$ is called an universal unfolding if every unfolding
$V({\mathcal{I}},b)$, $b\in\mathbb{R}^{m}$, of $w$ factors through $W$. That
is,

$V({\mathcal{I}},b)=W(\Lambda({\mathcal{I}},b),A(b))+K(b)$ (2.15)

where $A,K:\mathbb{R}^{m}\to\mathbb{R}^{\ell}$ and
$\Lambda:\mathbb{R}^{N}\times\mathbb{R}^{m}\to\mathbb{R}^{N}$.  $\Diamond$

It follows that every small perturbation $V(\cdot,b)$ is equivalent to a
perturbation $W(\cdot,A(b))$ of $w$ in the $W$ family.

######  Example 2.11 (Single Input Parameter).

Let $z:\mathbb{R}\to\mathbb{R}$ be singular at the origin. Because
$z({\mathcal{I}})$ is $1$-dimensional, we consider singularity types near the
origin of a $1$-variable function $w({\mathcal{I}})$. Such singularities are
determined by the first non-vanishing ${\mathcal{I}}$-derivative $w^{(k)}(0)$
(unless all derivatives vanish, which is an ‘infinite codimension’ phenomenon
that we do not discuss further). If such $k$ exists, the normal form is
$\pm{\mathcal{I}}^{k}$. When $k\geq 3$ the universal unfolding for catastrophe
theory equivalence is

$\pm{\mathcal{I}}^{k}+a_{k-2}{\mathcal{I}}^{k-2}+a_{k-3}{\mathcal{I}}^{k-3}+\cdots+a_{1}{\mathcal{I}}$

for parameters $a=(a_{j})$ and when $k=2$ (a nondegenerate singularity) the
universal unfolding is $\pm{\mathcal{I}}^{2}$. The codimension $\ell$ in this
setting is therefore $\ell=k-2$. To summarize: the normal form of the input-
output function near a nondegenerate singularity (codimension $0$) is

$z({\mathcal{I}})=\pm{\mathcal{I}}^{2}$ (2.16)

and no unfolding parameter is required. Here, ${\mathcal{I}}_{0}=0$ is a
simple homeostasis point. Similarly,

$z({\mathcal{I}})=\pm{\mathcal{I}}^{3}$ (2.17)

is the normal form of the input-output function near a least degenerate
singular point (codimension $1$), and

$z({\mathcal{I}};a)=\pm{\mathcal{I}}^{3}+a{\mathcal{I}}$ (2.18)

is a universal unfolding. The geometry of this universal unfolding is shown in
Figure 1. Here, the singularity occurs when $({\mathcal{I}}_{0},a_{0})=(0,0)$
and it is a chair homeostasis point, Figure 1(b). When $a\neq 0$ there is no
singularity at ${\mathcal{I}}_{0}=0$ but the plateau persists all $a$ near
$a_{0}=0$. The difference is that when $a<0$ the plateau has two nondegenerate
singular points, Figure 1(a) and when $a>0$ the plateau has no singular
points, Figure 1(c). See [26] for details.  $\Diamond$

(a) $a<0$

(b) $a=0$

(c) $a>0$

Figure 1:  Input-output function exhibiting chair homeostasis at some point
${\mathcal{I}}_{0}$. The corresponding to the universal unfolding is
$z({\mathcal{I}};a)=-{\mathcal{I}}^{3}+a{\mathcal{I}}$, up to right
equivalence. (a) the unfolding parameter $a<0$ and the degenerate singularity
unfolds into two nondegenerate singularities; (b) the unfolding parameter
$a=0$ and the singularity is degenerate; (c) the unfolding parameter $a>0$ and
the degenerate singularity unfolds into no singular points.

######  Example 2.12 (Multiple Input Parameters).

We will consider the $N=2$ case. Table 1 summarizes the classification when
$N=2$, so
${\mathcal{I}}=({\mathcal{I}}_{1},{\mathcal{I}}_{2})\in\mathbb{R}^{2}$. Here
the list is restricted to codimension $\leq 3$. The associated geometry,
especially for universal unfoldings, is described in [22, 49, 112] up to
codimension 4. Singularities of much higher codimension have also been
classified, but the complexities increase considerably. For example, Arnold
[12] provides an extensive classification up to codimension 10 (for the
complex analog). Note that in the $N=2$ case the normal forms for $N=1$ appear
again, but now there is an extra quadratic term $\pm{\mathcal{I}}_{2}^{2}$.
This term is a consequence of the _s_ plitting lemma in singularity theory,
arising here when the Hessian of $z$ has rank 1 rather than rank 0 (corank 1
rather than corank [22, 112, 142]. The presence of the
$\pm{\mathcal{I}}_{2}^{2}$ term affects the range over which
$z({\mathcal{I}})$ changes when ${\mathcal{I}}_{2}$ varies, but not when
${\mathcal{I}}_{1}$ varies.  $\Diamond$

The standard geometric features considered in catastrophe theory focus on the
gradient of the function $z({\mathcal{I}})$ in normal form. In contrast, what
matters here is the function itself. Specifically, we are interested in the
region in the ${\mathcal{I}}$-space where the function $z$ is approximately
constant, i.e. the plateau of $z$.

More specifically, for each normal form $z({\mathcal{I}})$ we choose a small
$\delta>0$ and form the set

$P_{\delta}=\\{{\mathcal{I}}\in\Omega:\|z({\mathcal{I}})\|\leq\delta\\}.$
(2.19)

This is the $\delta$-_plateau_ region on which $z({\mathcal{I}})$ is
approximately constant, where $\delta$ specifies how good the approximation
is. Even though the definition of infinitesimal homeostasis is qualitative it
is possible to extract quantitative bounds on the size of the plateau using
Taylor’s theorem with remainder. See [56, Chap. 5] for details.

As we mentioned before, if $z$ is perturbed slightly (by variation of the
suppressed parameters), the plateau $P_{\delta}$ varies continuously.
Therefore, one can compute the approximate plateau by focusing on the
singularity, rather than on its universal unfolding. In other words, for
sufficiently small perturbations plateaus of singularities depend mainly on
the singularity itself and not on its universal unfolding.

This observation is important because the universal unfolding has many zeros
of the gradient of $z({\mathcal{I}})$, hence ‘homeostasis points’ near which
the value of $z({\mathcal{I}})$ varies more slowly than linear. However, this
structure seems less important when considering the relationship of
infinitesimal homeostasis with homeostasis. See the discussion of the
unfolding of the chair summarized in [52, Figure 3].

The ‘qualitative’ geometry of the plateau – that is, its differential topology
and associated invariants – is characteristic of the singularity. This offers
one way to infer the probable type of singularity from numerical data; it also
provides information about the region in which the system concerned is
behaving homeostatically. We do not develop a formal list of invariants here,
but we indicate a few possibilities.

The main features of the plateaus associated with the six normal forms of
Table 1 are illustrated in Figure 2. Figure plots, for each normal form, a
sequence of contours from $-\delta$ to $\delta$; the union is a picture of the
plateaus. By unfolding theory, these features are preserved by small
perturbations of the model, and by the choice of $\delta$ in (2.19) provided
it is sufficiently small. Graphical plots of such perturbations (not shown)
confirm this assertion.

Table 1: Classification of singularities of input-output maps $\mathbb{R}^{2}\rightarrow\mathbb{R}$ of codimension $\leq 3$. name | normal form | codim | universal unfolding
---|---|---|---
Morse (simple±) | $\pm{\mathcal{I}}_{1}^{2}\pm{\mathcal{I}}_{2}^{2}$ | 0 | $\pm{\mathcal{I}}_{1}^{2}\pm{\mathcal{I}}_{2}^{2}$
fold (chair) | ${\mathcal{I}}_{1}^{3}\pm{\mathcal{I}}_{2}^{2}$ | 1 | ${\mathcal{I}}_{1}^{3}+a{\mathcal{I}}_{1}\pm{\mathcal{I}}_{2}^{2}$
cusp± | $\pm{\mathcal{I}}_{1}^{4}\pm{\mathcal{I}}_{2}^{2}$ | 2 | $\pm{\mathcal{I}}_{1}^{4}+a{\mathcal{I}}_{1}^{2}+b{\mathcal{I}}_{1}\pm{\mathcal{I}}_{2}^{2}$
swallowtail | ${\mathcal{I}}_{1}^{5}\pm{\mathcal{I}}_{2}^{2}$ | 3 | ${\mathcal{I}}_{1}^{5}+a{\mathcal{I}}_{1}^{3}+b{\mathcal{I}}_{1}^{2}+c{\mathcal{I}}_{1}\pm{\mathcal{I}}_{2}^{2}$
hyperbolic umbilic | ${\mathcal{I}}_{1}^{3}+{\mathcal{I}}_{2}^{3}$ | 3 | ${\mathcal{I}}_{1}^{3}+{\mathcal{I}}_{2}^{3}+a{\mathcal{I}}_{1}{\mathcal{I}}_{2}+b{\mathcal{I}}_{1}+c{\mathcal{I}}_{2}$
elliptic umbilic | ${\mathcal{I}}_{1}^{3}-3{\mathcal{I}}_{1}{\mathcal{I}}_{2}^{2}$ | 3 | ${\mathcal{I}}_{1}^{3}-3{\mathcal{I}}_{1}{\mathcal{I}}_{2}^{2}+a({\mathcal{I}}_{1}^{2}+{\mathcal{I}}_{2}^{2})+b{\mathcal{I}}_{1}+c{\mathcal{I}}_{2}$

simple+ simple- cusp+ cusp-

chair swallowtail hyperbolic elliptic

Figure 2: Plateaus shown by contour plots for each singularity in Table 1.
Reproduced from [55, Fig. 4]. Each plot has 200 equally spaced contour levels
for $\delta$ from $-0.2$ to $0.2$.

The notion of infinitesimal homeostasis point can be generalized in several
directions. For instance, Andrade et al. [3] consider the situation where the
infinitesimal homeostasis point ${\mathcal{I}}_{0}$ is at the boundary of the
domain of the input-output function, which may be finite or infinity. They
introduce the notion of _asymptotic infinitesimal homeostasis_ for input-
output functions $z$ defined on a domain of the form
$D=[{\mathcal{I}}^{-},+\infty)$, for example. In this case we say that $z$
exhibits _asymptotic infinitesimal homeostasis_ at $+\infty$ if

$\lim_{\mathcal{I}\to\infty}z^{\prime}(\mathcal{I})=0.$

In [3] the authors use some simple criteria for the existence of asymptotic
infinitesimal homeostasis to study a model for copper homeostasis (see
subsection 4.2).

######  Theorem 2.13.

Let $z:D\to\mathbb{R}$, with $D=(\mathcal{I}^{-},+\infty)$ be a smooth
function.

  1. [(1)]

  2. 1.

If $z$ exhibits near-perfect homeostasis on $D$ (i.e., it is upper and lower
bounded on $D$) and is monotonic then $z$ exhibits asymptotic infinitesimal
homeostasis.

  3. 2.

(Hadamard’s Lemma [61]) If the first $n$ derivatives
$z^{\prime},\ldots,z^{(n)}$ of $z$ are bounded, then
$\displaystyle\lim_{{\mathcal{I}}\to+\infty}z({\mathcal{I}})=c$ with
$c\in\mathbb{R}$ implies that

$\lim_{{\mathcal{I}}\to+\infty}z^{(j)}({\mathcal{I}})=0$

for all $1\leq j\leq n$.

It is interesting to mention that asymptotic infinitesimal homeostasis is
related to asymptotic uniformity [48] and singularities at infinity [129].

###  2.3 Homeostasis in Input-Output Networks

In this subsection we introduce input-output networks. This concept is crucial
for the development of a rich combinatorial theory.

######  Definition 2.14.

An _input-output network_ is a directed graph $\mathcal{G}$ with a
distinguished input node $\iota$ and a distinguished output node $o$. The
network $\mathcal{G}$ is a _core network_ if every node in $\mathcal{G}$ is
downstream from the input node $\iota$ and upstream from the output $o$.
$\Diamond$

An admissible system of differential equations associated with $\mathcal{G}$
has the form

$\displaystyle\begin{split}\dot{x}_{\iota}&=f_{\iota}(x_{\iota},x_{\rho},x_{o},{\mathcal{I}})\\\
\dot{x}_{\rho}&=f_{\rho}(x_{\iota},x_{\rho},x_{o})\\\
\dot{x}_{o}&=f_{o}(x_{\iota},x_{\rho},x_{o})\end{split}$ (2.20)

where ${\mathcal{I}}\in\mathbb{R}$ is the input parameter,
$X=(x_{\iota},x_{\rho},x_{o})\in\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}$
is the vector of state variables associated to nodes in $\mathcal{G}$, and
$F(X,{\mathcal{I}})=(f_{\iota}(X,{\mathcal{I}}),f_{\rho}(X),f_{o}(X))$ is a
smooth family of mappings on the state space
$\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}$. Note that ${\mathcal{I}}$
appears only in the equation of system (2.20) corresponding to the input node.

We can write (2.20) as

$\dot{X}=F(X,{\mathcal{I}}).$

We denote the partial derivative of the function associated to node $j$ with
respect to the state variable associated to node $k$, also called _linearized
couplings_ , by

$f_{j,x_{k}}=\frac{\partial}{\partial x_{k}}f_{j}$

We assume $f_{j,k}\equiv 0$ precisely when no arrow connects node $k$ to node
$j$. That is, $f_{j}$ is independent of $x_{k}$ when there is no arrow $k\to
j$. This is a modeling assumption made in $\mathcal{G}$. Finally, we assume
the _genericity condition_ (2.11), which in this case reads

$f_{\iota,{\mathcal{I}}}\neq 0$

generically.

Suppose $\dot{X}=F(X,{\mathcal{I}}_{0})$ has a hyperbolic equilibrium at
$X^{*}$ at ${\mathcal{I}}^{*}$. As before, we can define an input-output
function ${\mathcal{I}}\mapsto x_{o}({\mathcal{I}})$ defined on a neighborhood
of ${\mathcal{I}}^{*}$.

######  Remark 2.15 (SISO Systems).

Each node in an input-output network $\mathcal{G}$ corresponds to a one-
dimensional state variable of the admissible system. In particular, the output
node corresponds to a scalar quantity and the input parameter is a scalar
quantity. The class of systems described by input-output networks can be seen
as the class of single input, single output (SISO) systems. As we will see
later, it is possible to extend this definition to include input-output
networks with multiple input nodes, but a single input parameter and single
output, and multiple inputs, single output (MISO).  $\Diamond$

######  Remark 2.16 (Input Node $\neq$ Output Node).

In Definition 2.14 we explicitly exclude the possibility that the output node
is one of the input nodes. For now, this assumption is taken into account
purely for the sake of convenience. In fact, all the results should be valid
when the input $=$ output, but then all the theorems and proofs should be
properly adapted to take this particular case into account. This possibility
will be considered in greater detail in Subsection 4.2 $\Diamond$

######  Remark 2.17 (Multiple Outputs).

It is important to clarify why the input-output function $z$ should be always
real-valued ($1$-dimensional). This seems to be the correct way to formulate
the notion of simultaneously tracking ‘multiple outputs’. Even in the control
theoretic context of single input, multiple output (SIMO) and multiple input,
multiple output (MIMO) the observable $\phi$ is real-valued (typically it is a
linear functional). Actually, there is a singularity theoretic reason for this
requirement, as well.

Let us try to define a two-output input-output function by considering a
network with one input node $x_{\iota}$ and two output nodes $x_{o_{1}}$ and
$x_{o_{2}}$. Then the input-output function would be a mapping
$\mathbb{R}\to\mathbb{R}\times\mathbb{R}$ given by
${\mathcal{I}}\to\big{(}x_{o_{1}}({\mathcal{I}}),x_{o_{2}}({\mathcal{I}})\big{)}$.
This, in turn, would force us to consider more general coordinate changes on
the target space $\mathbb{R}\times\mathbb{R}$ (recall that in the real-valued
situation the action of the right equivalence on the target is by a
translation by a constant). Even in this simple two-output case one is led to
consider general _network preserving coordinate changes_ [53, 8]. Typically,
the most natural coordinate change on $\mathbb{R}\times\mathbb{R}$ that
preserves the network structure is by a diagonal diffeomorphism
$\varphi:\mathbb{R}\times\mathbb{R}\to\mathbb{R}\times\mathbb{R}$, given by
$\varphi(x_{1},x_{2})=(\varphi_{1}(x_{1}),\varphi_{2}(x_{2}))$, where
$\varphi_{j}:\mathbb{R}\to\mathbb{R}$ are diffeomorphisms. The set of smooth
mappings of the form $\mathbb{R}\to\mathbb{R}\times\mathbb{R}$, with the
equivalence relation generated by a pair of diffeomorphisms $(\psi,\varphi)$,
where $\psi$ acts on the source $\mathbb{R}$ and $\varphi$ is a diagonal
diffeomorphism acting on the target $\mathbb{R}\times\mathbb{R}$ is called the
set of _divergent diagrams_ [95]. The classification of singularities of
divergent diagrams have been obtained and the conclusion is that all
‘degenerate’ singularities have _infinite codimension_ [36, 37, 38, 95].
Therefore, we would not be able to define an analogue of chair homeostasis
(nor any other higher order singularity) in this setting.  $\Diamond$

Wang et al. [134] obtain a formula for the derivative of the input output
function $x_{o}({\mathcal{I}})$, with respect to ${\mathcal{I}}$, that is a
particular case of formula (2.5) for input-output networks. Here, the
observable is the projection onto the output coordinate. However, in this case
the formula for $x_{o}^{\prime}({\mathcal{I}})$ is much simpler than the
general case and provides a straightforward criterion for the occurrence of
homeostasis in an input-output network.

The _homeostasis matrix_ of (2.20) is the $(n+1)\times(n+1)$ matrix $H$ is
obtained from the $(n+2)\times(n+2)$ Jacobian matrix $J$ of (2.20) by deleting
its first row and last column. Indeed

$\displaystyle J=\begin{pmatrix}f_{\iota,\iota}&f_{\iota,\rho}&f_{\iota,o}\\\
f_{\rho,\iota}&f_{\rho,\rho}&f_{\rho,o}\\\
f_{o,\iota}&f_{o,\rho}&f_{o,o}\end{pmatrix}\Longrightarrow
H=\begin{pmatrix}f_{\rho,\iota}&f_{\rho,\rho}\\\
f_{o,\iota}&f_{o,\rho}\end{pmatrix}$

where $J$ and $H$ are both functions of $(X({\mathcal{I}}),{\mathcal{I}})$ as
in (2.20). More precisely:

######  Theorem 2.18 ([134, lemma 1.5]).

The derivative of the input-output function $x_{o}({\mathcal{I}})$, with
respect to ${\mathcal{I}}$ is given by

$x_{o}^{\prime}({\mathcal{I}})=\pm
f_{\iota,{\mathcal{I}}}\frac{\det(H)}{\det(J)}$ (2.21)

Therefore, (2.20) undergoes infinitesimal homeostasis at ${\mathcal{I}}_{0}$
if and only if $\det(H)=0$, evaluated at $(X_{0},{\mathcal{I}}_{0})$.

Homeostasis in a given network $\mathcal{G}$ can be determined by analyzing a
simpler network that is obtained by eliminating certain nodes and arrows from
$\mathcal{G}$. We call the network formed by the remaining nodes and arrows
the core subnetwork.

######  Definition 2.19.

A node $\tau$ in a network $\mathcal{G}$ is _downstream_ from a node $\rho$ in
$\mathcal{G}$ if there exists a path in $\mathcal{G}$ from $\rho$ to $\tau$.
Node $\rho$ is _upstream_ from node $\tau$ if $\tau$ is downstream from
$\rho$.

These relationships are important when trying to classify infinitesimal
homeostasis. For example, if the output node $o$ is not downstream from the
input node $\iota$, then the input-output function $x_{o}({\mathcal{I}})$ is
identically constant in ${\mathcal{I}}$. Although technically this is a form
of infinitesimal homeostasis, it is an uninteresting form.

######  Definition 2.20.

Let $\mathcal{G}$ be an input-output network.

  1. [(a)]

  2. 1.

The input-output network is a core network if every node is both upstream from
the output node and downstream from the input node.

  3. 2.

Every input-output network $\mathcal{G}$ has a core subnetwork
$\mathcal{G}_{c}$ whose nodes are the nodes in $\mathcal{G}$ that are both
upstream from the output node and downstream from the input node and whose
arrows are the arrows in $\mathcal{G}$ whose head and tail nodes are both
nodes in $\mathcal{G}_{c}$.

######  Theorem 2.21 ([134, Thm. 2.4]).

Let $\mathcal{G}$ be an input-output network and let $\mathcal{G}_{c}$ be the
associated core subnetwork. The input-output function associated with
$\mathcal{G}_{c}$ has a point of infinitesimal homeostasis at
${\mathcal{I}}_{0}$ if and only if the input-output function associated with
$\mathcal{G}$ has a point of infinitesimal homeostasis at ${\mathcal{I}}_{0}$.

It follows from Corollary 2.21 that classifying infinitesimal homeostasis for
networks $\mathcal{G}$ is equivalent to classifying infinitesimal homeostasis
for the core subnetwork $\mathcal{G}_{c}$.

######  Definition 2.22.

Let $\mathcal{G}_{1}$ and $\mathcal{G}_{2}$ be two input-output networks with
the same number of nodes.

  1. [(a)]

  2. 1.

We say that $\mathcal{G}_{1}$ and $\mathcal{G}_{2}$ are _core equivalent_ if
the determinants of their homeostasis matrices are identical as polynomials in
the variables $f_{j,x_{k}}$.

  3. 2.

A backward arrow is an arrow whose head is the input node $\iota$ or whose
tail is the output node $o$.

######  Corollary 2.23.

If two core networks differ from each other by the presence or absence of
backward arrows, then the core networks are core equivalent.

######  Proof 2.24.

The linearized couplings associated to backward arrows are of the form
$f_{\iota,x_{k}}$ and $f_{k,x_{o}}$, which do not appear in the homeostasis
matrix $H$. ∎

Therefore, backward arrows can be ignored when computing infinitesimal
homeostasis with the homeostasis matrix $H$. However, backward arrows cannot
be totally ignored, since they are involved in the determination of both the
equilibria of (2.20) and their stability.

Corollary 2.23 can be generalized to a theorem giving necessary and sufficient
graph theoretic conditions for core equivalence [134, Theorem 3.3].

Before proceeding to the classification results, we provide context for our
results by looking at some of the biochemical models discussed by Reed in
[115]. In doing so we show that input-output networks form a natural category
in which homeostasis may be explored.

###  2.4 $3$-node Input-Output Networks: Biochemical Networks

There are many examples of biochemical networks in the literature. In
particular examples, modelers decide which substrates are important and how
the various substrates interact. The network resulting from the detailed
modeling of the production of extracellular dopamine (eDA) by Best et al. [18]
and Nijhout et al. [103]. These authors derive a differential equation model
for this biochemical network and use the results to study homeostasis of eDA
with respect to variation of the enzyme tyrosine hydroxylase (TH) and the
dopamine transporters (DAT).

In another direction, relatively small biochemical network models are often
derived to help analyze a particular biochemical phenomenon. We present four
examples; three are discussed in Reed et al. [115] and one in Ma et al. [89].
These examples belong to a class that we call biochemical input-output
networks and will help to interpret the mathematical results.

######  Example 2.25 (Feedforward Excitation).

The input-output network corresponding to feedforward excitation is in Figure
3. This motif occurs in a biochemical network when a substrate activates the
enzyme that removes a product. The standard biochemical network diagram for
this process is shown in Figure 3(a). Here $\mathbf{X}$, $\mathbf{Y}$,
$\mathbf{Z}$ are the names of chemical substrates and their concentrations are
denoted by lower case $x$, $y$, $z$. Each straight arrow represents a flux
coming into or going away from a substrate. The differential equations for
each substrate simply state that the rate of change of the concentration is
the sum of the arrows going towards the substrate minus the arrows going away
(conservation of mass). The curved line indicates that substrate is activating
an enzyme.

(a) FE motif from Reed et al. [115]

(b) Input-output network derived from (a)

Figure 3: Feedforward excitation (FE).

Both diagrams in Figure 3 represent the same information, but in different
ways. Figure 3(b) uses nodes to represent variables, and arrows to represent
couplings. In other areas, conventions can differ, so it is necessary to
translate between the two representations. The simplest method is to write
down the model ODEs.

The equations corresponding to the biochemical network in Figure 3 (a) are as
follows

$\begin{split}\dot{x}&={\mathcal{I}}-g_{1}(x)-g_{4}(x)\\\
\dot{y}&=g_{1}(x)-g_{2}(y)-g_{5}(y)\\\
\dot{z}&=g_{2}(y)-f(x)g_{3}(z)\end{split}$ (2.22)

In this motif, one path consists of two excitatory couplings: $g_{1}(x)>0$
from $\bf X$ to $\bf Y$ and $g_{2}(y)>0$ from $\bf Y$ to $\bf Z$. The other
path is an excitatory coupling $f(x)>0$ from $\bf X$ to the synthesis or
degradation $g_{3}(z)$ of $\bf Z$ and hence is an inhibitory path from $\bf X$
to $\bf Z$ having a negative sign.

It is shown in [115] (and reproduced using the abstract theory in [58]) that
the model equations (2.22) for feedforward excitation leads to infinitesimal
homeostasis at $X_{0}$ if

$f_{x}(x_{0})=\frac{g_{1}^{\prime}(x_{0})g_{2}^{\prime}(y_{0})}{g_{3}(z_{0})(g_{2}^{\prime}(y_{0})+g_{5}^{\prime}(y_{0}))}$
(2.23)

where $X_{0}=(x_{0},y_{0},z_{0})$ is a stable equilibrium.

Figure 3(b) redraws the diagram in Figure 3(a) using the abstract definition
of input-output network, together with some extra features (the signs along
the arrows) that are crucial to this particular application. We consider $x$
to be a distinguished input variable, with $z$ as a distinguished output
variable, while $y$ is an intermediate regulatory variable. Accordingly we
change notation and write $x_{\iota}=x$, $x_{\rho}=y$ and $x_{o}=z$.

The general form of the equations associated to the diagram of Figure 3(b) in
the state variables $x_{\iota},x_{\rho},x_{o}$ is

$\begin{split}\dot{x}_{\iota}&=f_{\iota}(x_{\iota},{\mathcal{I}})\\\
\dot{x}_{\rho}&=f_{\rho}(x_{\iota},x_{\rho})\\\
\dot{x}_{o}&=f_{o}(x_{\iota},x_{\rho},x_{o})\end{split}$ (2.24)

In Figure 3(b) these variables are associated with three nodes $\iota,\rho,o$.
Each node has its own symbol: a square for $\iota$, circle for $\rho$, and
triangle for $o$. Here these symbols are convenient ways to show which type of
variable (input, regulatory, output) the node corresponds to. Arrows indicate
that the variables corresponding to the tail node occur in the component of
the ODE corresponding to the head node. For example, the component for
$\dot{x}_{o}$ is a function of $x_{\iota}$, $x_{\rho}$, and $x_{o}$. We
therefore draw an arrow from $\iota$ to $o$ and an arrow from $\rho$ to $o$.
We do not draw an arrow from $o$ to itself, however: by convention, every node
variable can appear in the component for that node. In a sense, the node
symbol (circle) represents this ‘internal’ arrow.

The Jacobian of (2.24) is

$J=\left(\begin{array}[]{cc;{2pt/2pt}c}f_{\iota,x_{\iota}}&0&0\\\
\hdashline[2pt/2pt]f_{\rho,x_{\iota}}&f_{\rho,x_{\rho}}&0\\\
f_{o,x_{\iota}}&f_{o,x_{\rho}}&f_{o, x_o}\end{array}\right)$

The bottom-left block of $J$ is the homeostasis matrix $H$ of of (2.24). Since
$J$ is lower triangular, it follows that linear stability occurs when the
linearized self-couplings are all negative.

The mathematics described here shows that infinitesimal homeostasis occurs in
the system in the second column of (2.22) if and only if

$\det(H)=f_{\rho,x_{\iota}}f_{o,x_{\rho}}-f_{\rho,x_{\rho}}f_{o,x_{\iota}}=0$
(2.25)

at the stable equilibrium $X_{0}$. It is easy to see that (2.23) is a
particular case of (2.25).

Figure 3(b) incorporates some additional information. The arrow from
${\mathcal{I}}$ to node $\iota$ indicates that ${\mathcal{I}}$ occurs in the
equation for $\dot{x}_{\iota}$ as the input parameter. Similarly the arrow
from node $o$ to the symbol $\mathcal{O}$ indicates that node $o$ is the
output node. Finally, the $\pm$ signs indicate which arrows are excitatory or
inhibitory. This extra information is special to biochemical networks and does
not appear as such in the general theory.

In Figure 7(a) we show the graph of the input-output function $x_{o}$, as a
function of the input parameter ${\mathcal{I}}\in[60,140]$, for the model
equations (2.22). Here, $g_{1}(x)=g_{2}(x)=g_{3}(x)=g_{4}(x)=x$, $g_{5}(x)=4x$
and

$f(x)=1+\frac{1}{1+\exp\left(\frac{50-x}{a}\right)}$

where $a=8.33$ and the infinitesimal homeostasis point occurs at
${\mathcal{I}}_{0}=100$. In this case, it can be shown that the infinitesimal
homeostasis occurring at ${\mathcal{I}}_{0}=100$ and $a\approx 8.33$ is in
fact a chair singularity and $\sigma$ is the unfolding parameter, see [115]
for details.  $\Diamond$

######  Example 2.26 (Product inhibition).

Also called _feedback inhibition_ , this is probably one of the simplest and
best known homeostatic mechanisms in biochemistry. In its simplest form,
feedback inhibition means that the product of a biochemical chain inhibits one
or more of the enzymes involved in its own synthesis. Thus if the
concentration of the end product goes up, synthesis is slowed, and if the
concentration goes down, the inhibition is partially withdrawn and the
synthesis goes faster. More specifically, suppose that substrate $\bf X$
influences $\bf Y$, which influences $\bf Z$, and $\bf Z$ inhibits the flux
$g_{1}$ from $\bf X$ to $\bf Y$. The biochemical network for this process is
shown in Figure 4(a).

(a) PI motif from Reed et al. [115]

(b) Input-output network derived from (a)

Figure 4: Product inhibition (PI).

This time the model equations for Figure 4(a) are

$\begin{split}\dot{x}&={\mathcal{I}}-g_{4}(x)-f(z)g_{1}(x)\\\
\dot{y}&=f(z)g_{1}(x)-g_{2}(y)-g_{5}(y)\\\
\dot{z}&=g_{2}(y)-g_{3}(z)\end{split}$ (2.26)

and the input-output equations associated to (2.26) can be read directly from
Figure 4(b)

$\begin{split}\dot{x}_{\iota}&=f_{\iota}(x_{\iota},x_{o},{\mathcal{I}})\\\
\dot{x}_{\rho}&=f_{\rho}(x_{\iota},x_{\rho},x_{o})\\\
\dot{x}_{o}&=f_{o}(x_{\rho},x_{o})\end{split}$ (2.27)

Reed et al. [115] explain that, due to biochemical reasons, the model equation
2.26 for product inhibition must satisfy some additional constraints

$f>0\quad g_{1}^{\prime}<0\quad g_{2}^{\prime}<0$ (2.28)

The Jacobian of (2.27) is

$J=\left(\begin{array}[]{cc;{2pt/2pt}c}f_{\iota,x_{\iota}}&0&f_{\iota, x_o}\\\
\hdashline[2pt/2pt]f_{\rho,x_{\iota}}&f_{\rho,x_{\rho}}&f_{\rho, x_o}\\\
0&f_{o,x_{\rho}}&f_{o, x_o}\end{array}\right)$

The bottom-left block of $J$ is the homeostasis matrix $H$ of of (2.27).

Our general mathematical results show that the model system (2.26) exhibits
infinitesimal homeostasis at a stable equilibrium $X_{0}$ if and only if

$\det(H)=f_{\rho,x_{\iota}}f_{o,x_{\rho}}=0$ (2.29)

That is, either

$f_{\rho,x_{\iota}}=f(z_{0})g_{1}^{\prime}(x_{0})=0\quad\text{or}\quad
f_{o,x_{\rho}}=g_{2}^{\prime}(y_{0})=0$ (2.30)

It follows from (2.28) and (2.30) that the model equation (2.26) cannot
exhibit infinitesimal homeostasis. Nevertheless, Reed et al. [115] show that
this biochemical network equations do exhibit near-prefect homeostasis; that
is, the output $z$ is almost constant for a broad range of input values
${\mathcal{I}}$. In the general admissible equations 2.27 infinitesimal
homeostasis can occur generically. However, due to the special form of the
model equation (2.26) and the additional constraints (2.28), infinitesimal
homeostasis is forced to occur at the ‘boundary of the universal unfolding’,
where $g_{1}^{\prime}=0$ or $g_{2}^{\prime}=0$.

In Figure 7(b) we show the graph of the input-output function $x_{o}$, as a
function of the input parameter ${\mathcal{I}}\in[0,200]$, for the model
equations (2.26). Here, $g_{1}(x)=g_{2}(x)=g_{3}(x)=g_{4}(x)=g_{5}(x)=x$ and

$f(x)=1+\frac{b}{1+\exp\left(\frac{x-50}{a}\right)}$

where $a=0.5$ and $b=200$. As explained above infinitesimal homeostasis cannot
occur in this model equation. The almost flat region occurs from
${\mathcal{I}}=50$ to ${\mathcal{I}}=110$, see [115] for details.  $\Diamond$

######  Example 2.27 (Substrate Inhibition).

The biochemical network model for substrate inhibition is given in Figure
5(a), and the associated model system is given in the first column of (2.31).
This biochemical network and the model system are discussed in Reed et al.
[115].

(a) SI motif from [115]

(b) Input-output network derived from (a)

Figure 5: Substrate inhibition (SI).

The equations associated to the diagram for substrate inhibition from Figure
5(a) are

$\begin{split}\dot{x}&={\mathcal{I}}-g_{1}(x)-g_{4}(x)\\\
\dot{y}&=g_{1}(x)-g_{5}(y)-g_{2}(y)\\\ \dot{z}&=g_{2}(y)-g_{3}(z)\end{split}$
(2.31)

Reed et al. [115] provides a biochemical justification for taking
$g_{1}^{\prime}>0$ for $x>0$, whereas the coupling (or kinetics term)
$g_{2}^{\prime}$ can change sign.

That model system of ODEs can be easily translated to the input-output system

$\begin{split}\dot{x}_{\iota}&=f_{\iota}(x_{\iota},{\mathcal{I}})\\\
\dot{x}_{\rho}&=f_{\rho}(x_{\iota},x_{\rho})\\\
\dot{x}_{o}&=f_{o}(x_{\rho},x_{o})\end{split}$ (2.32)

The Jacobian of (2.32) is

$J=\left(\begin{array}[]{cc;{2pt/2pt}c}f_{\iota,x_{\iota}}&0&0\\\
\hdashline[2pt/2pt]f_{\rho,x_{\iota}}&f_{\rho,x_{\rho}}&0\\\
0&f_{o,x_{\rho}}&f_{o, x_o}\end{array}\right)$

The bottom-left block of $J$ is the homeostasis matrix $H$ of of (2.32). Since
$J$ is lower triangular, it follows that linear stability occurs when the
linearized self-couplings are all negative.

The mathematics described here shows that the condition for infinitesimal
homeostasis is

$\det(H)=f_{\rho,x_{\iota}}f_{o,x_{\rho}}=0$ (2.33)

Note that this is the same condition as in the case of product inhibition
(2.29). Therefore, the network for product inhibition and the network for
substrate inhibition are _core equivalent_. In fact, the network for substrate
inhibition form Figure 5 can be obtained from the network from for product
inhibition form Figure 4 by removing two backward arrows: $o\to\iota$ and
$o\to\rho$. These two backward arrows correspond to two off-diagonal entries
in the last column of the Jacobian matrix that are absent in the substrate
inhibition. As explained before, they do not affect the homeostasis matrix but
can influence the linear stability of the system. In fact, as noted above, the
absence of these two arrows in the substrate inhibition network renders the
Jacobian Matrix lower triangular and so the linear stability depends only on
the linearized self-couplings of the three nodes.

On the assumption that $g_{1}^{\prime}=f_{\rho,x_{\iota}}>0$, infinitesimal
homeostasis is possible only if the coupling $\rho\to o$ is neutral, that is,
if $f_{o,x_{\rho}}=g_{2}^{\prime}=0$ at the equilibrium point. This conclusion
agrees with the observation in [115] that $\bf Z$ can exhibit infinitesimal
homeostasis in the substrate inhibition motif if the infinitesimal homeostasis
is built into the kinetics for the function $g_{2}$ which controls the
coupling between $\bf Y$ and $\bf Z$.

Reed et al. [115] note that neutral coupling can arise from substrate
inhibition of enzymes, enzymes that are inhibited by their own substrates. See
the discussion in [116]. This inhibition leads to reaction velocity curves
that rise to a maximum (the coupling is excitatory) and then descend (the
coupling is inhibitory) as the substrate concentration increases.
Infinitesimal homeostasis with neutral couplings arising from substrate
inhibition often has important biological functions and has been estimated to
occur in about $20\%$ of enzymes [116]. In the 1930s, Haldane [63] introduced
the concept of substrate inhibition in which the substrate of the reaction
itself inhibits the enzyme that catalyzes the reaction. Golubitsky and Wang
[58] have called homeostasis similar to the one occurring in substrate
inhibition _Haldane homeostasis_ , since it arise from neutral coupling, that
is, the linearized coupling between to nodes changes sign as the input
parameter is varied.

In Figure 7(c) we show the graph of the input-output function $x_{o}$, as a
function of the input parameter ${\mathcal{I}}\in[0,350]$, for the model
equations (2.31). Here, $g_{1}(x)=g_{2}(x)=g_{3}(x)=g_{4}(x)=g_{5}(x)=x$ and

$f(x)=1+\frac{1}{1+\exp\left(\frac{50-x}{a}\right)}$

where $a=8$. Here, unlike in the case of product inhibition, infinitesimal
homeostasis can occur by a Haldane type mechanism. After an initial growth
from $1$ to $2$ over the ${\mathcal{I}}$-range $[0,100]$ the curve becomes
flat and the system exhibits perfect homeostasis for ${\mathcal{I}}>100$.
$\Diamond$

######  Example 2.28 (Negative Feedback Loop).

Here, each enzyme $\mathbf{X}$, $\mathbf{Y}$, $\mathbf{Z}$ in the feedback
loop motif, see Figure 6(a) can have active and inactive forms. In the kinetic
equations (2.34) the coupling from $\bf X$ to $\bf Z$ is non-neutral according
to Ma et al. [89]. Hence, in this model only null-degradation homeostasis is
possible. In the kinetic equations in the model the $\dot{y}$ equation does
not depend on $y$ and homeostasis can only be perfect homeostasis. However,
this model is a simplification based on saturation in $y$ [89]. In the
original system $\dot{y}$ does depend on $y$ and we expect standard null-
degradation homeostasis to be possible in that system. The corresponding
input-output network to the negative feedback loop motif in Figure 6(a) is
shown in Figure 6(b).

(a) NF motif adapted from Ma et al. [89]

(b) Input-output network derived from (a)

Figure 6: Negative feedback loop (NF) Unlike the arrows in Figures 3, 5 and 4
that represent mass transfer between substrates, in this negative feedback
loop motif positive or negative arrows between enzymes indicate the activation
or inactivation of an enzyme by a different enzyme.

The equations associated to the negative feedback loop diagram of Figure 6(a)
are

$\begin{split}\dot{x}&={\mathcal{I}}k_{{\mathcal{I}}x}\frac{1-x}{(1-x)+K_{{\mathcal{I}}x}}-F_{x}k^{\prime}_{F_{x}}\frac{x}{x+K^{\prime}_{F_{x}}}\\\
\dot{y}&=zk_{zy}-F_{y}k^{\prime}_{F_{y}}\\\
\dot{z}&=xk_{xz}\frac{1-z}{(1-z)+K_{xz}}-yk^{\prime}_{yz}\frac{z}{z+K^{\prime}_{yz}}\end{split}$
(2.34)

where $k_{{\mathcal{I}}x}$, $K_{{\mathcal{I}}x}$, $F_{x}$,
$k^{\prime}_{F_{x}}$, $K^{\prime}_{F_{x}}$, $k_{zy}$, $F_{y}$,
$k^{\prime}_{F_{y}}$, $k_{xz}$, $K_{xz}$, $k^{\prime}_{yz}$, $K^{\prime}_{yz}$
are 12 constants. That model system of ODEs is easily translated to the input-
output system

$\begin{split}\dot{x}_{\iota}&=f_{\iota}(x_{\iota},{\mathcal{I}})\\\
\dot{x}_{\rho}&=f_{\rho}(x_{\rho},x_{o})\\\
\dot{x}_{o}&=f_{o}(x_{\iota},x_{\rho},x_{o})\end{split}$ (2.35)

The Jacobian of (2.35) is

$J=\left(\begin{array}[]{cc;{2pt/2pt}c}f_{\iota,x_{\iota}}&0&0\\\
\hdashline[2pt/2pt]0&f_{\rho,x_{\rho}}&f_{\rho, x_o}\\\
f_{o,x_{\iota}}&f_{o,x_{\rho}}&f_{o, x_o}\end{array}\right)$ (2.36)

The bottom-left block of $J$ is the homeostasis matrix $H$ of of (2.35).

The mathematics described here shows that the model system (2.35) exhibits
infinitesimal homeostasis if and only

$\det(H)=f_{o,x_{\iota}}f_{\rho,x_{\rho}}=0$ (2.37)

That is, either

$f_{\rho,x_{\iota}}=0\quad\text{or}\quad f_{\rho,x_{\rho}}=0$ (2.38)

The first case is the Haldane homeostasis type discussed before. The second
case is called _null-degradation_ , since it arises when the degradation
constant (i.e., the linearized self-coupling) of the regulatory node changes
sign when the input parameter is varied.

Stability of the equilibrium in this motif implies negative feedback between
$\rho$ and $o$. From the Jacobian matrix 2.36 of (2.35) it follows that, at
null-degradation homeostasis ($f_{\rho,x_{\rho}}=0$), linear stability implies
that

$f_{\iota,x_{\iota}}<0\;,\qquad f_{o,x_{o}}<0\;,\qquad
f_{\rho,x_{o}}f_{o,x_{\rho}}<0$ (2.39)

Conditions (2.39) imply that both the input node and the output node need to
degrade and the couplings $\rho\to o$ and $o\to\rho$ must have opposite signs.
This observation agrees with [89] that homeostasis is possible in the network
motif Figure 6(a) if there is a negative loop between $\bf Y$ and $\bf Z$ and
when the linearized internal dynamics of $\bf Y$ is zero. Therefore, the
negative feedback is ‘forced’ by the condition for null-degradation and
stability.

In Ferrell [46] the author reviews some motifs that are capable of displaying
perfect or near-perfect homeostasis and presents a simplification of
biochemical example of Ma et al. [89], that exhibits null-degradation
homeostasis. In Figure 7(d) we show the graph of the input-output function
$x_{o}$, as a function of the input parameter ${\mathcal{I}}\in[0,100]$, for
the model equations of [46, Fig 2]. After an almost instantaneous jump from
$0$ to $1$ the curve becomes flat and the system exhibits, in fact perfect
homeostasis for ${\mathcal{I}}>2$.  $\Diamond$

(a) FE motif from Reed et al. [115]

(b) PI motif from Reed et al. [115]

(c) SI motif from Reed et al. [115]

(d) NF motif from Ferrell [46]

Figure 7:  Input-output functions of the four $3$-node biochemical networks.
Input-output functions were computed by numerical continuation of an
equilibrium using Auto from XPPAut [44].

As we have seen in these four examples of $3$-node biochemical networks there
are plenty of possibilities of homeostatic behaviors already in very small
networks. Golubitsky and Wang [58] started a systematic investigation of the
$3$-node input-output networks in an attempt to organize the understanding of
the biochemical examples. Their success was the motivation for Wang et al.
[134] to undertake the general case of an $N+2$-node input-output network and
eventually led to a general theory. In hindsight, their theorem for
classification of the $3$-node input-output networks can be simply stated as

######  Theorem 2.29 ([58, 134]).

There are $3$ core equivalence classes of $3$-node input-output core networks
with input node $\neq$ output node. Representatives are given by: (i) the
Feedforward Excitation motif, (ii) the Substrate Inhibition motif and (iii)
the Negative Feedback motif. Any other $3$-node input-output network with two
distinguished nodes can be obtained from these $3$ representatives by adding
or removing backward arrows.

As we shall see later, if we allow input node $=$ output node then there are
two more core equivalence classes, thus giving a complete classification of
all $3$-node input-output networks (see subsection 4.2).

As one goes to four-node input-output networks the situation is much more
complicated. Huang and Golubitsky [70] employs the theory developed in Wang et
al. [134] to show that there are $20$ core equivalence classes of four-node
input-output core networks with input node $\neq$ output node.

###  2.5 Homeostasis Patterns in Input-Output Networks

Now we consider the notion of a homeostasis pattern on a given input-output
network $\mathcal{G}$. A homeostasis pattern is the set of nodes $j$ in
$\mathcal{G}$ (including the output node $o$) such that the node coordinate
$x_{j}$, as a function of ${\mathcal{I}}$, satisfies
$x^{\prime}_{j}({\mathcal{I}}_{0})=0$. In other words, a homeostasis pattern
is a set of nodes $S$ of $\mathcal{G}$, that includes the output node $o$, and
all nodes in $S$ are simultaneously (infinitesimally) homeostatic at a given
parameter value ${\mathcal{I}}_{0}$.

However, before going into the details of this theory, we use such
calculations to give some indication of why the input-output network in Figure
8 has exactly the $4$ homeostasis patterns exhibited in Figure 9.

Consider for example the $6$ node input-output network $\mathcal{G}$ shown in
Figure 8.

Figure 8: A $6$-node input-output network.

Although there are exactly $31$ subsets of nodes of $\mathcal{G}$ including
the output node $o$, only $4$ subsets define homeostasis patterns: $\\{o\\}$,
$\\{o,\tau_{3}\\}$, $\\{o,\tau_{2},\tau_{3}\\}$ and
$\\{o,\tau_{2},\tau_{3},\sigma,\iota\\}$. These homeostasis patterns can be
graphically represented by coloring the nodes of $\mathcal{G}$ that are
homeostatic (see Figure 9).

(a) $f_{\tau_{3},\tau_{3}}$

(b) $f_{\tau_{1},\tau_{1}}$

(c) $f_{\tau_{2},\tau_{2}}$

(d) $(f_{\sigma,\iota}f_{o,\sigma}-f_{\sigma,\sigma}f_{o,\iota})$

Figure 9: The four infinitesimal homeostasis patterns of (2.40). Cyan nodes
are homeostatic.

The admissible system of parametrized equations for the network in Figure 8,
in coordinates $X=(\iota,\sigma,\tau_{1},\tau_{2},\tau_{3},o)$, is:

$\begin{array}[]{ccl}\dot{\iota}&=&f_{\iota}(\iota,\tau_{1},{\mathcal{I}})\\\
\dot{\sigma}&=&f_{\sigma}(\iota,\sigma,\tau_{2})\\\
\dot{\tau}_{1}&=&f_{\tau_{1}}(\sigma,\tau_{1})\\\
\dot{\tau}_{2}&=&f_{\tau_{2}}(\tau_{2},\tau_{3},o)\\\
\dot{\tau}_{3}&=&f_{\tau_{3}}(\tau_{3},o)\\\
\dot{o}&=&f_{o}(\iota,\sigma,o)\end{array}$ (2.40)

The homeostasis matrix $H$ is obtained from the Jacobian matrix $J$ of (2.40)
and is given by

$H=\begin{pmatrix}f_{\sigma,\iota}&f_{\sigma,\sigma}&0&f_{\sigma,\tau_{2}}&0\\\
0&f_{\tau_{1},\sigma}&f_{\tau_{1},\tau_{1}}&0&0&\\\
0&0&0&f_{\tau_{2},\tau_{2}}&f_{\tau_{2},\tau_{3}}\\\
0&0&0&0&f_{\tau_{3},\tau_{3}}\\\ f_{o,\iota}&f_{o,\sigma}&0&0&0\end{pmatrix}$

Using row and column expansion it is straightforward to calculate

$\det(H)=f_{\tau_{3},\tau_{3}}\,f_{\tau_{1},\tau_{1}}\,f_{\tau_{2},\tau_{2}}\,(f_{\sigma,\iota}f_{o,\sigma}-f_{\sigma,\sigma}f_{o,\iota})$
(2.41)

Theorem 2.18says that the input-output function $x_{o}({\mathcal{I}})$
undergoes infinitesimal homeostasis at ${\mathcal{I}}_{0}$ if and only if
$\det(H)=0$, evaluated at $(X({\mathcal{I}}_{0}),{\mathcal{I}}_{0})$. Here
$X({\mathcal{I}}_{0})$ is the equilibrium used to construct the input-output
function. The expression $\det(H)$ is a multivariate polynomial in the partial
derivatives $f_{j,\ell}$ of the components of the admissible vector field. As
a polynomial, $\det(H)$ is reducible with $4$ irreducible factors, so
$\det(H)=0$ if and only if one of its irreducible factors vanishes.

According to [134] these irreducible factors determine the homeostasis types
(see section 3). Hence, (2.40) has $4$ homeostasis types. One of the main
results of this paper says that each homeostasis type determines a unique
homeostasis pattern. Moreover, our theory gives a purely combinatorial
procedure to find the set of nodes that belong to each homeostasis pattern.
When applied to Figure 8 it yields the four patterns in Figure 9.

In a simple example such as Figure 8 it is also possible to find the sets of
nodes in each homeostasis pattern by a ‘bare hands’ calculation based on the
admissible ODEs. Here this calculation serves as a check on the results.
Direct calculations with the ODEs can, of course, be used instead of the
combinatorial approach in sufficiently simple cases.

Finally, we determine the four homeostasis patterns by direct calculation. We
do this by assuming that there is a one-parameter family of stable equilibria
$X({\mathcal{I}})$ where $X({\mathcal{I}}_{0})=X_{0}$ using implicit
differentiation with respect to ${\mathcal{I}}$ (indicated by ′), and
expanding (2.40) to first order at ${\mathcal{I}}_{0}$. The linearized system
of equations is

$\begin{array}[]{rcl}0&=&f_{\iota,\iota}\iota^{\prime}+f_{\iota,\tau_{1}}\tau_{1}^{\prime}+f_{\iota,{\mathcal{I}}}\\\
0&=&f_{\sigma,\iota}\iota^{\prime}+f_{\sigma,\sigma}\sigma^{\prime}+f_{\sigma,\tau_{2}}\tau_{2}^{\prime}\\\
0&=&f_{\tau_{1},\sigma}\sigma^{\prime}+f_{\tau_{1},\tau_{1}}\tau_{1}^{\prime}\\\
0&=&f_{\tau_{2},\tau_{2}}\tau_{2}^{\prime}+f_{\tau_{2},\tau_{3}}\tau_{3}^{\prime}\\\
0&=&f_{\tau_{3},\tau_{3}}\tau_{3}^{\prime}\\\
0&=&f_{o,\iota}\iota^{\prime}+f_{o,\sigma}\sigma^{\prime}\end{array}$ (2.42)

Next we compute the homeostasis patterns corresponding to the $4$ homeostasis
types of (2.40).

(a) Homesotasis Type: $f_{\tau_{3},\tau_{3}}$; Homeostatic nodes: $\\{o\\}$
Equation (2.42) becomes

$\begin{array}[]{rcl}0&=&f_{\iota,\iota}\iota^{\prime}+f_{\iota,\tau_{1}}\tau_{1}^{\prime}+f_{\iota,{\mathcal{I}}}\\\
0&=&f_{\sigma,\iota}\iota^{\prime}+f_{\sigma,\sigma}\sigma^{\prime}+f_{\sigma,\tau_{2}}\tau_{2}^{\prime}\\\
0&=&f_{\tau_{1},\sigma}\sigma^{\prime}+f_{\tau_{1},\tau_{1}}\tau_{1}^{\prime}\\\
0&=&f_{\tau_{2},\tau_{2}}\tau_{2}^{\prime}+f_{\tau_{2},\tau_{3}}\tau_{3}^{\prime}\\\
0&=&0\\\ 0&=&f_{o,\iota}\iota^{\prime}+f_{o,\sigma}\sigma^{\prime}\end{array}$
(2.43)

Since $f_{\tau_{2},\tau_{2}}$ and $f_{\tau_{2},\tau_{3}}$ are generically
nonzero at homeostasis, the fourth equation implies that generically
$\tau_{2}^{\prime}$ and $\tau_{3}^{\prime}$ are nonzero. The second and sixth
equations can be rewritten as

$\begin{pmatrix}f_{\sigma,\iota}&f_{\sigma,\sigma}\\\
f_{o,\iota}&f_{o,\sigma}\end{pmatrix}\begin{pmatrix}\iota^{\prime}\\\
\sigma^{\prime}\end{pmatrix}=-\begin{pmatrix}f_{\sigma,\tau_{2}}\tau_{2}^{\prime}\\\
0\end{pmatrix}$

Generically the right hand side of this matrix equation at ${\mathcal{I}}_{0}$
is nonzero; hence generically $\iota^{\prime}$ and $\sigma^{\prime}$ are also
nonzero. The third equation implies that generically $\tau_{1}^{\prime}$ is
nonzero. Therefore, in this case, the only homeostatic node is $o$.

(b) Homeostasis Type: $f_{\tau_{1},\tau_{1}}$; Homeostatic nodes:
$\\{\iota,\tau_{2},\tau_{3},\sigma,o\\}$
In this case (2.42) becomes

$\begin{array}[]{rcl}0&=&f_{\iota,\iota}\iota^{\prime}+f_{\iota,\tau_{1}}\tau_{1}^{\prime}+f_{\iota,{\mathcal{I}}}\\\
0&=&f_{\sigma,\iota}\iota^{\prime}+f_{\sigma,\sigma}\sigma^{\prime}+f_{\sigma,\tau_{2}}\tau_{2}^{\prime}\\\
0&=&f_{\tau_{1},\sigma}\sigma^{\prime}\\\
0&=&f_{\tau_{2},\tau_{2}}\tau_{2}^{\prime}+f_{\tau_{2},\tau_{3}}\tau_{3}^{\prime}\\\
0&=&f_{\tau_{3},\tau_{3}}\tau_{3}^{\prime}\\\
0&=&f_{o,\iota}\iota^{\prime}+f_{o,\sigma}\sigma^{\prime}\end{array}$ (2.44)

The fifth equation implies that generically $\tau_{3}^{\prime}=0$. The fourth
equation implies that generically $\tau_{2}^{\prime}=0$. The third equation
implies that generically $\sigma^{\prime}=0$ and the sixth equation implies
that generically $\iota^{\prime}$ is zero. It follows that the infinitesimal
homeostasis pattern is
$\iota^{\prime}=\tau_{2}^{\prime}=\tau_{3}^{\prime}=\sigma^{\prime}=o^{\prime}=0$.

(c) Homeostasis Type: $f_{\tau_{2},\tau_{2}}$; Homeostatic nodes:
$\\{\tau_{3},o\\}$
Equation (2.42) becomes

$\begin{array}[]{rcl}0&=&f_{\iota,\iota}\iota^{\prime}+f_{\iota,\tau_{1}}\tau_{1}^{\prime}+f_{\iota,{\mathcal{I}}}\\\
0&=&f_{\sigma,\iota}\iota^{\prime}+f_{\sigma,\sigma}\sigma^{\prime}+f_{\sigma,\tau_{2}}\tau_{2}^{\prime}\\\
0&=&f_{\tau_{1},\sigma}\sigma^{\prime}+f_{\tau_{1},\tau_{1}}\tau_{1}^{\prime}\\\
0&=&f_{\tau_{2},\tau_{3}}\tau_{3}^{\prime}\\\
0&=&f_{\tau_{3},\tau_{3}}\tau_{3}^{\prime}\\\
0&=&f_{o,\iota}\iota^{\prime}+f_{o,\sigma}\sigma^{\prime}\end{array}$ (2.45)

The fourth or fifth equation implies that $\tau_{3}^{\prime}=0$. The first and
sixth equations imply that $\iota^{\prime}$, $\sigma^{\prime}$, and
$\tau_{2}^{\prime}$ are nonzero. The third equation implies that generically
$\tau_{1}^{\prime}$ is nonzero. Hence the infinitesimal homeostasis pattern is
$\\{\tau_{3},o\\}$.

(d) Homeostasis Type:
$f_{\sigma,\iota}f_{o,\sigma}-f_{\sigma,\sigma}f_{o,\iota}$; Homeostatic nodes
$\\{\tau_{2},\tau_{3},o\\}$
To repeat, Equation (2.42) is

$\begin{array}[]{rcl}0&=&f_{\iota,\iota}\iota^{\prime}+f_{\iota,\tau_{1}}\tau_{1}^{\prime}+f_{\iota,{\mathcal{I}}}\\\
0&=&f_{\sigma,\iota}\iota^{\prime}+f_{\sigma,\sigma}\sigma^{\prime}+f_{\sigma,\tau_{2}}\tau_{2}^{\prime}\\\
0&=&f_{\tau_{1},\sigma}\sigma^{\prime}+f_{\tau_{1},\tau_{1}}\tau_{1}^{\prime}\\\
0&=&f_{\tau_{2},\tau_{2}}\tau_{2}^{\prime}+f_{\tau_{2},\tau_{3}}\tau_{3}^{\prime}\\\
0&=&f_{\tau_{3},\tau_{3}}\tau_{3}^{\prime}\\\
0&=&f_{o,\iota}\iota^{\prime}+f_{o,\sigma}\sigma^{\prime}\end{array}$ (2.46)

The fifth equation implies generically that $\tau_{3}^{\prime}=0$ and the
fourth equation implies generically that $\tau_{2}^{\prime}=0$. Again, the
second and sixth equations can be rewritten in matrix form as

$\begin{pmatrix}f_{\sigma,\iota}&f_{\sigma,\sigma}\\\
f_{o,\iota}&f_{o,\sigma}\end{pmatrix}\begin{pmatrix}\iota^{\prime}\\\
\sigma^{\prime}\end{pmatrix}=-\begin{pmatrix}f_{\sigma,\tau_{2}}\tau_{2}^{\prime}\\\
0\end{pmatrix}=0$

Hence generically $\iota^{\prime}$ and $\sigma^{\prime}$ are nonzero. The
third equation implies that $\tau_{1}$ is nonzero. Hence the homeostatic nodes
are $\tau_{2},\tau_{3},o$.

The homeostasis types of (2.40) with the corresponding homeostasis patterns
are summarized in Table 2.

Table 2: Infinitesimal homeostasis patterns for admissible systems in (2.40). Homeostasis Type | Homeostasis Pattern | Figure 9
---|---|---
$f_{\tau_{3},\tau_{3}}$ | $\\{o\\}$ | (a)
$f_{\tau_{1},\tau_{1}}$ | $\\{\iota,\tau_{2},\tau_{3},\sigma,o\\}$ | (b)
$f_{\tau_{2},\tau_{2}}$ | $\\{\tau_{3},o\\}$ | (c)
$f_{\sigma,\iota}f_{o,\sigma}-f_{\sigma,\sigma}f_{o,\iota}$ | $\\{\tau_{2},\tau_{3},o\\}$ | (d)

In principle the homeostasis patterns can be computed in the manner shown
above; but in practice this becomes very complicated.

##  3 Classification

In this section we describe the theory of Wang et al. [134] and Duncan et al.
[40] about the classification of homeostasis types and homeostasis patterns in
an input-output network, respectively.

###  3.1 Homeostasis Subnetworks

In [134] the authors show that the determination of infinitesimal homeostasis
in an input-output networks reduces to the study of core networks. We assume
throughout that the input-output networks are core networks.

The first main result of [134] is the observation that one can apply the
Frobenius-König [118, 23] theorem to the homeostasis matrix $H$ and obtain a
decomposition of $\det(H)$ into its irreducible factors.

######  Theorem 3.1 ([134, Theorem 1.11]).

Assume (2.20) has a hyperbolic equilibrium at $(X^{*},{\mathcal{I}}^{*})$.
Then there are permutation matrices $P$ and $Q$ such that $PHQ$ is block upper
triangular with square diagonal blocks $B_{1},\ldots,B_{m}$. The blocks
$B_{j}$ are irreducible in the sense that each $B_{j}$ cannot be further block
triangularized. It follows that

$\displaystyle\det(H)=\det(B_{1})\cdots\det(B_{m})$ (3.47)

is an irreducible factorization of $\det(H)$.

######  Definition 3.2.

Let $\mathcal{G}$ be an input-output network and $H$ its homeostasis matrix.
Each irreducible square block $B_{\eta}$ in (3.47) is called a homeostasis
block. Further, we say that infinitesimal homeostasis in $\mathcal{G}$ is of
homeostasis type $B_{\eta}$ if for all $\xi\neq\eta$, $\det(B_{\eta})=0$ and
$\det(B_{\xi})\neq 0$.  $\Diamond$

######  Remark 3.3.

Let $B_{\eta}$ be a homeostasis type and let

$h_{\eta}({\mathcal{I}})\equiv\det B_{\eta}(X({\mathcal{I}}),{\mathcal{I}})$

A chair point of type $\eta$ occurs at ${\mathcal{I}}_{0}$ if
$h_{\eta}({\mathcal{I}}_{0})=h_{\eta}^{\prime}({\mathcal{I}}_{0})=0$ and
$h_{\eta}^{\prime\prime}({\mathcal{I}}_{0})\neq 0$.  $\Diamond$

In principle every homeostasis type can lead to infinitesimal homeostasis,
that is, $h_{\eta}({\mathcal{I}}_{0})=0$ for some input value
${\mathcal{I}}_{0}$. For simplicity, we say that node $x_{o}$ is _homeostatic
at ${\mathcal{I}}_{0}$_.

Let $B_{\eta}$ be a $k\times k$ diagonal block of $H$, hence $\det(B_{\eta})$
has degree $k$. Since the entries of $B_{\eta}$ are entries of $H$, these
entries have the form $f_{\rho,x_{\tau}}$; that is, the entries are either $0$
(if $\tau\to\rho$ is not an arrow in $\mathcal{G}$), self-coupling (if
$\tau=\rho$), or coupling (if $\tau\to\rho$ is an arrow in $\mathcal{G}$).

Since $P$ and $Q$ in Theorem 3.1 are constant permutation matrices, all
entries in each row (resp. column) of $B_{\eta}$ must lie in a single row
(resp. column) of $H$. Hence, $B_{\eta}$ has the form

$B_{\eta}=\begin{pmatrix}f_{\rho_{1},x_{\tau_{1}}}&\cdots&f_{\rho_{1},x_{\tau_{k}}}\\\
\vdots&\ddots&\vdots\\\
f_{\rho_{k},x_{\tau_{1}}}&\cdots&f_{\rho_{k},x_{\tau_{k}}}\end{pmatrix}$
(3.48)

It follows that the number of self-coupling entries of $B_{\eta}$ are the same
no matter which permutation matrices $P$ and $Q$ are used to bring $H$ to the
block-diagonal form.

######  Theorem 3.4 ([134, Theorem 4.4]).

Each $k\times k$ submatrix $B_{\eta}$ has either $k$ or $k-1$ self-coupling
entries.

######  Definition 3.5.

The homeostasis class of an irreducible component $B_{\eta}$ of degree $k$ is
_appendage_ if $B_{\eta}$ has $k$ self-couplings and _structural_ if
$B_{\eta}$ has $k-1$ self-couplings.  $\Diamond$

Now we can associate a _homeostasis subnetwork_
$\mathcal{K}_{\eta}\subset\mathcal{G}$ with each homeostasis block $B_{\eta}$
and give a graph-theoretic description of each $\mathcal{K}_{\eta}$.

######  Definition 3.6.

The _homeostasis subnetwork_ $\mathcal{K}_{\eta}$ of $\mathcal{G}$ associated
with the homeostasis block $B_{\eta}$ is defined as follows. The _nodes_ in
$\mathcal{K}_{\eta}$ are the union of nodes $p$ and $q$ where $f_{p,x_{q}}$ is
a nonzero entry in $B_{\eta}$ and the _arrows_ of $\mathcal{K}_{\eta}$ are the
union of arrows $q\to p$ where $p\neq q$.  $\Diamond$

In order to give a combinatorial characterization of the structural and
appendage homeostasis subnetworks we need some simple concepts from graph
theory.

######  Definition 3.7.

Let $\mathcal{G}$ be a core input-output network.

  1. [(a)]

  2. 1.

A _simple path_ from node $\kappa_{1}$ to node $\kappa_{2}$ in $\mathcal{G}$
is a directed path that starts at $\kappa_{1}$, ends at $\kappa_{2}$, and
visits each node on the path exactly once. We denote the existence of a simple
path from $\kappa_{1}$ to $\kappa_{2}$ by
$\kappa_{1}\rightsquigarrow\kappa_{2}$. A _simple cycle_ is a simple path
whose first and last nodes are identical.

  3. 2.

An _$\iota o$ -simple path_ is a simple path from the input node $\iota$ to
the output node $o$.

  4. 3.

A node $\sigma$ is _simple_ if it lies on an $\iota o$-simple path. A node
$\tau$ is _appendage_ if it is not simple.

  5. 4.

A simple node $\rho$ is _super-simple_ if it lies on every $\iota o$-simple
path.  $\Diamond$

We typically use $\sigma$ to denote a simple node, $\rho$ to denote a super-
simple node, and $\tau$ to denote an appendage node when the type of the node
is assumed a priori. Otherwise, we use $\kappa$ to denote an arbitrary node.
Note that $\iota$ and $o$ are super-simple nodes.

Let $\rho_{0},\rho_{1},\ldots,\rho_{q},\rho_{q+1}$ be the super-simple nodes,
where $\rho_{0}=\iota$ and $\rho_{q+1}=o$. The super-simple nodes are totally
ordered by the order of their appearance on any $\iota o$-simple path, and
this ordering is independent of the $\iota o$-simple path. We denote an $\iota
o$-simple path by

$\iota\rightsquigarrow\rho_{1}\rightsquigarrow\cdots\rightsquigarrow\rho_{q}\rightsquigarrow
o$

where $\rho_{j}\rightsquigarrow\rho_{j+1}$ indicates a simple path from
$\rho_{j}$ to $\rho_{j+1}$. The ordering of the super-simple nodes is denoted
by

$\rho_{0}\prec\rho_{1}\prec\cdots\prec\rho_{q}\prec\rho_{q+1},$

and $\prec$ is a total ordering. The ordering $\prec$ extends to a partial
ordering of simple nodes, as follows. If there exists a super-simple node
$\rho$ and an $\iota o$-simple path such that

$\iota\rightsquigarrow\sigma_{1}\rightsquigarrow\rho\rightsquigarrow\sigma_{2}\rightsquigarrow
o$

then the partial orderings

$\sigma_{1}\prec\rho\qquad\rho\prec\sigma_{2}\qquad\sigma_{1}\prec\sigma_{2}$

are valid. In this partial ordering every simple node is comparable to every
super-simple node but two simple nodes that lie between the same adjacent
super-simple nodes need not be comparable.

We recall the definition of transitive (or strong) components of a network.
Two nodes are equivalent if there is a path from one to the other and back. A
transitive component is an equivalence class for this equivalence relation.

######  Definition 3.8.

  1. [(a)]

  2. 1.

Let $S$ be an $\iota o$-simple path. The complementary subnetwork of $S$ is
the network $C_{S}$ whose nodes are nodes that are not in $S$ and whose arrows
are those that connect nodes in $C_{S}$.

  3. 2.

An appendage node $\tau$ is super-appendage if for each $C_{S}$ containing
$\tau$, the transitive component of $\tau$ in $C_{S}$ consists only of
appendage nodes.  $\Diamond$

Note that this definition of super-appendage leads to a slightly different,
but equivalent, definition of homeostasis subnetwork to the one given in
[134]. However, this change enables us to define pattern networks in a more
straightforward way (see Remark 3.11).

Now we start with the definition of structural subnetworks.

######  Definition 3.9.

Let $1\leq j\leq q+1$ and $\rho_{j-1}\prec\rho_{j}$ be two consecutive super-
simple nodes. Then the $j^{th}$ _simple subnetwork_
$\mathcal{L}_{j}^{\prime\prime}$, the $j^{th}$ augmented simple subnetwork
$\mathcal{L}_{j}^{\prime}$, and the $j^{th}$ _structural subnetwork_
$\mathcal{L}_{j}$ are defined in four steps as follows.

  1. [(a)]

  2. 1.

The $j^{th}$ _simple subnetwork_ $\mathcal{L}_{j}^{\prime\prime}$ consists of
simple nodes $\sigma$ where

$\rho_{j-1}\prec\sigma\prec\rho_{j}$

and all arrows connecting these nodes. Note that
$\mathcal{L}_{j}^{\prime\prime}$ does not contain the super-simple nodes
$\rho_{j-1}$ and $\rho_{j}$, and $\mathcal{L}_{j}^{\prime\prime}$ can be the
empty set.

  3. 2.

An appendage but not super-appendage node $\tau$ is linked to
$\mathcal{L}_{j}^{\prime\prime}$ if for some complementary subnetwork $C_{S}$
the transitive component of $\tau$ in $C_{S}$ is the union of $\tau$, nodes in
$\mathcal{L}_{j}^{\prime\prime}$, and non-super-appendage nodes. The set of
$j^{th}$-linked appendage nodes $T_{j}$ is the set of non-super-appendage
nodes that are linked to $\mathcal{L}_{j}^{\prime\prime}$.

  4. 3.

The $j^{th}$ augmented simple subnetwork $\mathcal{L}_{j}^{\prime}$ is

$\mathcal{L}_{j}^{\prime}=\mathcal{L}_{j}^{\prime\prime}\cup T_{j}$

and all arrows connecting these nodes.

  5. 4.

The $j^{th}$ _structural subnetwork_ $\mathcal{L}_{j}$ consists of the
augmented simple subnetwork $\tilde{\mathcal{L}}_{j}$ and adjacent super-
simple nodes. That is

$\mathcal{L}_{j}=\\{\rho_{j-1}\\}\cup\mathcal{L}_{j}^{\prime}\cup\\{\rho_{j}\\}$

and all arrows connecting these nodes.  $\Diamond$

Note that an structural subnetwork $\mathcal{L}_{j}$ is itself an input-output
network, with input node $\rho_{j-1}$ and output node $\rho_{j}$.

Next we define the appendage subnetworks, which were defined in Section 1.7.2
of [134] as any transitive component of the subnetwork consisting only of
appendage nodes and the arrows between them.

######  Definition 3.10.

An _appendage subnetwork_ $\mathcal{A}_{k}$, $1\leq k\leq r$ is a transitive
component of the subnetwork of super-appendage nodes.  $\Diamond$

######  Remark 3.11.

In [134] the authors define an appendage subnetwork as a transitive component
of appendage nodes $\mathcal{A}$ that satisfy the no cycle condition. The no-
cycle condition is formulated in terms of the non-existence of a cycle between
appendage nodes in $\mathcal{A}$ and the simple nodes in $C_{S}$ for all
simple $\iota o$-simple paths $S$. Here, we define an appendage subnetwork as
a transitive component of super-appendage nodes, which are defined in terms of
transitive components with respect to $C_{S}$ for all simple $\iota o$-simple
paths $S$. These two definitions are equivalent because two nodes belong to
the same transitive component if and only if both nodes lie on a (simple)
cycle.  $\Diamond$

######  Theorem 3.12 ([134]).

Let $\mathcal{K_{\eta}}$ be a homeostasis subnetwork of $\mathcal{G}$ with
corresponding diagonal block $B_{\eta}$.

  1. [(1)]

  2. 1.

If $B_{\eta}$ is structural then $\mathcal{K_{\eta}}=\mathcal{L}_{j}$ for some
$1\leq k\leq q$. Moreover, $\mathcal{L}_{j}$ is an input-output network, hence
$B_{\eta}$ is an homeostasis matrix and

$\det(B_{\eta})=\det(H(\mathcal{L}_{j}))$

where $H(\mathcal{L}_{j})$ is the homeostasis matrix of $\mathcal{L}_{j}$.

  3. 2.

If $B_{\eta}$ is appendage then $\mathcal{K_{\eta}}=\mathcal{A}_{k}$ for some
$1\leq k\leq r$. Moreover, $\mathcal{A}_{k}$ is a network without any
distinguished node, hence $B_{\eta}$ is a jacobian matrix and

$\det(B_{\eta})=\det(J(\mathcal{A}_{k}))$

where $J(\mathcal{A}_{j})$ is the jacobian matrix of $\mathcal{A}_{k}$.

Theorem 3.12 says that the network $\mathcal{G}$ decomposes into subnetworks
$\mathcal{K}_{\eta}$, each of which is of two types: structural or appendage.
Two of these subnetworks may have common nodes only when they are consecutive
structural subnetworks. Furthermore, the irreducible factors of $\det(H)$ can
be obtained directly from the set of homeostasis subnetworks.

###  3.2 Homeostasis Patterns

######  Definition 3.13.

A homeostasis pattern corresponding to the homeostasis block $B_{\eta}$ at
${\mathcal{I}}_{0}$ is the collection of all nodes, including the output node
$o$, that are simultaneously forced to be homeostatic at $\mathcal{I}_{0}$.
$\Diamond$

######  Definition 3.14.

Assume that the output node $o$ is homeostatic at $(X_{0},{\mathcal{I}}_{0})$,
that is, $x_{o}^{\prime}({\mathcal{I}}_{0})=0$ for the input value
${\mathcal{I}}_{0}$.

  1. [(a)]

  2. 1.

We call the homeostasis subnetwork $\mathcal{K}_{\eta}$ _homeostasis inducing_
if $h_{\mathcal{K}_{\eta}}\equiv\det(B_{\eta})=0$ at
$(X_{0},{\mathcal{I}}_{0})$.

  3. 2.

Homeostasis of a node $\kappa\in\mathcal{G}$ is _induced_ by a homeostasis
subnetwork $\mathcal{K}$, denoted $\mathcal{K}\Rightarrow\kappa$, if $\kappa$
is generically homeostatic whenever $\mathcal{K}$ is homeostasis inducing.

  4. 3.

A homeostasis subnetwork $\mathcal{K}$ _induces_ a subset of nodes
$\mathcal{N}$ ($\mathcal{K}\Rightarrow\mathcal{N}$), if
$\mathcal{K}\Rightarrow\kappa$ for each node
$\kappa\in\mathcal{N}\subset\mathcal{G}$.  $\Diamond$

By definition, every homeostasis subnetwork $\mathcal{K}$ induces homeostasis
in the output node $o$, that is, $\mathcal{K}\Rightarrow o$.

In this subsection we construct the _homeostasis pattern network_
$\mathcal{P}$ associated with $\mathcal{G}$ (see Definition 3.21), which
serves to organize the homeostasis subnetworks and to clarify how each
homeostasis subnetwork connects to the others.

The homeostasis pattern network is defined in the following steps. First, we
define the structural pattern network $\mathcal{P}_{\mathcal{S}}$ in terms of
the structural subnetworks of $\mathcal{G}$. Second, we define the appendage
pattern network $\mathcal{P}_{\mathcal{A}}$ in terms of the appendage
subnetworks of $\mathcal{G}$. Finally, we define how the nodes in
$\mathcal{P}_{\mathcal{S}}$ connect to nodes in $\mathcal{P}_{\mathcal{A}}$
and conversely.

######  Definition 3.15.

The structural pattern network $\mathcal{P}_{\mathcal{S}}$ is the feedforward
network whose nodes are the super-simple nodes $\rho_{j}$ and the _backbone_
nodes $\widetilde{\mathcal{L}}_{j}$, where $\widetilde{\mathcal{L}}_{j}$ is
the augmented structural subnetwork $\mathcal{L}_{j}^{\prime}$ treated as a
single node. The nodes and arrows of $\mathcal{P}_{\mathcal{S}}$ are given as
follows.

$\iota=\rho_{0}\to\widetilde{\mathcal{L}}_{1}\to\rho_{1}\to\cdots\to\widetilde{\mathcal{L}}_{q+1}\to\rho_{q+1}=o$
(3.49)

If a structural subnetwork ${\mathcal{L}}$ consists of an arrow between two
adjacent super-simple nodes (Haldane type) then the corresponding augmented
structural subnetwork ${\mathcal{L}^{\prime}}$ is the empty network;
nevertheless the corresponding backbone node $\widetilde{\mathcal{L}}$ must be
included in the structural pattern network $\mathcal{P}_{\mathcal{S}}$.
$\Diamond$

######  Definition 3.16.

The appendage pattern network $\mathcal{P}_{\mathcal{A}}$ is the network whose
nodes are the components $\widetilde{\mathcal{A}}$ in the condensation of the
subnetwork of super-appendage nodes. Such a node $\widetilde{\mathcal{A}}$ is
called an appendage component. An arrow connects nodes
$\widetilde{\mathcal{A}}_{1}$ and $\widetilde{\mathcal{A}}_{2}$ if and only if
there are super-appendage nodes $\tau_{1}\in\widetilde{\mathcal{A}}_{1}$ and
$\tau_{2}\in\widetilde{\mathcal{A}}_{2}$ such that $\tau_{1}\to\tau_{2}$ in
$\mathcal{G}$.  $\Diamond$

######  Definition 3.17.

Let $\mathcal{G}$ be a directed graph. Consider the equivalence relation
$\sim$ given by path equivalence on the set of nodes of $\mathcal{G}$. The
equivalence classes are the path components of $\mathcal{G}$. Let
$\mathcal{G}^{c}=\mathcal{G}/\\!\\!\sim$ be the quotient graph called the
_condensation of $\mathcal{G}$_. That is, the nodes of $\bar{\mathcal{G}}$ are
the path components of $\mathcal{G}$ and the arrows of $\bar{\mathcal{G}}$ are
defined as follows. There is a directed arrow form path component
$\mathcal{A}_{1}$ to path component $\mathcal{A}_{2}$ if there are nodes
$\tau_{1}\in\mathcal{A}_{1}$, $\tau_{2}\in\mathcal{A}_{2}$ and a direct arrow
form $\tau_{1}$ to $\tau_{2}$.  $\Diamond$

To complete the homeostasis pattern network, we describe how the nodes in
$\mathcal{P}_{\mathcal{A}}$ and the nodes in $\mathcal{P}_{\mathcal{S}}$ are
connected. To do so, we take advantage of the feedforward ordering of the
nodes in $\mathcal{P}_{\mathcal{S}}$ and the feedback ordering of the nodes in
$\mathcal{P}_{\mathcal{A}}$.

######  Definition 3.18.

A simple path from $\kappa_{1}$ to $\kappa_{2}$ is an _appendage path_ if some
node on this path is an appendage node and every node on this path, except
perhaps for $\kappa_{1}$ and $\kappa_{2}$, is an appendage node.  $\Diamond$

######  Definition 3.19 (How $\mathcal{P}_{\mathcal{A}}$ connects to
$\mathcal{P}_{\mathcal{S}}$).

Given a node $\widetilde{\mathcal{A}}\in\mathcal{P}_{\mathcal{A}}$, we
construct a unique arrow from $\widetilde{\mathcal{A}}$ to the structural
pattern network $\mathcal{P}_{\mathcal{S}}$ in two steps:

  1. [(a)]

  2. 1.

Consider the collection of nodes $\mathcal{V}$ in $\mathcal{P}_{\mathcal{S}}$
for which there exists a simple node $\sigma\in\mathcal{V}$ and appendage node
$\tau\in\widetilde{\mathcal{A}}$, such that there is an appendage path from
$\tau$ to $\sigma$.

  3. 2.

Let $\mathcal{V}_{max}(\widetilde{\mathcal{A}})$ be a maximal node in this
collection, that is, the most downstream in $\mathcal{P}_{\mathcal{S}}$. It
follows from (3.49) that $\mathcal{V}_{max}$ is either a super-simple node
$\rho_{j}$ or a backbone node $\widetilde{\mathcal{L}}_{j}$. Maximality
implies that $\mathcal{V}_{max}$ is uniquely defined. We then say that there
is an arrow from $\widetilde{\mathcal{A}}$ to
$\mathcal{V}_{max}\in\mathcal{P}_{\mathcal{S}}$.  $\Diamond$

######  Definition 3.20 (How $\mathcal{P}_{\mathcal{A}}$ is connected from
$\mathcal{P}_{\mathcal{S}}$).

Given a node $\widetilde{\mathcal{A}}\in\mathcal{P}_{\mathcal{A}}$ we choose
uniquely an arrow from the structural pattern network
$\mathcal{P}_{\mathcal{S}}$ to $\widetilde{\mathcal{A}}$ in two steps:

  1. [(a)]

  2. 1.

Consider the collection of nodes $\mathcal{V}$ in $\mathcal{P}_{\mathcal{S}}$
for which there exists a simple node $\sigma\in\mathcal{V}$ and appendage node
$\tau\in\widetilde{\mathcal{A}}$, such that there is an appendage path from
$\sigma$ to $\tau$.

  3. 2.

Let $\mathcal{V}_{min}(\widetilde{\mathcal{A}})$ be a minimal node in this
collection, that is, the most upstream node in $\mathcal{P}_{\mathcal{S}}$.
Then $\mathcal{V}_{min}$ is either a super-simple node $\rho_{j}$ or a
backbone node $\widetilde{\mathcal{L}}_{j}$, and the minimality implies
uniqueness of $\mathcal{V}_{min}$. We then say that there is an arrow from
$\mathcal{V}_{min}\in\mathcal{P}_{\mathcal{S}}$ to $\widetilde{\mathcal{A}}$.
$\Diamond$

Since we consider only core input-output networks, all appendage nodes are
downstream from $\iota$ and upstream from $o$. Hence, for any node
$\widetilde{\mathcal{A}}\in\mathcal{P}_{\mathcal{A}}$, there always exist
nodes $\mathcal{V}_{min},\mathcal{V}_{max}\in\mathcal{P}_{\mathcal{S}}$ as
mentioned above.

######  Definition 3.21.

The _homeostasis pattern network_ $\mathcal{P}$ is the network whose nodes are
the union of the nodes of the structural pattern network
$\mathcal{P}_{\mathcal{S}}$ and the appendage pattern network
$\mathcal{P}_{\mathcal{A}}$. The arrows of $\mathcal{P}$ are the arrows of
$\mathcal{P}_{\mathcal{S}}$, the arrows of $\mathcal{P}_{\mathcal{A}}$, and
the arrows between $\mathcal{P}_{\mathcal{S}}$ and $\mathcal{P}_{\mathcal{A}}$
as described above.  $\Diamond$

######  Remark 3.22.

The super-simple nodes in $\mathcal{P}$ correspond to the super-simple nodes
of $\mathcal{G}$. Each super-simple node $\rho_{j}\in\mathcal{G}$ (for $1\leq
j\leq q$) belongs to exactly two structural subnetworks $\mathcal{L}_{j-1}$
and $\mathcal{L}_{j}$. Thus they are not associated to a single homeostasis
subnetwork of $\mathcal{G}$.  $\Diamond$

It follows from Remark 3.22 that there is a correspondence between the
homeostasis subnetworks of $\mathcal{G}$ and the non-super-simple nodes of
$\mathcal{P}$.

######  Remark 3.23.

  1. [(a)]

  2. 1.

Each structural subnetwork $\mathcal{L}\subseteq\mathcal{G}$ corresponds to
the backbone node $\widetilde{\mathcal{L}}\in\mathcal{P}_{\mathcal{S}}$. Note
that the augmented structural subnetworks
$\mathcal{L}^{\prime}\subsetneq\mathcal{L}$ are not homeostasis subnetworks.

  3. 2.

Each appendage subnetworks $\mathcal{A}\subset\mathcal{G}$ corresponds to a
appendage component $\widetilde{\mathcal{A}}\in\mathcal{P}_{\mathcal{A}}$.

  4. 3.

For simplicity in notation we let $\mathcal{V}_{\mathcal{S}}$ denote a node in
$\mathcal{P}_{\mathcal{S}}$. Further we let $\widetilde{\mathcal{V}}$ denote a
non-super-simple node of $\mathcal{P}$ and $\mathcal{V}$ denote its
corresponding homeostasis subnetwork.  $\Diamond$

The main point of introducing the homeostasis pattern network $\mathcal{P}$ is
to relate homeostatic induction between the set of homeostasis subnetworks of
$\mathcal{G}$ to induction between nodes in $\mathcal{P}$. In Definition 3.24
bellow we formalize this notion. Hence, every node in a homeostasis pattern
(which can be backbone or appendage) is induced by either a backbone node or
an appendage node in the homeostasis pattern network $\mathcal{P}$.

######  Definition 3.24.

Let $\widetilde{\mathcal{V}}_{1},\widetilde{\mathcal{V}}_{2}\in\mathcal{P}$ be
non-super-simple nodes and $\rho\in\mathcal{P}$ be a super-simple node. Let
$\mathcal{V}_{1},\mathcal{V}_{2}\subset\mathcal{G}$ be the corresponding
homeostasis subnetworks to
$\widetilde{\mathcal{V}}_{1},\widetilde{\mathcal{V}}_{2}\in\mathcal{P}$. We
say that $\widetilde{\mathcal{V}}_{1}$ induces $\widetilde{\mathcal{V}}_{2}$,
denoted by
$\widetilde{\mathcal{V}}_{1}\Rightarrow\widetilde{\mathcal{V}}_{2}$, if and
only if $\mathcal{V}_{1}\Rightarrow\mathcal{V}_{2}$. We say that
$\widetilde{\mathcal{V}}_{1}$ induces $\rho$, denoted by
$\widetilde{\mathcal{V}}_{1}\Rightarrow\rho$, if and only if
$\mathcal{V}_{1}\Rightarrow\rho$.  $\Diamond$

We exclude super-simple nodes of $\mathcal{P}$ from being ‘homeostasis
inducing’ because they are not associated to a homeostasis subnetwork of
$\mathcal{G}$ (see Remark 3.22). However, when a backbone node
$\widetilde{\mathcal{L}}_{j}\in\mathcal{P}$ induces homeostasis on other nodes
of $\mathcal{P}$, it is the corresponding structural subnetwork
$\mathcal{L}_{j}$, with its two super-simple nodes $\rho_{j-1},\rho_{j}$ that
induce homeostasis.

As explained before, the homeostasis pattern network $\mathcal{P}$ allows us
to characterize homeostasis patterns by reducing to four possibilities that
are covered by Theorems 3.25 \- 3.28.

Structural homeostasis patterns are given by the following two theorems.

######  Theorem 3.25 (Structural $\Rightarrow$ Structural).

A backbone node $\widetilde{\mathcal{L}}_{j}\in\mathcal{P}_{\mathcal{S}}$
induces every node of the structural pattern network $\mathcal{P}_{S}$
strictly downstream from $\widetilde{\mathcal{L}}_{j}$, but no other nodes of
$\mathcal{P}_{S}$.

See Figure 10 for an application of Theorem 3.25.

Figure 10: An example of structural network induced by structural homeostasis.
Suppose the backbone node $\widetilde{\mathcal{L}}_{j}$ in red is homeostasis
inducing. Then Theorem 3.25 implies that the blue nodes in the structural
pattern network are all homeostatic.

######  Theorem 3.26 (Structural $\Rightarrow$ Appendage).

A backbone node $\widetilde{\mathcal{L}}_{j}\in\mathcal{P}_{\mathcal{S}}$
induces every appendage component of $\mathcal{P}_{\mathcal{A}}$ whose
$\mathcal{V}_{min}$ (see Definition 3.20) is strictly downstream, but no other
nodes of $\mathcal{P}_{\mathcal{A}}$.

See Figure 11 for an application of Theorem 3.26.

Figure 11: An example of appendage subnetworks induced by structural
homeostasis. Suppose the backbone node $\widetilde{\mathcal{L}}_{j}$ in red is
homeostasis inducing. Then by Theorem 3.26 and the fact that the super-simple
node $\rho_{j}$ is strictly downstream from $\widetilde{\mathcal{L}}_{j}$ the
blue appendage components downstream from $\rho_{j}$ are homeostatic.

The appendage homeostasis patterns are characterized by the following two
theorems.

######  Theorem 3.27 (Appendage $\Rightarrow$ Structural).

An appendage component $\widetilde{\mathcal{A}}\in\mathcal{P}_{\mathcal{A}}$
induces every super-simple node of $\mathcal{P}_{S}$ downstream from
$\mathcal{V}_{max}(\widetilde{\mathcal{A}})$ (see Definition 3.19), but no
other super-simple nodes. Further, an appendage component
$\widetilde{\mathcal{A}}\in\mathcal{P}_{\mathcal{A}}$ induces a backbone node
$\widetilde{\mathcal{L}}_{j}$ if and only if $\widetilde{\mathcal{L}}_{j}$ is
strictly downstream from $\mathcal{V}_{max}(\widetilde{\mathcal{A}})$.

See Figure 12 for an application of Theorem 3.25.

Figure 12: An example of structural subnetworks induced by appendage
homeostasis. Suppose the appendage component $\widetilde{\mathcal{A}}$ in red
is homeostasis inducing. Since $\widetilde{\mathcal{A}}$ connects to the
super-simple node $\rho_{j-1}$, then by Theorem 3.27 the blue nodes in the
structural pattern network are homeostatic.

######  Theorem 3.28 (Appendage $\Rightarrow$ Appendage).

An appendage component
$\widetilde{\mathcal{A}}_{i}\in\mathcal{P}_{\mathcal{A}}$ induces an appendage
component $\widetilde{\mathcal{A}}_{j}\in\mathcal{P}_{\mathcal{A}}$ if and
only if $\widetilde{\mathcal{A}}_{i}$ is strictly upstream from
$\widetilde{\mathcal{A}}_{j}$ and every path from
$\widetilde{\mathcal{A}}_{i}$ to $\widetilde{\mathcal{A}}_{j}$ in
$\mathcal{P}$ contains a super-simple node $\rho$ satisfying
$\mathcal{V}_{max}(\widetilde{\mathcal{A}}_{i})\preceq\rho\preceq\mathcal{V}_{min}(\widetilde{\mathcal{A}}_{j})$.

See Figure 13 for an application of Theorem 3.28.

Figure 13:  An example of appendage subnetworks induced by appendage
homeostasis. Suppose the appendage component $\widetilde{\mathcal{A}}_{3}$ in
red is homeostasis inducing. Since $\widetilde{\mathcal{A}}_{3}$ has only one
path to the blue appendage component $\widetilde{\mathcal{A}}_{5}$ containing
the super-simple node $\rho_{j}$, then by Theorem 3.28
$\widetilde{\mathcal{A}}_{5}$ is homeostatic, but no other appendage
subnetwork is homeostatic.

######  Example 3.29.

We consider homeostasis patterns for the admissible systems in (2.40) obtained
from Figure 8. Specifically, we show how the theorems in this section lead to
the determination of the homeostasis patterns that were derived by direct
calculation from the equations (2.42). The corresponding homeostasis pattern
network $\mathcal{P}$ is shown in Figure 14. The answer is listed in Table 2.

Case (a) $f_{\tau_{3},\tau_{3}}$: homeostasis is induced by the node
$\widetilde{\mathcal{A}}_{3}$ of $\mathcal{P}_{\mathcal{A}}$. Theorem 3.27
shows $\widetilde{\mathcal{A}}_{3}$ induces $\\{o\\}$, which is the only
super-simple node of $\mathcal{P}_{\mathcal{S}}$ downstream from
$\mathcal{V}_{max}(\widetilde{\mathcal{A}}_{3})$. And
$\widetilde{\mathcal{A}}_{3}$ induces no backbone node. Theorem 3.28 shows
$\widetilde{\mathcal{A}}_{3}$ induces no appendage component. Therefore, in
this case the homeostasis pattern is $\\{o\\}$.

Case (b) $f_{\tau_{1},\tau_{1}}$: homeostasis is induced by the node
$\widetilde{\mathcal{A}}_{1}$ of $\mathcal{P}_{\mathcal{A}}$. Theorem 3.27
shows $\widetilde{\mathcal{A}}_{1}$ induces $\\{\iota,o\\}$, which are the
super-simple nodes of $\mathcal{P}_{\mathcal{S}}$ downstream from
$\mathcal{V}_{max}(\widetilde{\mathcal{A}}_{1})$. Also
$\widetilde{\mathcal{A}}_{1}$ induces $\\{\widetilde{\mathcal{L}}_{1}\\}$,
which is the backbone node of $\mathcal{P}_{\mathcal{S}}$ downstream from
$\iota$. Theorem 3.28 shows $\widetilde{\mathcal{A}}_{1}$ induces
$\\{\widetilde{\mathcal{A}}_{2},\widetilde{\mathcal{A}}_{3}\\}$, which are the
nodes of $\mathcal{P}_{\mathcal{A}}$ downstream from
$\widetilde{\mathcal{A}}_{1}$ and each path contains the super-simple node $o$
with $\mathcal{V}_{max}(\widetilde{\mathcal{A}}_{1})\preceq
o\preceq\mathcal{V}_{min}(\widetilde{\mathcal{A}}_{2})$ or
$\mathcal{V}_{min}(\widetilde{\mathcal{A}}_{3})$. Therefore, in this case the
homeostasis pattern is $\\{\iota,\ \sigma,\ o,\ \tau_{2},\ \tau_{3}\\}$.

Case (c) $f_{\tau_{2},\tau_{2}}$: homeostasis is induced by the node
$\widetilde{\mathcal{A}}_{2}$ of $\mathcal{P}_{\mathcal{A}}$. Theorem 3.27
shows $\widetilde{\mathcal{A}}_{2}$ induces $\\{o\\}$, which is the super-
simple nodes of $\mathcal{P}_{\mathcal{S}}$ downstream from
$\mathcal{V}_{max}(\widetilde{\mathcal{A}}_{2})$. And
$\widetilde{\mathcal{A}}_{2}$ induces no backbone node. Theorem 3.28 shows
$\widetilde{\mathcal{A}}_{2}$ induces $\\{\widetilde{\mathcal{A}}_{3}\\}$,
which is the node of $\mathcal{P}_{\mathcal{A}}$ downstream from
$\widetilde{\mathcal{A}}_{2}$ and each path contains the super-simple node $o$
with $\mathcal{V}_{max}(\widetilde{\mathcal{A}}_{2})\preceq
o\preceq\mathcal{V}_{min}(\widetilde{\mathcal{A}}_{3})$. Therefore, in this
case the homeostasis pattern is $\\{o,\ \tau_{3}\\}$.

Case (d): $\det\left(\begin{smallmatrix}f_{\sigma,\iota}&f_{\sigma,\sigma}\\\
f_{o,\iota}&f_{o,\sigma}\end{smallmatrix}\right)$: homeostasis is induced by
node $\widetilde{\mathcal{L}}_{1}$ of $\mathcal{P}_{\mathcal{S}}$. Theorem
3.25 shows $\widetilde{\mathcal{L}}_{1}$ induces $\\{o\\}$, which is the node
of $\mathcal{P}_{\mathcal{S}}$ downstream from $\widetilde{\mathcal{L}}_{1}$.
Theorem 3.26 shows $\widetilde{\mathcal{L}}_{1}$ induces
$\\{\widetilde{\mathcal{A}}_{2},\widetilde{\mathcal{A}}_{3}\\}$, whose
$\mathcal{V}_{min}$ are strictly downstream from
$\widetilde{\mathcal{L}}_{1}$. Therefore, in this case the homeostasis pattern
is $\\{o,\ \tau_{2},\ \tau_{3}\\}$.  $\Diamond$

Figure 14: Pattern network $\mathcal{P}$ obtained from input-output network
$\mathcal{G}$ shown in Figure 8.

Finally, we give three general results about the induction relation. First, in
Theorem 3.30 we state that the induction relation is characterized by its
behavior on homeostasis subnetworks. Next, in Theorems 3.31 we state that
induction applies in at least one direction for distinct homeostasis
subnetworks and that no subnetwork induces itself. Finally, in Theorem 3.32 we
state that distinct subnetworks have distinct homeostasis patterns. That is,
the set of nodes induced by a homeostasis subnetwork is unique among all
homeostasis subnetworks. For proofs see [40].

######  Theorem 3.30.

Suppose $\mathcal{K}_{1}$ and $\mathcal{K}_{2}$ are distinct homeostasis
subnetworks of $\mathcal{G}$. Let $\kappa$ be a node of $\mathcal{K}_{2}$
where $\kappa\neq o$. If $\mathcal{K}_{1}\Rightarrow\kappa$, then
$\mathcal{K}_{1}\Rightarrow\mathcal{K}_{2}$.

######  Theorem 3.31.

Let $\mathcal{K}_{1}$ be a homeostasis subnetwork of $\mathcal{G}$. Then
generically $\mathcal{K}_{1}\not\Rightarrow\mathcal{K}_{1}$. Moreover, let
$\mathcal{K}_{2}$ be some other homeostasis subnetwork of $\mathcal{G}$. Then
one of the following relations holds:

  1. [(a)]

  2. 1.

$\mathcal{K}_{1}\Rightarrow\mathcal{K}_{2}$ and
$\mathcal{K}_{2}\not\Rightarrow\mathcal{K}_{1}$,

  3. 2.

$\mathcal{K}_{2}\Rightarrow\mathcal{K}_{1}$ and
$\mathcal{K}_{1}\not\Rightarrow\mathcal{K}_{2}$,

  4. 3.

$\mathcal{K}_{1}\Rightarrow\mathcal{K}_{2}$ and
$\mathcal{K}_{2}\Rightarrow\mathcal{K}_{1}$.

Theorems 3.30 and 3.31 imply that each homeostasis subnetwork has a unique
homeostasis pattern associated to it. Specifically:

######  Theorem 3.32.

Let $\mathcal{K}_{1}$ and $\mathcal{K}_{2}$ be two distinct homeostasis
subnetworks of $\mathcal{G}$. Then the set of subnetworks induced by
$\mathcal{K}_{1}$ and the set of subnetworks induced by $\mathcal{K}_{2}$ are
distinct.

Theorem 3.32 implies that there is a bijective correspondence between
homeostasis subnetworks of $\mathcal{G}$ and homeostasis patterns supported by
$\mathcal{G}$. Recall that theorem 3.12 implies that there is a bijective
correspondence between the homeostasis types of $G$, that is, the irreducible
factors of $\det(H)$ (or the diagonal blocks of $H$) and homeostasis
subnetworks of $\mathcal{G}$. Therefore we can summarize the main result of
our theory as bijective correspondence between three sets defined for any
input-output network $\mathcal{G}$:

$\begin{split}\\{\text{Homeostasis types}\\}&\simeq\\{\text{Homeostasis
subnetworks}\\}\\\ &\simeq\\{\text{Homeostasis patterns}\\}\end{split}$

##  4 Applications

In this section we consider four examples of applications of the theory of
[134, 40] presented in the previous section. The first application is to gene
regulatory networks, following the approach in [7, 6]. The second application
is to intracellular metal ion regulation, following the approach of [5]. The
third application is to bacterial chemotaxis. Here, the theory of [134] needs
to be adapted to the case of multiple input nodes (with one input parameter)
and this was accomplished in [92]. In the last application we consider a
generalization of the notion of homeostasis to the periodic case and apply
this to the modeling of period homeostasis in the circadian rhythm.

###  4.1 Genetic Regulatory Networks

Gene expression is the process by which the information encoded in a gene is
turned into a biological function, which ultimately manifests itself as a
phenotype effect. This is accomplished by a complex series of enzymatic
chemical reactions within the cell leading to the synthesis of specific macro-
molecules called the gene product. The process of gene expression is used by
all known life – eukaryotes (including multicellular organisms), prokaryotes
(bacteria and archaea), and even viruses – to generate the molecular machinery
for life.

There are, basically, two types of gene products: (i) for protein-coding genes
the gene product is a protein; (ii) non-coding genes, such as transfer RNA
(tRNA) and small nuclear RNA (snRNA), the gene product is a functional non-
coding RNA (ncRNA).

Regulation of gene expression, or simply gene regulation, is the range of
mechanisms that are used by cells to increase or decrease the amount of
specific gene products. Sophisticated schemes of gene expression are widely
observed in biology, going from triggering developmental pathways, to
responding to environmental stimuli.

A gene (or genetic) regulatory network (GRN) is a collection of molecular
regulators that interact with each other and with other substances in the cell
to govern the gene expression levels of mRNA and proteins. The molecular
regulators can be DNA, RNA, protein or any combination of two, or more of
these three that form a complex, such as a specific sequence of DNA and a
transcription factor to activate that sequence. The interaction can be direct
or indirect (through transcribed RNA or translated protein). When a protein
acts as a regulator it is called a transcription factor, which is one of the
main players in regulatory networks. By binding to the promoter region of a
coding gene they turn them on, initiating the production of another protein,
and so on. Transcription factors can be excitatory (activators) or inhibitory
(repressors).

The development of advanced experimental techniques in molecular biology is
producing increasingly large amounts of experimental data on gene regulation.
This, in turn, demands the development of mathematical modeling methods for
the study and analysis of gene regulation. Mathematical models of GRNs
describe both gene expression and regulation, and in some cases generate
predictions that support experimental observations. Formally, a GRN is
represented by a directed graph. Nodes represent the variables associated to
genes (e.g., mRNA and/or protein concentration) and directed links represent
couplings between genes (e.g., effect of one gene product on other genes). In
this review we will focus on the mathematical modeling of GRNs using coupled
systems of ordinary differential equations. As we will see in a moment, in the
coupled ODE setting there is an important issue concerning the number of
variables / equations associated to each node.

Even though the notion of homeostasis is often associated with regulating
global physiological parameters, such as temperature, hormone levels, or
concentrations of molecules in the bloodstream in complex multicellular
organisms, it also can be applied to unicellular organisms. However, the issue
here is how some internal cell state of interest (such as the concentration of
some gene product) responds to changes in the intra-cellular and/or extra-
cellular environment [89, 125, 120]. For instance, Antoneli et al. [7] study
the occurrence of homeostasis in a _feedforward loop_ motif from the GRN of
_S. cerevisiae_ (see Figure 15).

In [7] the authors use theorem 2.18 and bare hands calculations to find
infinitesimal homeostasis in a small $3$-node GRN called ‘feedforward loop
motif’ (see Figure 15). Assuming that the regulation of the three genes is
inhibitory (repression) and that only gene SPF1 is regulated by upstream
transcription factors (the input parameter), they show that the protein
concentration of gene GAP1 (the output node) robustly exhibits infinitesimal
homeostasis with respect to variation on the regulation level of SPF1 over a
wide range. Moreover, SPF1 and GZF3 (i.e., their protein concentrations) are
not homeostatic for any value of the input parameter. Here, ‘robustly’ means
that the occurrence of infinitesimal homeostasis (or not) on the protein
concentrations of the three genes described above is persistent under
variation of kinetic parameters of the defining differential equations (the
rates of synthesis and degradation of the mRNA and protein concentrations). In
order to obtain compatibility of the infinitesimal homeostasis formalism with
the GRN structure they use the protein-mRNA network (PRN) representation.
Moreover, in this particular setting, they were able to explicitly compute the
homeostasis point by assuming a special functional form for the equations.

Figure 15: Feedforward loop motif from the GRN of _S. cerevisiae_ (see
Antoneli et al.[7]). Note that autoregulation implies that the protein
concentration associated with the gene SFP1 affects directly the mRNA
concentration associated with that gene.

Antoneli et al. [6] employ the theory for input-output networks explained
before to deal with infinitesimal homeostasis in general gene regulatory
networks (GRN). As we will explain in a moment a GRN is not exactly a
biochemical network as in [115]. There is a ‘mismatch’ between the number of
nodes in the network and the number of state variables of the underlying
system of ODEs. In order to resolve this ‘mismatch’, we generalize the
approach of [7], used to analyze feedforward loops, to arbitrary GRNs. The
idea there was to replace each ‘gene’ node of a GRN by a pair of ‘protein’ and
‘mRNA’ nodes to obtain a protein-mRNA network (PRN). Now, PRNs have a
mathematical structure similar to that of biochemical networks and hence the
theory of [134, 40] can be readily applied. More importantly, since the
classification results have a purely combinatorial side, we can consider the
GRN and its associated PRN simultaneously, and work out a correspondence
between the classifications of homeostasis types and homeostasis patterns in
both of them. Even though infinitesimal homeostasis only makes sense on the
PRN (dynamical) level, its purely combinatorial aspects can be transferred to
the GRN.

The main result of [6] is a complete characterization of homeostasis types and
homeostasis patterns on the PRN that have correspondent on the GRN. A
byproduct of this characterization is the discovery of homeostasis types and
homeostasis patterns on the PRN that do not have counterparts on the GRN. The
novelty of the approach of [6] is the simultaneous use of two networks, the
GRN and the PRN, in the analysis of gene expression homeostasis, and the lack
of assumptions about the functional form of the differential equations.

At the abstract level a GRN is a directed graph whose nodes are the genes and
a directed link from a source gene to a target gene indicates that the gene
product of the source gene acts as a molecular regulator of the target gene.
Autoregulation occurs when a link connects a gene to itself, that is, the gene
product is a molecular regulator of the gene itself.

At the dynamical level, a gene (a node in the GRN) represents the collection
of processes that ultimately lead to the making of the gene product. A
protein-coding gene should have at least two processes: (i) transcription,
that is, the synthesis of a mRNA copy from a DNA strand and (ii) translation,
that is, the synthesis of protein from a mRNA. The ‘output’ of the
corresponding node in the GRN is the protein concentration at a given time.

Lut us assume the simplest scenario, namely, a GRN containing only protein-
coding genes. Then, each gene (node) represents two concentrations (the
concentration of mRNA and the concentration of protein). That is, there is a
mismatch between the number of nodes $N$ and the number of state variables