A Model of the Cellular Iron Homeostasis Network Using
Semi-Formal Methods for Parameter Space Exploration
Nicolas Mobilia

Alexandre Donzé

UJF-Grenoble 1 / CNRS
TIMC-IMAG UMR 5525,
Grenoble, F-38041, France
nicolas.mobilia@imag.fr

EECS Department
University of California Berkeley
Berkeley, CA 94720 USA
donze@eecs.berkeley.com

Jean Marc Moulis

Éric Fanchon

UJF-Grenoble 1 / CNRS UMR 4952,
Institut de Recherches en Technologies et Sciences pour le Vivant,
17 rue des Martyrs,
38054 Grenoble
jean-marc.moulis@cea.fr

UJF-Grenoble 1 / CNRS
TIMC-IMAG UMR 5525,
Grenoble, F-38041, France
eric.fanchon@imag.fr

This paper presents a novel framework for the modeling of biological networks. It makes use of
recent tools analyzing the robust satisfaction of properties of (hybrid) dynamical systems. The main
challenge of this approach as applied to biological systems is to get access to the relevant parameter
sets despite gaps in the available knowledge. An initial estimate of useful parameters was sought by
formalizing the known behavior of the biological network in the STL logic using the tool Breach.
Then, once a set of parameter values consistent with known biological properties was found, we tried
to locally expand it into the largest possible valid region. We applied this methodology in an effort to
model and better understand the complex network regulating iron homeostasis in mammalian cells.
This system plays an important role in many biological functions, including erythropoiesis, resistance
against infections, and proliferation of cancer cells.

1

Introduction

In cellular biology, the value of the parameters are uncertain (when a measurement does exist) or not
known at all (unmeasured or not computable from first principles). It is well known that values of kinetic
constants obtained from in vitro measurements on purified enzymes can be significantly different from
the in vivo values, due to interactions with other cell components, to regulatory modifications of enzymes,
to sequestration, to anomalous diffusion or still other phenomena. In addition, the values available from
the literature are heterogeneous (measurement performed in different species, or on different cell types, or
in different conditions). Two sets of measurements performed on identical cell types placed in supposedly
identical conditions can also be qualitatively different for multiple reasons (undetected heterogeneity
of cell populations, different batches of antibodies used for detection, for instance). The fact that the
available data and pieces of information are generally too scarce with respect to the size of dynamical cell
models forces the modeler to include every bit of available data, at the cost of increased heterogeneity
and thus increased uncertainty. But even with such an inclusive strategy the amount of data remains
inadequate to identify a fully instantiated model.
The complexities and heterogeneity of both the parameter values and the analytical expressions of the
terms entering the differential equations describing a biological system implies such a system results from
a family of dynamical systems, instead of a single one. The properties which formalize experimentally
observed behaviors should thus be robust with respect to the uncertainties mentioned above.
Ezio Bartocci and Luca Bortolussi (Eds.): HSB 2012
EPTCS 92, 2012, pp. 42–57, doi:10.4204/EPTCS.92.4

c N. Mobilia, A. Donzé, J.M. Moulis, É. Fanchon

N. Mobilia, A. Donzé, J.M. Moulis, É. Fanchon

43

These issues are considered in the present paper and the semi-formal method implemented to explore
the parameter space is applied to an important question in biology: how iron homeostasis is maintained
in mammalian cells? Iron is the most abundant transition metal required by most forms of life. The
largest part needed by animals ends up in red blood cells, but all cell types do use the metal [2]. Iron
is handled by specific molecules in animal physiology, so as to avoid the deleterious reactions the iron
ions catalyze, conversion of oxygen derivatives into reactive species (the reactive oxygen species, ROS)
prominent among them [15]. The mandatory requirement for iron and the negative impact the metal may
have on cellular function call for a tight regulation of its use [7]. Deregulation translates into cellular
oxidative damage: this links iron homeostasis to the redox balance in animal, and other, cells. These
connections underlie scores of (human) pathological situations in which these homeostatic systems (iron
and redox homeostasis) are perturbed.
Our approach rests on the following points: we deal explicitly with the parameter space, and we do
not focus on a single instantiation of the model parameters. When the parameter values are uncertain
due to the fact that they have been measured in different cell types or in different conditions, the data are
formally represented by inequalities. In other words, we do not assign the measured value directly to the
parameter but rather impose inequality constraints on that parameter. If there is no measurement at all,
bounds are chosen to define a physiologically reasonable domain. Observed behavioral properties can be
formalized using a temporal logic formalism. Then the behaviors of the parametrized system are explored
by sampling predefined volumes of parameter space (including initial conditions and model parameters),
and the instantiated models that satisfy the temporal logic formula are identified. As illustrated below
this very basic scheme needs some elaboration to be applied to the modeling of a complex, yet simplified,
biological phenomenon.
The biological network we have built for iron homeostasis is presented in section 2. The system of
differential equations is introduced in section 3 and its dynamical properties explored in the following
section.

2

Biological network

We are considering regulation of iron handling in cell types which acquire iron by endocytosis of the
transferrin-transferrin receptor complex (most of them) and express an iron exporter in the form of ferroportin. Proper supply of iron to animal cells is required to avoid deleterious effects of too much or too
little iron as such imbalance jeopardizes the proper function of cells and may lead to death. Regulation
of iron supply occurs at two levels, a systemic one and a local one organized around Iron Regulatory
Proteins. The latter are the main focus of this study since they directly act on the molecular nodes of the
considered network.

2.1

Species

The considered network is shown in Figure 1 A. It contains five species: Iron Regulatory Protein (IRP,
no difference is made between the two known proteins displaying this function), iron (Fe), ferritin (Ft),
the transferrin receptor (TfR1) and ferroportin (FPN1a). Moreover, transferrin (Tf) is a protein which
carries iron in the blood stream throughout the body, and which can bind zero, one, or two iron atoms. It
is the main provider of iron to most cells and, as such, it is an input of the network. The circle labeled by
Fe out is the iron exported out of the cell. The circle labeled by Fe cons represents the iron consumed
for the cell needs.

44

A Model of the Cellular Iron Homeostasis Network
-

-

TfR1

+

+

-

Fe

-

Ft
+

-

FPN1a

+

IRP

A)

B)

Figure 1: A) Schematic representation of the main biological processes involved in the cellular homeostatic control of iron.
The dashed arrows represent translation of mRNA into proteins. The arrow ending in an unfilled white
arrow represents amplification of translation, while the lines with a broken ending arrowhead represent
inhibition of translation. The lines ending with a combined perpendicular stroke and arrow represent iron
transport through membranes. The arrow leading to a pink circle indicates degradation. Finally, the two
regular arrows represent the loading/unloading of iron into/from the ferritins and the consumption of iron
for the cell needs. The green rounded rectangles represent proteins, the green parallelograms represent
mRNA, and the blue circles represent atoms. The yellow concave hexagon represents the transferrin
receptor. This diagram was drawn with the software CellDesigner [14].
B) Interaction graph defined by the equations. The signs on the arrows indicate whether the interaction
is positive or negative.
The IRP are proteins which can bind to a specific hairpin-like structure on mRNA, called Iron Responsive Element (IRE). Depending on the position of the IRE on the non-coding sequence of mRNA,
two main outputs may result:
• If the IRE is located at the beginning of the mRNA, in the 5’ untranslated region (UTR), the
binding of the IRP prevents the translation machinery from binding to the mRNA and translating
it, so the binding of the IRP leads to a reduced quantity of the corresponding protein;
• If the IRE is located at the end of the mRNA (3’-UTR), the binding of the IRP inhibits the degradation of this mRNA, leading to an increase in the quantity of translated protein.
In the presence of enough iron, the regulatory activity of the IRP is decreased by different molecular
mechanisms [24].
The iron species (Fe) represents the cellular iron which is available for cellular processes. The actual
form of this iron in biological cells is still a topic of discussion. A virtual species, which we call “Fe”, is
introduced in order to model the usable iron in the cell. If this amount of intracellular iron becomes too
low, the cell cannot maintain vital functions and dies.

N. Mobilia, A. Donzé, J.M. Moulis, É. Fanchon

45

Ferritin is a protein which stores iron [3]. The complex mechanisms of iron loading into and release
from ferritin are not considered here. The mRNA of both ferritin subunits contain an IRE on their 5’-UTR
regions, and the IRP inhibit their translation.
Transferrin receptors are located on the surface of the cell [1]. The fixation of transferrin on its
receptor leads to the endocytosis of the complex. This process involves several steps and releases iron
into the cell. Transferrin receptor mRNA contain five IRE in the 3’-UTR region, which is stabilized upon
IRP binding. The density of transferrin receptor at the cell surface is thus correlated with the IRP activity.
Finally, ferroportin (FPN1a) is an iron exporter located at the surface of some cells [20]. A form of
FPN1a contains an IRE on its 5’-UTR mRNA region, so the translation is inhibited by the presence of
IRP.

2.2

Behavior

A qualitative description of the system has been obtained through a large body of biological experiments
and we selectively summarize them here.
If the amount of iron is sufficient (iron-replete situation), the cell is in a stationary state. Assuming
that the iron provision stays constant, the concentration of each protein described in Figure 1 A does not
vary. In this state, the IRP activity is expected to be relatively low. In the case of too much iron (which
is not studied here), the saturation level of transferrin increases, the IRP activity reaches a minimum and
inhibits translation of the iron importer (TfR1) to minimize cellular iron input.
From the iron-replete stationary state, if the cells become iron-depleted, the inhibition of IRP activity
by iron decreases. Thus, the increasing IRP activity leads to an increase in the transferrin receptor
concentration and to a decrease in ferroportin concentration. This, in turn, leads to an augmentation in
the amount of iron entering the cell and a reduction in the amount of iron leaving the cell. Moreover, the
iron loaded in ferritin is released to convert stored iron into a form available for biosynthetic purposes.
This release allows the cell to face a temporary reduction in the level of iron (typically because of a
reduction in the iron input) for a certain amount of time. When ferritin is depleted, and if there is still
no iron input mediated by the receptor TfR1, the cell runs out of iron and enters in a cellular death
process. The IRP regulating system is effective as long as the iron level remains low. When the iron level
regains an acceptable level again, the IRP activity decreases, and the amount of filled ferritin goes up,
regenerating the stock of stored iron.

3

ODE modeling

In this section we describe the differential equations that we derived to specifically describe the evolution
of the concentration of each species over time. Some processes in the iron homeostasis network are in
fact a composition of many reactions, and the choice of an analytical expression to represent them is not
straightforward. This is the case of the iron entry mediated by TfR1, which involves endocytosis, and of
the loading of iron in ferritin, for example.

3.1

Iron equation

In the cell, when the level of iron becomes low, some iron-proteins release their iron content, by protein
degradation for instance, which contributes to refill the iron pool; however, only iron release from ferritin
is considered in our model without affecting the general validity of the reasoning. The equation defining

46

A Model of the Cellular Iron Homeostasis Network

the iron concentration is:
dFe
dFt
= kFe input · TfR1 · Tfsat − nFt ·
− kFe export · Fe · FPN1a − kFe cons · Fe
(1)
dt
dt
The import of iron into the cell is proportional to the concentration of transferrin receptor and to
the average amount of iron bound to transferrin (called transferrin saturation): Tfsat . The second term
describes the storage or release of iron due to the synthesis or degradation of ferritin. This variation
is equal to the number of iron atoms per ferritin molecule (aka nFt ) times the variation in the ferritin
concentration (we consider that (i) iron release occurs only when ferritin is degraded (ii) synthesized
ferritins are very quicky filled with iron). The third term describes the export of iron out of the cell as
proportional to the iron and FPN1a concentration. The export parameter is kFe export . Finally, within
the cell, iron is used for many purposes; this consumption of iron is represented by the last term in our
equation (kFe cons · Fe). For the sake of simplicity, we consider this consumption term (which represents a
whole set of cellular processes) to be proportional to the iron concentration. The export of iron by FPN1
forms that are not regulated by IRP can be accounted for by this term.

3.2

IRP equation

The equation describing the IRP concentration is:
dIRP
= kIRP prod − kFe→IRP · sig+ (Fe, θFe→IRP ) · IRP − kIRP deg · IRP
(2)
dt
We define a constant rate of production of IRP by kIRP prod . The concentration of iron regulates
the activity of IRP. The degradation of IRP is described by a constant basal term (kIRP deg · IRP) and an
iron-related term (kFe→IRP · sig+ (Fe, θFe→IRP ) · IRP). Then, if the iron level is significantly below the
threshold θFe→IRP , the degradation rate is kIRP deg · IRP. If the iron concentration is significantly above
this threshold, the degradation rate is (kFe→IRP + kIRP deg ) · IRP, where kFe→IRP is the parameter describing
the inhibition of IRP by iron.
The sigmoid function is defined such that:
sig+ (x, θ ) =

xn
xn + θ n

where x is a concentration variable, θ is the threshold and n defines the steepness of the sigmoid. In our
model, we use n = 30, which corresponds to a very steep sigmoid.

3.3

Ferritin equation

To model the temporal evolution of ferritin, we assume that each ferritin protein is rapidly filled with iron
atoms, right after its synthesis is completed. In other words the amount of empty ferritin is negligible.
The variable Ft thus represents the concentration of filled ferritin proteins within the cell. The equation
describing the evolution of the concentration of ferritin is:
dFt
= kFt prod − kIRP→Ft · sig+ (IRP, θIRP→Ft ) − kFt deg · Ft
(3)
dt
We consider a basal production rate described by kFt prod . The second term describes the regulation by
IRP. We use a sigmoidal regulation here, but we could consider other dynamics. If the activity of IRP
is above the threshold θIRP→Ft , the production rate is lowered by kIRP→Ft . It follows that kIRP→Ft has to
be lower or equal to kFt prod : as IRP inhibits translation, this inhibition should never lead to a negative
production rate. Finally, the third term (kFt deg · Ft) describes the spontaneous degradation of ferritin.

N. Mobilia, A. Donzé, J.M. Moulis, É. Fanchon

3.4

47

Ferroportin equation

The equation describing the ferroportin concentration is very similar to that of ferritin. The equation is:
dFPN1a
= kFPN1a prod − kIRP→FPN1a · sig+ (IRP, θIRP→FPN1a ) − kFPN1a deg · FPN1a
(4)
dt
The first term describes the basal production of FPN1a. The second term expresses the regulation
of IRP on the translation of FPN1a. With the same reasoning as previously, the regulation parameter
kIRP→FPN1a has to be lower or equal to the production parameter kFPN1a prod . The dynamics of IRP
regulation on Ft and FPN1a are considered similar because both FPN1a and Ft mRNAs have a single
IRE located in the 5’-UTR region without experimental evidence of major differences between the two.

3.5

Transferrin receptor equation

The concentration of the transferrin receptor evolves according to the following equation:
dTfR1
= kTfR1 prod + kIRP→TfR1 · IRP − kTfR1 deg · TfR1
(5)
dt
The production of TfR1 includes a basal rate kTfR1 prod and it is increased proportionally to the IRP
concentration. This is taken into account in the term kIRP→TfR1 · IRP. This term is different from those
describing the regulation of ferritin and ferroportin in an effort to consider the binding of several IRPs to
a single TfR mRNA molecule (since it contains several IREs). The half-life of TfR1 mRNA is considered
proportional to the IRP concentration, although this is a rough approximation [13]. The last term of the
equation represents the spontaneous degradation.

3.6

Interaction graph

The interaction graph is a useful tool to visualize the pattern of interactions in a system of differential
equations. It is defined from the Jacobian matrix of the system. Each non-zero element in this matrix
indicates that a species influences another species. The nodes of the interaction graph are the dynamical
variables of the system (concentration of biological species), and each arc represents the influence of a
node on another, i.e. a non-zero element of the Jacobian matrix. The sign of this element determines
whether the influence is positive or negative. Of course it is possible that the sign of the elements depends
on the point of phase space at which the Jacobian is evaluated, but in our model the signs are constant.
The interaction graph we obtain is shown in Figure 1 B.
We define a negative (resp. positive) circuit (closed path) as a succession of arrows in the interaction
graph such that the product of the signs labelling these arrows is negative (resp. positive). This graph
contains three negative circuits of length greater than one (IRP→Fe→IRP ; IRP→FPN1a→Fe→IRP and
IRP→TfR1→Fe→IRP), and five negative loops corresponding to spontaneous degradation or consumption. It also contains one positive circuit ( IRP→Ft→Fe→IRP ). The presence of these circuits, together
with the presence of non-linear terms in the equations, generate a complex dynamical behavior.

4

Methodology

In this section we present our strategy to study the system. We first translate all available biological data
into inequalities on parameters and variables. Then, we formally define the expected behavior for the
two modes described in Section 2.2: the iron-replete mode and the iron-depleted mode. Finally we look
for a subspace which robustly verifies this behavior.

48

A Model of the Cellular Iron Homeostasis Network

4.1

Parameter search space

The parameter space of the system is twenty-dimensional. To define an initial parameter search space, we
combined published data from different sources as shown in table 1. Whenever possible, the available
data obtained with erythroleukemic cells have been used, as these cells are relatively immature and
proliferate, as considered for the initial steady state.
In mouse macrophages, half-life of FPN1a
seems to be more than 20 hours [18]
The production rate of TfR1 was estimated at
7×10−6 pg/cell/min as in [22]
The half-life of IRP is 12-15 hours in H1299 human lung cancer cells [11]
In human erythroleukemia cell line K562, halflife of TfR1 is 8 hours [25]
The TfR1 production rate enhanced by IRP is
6×10−5 pg/cell/min [22]
In human erythroleukemia cell line K562,
iron input is up to 7.2±2.4 × 104 Fe
atoms/cell/min [23]
In mammalian cells, the number of iron atoms
in ferritin up to 4500 [3]
Human normal saturation level of transferrin is
within 25% and 45% [21], we take 30%
In human erythroid cells, a fraction of the iron
exporter (FPN1) is not IRP-regulated [6]. The
amount of iron needed for biosynthesis (“consumed”) is higher than the exported one.
A higher limit for the iron pool in erythroleukemic cells is set at 2 µmol/L [12]
In rats liver, the IRP concentration in around
0.11 pmol/mg of cytosolic protein under normal
conditions [4]
In human erythroleukemic (K562) cells, there
are 140,000 transferrin receptors by cell [17]

0 < kFPN1a deg ≤ 9.6 × 10−6 s−1
1.0 × 10−13 ≤ kTfR1 prod ≤ 2.0 × 10−13 mol/L/s
1.28 × 10−5 ≤ kIRP deg ≤ 1.6 × 10−5 s−1
2.0 × 10−5 ≤ kTfR1 deg ≤ 3.0 × 10−5 s−1
4.2 × 10−5 ≤ kIRP→TfR1 ≤ 14.4 × 10−5 s−1
2 × 10−2 ≤ kFe input ≤ 3.9 × 10−2 s−1
0 < nFt ≤ 4500
Tfsat = 0, 3 (in iron-replete situation)
kFe cons ≥ kFe export · FPN1a (in iron-replete situation)
Fe ≤ 2 × 10−6 mol/L
3 × 10−9 ≤ IRP ≤ 10.7 × 10−9 mol/L (in ironreplete situation)
1.0 × 10−8 ≤ TfR1 ≤ 10.0 × 10−8 mol/L (in
iron-replete situation)

Table 1: The left column shows the data collected in the literature; the right column describes the corresponding interval on parameters and variables, and the inequalities deduced.
We then obtained intervals for eight parameters, including one in the iron-replete situation. Moreover,
we derived an inequality linking kFe cons , kFe export and FPN1a. Finally, we found data on the value of
three variables (Fe, IRP and TfR1) among which two are related to the iron-replete situation.
In order to evaluate the search space for the twelve other parameters, we proceeded by analogy. For
example, we had no information on the basal production rate of ferritin, ferroportin and IRP, but we
knew an interval for the production rate of the transferrin receptor. We then assumed that the basal production rate of ferroportin was neither more than 1000 times higher than the highest possible transferrin

N. Mobilia, A. Donzé, J.M. Moulis, É. Fanchon

49

receptor basal production rate, nor 1000 times lower than the lowest. Concerning the thresholds parameters, we simply assumed that their possible intervals were the same as for the corresponding protein
concentrations.

4.2

The iron-replete steady state

The iron-replete situation is modeled by a steady state, so we do not consider evolutions of the variables
over time. In this state, the value of the variable representing the iron concentration is expected to be significantly higher than the threshold θFe→IRP . We can thus infer that the sigmoidal term sig+ (Fe, θFe→IRP )
is very close to one, and we replace it. The value of the variable representing the IRP activity is also
expected to be largely below the thresholds θIRP→Ft and θIRP→FPN1a . We then consider that the two
sigmoidal functions sig+ (IRP, θIRP→Ft ) and sig+ (IRP, θIRP→FPN1a ) are very close to zero. These simplifications contribute to reducing the complexity of the system. Note that this amounts to approximating
the sigmoid by step functions, and consequently the system could be viewed as an hybrid system.
We developed two complementary approaches for studying the system in the iron-replete steady
state. The first consists in a mathematical analysis of the system. The second consists in performing
deductions from the quantitative data.
Looking for a steady state, we set the derivatives in the model to zero and solved the resulting
algebraic system using the symbolic solver Maxima [19]. This way, we could prove the existence of one
unique stationary point, given by Equations (6) to (10).

Fe =
Ft =
FPN1a =
IRP =
T f R1 =



T fsat · kFPN1a deg · kFe input · kT f R1 prod · kIRP deg + kFe→IRP + kIRP→T f R1 · kIRP prod


(6)
kFe export · kFPN1a prod + kFe cons · kFPN1a deg · kIRP deg + kFe→IRP · kT f R1 deg
kFt prod
(7)
kFt deg
kFPN1a prod
(8)
kFPN1a deg
kIRP prod
(9)
kIRP deg + kFe→IRP

kIRP deg + kFe→IRP · kT f R1 prod + kIRP→T f R1 · kIRP prod

(10)
kIRP deg + kFe→IRP · kT f R1 deg

We then considered the stability of this stationary point by computing its Jacobian matrix since an
unstable stationary point cannot obviously represent an observed cellular state. The obtained Jacobian
matrix is shown in equation (11).


−FPN1a · kFe export − kFe cons Tfsat · kFe input −Fe · kFe export kFt deg · nFt
0


0
−kTfR1 deg
0
0
kIRP→TfR1




0
0
−k
0
0
FPN1a
deg




0
0
0
−kFt deg
0
0
0
0
0
−kIRP deg − kFe→IRP
(11)
As this matrix is upper triangular, its eigenvalues are equal to its diagonal elements. It follows that all
eigenvalues are real and negative, which proves the stability of the stationary point.

50

A Model of the Cellular Iron Homeostasis Network

Another way to study the steady state consists in using the specific data we have on the variables at
the stationary state and propagate them on the parameter space. For this purpose, we use Realpaver [16],
an interval solver. The solver uses the algebraic equations between parameters and variables, as well as
defined intervals and constraints, as inputs for the calculation. The aim of Realpaver is to reduce the
provided intervals and consequently reduce the parameter space. The results are the following:
• kFPN1a deg ≥ 1.0 × 10−13 s−1
(initial interval: [0.0, 9.6 × 10−6 ] )
• kFPN1a prod ≤ 9.6 × 10−11 mol/L/s

(initial interval: [1.0 × 10−18 , 1.0 × 10−10 ] )

• kIRP prod ≥ 3.84 × 10−14 mol/L/s

(initial interval: [1.0 × 10−18 , 1.0 × 10−10 ] )

• kFe cons ≤ 3.39 × 10−4 s−1

(initial interval: [0.0, 1.0] )

• kFe→IRP ≤ 3.33 × 10−2 s−1
(initial interval: [1.0 × 10−9 , 1.0] )
Realpaver returns intervals which are reduced by many orders of magnitude for four of the parameters. These results prove that some subspaces of the parameter space do not contain valid parameter sets,
and then, need not be considered when searching for an parameter set producing an adequate iron-replete
steady state. This will reduce the search of the parameter space in the next step (analysis of the dynamical
behavior).
Besides, Realpaver makes no difference between variables and parameters considered in the algebraic
system, they are all variables (unknowns). Consequently it also propagates the constraints from the
parameters to the variables, and provides deductions on three variables in the iron-replete situation:
• FPN1a ≥ 1.04 × 10−13 mol/L
(initial interval: [1.0 × 10−13 , 1.0 × 10−5 ] )
• Fe ≥ 3.5 × 10−8 mol/L

(initial interval: [0.0, 2.0 × 10−6 ] )

• TfR1 ≤ 8.7 × 10−8 mol/L
(initial interval: [1.0 × 10−8 , 1.0 × 10−7 ] )
These results, provided by Realpaver in less than one second on a Core 2 Duo 3 GHz with 8 Gb of
RAM, can be considered as deductions made from the initial data. Indeed, the difficulty in finding a valid
parameter set is that steady state values must be within their intervals. These two approaches facilitate
the manual search of an initial parameter space. The analytical solution of the equations indicates, for
each parameter, the variables impacted by a change in the parameter value. We use this information to
define an order for checking variables, such that there are parameters allowing us to adjust a variable
value but not changing the value of the already correct variables. Consequently, to find an initial valid
parameter set, one has only to tune some parameters to fix the value of the first variable, then tune some
other parameters to fix the value of the second variable, and so on. If the value of one variable cannot
be set in its interval, one has to step back and change the value of parameters used to fix the value of
the previous variable. The variable order we use is: IRP, TfR1 (through the adjustment of kIRP→TfR1 ,
kTfR1 prod and kTfR1 deg ), FPN1a, Fe (through kFe cons , kFe export and kFe input ) and finally Ft. Using this
method, we found the following initial parameter set, denoted as P0 :
kTfR1 prod = 1.7 × 10−13 mol/L/s
kFPN1a deg = 5.0 × 10−6
s−1
kTfR1 deg = 2.4 × 10−5
s−1
kFe input = 3.0 × 10−2
s−1
−11
nFt = 400
kFt prod = 7 × 10
mol/L/s
kIRP prod = 8.0 × 10−12 mol/L/s
kFt deg = 5.0 × 10−3
s−1
kFe export = 300
L/mol/s
kFe→IRP = 1.0 × 10−3
s−1
−8
−7
θIRP→FPN1a = 3.0 × 10
mol/L
θFe→IRP = 1.5 × 10
mol/L
−13
−5
kIRP→FPN1a = 5.0 × 10
mol/L/s
kIRP deg = 1.4 × 10
s−1
−4
−1
−12
kIRP→TfR1 = 1.4 × 10
s
kFPN1a prod = 2.5 × 10
mol/L/s
kFe cons = 3.0 × 10−4
s−1
θIRP→Ft = 3.0 × 10−8
mol/L
−11
kIRP→Ft = 7 × 10
mol/L/s

N. Mobilia, A. Donzé, J.M. Moulis, É. Fanchon

51

It will be used as the initial parameter set for the exploration of the parameter space, as explained in
the following section.

4.3

Dynamics in iron-depleted situation

We are now interested in the behavior of the network when the input of iron is cut. As a consequence,
we need to compute the evolution of the system over time. For that, we use Breach [8, 9], a tool based
on Matlab/C++. This tool allows not only to simulate differential systems, but also to express temporal
logic formula and to check whether a trajectory satisfies or not a given formula.
We introduce briefly STL (Signal Temporal Logic). A more complete description can be found in [9].
Initially created for verifying the correctness of reactive programs, temporal logics provide languages to
express properties over time-dependent phenomena for any kind of dynamical system. STL and its
quantitative semantics [10] specializes into systems with continuous time and with real-valued states.
An STL formula has the following syntax:
ϕ = µ | (ϕ1 ) and (ϕ2 ) | ev_[a,b] (ϕ) | alw_[a,b] (ϕ)

(12)

where µ is an inequality constraint between expressions involving variables (written var[t]), constants
or derivatives of a variable over time (written ddt{var}[t]). A formula is evaluated at each requested
time (thereafter called “the current time”). A formula with a and statement is true if both ϕ1 and ϕ2 are
true at the current time. The ev_[a,b] statement is true if ϕ holds at least once between the current time
plus a and the current time plus b. If the interval is omitted, then a is 0 and b is inf, that is, the formula is
true if ϕ holds at least once, anytime after the current time. The statement alw_[a,b] is true if ϕ holds
between current time plus a and current time plus b. Again, if the interval is omitted then a is 0 and b is
inf meaning in that case that the formula must hold all the time.
The experiment we want to simulate is the cut-off of iron entry for cells initially in the iron-replete
situation. The total duration of the simulation is 48 hours, as the qualitative description of the irondepleted situation tells us that significant changes can be seen in 24 hours. First we want the simulations
to stabilize in a steady state corresponding to the iron-replete situation. Then, at time t=6 hours, we cut
the entry of iron by setting Tfsat to zero. We expect the system to turn on the IRP regulatory system in
order to maintain a correct iron level for at least ten hours. To identify parameter values for which the
simulations satisfy all theses properties, we use temporal logic formula.
We first want the system to be in a steady state corresponding to the iron-replete situation. To avoid an
inappropriately long stabilization phase, we specify initial conditions close to the value of the stationary
state. To ensure that the system reaches a steady state, we use the following temporal logic formula:
ϕS1 = abs(ddt{Fe}[t]

/ Fe[t])

< 1.0e-4

ϕS2 = abs(ddt{Ft}[t]

/ Ft[t])

< 1.0e-4

ϕS3 = abs(ddt{FPN1a}[t] / FPN1a[t]) < 1.0e-4
ϕS4 = abs(ddt{IRP}[t]

/ IRP[t])

< 1.0e-4

ϕS5 = abs(ddt{TfR1}[t]

/ TfR1[t])

< 1.0e-4

(13)

(To make the reading easier, all formula dealing with the iron-replete steady state will begin by “ϕS ”,
formula related to parameters will begin by “ϕP ” and formula related to the behavior will begin by “ϕB ”.)
To check that the value of each variable at the steady state is in the intervals shown in table 1, or deduced

52

A Model of the Cellular Iron Homeostasis Network

in 4.2, we use the following formula:
ϕS6 = (3.5e-8 < Fe[t])
ϕS7 = (3.0e-9 < IRP[t]) and (IRP[t] < 1.07e-8)
ϕS8 = (1.0e-8 < TfR1[t]) and (TfR1[t] < 8.7e-8)

(14)

ϕS9 = (1.04e-13 < FPN1a[t]) and (FPN1a[t] < 1.0e-5)
ϕS10 = (1.0e-13 < Ft[t]) and (Ft[t] < 1.0e-5)
As the maximal concentration of iron is 2 × 10−6 mol/L all the time, and not only at the steady state, we
do not enforce this property here.
Then, when the iron input is switched off, and the iron level becomes too low, the ferritins release
iron. This release is due to a lower ferritin production rate, which provokes a decrease of ferritin concentration. This can only be triggered by crossing the threshold θIRP→Ft . This agrees with the previous
statement that the variable describing IRP activity is lower than θIRP→Ft under iron-replete conditions.
We verify that property with the following formula:
ϕS11 = IRP[t] < theta_IRP_Ft

(15)

This increase of IRP activity can only be triggered by a decrease in iron concentration important enough
to cross the threshold θFe→IRP . This is in line with the already indicated inequality Fe > θFe→IRP under
iron-replete conditions. The corresponding formula is:
ϕS12 = Fe[t] > theta_Fe_IRP

(16)

Moreover, when the iron level becomes too low, the FPN1a concentration is expected to become lower.
This, also, can only be triggered by a crossing by IRP of the threshold θIRP→FPN1a . That means that in
the iron-replete steady state the FPN1a concentration is above this threshold, which is expressed by this
formula:
ϕS13 = IRP[t] < theta_IRP_FPN1a
(17)
Nevertheless , as explained in section 3.4 (resp. section 3.3), the regulation by IRP on FPN1a (kIRP→FPN1a )
(resp. on Ft (kIRP→Ft )) cannot be higher that the basal production rate of FPN1a (kFPN1a prod ) (resp. Ft
(kFt prod )). The formula imposing that are:
ϕP1 = k_IRP_FPN1a <= k_FPN1a_prod

(18)

ϕP2 = k_IRP_Ft <= k_Ft_prod

(19)

and
The activation of the regulatory system leads to the conservation of the iron level. We express that by
enforcing the existence of a plateau of this level after shutting the iron import off. We want this plateau
to exist for at least ten hours. To avoid incorrect simulations for which the concentration of this plateau
is zero, we set it at, at least, one tenth of the concentration in the iron-replete situation. The following
property expresses that:
ϕB1 = ev_[6*3600, inf] (alw_[0, 10*3600]
( (ϕS1 )

and

(Fe[t] > 0.01*Fe[4*3600]) ))

(20)

N. Mobilia, A. Donzé, J.M. Moulis, É. Fanchon

53

When the iron stored in ferritin is exhausted, the iron concentration eventually decreases to zero.
Ferritin disappearance precedes iron depletion which is not rescued by replenishment from any internal
store, in contrast to the transient phase immediately following removal of iron input from transferrin
(Figure 2 A. Otherwise, it would mean that the cell has stored iron but is not using it, which contradicts
the conditions set in 3.1. This implies that the parameter describing the strength of the regulation by the
IRP (i.e. kIRP→Ft ) is slightly lower or equal to the production rate (i.e. kFt prod ). We input this constraint
with the following formula:
ϕP3 = k_Ft_prod*0.95 <= k_IRP_Ft

(21)

Next, in table 1, two properties are expressed: the first enforces that the iron concentration is never
higher than 2 × 10−6 mol/L, and the second enforces that the inequality kFe cons ≥ kFe export · FPN1a
holds in iron-replete situation. The corresponding formula are:
ϕB2 = alw (Fe[t] < 2e-6)

(22)

ϕS14 = k_Fe_cons > k_Fe_export*FPN1a[t]

(23)

and
Finally, we have to define formula linking all the previous formula. For this purpose, we first define
the formula ϕSall which aggregates all formula related to the steady state. As we want to be sure that,
before the iron input is switched off, the system reaches and stays at the steady state, we enforce that
the properties related to the steady state hold at least for one hour before the sixth hour. This formula
expresses that:
ϕSall = ev_[0,6*3600] (alw_[0,3600] (ϕS1 and (ϕS2 and (ϕS3 and ( . . . and ϕS14 ))))
(24)
We also define the formula ϕBPall enforcing that all formula related to parameters or behavior are satisfied:
ϕBPall = ϕP1 and (ϕP2 and (ϕP3 and (ϕB1 and ϕB2 )))
(25)
And finally, the formula ϕall , which enforces that all properties are satisfied:
ϕall = (ϕSall ) and (ϕBPall )

(26)

It is noticeable that no information on TfR1 in the iron-depleted situation nor on the parameters
acting on it appears in this analysis. This is because, as we totally cut the entry of iron into the cell,
variations in the concentration of TfR1 do not change the amount of incoming iron.

4.4

A robust valid region

Having defined the formula which expresses the relevant properties, we search for a subspace where
all parameter sets verify these properties, by first finding a valid instantiation of the parameters and
then exploring a neighborhood of this instantiation. In section 4.2, we found an initial parameter set,
P0 , satisfying constraints on the data for the steady state. It turns out that the simulation computed
with this parameter set, shown in Figure 2 A, also verifies ϕall . As the temporal formula describes the
expected qualitative behavior, this simulation is representative of the behavior of all trajectories in the
region that we are looking for. We look around this point (in parameter space) to find a region exhibiting
the expected behavior. To ensure that the formula ϕP2 and ϕP3 are satisfied, we enforce that kIRP→Ft is

54

A Model of the Cellular Iron Homeostasis Network

Tf_sat

0.4
0.2

Fe

0
0
−6
x 10
2

IRP

20

30

40

50

10

20

30

40

50

10

20

30

40

50

10

20

30

40

50

10

20

30

40

50

10

20
30
time (hours)

40

50

1
0
0
−7
x 10
4
2
0
0
−8
x 10
2

Ft

10

1

FPN1a

0
0
−7
x 10
5

TfR1

4.5
0
−6
x 10
2
1
0
0

A)

B)

Figure 2: A) Characteristic simulation of the expected qualitative behavior. It shows that the system
switches to a new mode when iron is lacking. The parameter set used is P0 . At time t=6 hours, Tfsat is
set to zero. The iron concentration consequently decreases, leading to an increase of IRP activity. This
increase leads to a decrease in Ft concentration, keeping the iron level constant. This also leads to a
decrease in FPN1a concentration and to an increase in TfR1 concentration. During this phase when iron
is kept constant, the concentration of iron is close to the threshold θFe→IRP and IRP activity is close to
the thresholds θIRP→FPN1a and θIRP→Ft . Once Ft concentration is zero, the iron concentration tends also
to zero.
B) Histogram showing the sensitivity of iron and ferritin concentration to variations of parameter value.
These results are normalized: for example, a variation of 1% in kFe cons value leads to a variation of 0.2%
in Fe value. To avoid incorrect sensitivity due do an initial transient phase of large amplitude, we do not
consider the first three hours when computing the sensitivity. Moreover, as the behavior of the cell (i.e.
death) changes when iron is totally exhausted, we consider the sensitivity only between t=3 hours and
t=20 hours.
computed equal to 0.97 × kft prod . Moreover, if the parameter kFPN1a prod is lower than kIRP→FPN1a , we
switch them to force that the formula ϕP1 is satisfied. So, the dimension of the search space is slightly
reduced.
To guide the search, we compute the sensitivity of the system to parameters. We show in Figure 2 B
the sensitivity of iron and ferritin concentration to a variation of parameters value. Sensitivity to ferro-

N. Mobilia, A. Donzé, J.M. Moulis, É. Fanchon

55

portin, transferrin receptor and IRP are also computed, but are not shown. We define the sensitivity to a
parameter as the maximal sensitivity of the five variables of this parameter. This sensitivity is used in the
following manner: if the system is very sensitive to a parameter, we will explore a close neighborhood
of this parameter. On the opposite, if the system sensitivity to a parameter is low, we will explore a large
neighborhood. Nevertheless, it has to be noted that the system can be very sensitive to a parameter, and
that even with large variations on this parameter the qualitative property is still satisfied. This happens
when the sensitivity which is detected affects an aspect of the trajectories which is not captured by the
property of interest.
To verify that all simulations corresponding to a parameter subspace exhibit the expected behavior,
we performed intensive trajectory computations. We randomly picked 10000 parameter sets within this
subspace and checked the truth value of the formula ϕall . By expanding as much as possible the parameter
interval, we found an hyper-rectangle centered on P0 such that all parameter sets in this region exhibit
the expected behavior. Some parameters of P0 could change up to 70% of their initial value, with an
average of 34% for all, while conserving the behavior of the system.
We observe that for the whole region of parameter space characterized in the iron-depleted situation,
the FPN1a concentration remains at a rather high level, implying the considered cell keeps its ability to
export iron despite shortage of the metal.

5

Conclusion

Two models of the core iron homeostasis system have been published previously. Omholt et al [22] aim
at building a “unifying meta-model”, that is a model which is not tied to a particular cell type. They
include explicitly the two different Iron Regulatory Proteins, IRP1 and IRP2. Both proteins have the
same mRNA targets, so they seem to be redundant, but they respond to the iron level in the cell through
different molecular mechanisms. The model of Omholt et al [22] does not explicitly include ferritin nor
ferroportin. In addition, all interactions in the network are switch-like (sigmoidal kinetics). A numerical
value was assigned to all parameters except the sigmoid steepness, and the stability of the unique steady
state was monitored as the sigmoid steepness is varied. The more recent paper of Chifman et al [5]
is closer to our approach. This model considers the same five biological species, but the analytical
expression of their differential system is different from ours. The major differences are for the IRP,
ferritin and ferroportin equations (note that the graph shown in figure 1 of [5] does not follow the formal
definition given here for an interaction graph). The mathematical analysis of the parametrized model
of [5] is general in the sense that it does not rely on parameter instantiation. They showed that the model
produces a single steady state, the stability of which was studied by performing simulations with sampled
parameter values.
Here we have chosen an approach intermediate between Omholt et al [22] and Chifman et al [5]: a
parametrized model has been built without imposing numerical values to the parameters. The available
experimental data allowed us to apply inequality constraints between the parameters, together with behavioral constraints expressed in a temporal logic formalism. We could thus study not only steady states,
but also dynamical responses to perturbations, as, for instance, the sharp decrease of iron input, starting
from an iron-replete cell state. The region of parameter space for which the model exhibits the observed
behavior has been characterized. It appeared that, in the iron-depleted situation, the FPN1a concentration
remains at a rather high level, implying the considered cell kept its ability to export iron despite shortage
of the metal. This is a significant outcome of modeling which justifies that another regulatory mechanism
than that of the IRP may be needed under these conditions. This additional regulation may be provided

56

A Model of the Cellular Iron Homeostasis Network

by FPN1a degradation by hepcidin.
We plan to extend the present model. For example the influence of iron on the inhibition of IRP
activity involves the iron cluster biosynthesis pathway which will be introduced in the model. Also, the
regulatory mechanisms of IRP could be described more precisely by integrating mRNA concentrations
in the model. These extensions will be performed in a controlled and stepwise manner. The qualitative
behavior will be compared with that of the core model presented here. This will give an insight into the
robustness of the dynamical property of interest.

Acknowledgments
This work was supported by Microsoft Research through its PhD Scholarship Programme.
The authors thank the Région Rhône-Alpes for financial support within the Cible 2010 Program.

References
[1] P. Aisen (2004): Transferrin receptor 1. The International Journal of Biochemistry &; Cell Biology 36(11),
pp. 2137–2143, doi:10.1016/j.biocel.2004.02.007.
[2] N. C. Andrews (2008): Forging a field: the golden age of iron biology. Blood 112(2), pp. 219–230,
doi:10.1182/blood-2007-12-077388.
[3] P. Arosio & S. Levi (2010): Cytosolic and mitochondrial ferritins in the regulation of cellular iron homeostasis and oxidative damage. Biochim. Biophys. Acta 1800(8), pp. 783–792, doi:10.1016/j.bbagen.2010.02.005.
[4] O. S. Chen, K. P. Blemings, K. L. Schalinske & R. S. Eisenstein (1998): Dietary iron intake rapidly influences
iron regulatory proteins, ferritin subunits and mitochondrial aconitase in rat liver. J. Nutr. 128(3), pp. 525–
535.
[5] J. Chifman, A. Kniss, P. Neupane, I. Williams, B. Leung, Z. Deng, P. Mendes, V. Hower, F. M. Torti, S. A.
Akman, S. V. Torti & R. Laubenbacher (2012): The core control system of intracellular iron homeostasis: a
mathematical model. J. Theor. Biol. 300, pp. 91–99, doi:10.1016/j.jtbi.2012.01.024.
[6] L. Cianetti, P. Segnalini, A. Calzolari, O. Morsilli, F. Felicetti, C. Ramoni, M. Gabbianelli, U. Testa & N. M.
Sposi (2005): Expression of alternative transcripts of ferroportin-1 during human erythroid differentiation.
Haematologica 90(12), pp. 1595–1606.
[7] I. De Domenico, D. McVey Ward & J. Kaplan (2008): Regulation of iron acquisition and storage: consequences for iron-linked disorders. Nat. Rev. Mol. Cell Biol. 9(1), pp. 72–81, doi:10.1038/nrm2295.
[8] A. Donzé (2010): Breach, A Toolbox for Verification and Parameter Synthesis of Hybrid Systems. In: CAV,
pp. 167–170, doi:10.1007/978-3-642-14295-6 17.
[9] A. Donzé, É. Fanchon, L. M. Gattepaille, O. Maler & P. Tracqui (2011): Robustness Analysis and Behavior Discrimination in Enzymatic Reaction Networks.
PLoS ONE 6(9), p. e24246,
doi:10.1371/journal.pone.0024246.
[10] A. Donzé & O. Maler (2010): Robust Satisfaction of Temporal Logic over Real-Valued Signals. In: FORMATS, pp. 92–106. Available at http://dx.doi.org/10.1007/978-3-642-15297-9_9.
[11] C. Dycke, P. Charbonnier, K. Pantopoulos & J. M. Moulis (2008): A role for lysosomes in the
turnover of human iron regulatory protein 2. Int. J. Biochem. Cell Biol. 40(12), pp. 2826–2832,
doi:10.1016/j.biocel.2008.05.015.
[12] S. Epsztejn, O. Kakhlon, H. Glickstein, W. Breuer & I. Cabantchik (1997): Fluorescence analysis of the
labile iron pool of mammalian cells. Anal. Biochem. 248(1), pp. 31–40.

N. Mobilia, A. Donzé, J.M. Moulis, É. Fanchon

57

[13] R. Erlitzki, J. C. Long & E. C. Theil (2002): Multiple, conserved iron-responsive elements in the 3’untranslated region of transferrin receptor mRNA enhance binding of iron regulatory protein 2. J. Biol.
Chem. 277(45), pp. 42579–42587, doi:10.1074/jbc.M207918200.
[14] . Funahashi, M. Morohashi, H. Kitano & N. Tanimura (2003): CellDesigner: a process diagram editor for gene-regulatory and biochemical networks. BIOSILICO 1(5), pp. 159–162, doi:10.1016/S14785382(03)02370-9.
[15] D. Galaris & K. Pantopoulos (2008): Oxidative stress and iron homeostasis: mechanistic and health aspects.
Crit Rev Clin Lab Sci 45(1), pp. 1–23, doi:10.1080/10408360701713104.
[16] L. Granvilliers & F. Benhamou (2006): Algorithm 852: RealPaver: an interval solver using constraint satisfaction techniques. ACM Transactions on Mathematical Software 32(1), pp. 138–156,
doi:10.1145/1132973.1132980.
[17] R. C. Hunt & L. Marshall-Carlson (1986): Internalization and recycling of transferrin and its receptor. Effect
of trifluoperazine on recycling in human erythroleukemic cells. J. Biol. Chem. 261(8), pp. 3681–3686.
[18] M. D. Knutson, M. Oukka, L. M. Koss, F. Aydemir & M. Wessling-Resnick (2005): Iron release from
macrophages after erythrophagocytosis is up-regulated by ferroportin 1 overexpression and down-regulated
by hepcidin. Proc. Natl. Acad. Sci. U.S.A. 102(5), pp. 1324–1328, doi:10.1073/pnas.0409409102.
[19] Maxima (2011): Maxima, a Computer Algebra System. Version 5.24.0. Available at http://maxima.
sourceforge.net/.
[20] R. Mayr, A. R. Janecke, . Schranz, W. J.H. Griffiths, W. Vogel, A. Pietrangelo & H. Zoller (2010): Ferroportin
disease: A systematic meta-analysis of clinical and molecular findings. Journal of Hepatology 53(5), pp.
941–949, doi:10.1016/j.jhep.2010.05.016.
[21] R. Moirand, A. M. Mortaji, O. Loréal, F. Paillard, P. Brissot & Y. Deugnier (1997): A new syndrome of liver
iron overload with normal transferrin saturation. The Lancet 349(9045), pp. 95–97, doi:10.1016/S01406736(96)06034-5.
[22] S. W. Omholt, X. Kefang, Ø. Andersen & E. Plahte (1998): Description and Analysis of Switchlike Regulatory
Networks Exemplified by a Model of Cellular Iron Homeostasis. Journal of Theoretical Biology 195(3), pp.
339–350, doi:10.1006/jtbi.1998.0800.
[23] M. Shvartsman, E. Fibach & Z. I. Cabantchik (2010): Transferrin-iron routing to the cytosol and
mitochondria as studied by live and real-time fluorescence.
Biochem. J. 429(1), pp. 185–193,
doi:10.1042/BJ20100213.
[24] J. Wang & K. Pantopoulos (2011): Regulation of cellular iron metabolism. Biochem. J. 434(3), pp. 365–381,
doi:10.1042/BJ20101825.
[25] A. M. Weissman, R. D. Klausner, K. Rao & J. B. Harford (1986): Exposure of K562 cells to anti-receptor
monoclonal antibody OKT9 results in rapid redistribution and enhanced degradation of the transferrin receptor. The Journal of Cell Biology 102(3), pp. 951–958, doi:10.1083/jcb.102.3.951.

