Modelling of cell choice between differentiation

and apoptosis on the basis of intracellular and

extracellular regulations and stochasticity

M. Banerjee1, V. Volpert2

1 Department of Mathematics and Statistics, Indian Institute of Technology
Kanpur

Kanpur 208016, India, E-mail<EMAIL_ADDRESS>
2 Institut Camille Jordan, UMR 5208 CNRS, University Lyon 1

69622 Villeurbanne, France, E-mail<EMAIL_ADDRESS>
Abstract. The work is devoted to the analysis of cell population dynamics
where cells make a choice between differentiation and apoptosis. This choice
is based on the values of intracellular proteins whose concentrations are
described by a system of ordinary differential equations with bistable
dynamics. Intracellular regulation and cell fate are controlled by the
extracellular regulation through the number of differentiated cells. Initial
intracellular protein concentrations are considered for each cell as random
variables with a given area of variation. Intracellular regulation,
extracellular regulation and random initial conditions work together to
produce differentiated cells and to control their number. The role of
intracellular regulation is to provide a possible choice between
differentiation and apoptosis, extracellular regulation controls the number of
differentiated cells, stochastic initial conditions can suppress oscillations
and provide stability of the system.

Key words: Cell differentiation, apoptosis, extra cellular regulation, cell
fate.

##  1 Introduction

###  1.1 Cell fate and multi-scale modelling

Cell population dynamics is determined by the equilibrium between self-
renewal, differentiation and apoptosis. Cell fate depends on complex
intracellular and extracellular regulations which can act on the level of the
whole cell population and not only on individual cells [28], [29], [30], [35].
Multi-scale modelling provides an appropriate tool to study cell population
dynamics. There are many different interpretations of multi-scale modelling in
biology, as well as many models and applications (see [2], [4], [16], [36],
[42] and the references therein). In the case of physiological processes, by
multi-scale models we will understand the models which consider some
populations of cells (tissue), intracellular regulation of cell fate and
possibly of some its properties, extracellular regulation by the surrounding
cells and/or by other organs and tissues.

We will distinguish local and global extracellular regulation. By local
regulation we will understand regulation from the surrounding cells of the
same tissue. It can occur by the direct cell-cell contact or by intermediate
of various molecules produced by cells of the considered tissue and diffusing
in the extracellular matrix. In particular local regulation occurs in stem
cell niche. Stem cells can self-renew only being in contact with other cells
in the niche (e.g., stromal cells). When they lose this contact, they
differentiate and cannot self-renew any more. This local control of self-
renewal is important to prevent uncontrolled self-renewal leading to cancer.

Global extracellular regulation is effectuated by other organs and tissues. In
particular, this can be hormones, nutrients, growth factors. Global regulation
takes into account some overall information about the tissue and ignores the
details about local cell distribution inside it. This information can be based
on the number of cells of certain type. For example, production of the hormone
erythropoietin in the kidney depends on the number of erythrocytes in blood
(more precisely, on hemoglobin). Erythropoietin arrives to the bone marrow and
regulates production of erythrocytes [25]. We will return to this example
below.

Hybrid discrete-continuous models are well adapted for multi-scale modelling
in biology. They consider cells as individual objects, intracellular
regulatory networks can be described by ordinary differential equation and
biochemical substances in the extracellular matrix by partial differential
equations. They can be completed by various models of global extracellular
regulation. Such models are applied to study various physiological processes
[21] and dynamics of cell populations [7], [8], [43].

Among other applications, hybrid models are developed to study hematopoiesis
(blood cell production) [7], [9]. Hematopoiesis represents an interesting
example to study dynamics of cell populations. It is a complex process which
begins with hematopoietic stem cells and results in production of red blood
cells, platelets and leucocytes. During this process, cells undergo many
consecutive steps of transformation controlled by numerous local and global
regulations. The model developed and studied in this work is inspired by
hematopoiesis but it admits also some other interpretations.

Traditionally mathematical models of hematopoesis are based on ordinary and
delay differential equations and on structured cell dynamics (transport
equations). Dynamics of hematopoietic stem cells was studied in [31], [32]
(see also [38] and references therein). Periodic oscillations in blood cell
count and various related blood diseases were studied in [3], [5], [13], [39].
A model of platelet production was developed in [18], [44]. Mathematical
modelling of erythropoiesis and anemia began in [45], [47], [46]. Age-
structured models of erythropoiesis which take into account the action of
erythropoietin are considered in [1], [15], [33]. Interaction of different
cell lineages was studied in [12], [13]. There are numerous works devoted to
modelling of leukemia development and treatment [6], [23], [24], [37], [41].

It is important to indicate that in all these works the rates of self-renewal,
differentiation and apoptosis are imposed as given constants or functions of
some concentrations. However usually these values are not known or, what is
more important, they even cannot be considered as given. Indeed, for the same
concentrations of hormones, drugs, nutrients and other extracellular
substances their action can depend on the concentrations of intracellular
proteins and, more general, on the cell state and history. Therefore we need
to take into account intracellular regulation and to introduce multi-scale
models.

Multi-scale hybrid models of hematopoiesis were developed in [10], [15], [17],
[20], [14], [26], [27]. These models descrbe self-renewal, differentiation and
apoptosis of erythroid progenitors on the basis of intracellular and
extracellular regulation. The intracellular regulation takes into account the
main proteins involved in this process. The extracellular regulation is based
on the interaction of Fas-lignad producing cells which stimulate
differentiation and apoptosis, on the interaction with macrophages which
promote self-renewal of erythroid progenitors and on the influence of
erythropoietin which downregulates apoptosis of erythroid progenitors. This
approach allowed a detailed description of erythropoiesis in normal and
pathological situations.

###  1.2 About this work

##### Biological motivation.

Proliferation and differentiation of hematopoietic stem cells give all
lineages of blood cells, erythrocytes, platelets and several lineages of white
blood cells. Each of them begins with progenitors, immature cells which can
differentiate, die by apoptosis and some of them self-renew. The equilibrium
between these three possibilities has a crucial importance for normal
functioning of hematopoiesis. If it is not preserved, various blood diseases
including leukemia can develop.

Let us consider this process in more detail with the example of
erythropoiesis, red blood cell production. This lineage begins with erythroid
progenitors (colony forming units) which can self-renew, differentiate or die
by apoptosis. They form erythroblastic islands, small units with several
dozens of cells around a macrophage [19], [20]. Their self-renewal occurs only
when they are close to macrophage, otherwise they differentiate or die. This
confirms the observation that self-renewal can take place only under a tight
local regulation. Only global regulation is not sufficient. In this case
erythroblastic islands can have unbounded growth or extinction [10].
Differentiation and apoptosis of erythroid progenitors is determined by the
concentration of Fas-ligand [34]. It is a growth factor produced either by the
progenitors themselves (mice) [40] or by more mature cells (humans) [11].

Successive differentiation of erythroid progenitors gives proerythroblast,
erythroblasts and finally reticulocytes which leave the bone marrow into blood
flow where they become mature erythrocytes. In normal case, these cells can
only differentiate or die by apoptosis, but not self-renew. The latter can
occur due to genetic mutations resulting in erythroleukemia [6]. The choice
between differentiation and apoptosis of these cells depends on the level of
the hormone erythropoietin produced in the kidney. The rate of its production
depends on the relation between the quantity of oxygen required by the
organism and the global quantity of hemoglobin in blood. In the first
approximation we can relate the quantity of hemoglobin to the number of
erythrocytes. Erythrocytes have a given life span in blood, after which they
die. It is about 40 days for mice and 120 days humans.

Thus, simplifying this description, we can summarize it as follows.
Erythroblasts can either differentiate and give mature erythrocytes or die by
apoptosis. The rates of differentiation and apoptosis depend on the number of
erythrocytes in blood. We do not consider here erythroid progenitors because
they can also self-renew, and we omit intermediate staged of erythroblast
maturation. Similar processes occur in other cell lineages of hematopoiesis
and in other tissues. We will use this schematic description as biological
motivation for the model considered in this work.

##### Model and objectives.

The purpose of this work is to suggest a relatively simple model of cell
population dynamics which would allow us to study how intracellular and
extracellular regulations and stochasticity influence cell fate. In order to
simplify the model, we will deliberately not take into account here spatial
heterogeneity and local cell-cell interactions. As discussed above, this is
appropriate for some stages of hematopoiesis and for other tissues.

Hence we will consider only two cell types. Cells of the first type can
differentiate or die by apoptosis. When they differentiate, they produce cells
of the second type whose number determines the rates of differentiation and
apoptosis of cells of the first type. Such models are often considered in
modelling of hematopoiesis (Section 1.1) with the only difference that
differentiation and apoptosis will be determined here by the intracellular
regulation and will not be considered as given constants or functions. This
difference is crucial. Hybrid models of hematopoiesis presented above take
into account intracellular regulation. However they are developed to study
functioning of erythroblastic islands with self-renewal of erythroid
progenitors, local cell interactions, cell motion and other factors. We need a
simpler model to study cell fate in a more pure and general setting. We will
present it in more detail in the next section.

Thus, we will consider in this work a hybrid multi-scale model with
intracellular regulation and global extracellular regulation. We will describe
the intracellular regulation with systems of ordinary differential equations
for the concentrations of intracellular proteins. Therefore we need to specify
the initial values of these concentrations for each cell. If these initial
concentrations are the same for all cells then, under constant external
conditions (which is the case for a steady state solution), the fate of all
cells will also be the same. Hence in order to describe cell choice between
differentiation and apoptosis, we need to introduce some random variations in
the initial conditions. This is one of the important features of the model. We
will see below that in some cases randomness of initial conditions influences
the behavior of the system. It can suppress oscillations and stabilize the
system.

Let us note that random perturbations can occur not only for the initial
conditions but also during the whole cell life. In this case, we need to
consider stochastic differential equations for the intracellular regulation
(see [22] and the references therein). We will restrict ourselves here to the
case of random initial conditions for two reasons. First of all, to simplify
the model and to separate the effects of random initial conditions from random
perturbations during whole cell evolution. Second, in bistable dynamics, when
the trajectory is already inside the basin of attraction of one of stable
stationary points, small random perturbations are not so important. They are
more essential in the beginning of this evolution when the trajectory is close
to the manifold separating basins of attraction. In this case, it is similar
to the introduction of random initial conditions.

Another key feature of the model will be the presence of a global regulation
of the number of differentiated cells. If the number of cells in the
population at time $t$ is $N(t)$ and the target number which should be
achieved by the system is $N_{0}$, then the system will adapt its parameters
depending on the difference $N(t)-N_{0}$. We will show that under some
conditions behavior of the system is determined by the global regulation and
its dependance on the intracellular regulation and on stochasticity is weak.
In this case we can model cellular systems with an uncomplete information
about intracellular regulation and about mechanisms and parameters of possible
stochasticity.

Thus, the main goal of this work is to show how intracellular regulation,
global extracellular regulation and random initial conditions determine
evolution of cell populations. We can summarize their action as follows.
Intracellular regulation provides a choice of cell fate, extracellular
regulation controls this choice, random initial conditions stabilize the
system. We will introduce the model in the next section. Numerical simulations
will be presented in Section 3 and an approximate analytical model in Section
4.

##  2 Model

We will consider cells of two types, $A$ and $B$. Every given time interval
$\tau$ a new cell $A$ appears due to differentiation of less mature cells into
cells $A$. Cells of the type $A$ can either differentiate into cells of the
type $B$ or, otherwise, they die by apoptosis. Hence we restrict ourselves
here to the cell choice between differentiation and apoptosis. We do not
consider here self-renewal of cells $A$ in order to study first the
functioning of the system in this simpler case. Self-renewal takes place for
more immature cells which provide a constant influx of cells $A$. This model
is shown schematically in Figure 1.

Figure 1: Schematic representation of the model. Cells of the type $A$ can
differentiate into cells $B$ or die by apoptosis. Cells $B$ have a fixed life
span after which they die. The number of differentiated cells $B$ influence
the choice of cells $A$ between differentiation and apoptosis. This choice is
determined by intracellular proteins described by a system of ordinary
differential equations.

Cells of the type $B$ have a finite and fixed life span $T$. After that they
die. The total number $N(t)$ of cells $B$ at time $t$ is determined by the
number of cells $A$ differentiated during the time interval $[t-T,t]$. Suppose
that the system should possess $N_{0}$ cells $B$ for the normal functioning of
the organism. Then the rate of differentiation of cells $A$ should depend on
$N(t)$ and $N_{0}$.

In general, the rates of differentiation and apoptosis cannot be considered as
given constants or functions. Indeed, cell response on the extracellular
regulation can depend on its intracellular regulation. Therefore we need to
introduce this regulation in the model. We will study the model with
intracellular regulation in Section 3. When we have this more complete model,
it can be possible to formulate some simplifying assumptions under which it
can be reduced to delay differential equations (Section 4).

##### Intracellular regulation.

We consider here the simplest regulation with only two proteins $p_{1}$ and
$p_{2}$ and, respectively, two equations for their concentrations:

$\frac{dp_{1}}{dt}=F_{1}(p_{1},p_{2}),\;\;\;\;\frac{dp_{2}}{dt}=F_{2}(p_{1},p_{2})\;.$
(2.1)

Each cell has its own values $p_{1}(t)$ and $p_{2}(t)$ described by this
system of equations. Obviously, these concentrations can be different in
different cells.

When a new cell $A$ appears, we prescribe it some initial values $p_{1}^{0}$
and $p_{2}^{0}$. Hence we can determine the evolution of these concentrations
inside each cell. We suppose that if the concentration $p_{1}(t)$ reaches some
critical level $p_{1}^{*}$, then the cell dies by apoptosis. If the
concentration $p_{2}(t)$ reaches some critical level $p_{2}^{*}$, then the
cell differentiates.

Figure 2: Trajectories of system (2.1).

Evolution of the intracellular concentrations depends on the functions $F_{1}$
and $F_{2}$. We will consider them in the following form:

$F_{1}(p_{1},p_{2})=k_{1}p_{1}(1-a_{11}p_{1}-a_{12}p_{2}),\;\;\;\;F_{2}(p_{1},p_{2})=k_{2}p_{2}(1-a_{21}p_{1}-a_{22}p_{2}),$
(2.2)

where $k_{1}$, $k_{2}$ and $a_{ij}$ are positive parameters. We consider the
functions $F_{1}$ and $F_{2}$ in this form in order to describe bistable
dynamics. They can be different depending on applications but their specific
form is not essential for what follows. This system can have up to four
stationary points with non-negative coordinates: $P_{0}=(0,0)$,
$P_{1}=(1/a_{11},0)$, $P_{2}=(0,1/a_{22})$,
$P_{3}=(\tilde{p}_{1},\tilde{p}_{2})>0$, where $\tilde{p}_{1}$ and
$\tilde{p}_{2}$ satisfy the system of equations

$a_{11}p_{1}+a_{12}p_{2}=1,\;\;\;\;a_{21}p_{1}+a_{22}p_{2}=1.$

We suppose that it has a positive solution, and the points $P_{1}$ and $P_{2}$
are stable, while the point $P_{3}$ is unstable. These conditions are
satisfied if $a_{21}>a_{11}>0$, $a_{12}>a_{22}>0$. The separatrix of the point
$P_{3}$ separates the basins of attraction of the points $P_{1}$ and $P_{2}$.
If the initial condition is below this line, then the trajectory converges to
the point $P_{1}$, if the initial condition is above this line, then the
trajectory converges to the point $P_{2}$. If the trajectory approaches the
point $P_{1}$, then the concentration $p_{1}$ becomes large, the concentration
$p_{2}$ small. In this case the cell will die by apoptosis. In the second
case, $p_{2}$ becomes large, $p_{1}$ small, the cell differentiates. Figure 2
shows behavior of trajectories of this system for random initial conditions
chosen in $[0.1,1]\times[0.1,1]$.

##### Cell number and extracellular regulation.

The number of cells $B$ at time $t$ is denoted by $N(t)$. It equals the number
of cells $A$ differentiated during the time interval from $t-T$ to $t$.
Suppose that the system should produce $N_{0}$ cells. If $N(t)$ is different
from $N_{0}$, then there is a feedback control which acts on the rate of
differentiation. Since the choice between differentiation and apoptosis of
cells $A$ is determined by the concentrations $p_{1}$ and $p_{2}$ described by
system (2.1), we will suppose that its coefficients depend on the difference
$N(t)-N_{0}$. Namely, we set

$a_{21}(N)=a_{21}^{0}+\alpha(N-N_{0}).$ (2.3)

In order to reduce the number of parameters, we consider a linear dependence
of the coefficient on $N$. We verify during the simulations that it remains
positive. All other coefficients are independent of $N$.

If $N(t)>N_{0}$, then the coefficient $a_{21}$ increases in comparison with
$a_{21}^{0}$. Therefore the function $F_{2}(p_{1},p_{2})$ in (2.2) decreases.
Hence growth of the concentration $p_{2}$ decelerates and, consequently, the
rate of differentiation also decreases. Thus the dependence of the
coefficients of the intracellular regulation on the cell number provides a
global control over the system with the purpose to get the required number of
differentiated cells.

##### Stochasticity in the initial conditions.

Let us recall that there is a constant influx of cells $A$. This means that
every time interval $\tau$ a new cell enters the compartment with cells $A$.
We need to prescribe the initial concentrations $p_{1}^{0}$ and $p_{2}^{0}$
for each new cell. Since these values are usually not known and they can
hardly be measured experimentally, this becomes an additional unknown
parameter which can influence behavior of the system. We will choose the
initial conditions from the square domains

$D=\\{a-k\leq p_{1}^{0}\leq a+k,\;\;\;b-k\leq p_{2}^{0}\leq b+k\\},$ (2.4)

where $a$, $b$ and $k$ are chosen in such a way that the whole domain $D$ has
non-negative coordinates. The initial condition will be chosen from the domain
$D$ as a random variable with a uniform distribution.

##  3 Numerical results

We present results of numerical simulations of the model described in Section
2. We set $\tau=0.1$, $T=20$, $a_{11}=a_{22}=0.5$, $a_{12}=a_{21}^{0}=1$,
$p_{1}^{*}=p_{2}^{*}=1.8$ and we will vary other parameters.

###  3.1 Deterministic initial conditions

We begin with the case where the initial conditions for the concentrations
$p_{1}^{0}$, $p_{2}^{0}$ are fixed. Therefore all cells $A$ have exactly the
same initial conditions in the intracellular regulation. Then this problem
does not have solutions with a constant value $N(t)\equiv N$. In this sense we
can say that the problem with deterministic initial conditions does not have
stationary solution. Indeed, suppose that such solution exists and consider
the phase plane of the system (2.1). Since $a_{21}(N)$ is constant (time
independent) then the separatrix which separates basins of attraction of the
points $P_{1}$ and $P_{2}$ is also fixed. If the initial condition
$(p_{1}^{0},p_{2}^{0})$ belongs to the basin of attraction of the point
$P_{1}$, then all cells will die. If it belongs to the basin of attraction of
the point $P_{2}$, then all cells will differentiate. Finally, if the initial
condition is exactly at the sepratrix, the corresponding trajectory will
converge to the unstable point $P_{3}$ and it will never reach the critical
values $p_{1}^{*}$ or $p_{2}^{*}$. Hence in all three cases the number of
differentiated cells cannot be equal $N$ if it is different from $0$ or
$T/\tau$.

Figure 3: The number of differentiated cells $N(t)$ for a constant initial
condition $p_{1}^{0}=p_{2}^{0}=0.1$. $N_{0}=150,\alpha=0.003$ (left),
$N_{0}=150,\alpha=0.009$ (middle), $N_{0}=200,\alpha=0.003$ (right).

Figure 3 shows numerical simulations in the case of a fixed (deterministic)
initial condition. The number of differentiated cells oscillates around the
target value $N_{0}$. Therefore the separatrix oscillates and the initial
condition $(p_{1}^{0},p_{2}^{0})$ belongs to different basins of attraction
during different time intervals. These oscillations do not result from the
instability of a stationary solution since the stationary solution does not
exist.

###  3.2 Random initial conditions

####  3.2.1 Existence of stationary solutions

We indicated in the previous section that the problem with fixed deterministic
initial conditions does not have stationary solutions. We will consider now
random initial conditions.

We consider system (2.1) with the functions $F_{1}$ and $F_{2}$ given by
equalities (2.2), (2.3). The initial conditions are taken in domain $D$
defined in (2.4) with a uniform distribution. For each given $N$, denote by
$R_{N}(t)$ the number of differentiated cells produced by the system during
the time interval $[t-T,t]$. Since the initial conditions are random, this
number can also have random oscillations. Let $\widehat{R}_{N}$ be its average
value with respect to asymptotically large time interval,

$\widehat{R}_{N}=\lim_{S\to\infty}\frac{1}{S}\int_{0}^{S}R_{N}(t)dt.$

We will consider that the problem has a stationary solution if the following
equality holds:

$\widehat{R}_{N}=N.$ (3.1)

We will show that this equation has an approximate solution.

Consider system (2.1) for this value of $N$. Separatrix $S$ of the point
$P_{3}$, that is the trajectory which approach this saddle point, separates
basins of attraction of the points $P_{1}$ and $P_{2}$. Denote by $D_{i}(N)$
the part of the domain $D$ which belongs to the basin of attraction of the
point $P_{i}$, $i=1,2$. Let $|D_{i}(N)|$ be measures of these subdomains. Then
$\gamma_{N}=|D_{2}(N)|/|D|$ is the proportion of differentiated cells for
asymptotically large time, and $\gamma_{N}T/\tau$ is their average number
during time interval $T$. Equation

$\gamma_{N}T/\tau=N$ (3.2)

can be considered as approximation of equation (3.1). We will indicate
conditions when this equation has a solution. Let us note first of all that
$\gamma_{N}$ is a decreasing function of $N$. The value of $N$ changes from
$0$ to $T/\tau$. If system (2.1) is well defined for all these values of $N$
and possesses the same structure (number and stability of stationary points),
then equation (3.2) has a unique solution.

####  3.2.2 Stable stationary solutions

In this section the initial conditions $p_{1}^{0}$ and $p_{2}^{0}$ will be
taken from the square domain $0\leq p_{1}^{0},p_{2}^{0}\leq k$ as a random
variable with a uniform distribution. During the time interval $T$ there are
$T/\tau$ new cells $A$. In the case of symmetric parameters and initial
conditions, a half of them will differentiate. For the values of parameters
considered here, we obtain $N=100$ differentiated cells $B$. Figure 4 shows
examples of numerical simulations for $N_{0}=100$ and different values of
$\alpha$. If $\alpha$ is sufficiently large, then there is an overshoot in
cell number. Its value for large time converges to $N_{0}$ after small
decaying oscillations.

Figure 4: The number of cells $N(t)$ as a function of time for $N_{0}=100$,
$\alpha=0$ (left) and $\alpha=0.009$ (middle). Four curves show different
simulations with the same values of parameters. The curves in the right figure
correspond to different values of $k$. Average value of $N(t)$ for 50
simulations, $N_{0}=100$, $\alpha=0.009$.

Figure 5: The number of differentiated cells $N(t)$ as a function of time for
$N_{0}=150$, $\alpha=0.009$ (left), $N_{0}=200$, $\alpha=0.003$ (middle) and
$\alpha=0.009$ (right) and different values of $k$. Average for 50
simulations.

Figures 4 and 5 show the evolution of the number of differentiated cells
$N(t)$ in time for different values of parameters. In all cases it converges
to some limiting value $N_{\infty}$. It is interesting to note that the
function $N(t)$ is practically independent of the value of $k$ which
determines the choice of initial conditions. In the wide range of variation of
this parameter, $0.1\leq k\leq 1.5$, the curves coincide. A slight difference
appears for $k=2$.

Let us recall that in the case of deterministic initial conditions (Section
3.1), $N(t)$ does not converge to a stationary value. Hence we can expect
appearance of oscillations for very small values of $k$. We will return to
this question in the next section. Therefore convergence to a stationary
solution is determined by randomness in the initial conditions. We illstrate
this mechanism in Figure 6. The domain $D$ of initial conditions is split by
the separatrix into two subdomains which belong to the basins of attraction of
the points $P_{1}$ and $P_{2}$. Depending on the position of the initial
condition, the corresponding trajectory will converge to one of these
stationary points. The areas of these two subdomains determine the proportion
of differentiated cells and the value $N_{\infty}$. Let us note that system of
equations (2.1)-(2.3) has time dependent coefficients. However for $t$
sufficiently large, $N(t)$ is close to a constant, and the curves in Figure 6
can be considered as trajectories of the autonomous system.

The global regulation by the number of differentiated cells $N(t)$ through the
coefficients of the system (see (2.3)) influences the limiting value
$N_{\infty}$. It increases with the increase of $N_{0}$ (Figures 5, 6) but in
general it is different from $N_{0}$. Clearly, if $N(t)\equiv N_{0}$, then
from (2.3) we have $a_{21}=a_{21}^{0}$. In this case, the system is symmetric
and exactly half of the total cell number differentiate. For the values of
parameters under consideration $N_{\infty}=100$. Therefore, if $N_{0}\neq
100$, then $N_{\infty}<N_{0}$. The value of $N_{\infty}$ increases and
approaches $N_{0}$ when we take greater values of $\alpha$ (Figure 5).

Figure 6: Trajectories of system (2.1) for the same simulations as in Figure 5
and for $95<t<100$. The domain of initial conditions (red) is split by the
separatrix of the saddle point $P_{3}$ in two subdomains. Depending on the
initial condition, trajectories approach either point $P_{1}$ or $P_{2}$;
$N_{0}=150$ (left), $N_{0}=200$ (right).

The simulations presented above are carried out for a fixed value $a_{21}^{0}$
in (2.3). We can also introduce its dependence on $N$. Then we can obtain a
better approximation of $N_{0}$ by $N_{\infty}$. We varied $a_{21}^{0}$ in
large limits and obtained convergence to a stationary value $N_{\infty}$ which
depends on the value of $a_{21}^{0}$ (not shown here).

Thus global regulation allows us to obtain a given stationary value of the
number of cells for a large interval of variation of the coefficients and of
the range of random initial conditions.

Figure 7: The number of differentiated cells $N(t)$ for $N_{0}=150$,
$\alpha=0.009$, $a=0.4,b=0.6$ and different values of $k$. When $k$ increases,
the amplitude of oscillations decay.

####  3.2.3 Oscillating solutions

In Section 3.2.1 we considered deterministic initial conditions and observed
periodic oscillations of the number of differentiated cells. In Section 3.2.2
we studied the case with random initial conditions from a sufficiently large
domain. In this case, the number of differentiated cells converges to a
constant value.

Let us now consider transition between these two cases. Initial conditions
will be uniformly distributed in the domain $[a-k,a+k]\times[b-k,b+k]$. Figure
7 shows the number of differentiated cells in time for different values of $k$
and all other parameters fixed. The amplitude of oscillations decreases when
$k$ increases, and they practically disappear for $k=0.4$. Hence randomness of
initial conditions can suppress oscillations. We will discuss a possible
mechanism of this stabilization in Section 4.

Figure 8: The number of differentiated cells $N(t)$ for $N_{0}=150$ (left),
$N_{0}=200$ (right), $\alpha=0.009$, $a=1,b=0.2$. The oscillations are
independent of the value of $k$ in the considered range.

Figure 9: Curves $p_{1}(t),p_{2}(t)$ on the $(p_{1},p_{2})$-plane for the same
values of parameters as in the previous figure. Initial conditions are shown
as red dots.

If we now take other values of $a$ and $b$ (Figure 8), then the oscillations
are practically independent of the value of $k$ (in the considered range).
Persistence of the oscillations can be related to large time delay for cells
to differentiate. Indeed, the time which the trajectory needs to reach the
critical value $p_{2}^{*}$ depends on the initial condition.

System (2.1) has time dependent coefficients (2.3). Therefore the
extracellular regulation through $N(t)$ acts on the cell during its all life
time. Indeed, some of the curves $(p_{1}(t),p_{2}(t))$ shown in Figure 9 cross
each other because the system is not autonomous. This extracellular regulation
can have more influence on the intracellular concentrations than the choice of
initial conditions. This is why the results shown in Figure 8 are practically
independent of the range where the initial conditions are chosen from.

Figure 10: The number of differentiated cells $N(t)$ for
$N_{0}=150,\alpha=0.009$, $k=0.01$ and different values of $\epsilon$ in
system (3.3), (3.4).

In order to analyze how time oscillations are related to the time delay of the
intracellular regulation, we introduce a small parameter $\epsilon$ in system
(2.1):

$\epsilon\frac{dp_{1}}{dt}=k_{1}p_{1}(1-a_{11}p_{1}-a_{12}p_{2}),$ (3.3)

$\epsilon\frac{dp_{2}}{dt}=k_{2}p_{2}(1-a_{21}p_{1}-a_{22}p_{2}).$ (3.4)

Decreasing $\epsilon$ we accelerate intracellular regulation. The amplitude of
oscillations also decreases (Figure 10). This confirms the hypothesis that
time oscillations are related to time delay as it is the case in delay
differential equations (Section 4).

##  4 Analytical approximation

The total number $N(t)$ of differentiated cells can be found from the equation

$\frac{dN}{dt}=H(t)-M(t),$ (4.1)

where $H(t)$ and $M(t)$ the rates of cell birth and death at time $t$. Since
cells have a fixed life time $T$, then

$M(t)=H(t-T),$

and equation (4.1) can be written as follows:

$\frac{dN}{dt}=H(t)-H(t-T).$ (4.2)

In the numerical model considered above, instead of the birth rate $H(t)$ we
consider cell choice between differentiation and apoptosis of cells $A$.
Therefore $H(t)$ is the rate of appearance of newly differentiated cells at
time $t$. It is determined by solutions of the system of differential
equations for the intracellular variables. So we cannot express it through $N$
in order to have a closed equation for this variable. However we can derive
such equation if we do some simplifying assumptions.

Let us recall that the choice between differentiation and apoptosis is
determined by the intracellular variables $p_{1}$ and $p_{2}$ described by
system (2.1). Its solution depends on the initial conditions and on the values
of the coefficients which depend on $N(t)$ due to the extracellular
regulation.

Let us fix some $t=t_{0}$. The value $N(t_{0})$ determines the values of the
coefficients of system (2.1). Basins of attraction of stable points $P_{1}$
and $P_{2}$ are separated by the stable manifold of the unstable point
$P_{3}$. Denote the subdomains of $D$ these basins of attractions by $D_{1}$
and $D_{2}$, respectively. If the initial condition is located in $D_{1}$, the
cell will go in the direction of the point $P_{1}$ and will die by apoptosis,
if it is in $D_{2}$, then it will go in the direction of $P_{2}$ and will
differentiate. If we know the ratio of initial conditions in $D_{1}$ and
$D_{2}$, then we can try to estimate the rate of differentiation. However, it
cannot be done directly for the model under consideration because: a) The
coefficients of the system depend on time through $N(t)$. Hence a trajectory,
which is initially in the basin of attraction of one of the two stable points,
can change it and go to another point, b) Even inside the same basin of
attraction and for fixed values of the coefficients, the time to reach the
critical values $p_{1}^{*}$ and $p_{2}^{*}$ depends on the initial condition.
This time increases when the initial condition approaches the stable manifold
of the point $P_{3}$. Hence the cell appeared at time $t=t_{0}$ can
differentiate or die at any other time $t>t_{0}$. c) The initial condition is
chosen randomly from some given domain $D$.

Therefore we cannot reduce the model with intracellular regulation described
by ordinary differential equations to a closed delay equation (4.2). Consider
now a simplified model where a cell appeared at time $t_{0}$ makes its choice
between differentiation and apoptosis exactly at the moment $t=t_{0}$ of its
appearance. This choice is not influenced by further change of $N(t)$. Then
for each cell we consider system (2.1) with constant coefficients. These
coefficients can be different for cells appeared at different moments of time.
In this case, the probability for the initial condition to be in domains of
attraction of each of two stable stationary points at time $t_{0}$ are given
by the formula

$\pi_{1}(t_{0})=|D_{1}|/|D|,\;\;\;\pi_{2}(t_{0})=|D_{2}|/|D|,$

where $|D|$ is the area of the domain $D$, $|D_{i}|$ is the area of the
corresponding subdomain. Since the coefficients of system (2.1) depend on
time, then the subdomains $D_{i}$ are also time dependent,
$D_{i}=D_{i}(t_{0})$.

Next, let us assume that the time to reach the critical values $p_{1}^{*}$ and
$p_{2}^{*}$, when the cell differentiates or dies, does not depend on the
initial condition. This means that instead of the ordinary differential system
of equations we consider a binary (or Boolean) mode with some given time delay
$\tau_{0}$. If the initial condition belongs to $D_{1}(t_{0})$, then the cell
will differentiate at time $t=t_{0}+\tau_{0}$, otherwise it will die. Then the
rates of differentiation and apoptosis are given by the relations

$d(t)=\gamma\pi_{1}(t-\tau_{0}),\;\;\;a(t)=\gamma\pi_{2}(t-\tau_{0}),$

where $\gamma$ is the rate of appearance of new cells. It remains to note that
$D_{i}(t)$ are some given functions of $N(t)$. Hence

$d(t)=\gamma f(N(t-\tau_{0})),\;\;\;a(t)=\gamma(1-f(N(t-\tau_{0})).$

Let us note that $f(N)$ is a decreasing function: if there are more
differentiated cells, then the rate of their production is less.

Thus, equation (4.2) can be written in the form

$\frac{1}{\gamma}\;\frac{dN}{dt}=f(N(t-\tau))-f(N(t-T-\tau)).$ (4.3)

This equation is obtained as an approximation of equation (4.2) under some
simplifying assumptions. It can also be considered independently, as a simple
phenomenological model of cell population under extracellular regulation. The
first term in the right-hand side of this equation determines the influx of
cells in the tissue and the second term their death. Equation (4.3) is
equivalent to the integral equation

$N(t)=\gamma\int_{t-T}^{t}f(N(s-\tau))ds+n_{0},$

where $n_{0}$ is an arbitrary constant. We note that $n_{0}=0$ since $N(t)$
equals exactly the number of cells differentiated from $t-T$ to $t$ given by
the integral in the right-hand side. The last equation can also be written in
the form:

$N(t)=\gamma\int_{0}^{T}f(N(t-y-\tau))dy.$

If we look for a constant solution $N(t)=N$, then we obtain the equation

$N=\gamma Tf(N).$

Since $f(N)$ is a positive decreasing function, then it has a unique solution
$N=N^{*}$ (cf. Section 3.2.1).

##### Linear stability analysis.

Let us linearize equation (4.3) about the constant solution $N^{*}$:

$\frac{1}{\gamma}\;\frac{du}{dt}=f^{\prime}(N^{*})(u(t-\tau))-u(t-T-\tau)).$
(4.4)

We look for its solution as $u(t)=\exp(\lambda t)$. Then we get

$\frac{1}{h}\;\lambda=e^{-\lambda\tau}-e^{-\lambda(T+\tau)},$ (4.5)

where $h=\gamma f^{\prime}(N^{*})$. We note that $\lambda=0$ is a solution of
this equation. It is related to the fact that equation (4.3) has a family of
solutions. We will look for the solution $\lambda=i\phi$. Then we get the
system of two equations

$\cos(\phi\tau)=\cos(\phi(T+\tau)),$

$-\frac{1}{h}\phi=-\sin(\phi\tau)+\sin(\phi(T+\tau)).$

From this system of equations we can find $\phi$ and the stability boundary.
From the first equation we get

$\phi(T+\tau)=2\pi n-\phi\tau$

(sinus of these two values should have opposite signs). Then from the second
equation we have

$\sin x=\frac{x}{2h\tau},$ (4.6)

where

$x=\frac{2\pi n\tau}{T+2\tau}\;.$

Equation (4.6) determines the stability boundary. It has a nonzero solution if
$\tau>1/(2h)$.

##### Implications for the original problem.

In the analytical model we obtain that oscillations appear if the time delay
$\tau$ of the intracellular regulation is sufficiently large. This conclusion
is in agreement with the result of the numerical simulations presented in
Section 3. Indeed, let the initial conditions be uniformly distributed in the
square $[a-k,a+k]\times[b-k,b+k]$. If the point $(a,b)$ is more close to the
point $P_{2}$, the time delay to reach the critical value $p_{1}^{*}$ is
greater than in the case where it is close to $P_{1}$. Accordingly, if
$a=1,b=0.2$ the oscillations are strong (Figure 8), for $a=0.4,b=0.6$ they are
weaker (Figure 7), for $a=0.1,b=1$ no oscillations (not shown).

We can also use this approximate analytical result to understand the influence
of stochasticity. When $k$ is small and $a=0.6,b=0.6$, there are oscillations.
When the initial conditions are taken from the square $[0,k]\times[0,k]$,
$k=1,1.5$, there are no oscillations. When we take a big square, then there
are different initial conditions. For some of them time delay is small (as in
the case where $(a,b)$ is close to $P_{1}$). Such solutions with small time
delay can suppress oscillations (see the next paragraph).

##### Two delays in the intracellular regulation.

In order to understand the role of variation of initial conditions in the
intracellular regulation, we can consider two production functions $f$ and $g$
with different delays:

$\frac{1}{\gamma}\;\frac{dN}{dt}=f(N(t-\tau_{1}))-f(N(t-T-\tau_{1}))+g(N(t-\tau_{2}))-g(N(t-T-\tau_{2})).$
(4.7)

This equation corresponds to the case of two different initial conditions with
two different time delays. We can expect that if the first delay is large
enough to produce oscillations considered along, and the second delay is
short, then it can remove the oscillations. Let us verify this conjecture.
After linearization we obtain the equation for the eigenvalues:

$\lambda=h_{1}\left(e^{-\lambda\tau_{1}}-e^{-\lambda(T+\tau_{1})}\right)+h_{2}\left(e^{-\lambda\tau_{2}}-e^{-\lambda(T+\tau_{2})}\right),$
(4.8)

where $h_{1}=\gamma f^{\prime}(N^{*})$, $h_{2}=\gamma g^{\prime}(N^{*})$,
$N^{*}$ is a solution of the equation

$N=\gamma T[f(N)+g(N)].$

We set $\tau_{1}=\tau$ as in equation (4.5) and substitute in the right-hand
side of equality (4.8) the solution $\lambda=i\phi$ of equation (4.5). Then we
get

$\lambda=\frac{h_{1}}{h}\;i\phi+h_{2}\;\left(e^{-i\phi\;\tau_{2}}-e^{-i\phi\;(T+\tau_{2})}\right)\;.$
(4.9)

Denote by $R$ the real part of the expression in the right-hand side of this
equality. Then

$R=h_{2}(\cos(\phi\tau_{2})-\cos(\phi(T+\tau_{2}))).$

If $\tau_{2}=0$, then $R\leq 0$ since $h_{2}<0$. Moreover, $R=0$ only if $\phi
T=2\pi n$ for some integer $n$. If we do not take into consideration these
exceptional cases, then $R<0$ for $\tau_{2}=0$ and, by continuity, $R$ remains
negative for sufficiently small $\tau_{2}$.

Thus, if we take the same values of parameters which determine the stability
boundary in the case of one time delay and introduce the second (short) time
delay, then we get in the stability region. Hence the second time delay
stabilizes the stationary solution. This model with two delays confirms the
results of the numerical simulations which show that random initial conditions
can suppress oscillations. Among random initial conditions, those which
correspond to large time delay, lead to the oscillations, while the initial
conditions, which correspond to small time delay, can suppress the
oscillations.

##  5 Discussion

The goal of this work is to study how intracellular regulation, extracellular
regulation and stochasticity in the initial intracellular concentrations
determine cell fate. The results of this work allow us to suggest that they
act in the following way. Intracellular regulation provides a choice between
two or several options. Global extracellular regulation controls the
realization of this choice. Stochasticity in the initial conditions stabilizes
the system.

Each of these three elements should satisfy certain conditions. Intracellular
regulation should possess some kind of bistability. If we describe it with a
system of ordinary differential equations, this can be two stable stationary
points which correspond to different cell fates. However bistabilty can be
understood here in a larger sense. The phase space should be split into two
invariant manifolds by a central manifold. The invariant manifolds correspond
to different cell fates. The central manifold should be controlled by the
global extracellular regulation through the coefficients of the ordinary
differential system of equations in order to adapt the volume of invariant
manifolds to get the target value of cell number.

Intracellular regulation is extremely complex and in many cases there is only
a partial information about it. An important methodological question is how to
take it into account if only a part of this network is known and the reaction
constants are not known. From the point of view of the modelling presented in
this work, we should take the known part of the intracellular regulation and
complete it to the system which described the required cell fate. This
additional part, since it is not known, can be constructed in different ways.
However, this difference may not be essential from the point of view of the
cell population because the global extracellular regulation will control it.
This control will be different for different intracellular regulations but it
will finally produce the same result.

We will finish this discussion with the last element of our regulatory system,
stochasticity in the initial conditions. We showed in numerical simulations
and in analytical models that it can stabilize the system. Oscillations is an
intrinsic property of cellular systems because of time delay in their
regulation. So stochasticity can be one of possible mechanisms to suppress
these oscillations. In this case an interesting question is about the origin
of stochasticity. It can result from random perturbations, from stochastic
dynamics in the case of small number of molecules or from some underlying
mechanism which can produce some variations in initial protein concentrations.
The question about possible mechanisms producing variation of initial
condition is completely open.

## References

  * [1] A.S Ackleh, K. Deng, K. Ito, J. Thibodeaux. A structured erythropoiesis model with nonlinear cell maturation velocity and hormone decay rate. Math. Bios., 204, (2006), 21–48.
  * [2] A. R. A. Anderson, M. Chaplain, K.A. Rejniak. Single cell based models in biology and medicine. Birkhäuser, Basel, 2007.
  * [3] R. Apostu, M.C. Mackey. Understanding cyclical thrombocytopenia: A mathematical modeling approach. J. Theor. Biol., 251 (2008), 297–316.
  * [4] S. Bernard. Modélisation multi-échelles en biologie. In: Le vivant discret et continu. N. Glade, A. Stephanou, Editeurs, Editions Materiologiques, 2013, 65-89.
  * [5] S. Bernard , J. Bélair, M.C. Mackey. Oscillations in cyclical neutropenia: New evidence based on mathematical modeling. J. Theor. Biol., 223, (2003), 283–298.
  * [6] N. Bessonov, F. Crauste, I. Demin, V. Volpert. Dynamics of erythroid progenitors and erythroleukemia. Math. Model. Nat. Phenom., 4 (2009), no. 3, 210-232.
  * [7] N. Bessonov, P. Kurbatova, V. Volpert. Particle dynamics modelling of cell populations. Proc. Conf. JANO, Mohhamadia, 2008. Math. Model. Nat. Phenom., 5 (2010), no. 7, 42–47.
  * [8] N. Bessonov, P. Kurbatova, V. Volpert. Pattern Formation in Hybrid Models of Cell Populations. In: Pattern formation in morphogenesis. V. Capasso, M. Gromov, A. Harel-Bellan, N. Morozova, L. Pritchard. Editors, Springer, 2012, 107-119.
  * [9] N. Bessonov , F. Crauste, S. Fischer, P. Kurbatova, V. Volpert. Application of Hybrid Models to Blood Cell Production in the Bone Marrow. Math. Model. Nat. Phenom., 6 (2011), no. 7, 2–12.
  * [10] N. Bessonov, N. Eymard, P. Kurbatova, V. Volpert. Mathematical modeling of erythropoiesis in vivo with multiple erythroblastic islands. Applied Mathematics Letters, 25 (2012), 1217–1221.
  * [11] J.A. Chasis, N. Mohandas. Erythroblastic islands: niches for erythropoiesis. Blood, 112 (2008).
  * [12] C. Colijn, M.C. Mackey. A mathematical model of hematopoiesis – I. Periodic chronic myelogenous leukemia. J. Theor. Biol., 237, (2005), 117–132.
  * [13] C. Colijn, M.C. Mackey. A mathematical model of hematopoiesis – II. Cyclical neutropenia. J. Theor. Biol., 237 (2005), 133–146.
  * [14] F. Crauste, I. Demin, O. Gandrillon, V. Volpert. Mathematical study of feedback control roles and relevance in stress erythropoiesis. J. Theor. Biology, 263 (2010), 303–316.
  * [15] F. Crauste, L. Pujo-Menjouet, S. Génieys, C. Molina, O. Gandrillon. Adding Self-Renewal in Committed Erythroid Progenitors Improves the Biological Relevance of a Mathematical Model of Erythropoiesis, J. Theor. Biology, 250 (2008), 322–338.
  * [16] V. Cristini, J. Lowengrub. Multiscale modeling of cancer: an integrated experimental and mathematical modeling approach. Cambridge University Press, Cambridge, 2010.
  * [17] I. Demin, F. Crauste, O. Gandrillon, V. Volpert. A multi-scale model of erythropoiesis. Journal of Biological Dynamics, 4 (2010), 59–70.
  * [18] J. Eller , I. Gyori, M. Zollei, F. Krizsa. Modelling Thrombopoiesis Regulation - I Model description and simulation results. Comput. Math. Appli, 14 (1987), 841–848.
  * [19] N. Eymard, N. Bessonov, O. Gandrillon, M.J. Koury, V. Volpert, The role of spatial organisation of cells in erythropoiesis. J. Math. Biol., (2014), in press.
  * [20] S. Fischer, P. Kurbatova, N. Bessonov, O. Gandrillon, V. Volpert, F. Crauste. Modelling erythroblastic islands : using a hybrid model to assess the function of central macrophage. J. Theor. Biol., 298 (2012), 92–106.
  * [21] N. Glade, A. Stephanou, Editeurs. Le vivant discret et continu. Editions Materiologiques, 2013.
  * [22] A. Golubev. Random discrete competing events vs. dynamic bistable switches in cell proliferation in differentiation. J. Theor. Biol., 267 (2010), no. 3, 341–354.
  * [23] A. Halanay. Periodic solutions in a mathematical model for the treatment of chronic myelogenous leukemia. Math. Model. Nat. Phenom., 7 (2012), no. 1, 235–244.
  * [24] A. Halanay, D. Candea, I. R. Radulescu. Existence and stability of limit cycles in a two-delays model of hematopoiesis including asymmetric division. Math. Model. Nat. Phenom., 9 (2014), no. 1, 58–78.
  * [25] M. J. Koury, M.C. Bondurant. Erythropoietin retards DNA breakdown and prevents programmed death in erythroid progenitor cells. Science, 248 (1990), 378–381.
  * [26] P. Kurbatova, S. Bernard, N. Bessonov, F. Crauste, I. Demin, C. Dumontet, S. Fischer, V. Volpert. Hybrid Model of Erythropoiesis and Leukemia Treatment with Cytosine Arabinoside. 2011, SIAM J. Appl. Math., 71 (2011), no. 6, 2246–2268.
  * [27] P. Kurbatova, N. Eymard, V. Volpert. Hybrid model of erythropoiesis. Acta Biotheoretica, 61 (2013), no. 3, 305-315.
  * [28] J. Lei, S.A. Levin, Q. Nie. Mathematical model of adult stem cell regeneration with cross-talk between genetic and epigenetic regulation. PNAS, 2014, E880-–E887.
  * [29] B.D. MacArthur.  Collective dynamics of stem cell populations. PNAS, 111 (2014), no. 10, 3653–-3654.
  * [30] B.D. MacArthur, I.R. Lemischka. Statistical mechanics of pluripotency. Cell 154 (20130, 484–489.
  * [31] M.C. Mackey. Unified hypothesis of the origin of aplastic anaemia and periodic hematopoïesis. Blood 51, (1978), 941–956.
  * [32] M.C Mackey. Dynamic hematological disorders of stem cell origin. In: G. Vassileva-Popova and E. V. Jensen, Editors. Biophysical and Biochemical Information Transfer in Recognition, Plenum Press, New York, (1979), 373–409.
  * [33] J.M. Mahaffy, J. Belair, M.C. Mackey. Hematopoietic model with moving boundary condition and state dependent delay: applications in erythropoiesis. J. Theor. Biol., 190 (1998) 135–146.
  * [34] R. De Maria, U. Testa, L. Luchetti, A. Zeuner, G. Stassi, E. Pelosi, R. Riccioni, N. Felli, P. Samoggia, C. Peschle. Apoptotic Role of Fas/Fas Ligand System in the Regulation of Erythropoiesis. Blood, 93 (1999), 796–803.
  * [35] R. Morris, I. Sancho-Martinez, T.O. Sharpee, J.C. Izpisua Belmonte. Mathematical approaches to modeling development and reprogramming. PNAS, 111 (2014), no. 14, 5076–-5082.
  * [36] J.M. Osborne, A. Walter, S.K. Kershaw, G.R. Mirams, A.G. Fletcher, P. Pathmanathan, D. Gavaghan, O.E. Jensen, P.K. Maini, H.M. Byrne. A hybrid approach to multi-scale modelling of cancer. Phil. Trans. R. Soc. A, 368 (2010), 5013–5028.
  * [37] H. Ozbay, C. Bonnet, H. Benjelloun, J. Clairambault. Stability analysis of cell dynamics in leukemia. Math. Model. Nat. Phenom., 7 (2012), no. 1, 203–234.
  * [38] I. Roeder. Quantitative stem cell biology: computational studies in the hematopoietic system. Curr. Opin. Hematol., 13 (2006), 222–228.
  * [39] M. Santillan, J.M. Mahaffy, J. Belair, M.C. Mackey. Regulation of platelet production: The normal response to perturbation and cyclical platelet disease. J. Theor. Biol., 206 (2000), 585–603.
  * [40] M. Socolovsky. Molecular insights into stress erythropoiesis. Current opinion in hematology, 14 (2007).
  * [41] T. Stiehl, A. Marciniak-Czochra. Mathematical modeling of leukemogenesis and cancer stem cell dynamics. Math. Model. Nat. Phenom., 7 (2012), no. 1, 166–202.
  * [42] V. Volpert. Elliptic partial differential equations. Volume 2. Reaction-diffusion equations. Birkhäuser, 2014.
  * [43] V. Volpert, N. Bessonov, N. Eymard, A. Tosenberger. Modèle multi-échelle de la dynamique cellulaire. In: Le vivant discret et continu. N. Glade, A. Stephanou, Editeurs, Editions Materiologiques, 2013, 91-111.
  * [44] H.E. Wichmann, M.D. Gerhardts, H. SPechtmeyer, R. Gross. A mathematical model of thrombopoiesis in rats. Cell Tissue Kinet., 12 (1979), 551–567.
  * [45] H.E. Wichmann, M. Loeffler. Mathematical Modeling of Cell Proliferation. Boca Raton, FL, CRC, (1985).
  * [46] H. Wulff, H.E. Wichmann, M. Loeffler, K. Pantel. A mathematical model of erythropoiesis in mice and rats. Part 3. Suppressed erythropoiesis. Cell Tissue Kinet., 22 (1989), 51–61.
  * [47] H.E. Wichmann, M. Loeffler, K. Pantel, H. Wulff. A mathematical model of erythropoiesis in mice and rats. Part 2. Stimulated erythropoiesis. Cell Tissue Kinet., 22 (1989), 31–49.


In this section we prove that the above constructed asymptotically flat spacetime $(\mathcal{M}^{(4+1)},g^{(4+1)})$ with a locally naked singularity is unstable in the $SO(3)\times U(1)$ symmetry class. We point out that showing the initial data is non-generic in a whole well-posedness class is beyond the purpose of this paper. On the other hand, we illustrate how the blue-shift instability mechanism of [6] applies in our case, by providing a two dimensional set of perturbations of our solution to spacetimes without locally naked singularities.

Given a spherically symmetric solution of the Einstein-Scalar field equations $(M,g,\phi)$ with a first singularity at the center $b_{\Gamma},$ we denote by $C_0^-$ the incoming null cone going into $b_{\Gamma}$ and by $C_0^+$ an outgoing null cone. We write the metric in double null coordinates $(u,v)$ such that $C_0^-=\{v=0\},$ $u=-2r$ on $C_0^-,$ and $v=2r$ on $C_0^+.$ Note that with this normalization we have $b_{\Gamma}=\{u=0,\ v=0\}.$ The blue-shift is defined:

\[\gamma(u)=\int_{u_0}^u\frac{1}{r}\cdot\frac{\mu}{1-\mu}(u',0)du'.\]

As pointed out before, given our normalization it suffices to prescribe $\phi$ on $C_0^+$ to obtain characteristic initial data. Using the two instability theorems of [6], one gets the following result:

Let $(M,g,\phi)$ be a bounded variation spherically symmetric solution of the Einstein-Scalar field equations. Assume that the spacetime has a locally naked singularity and that the blue-shift is unbounded on $C_0^-$. There exist functions $f_1,\ f_2$ on $C_0^+$ with $\partial_vf_1\in BV,\ f_2\in C^{1,1/2}$ and $f_1=f_2=0$ for $v<0$, with the following property: for any two constants $\lambda_1,\ \lambda_2$ such that at least one is nonzero, any exterior incoming null cone in the future development of the initial data $\phi+\lambda_1f_1+\lambda_2f_2$ reaches a spacelike singularity $\mathcal{B}.$

The instability theorem gives a 2-dimensional space of perturbations of solutions with locally naked singularities and unbounded blue-shift to spacetimes without locally naked singularities, where the boundary $\mathcal{B}$ is preceded by trapped surfaces. We represent the Penrose diagram of such Einstein-Scalar field solutions:

[dashed] (0,0) – (3,3);
[dashed] (3,3) – (2,4);
(0,0) – (0,4);
[very thin] (0,4) – (2,2);
[thick, dashdotted] (0,4) .. controls (1,3.8)  .. (2,4);
[color=black, fill=white] (0,4) circle (2pt) node[anchor=east] $b_{\Gamma}$;
(0,2.5)  node[anchor=east] $\Gamma$;
(1.05,3)  node[anchor=west] $C_0^-$;
(1.5,1.5)  node[anchor=north west] $\mathcal{I}^-$;
(1.3,3.9)  node[anchor=south east] $\mathcal{B}$;
[color=black, fill=white] (3,3) circle (2pt) node[anchor=north west] $i^0$;
[color=black, fill=white] (2,4) circle (2pt) node[anchor=south west] $i^+$;
[color=black, fill=white] (0,0) circle (2pt) node[anchor=north west] $i^-$;
(2.45,3.55)  node[anchor=south west] $\mathcal{I}^+$;

Given $(\mathcal{M}^{(4+1)},g^{(4+1)})$ our extension of a Lorentzian Hawking–Page solution to an asymptotically flat spacetime with a locally naked singularity, there exists in the $SO(3)\times U(1)$ symmetry class a two dimensional space of perturbations to globally hyperbolic spacetimes without locally naked singularities, where the boundary $\mathcal{B}$ is preceded by trapped surfaces.

We consider the Einstein-Scalar field solution $\big(\mathcal{Q}^{(1+1)}\times S^2,g,\phi\big)$ associated to the $(4+1)$-dimensional vacuum spacetime $(\mathcal{M}^{(4+1)},g^{(4+1)})$. This solution has a locally naked singularity, and $\mu=\frac{1}{4}$ on $C_0^-,$ so the blue-shift is unbounded. However, the center $\Gamma^{(3+1)}$ is singular, so we cannot apply the theorem directly. We consider the characteristic initial data $\phi+\lambda_1f_1+\lambda_2f_2$ on $C_0^+,$ and obtain a unique future development of the data with the same argument that we used in Section <ref>. We then notice that the proof of the above theorem depends only on the solution on $C_0^-$ and the cone exterior region, so we obtain that any exterior incoming null cone in the perturbed solution reaches the trapped region. Using the Kaluza–Klein reduction, we have the same behavior for the corresponding $(4+1)$-dimensional vacuum solution.

In [6], Christodoulou proves Conjecture <ref> for the Einstein-Scalar field system in spherical symmetry, which then implies weak cosmic censorship for this matter model. In addition to the instability theorem, Christodoulou also proves that for AC solutions with locally naked singularities the blue-shift is unbounded, and that the one-dimensional spaces of perturbations obtained using $f_2$ are pairwise disjoint. In the case of the $(4+1)$-dimensional vacuum equations with an $SO(3)\times U(1)$ symmetry, one could also prove weak cosmic censorship as a consequence of Conjecture <ref>. In this context, our above proof illustrates how the blue-shift effect implies the instability of the locally naked singularities that we constructed. However, to prove the conjecture we also need an extension principle at the center to conclude that given a first singularity at the center the blue-shift is unbounded.

We could also consider proving the instability of the solutions with locally naked singularities $(\mathcal{M}^{(4+1)},g^{(4+1)})$ outside of symmetry. A potential candidate for perturbing could be constructed starting from localizing the functions $f_1,\ f_2$ near a fixed angle $(\theta_0,\varphi_0,\psi_0).$ An essential step for this problem is a formation of trapped surfaces result analogue to [3], but without symmetry assumptions. We remark that the short-pulse solutions of [8] cannot be used, since we need to allow for small perturbations and use the unbounded blue-shift on $C_0^-.$ Since the incoming cone is far from Minkowski space, the situation is closer to that of [17]. Finally, we point out that the problem of perturbing spherically symmetric naked singularities of the Einstein-Scalar field system outside of symmetry was considered in [16].

§ APPENDIX

§.§ Expansion of the solutions near the center

In this section we compute in detail the expansions used in Section <ref>. The goal is to study the autonomous system of ODEs:

\begin{equation}\label{appendix system}
\begin{dcases}
\frac{d\alpha}{ds}=\alpha\big[(\theta+k)^2+(1-k^2)(1-\alpha)\big] \\
\frac{d\theta}{ds}=k\alpha\big(k\theta-1\big)+\theta\big[(\theta+k)^2-(1+k^2)\big]
\end{dcases}
\end{equation}

in a neighborhood of the critical point $(\alpha,\theta)=(0,\pm\sqrt{3}).$ We set $x_1=\alpha,\ x_2=\theta\mp\sqrt{3},$ to obtain the system:

\begin{equation}
\begin{dcases}
\frac{dx_1}{ds}=2x_1+x_1\bigg[x_2^2\pm\frac{4}{\sqrt{3}}x_2-\frac{2}{3}x_1\bigg] \\
\frac{dx_2}{ds}=-2kx_1+4x_2+\frac{1}{3}x_1x_2\pm\frac{7}{\sqrt{3}}x_2^2+x_2^3
\end{dcases}
\end{equation}

We diagonalize the linear part of this system by setting $y_1=x_1,\ y_2=-kx_1+x_2$. We obtain the system:

\begin{equation}\label{diagonal system}
\begin{dcases}
\frac{dy_1}{ds}=2y_1-2y_1^2\pm\frac{4}{\sqrt{3}}y_1y_2+U_1^{(3)}(y_1,y_2) \\
\frac{dy_2}{ds}=4y_2\pm\frac{7}{\sqrt{3}}y_2^2-3y_1y_2+U_2^{(3)}(y_1,y_2)
\end{dcases}
\end{equation}

where $U_1^{(3)}(y_1,y_2),\ U_2^{(3)}(y_1,y_2)$ are homogeneous polynomials of degree 3. The critical point $(0,0)$ is stable to the past. The eigenvalues $\lambda_1=2,\ \lambda_2=4$ are resonant, since $(2,0)\cdot(\lambda_1,\lambda_2)=\lambda_2.$ Using the Poincaré-Dulac theorem, we obtain that the solutions have the form:

\begin{equation}
\begin{dcases}
y_1=z_1+Q_1^{(2)}(z_1,z_2) \\
\end{dcases}
\end{equation}

where $Q_1^{(2)}(z_1,z_2),\ Q_2^{(2)}(z_1,z_2)$ are analytic function that vanish to order 2 at the origin, and $z_1,\ z_2$ are solutions to the system:

\begin{equation}
\begin{dcases}
\frac{dz_1}{ds}=2z_1\\
\frac{dz_2}{ds}=4z_2+C_0z_1^2
\end{dcases}
\end{equation}

We compute $z_1=a_1e^{2s},\ z_2=a_2e^{4s}+C_0a_1^2se^{4s}.$ We introduce the equivalence relation $(a_1,a_2)\sim(Aa_1,A^2a_2)$, for some $A>0,$ and we denote by $\mathfrak{a}$ the set of equivalence classes. There is a one to one correspondence between $\mathfrak{a}$ and the orbits of (<ref>) near the origin.  We write:

\[Q_1^{(2)}(z_1,z_2)=b^{ij}z_iz_j+Q_1^{(3)}(z_1,z_2),\ Q_2^{(2)}(z_1,z_2)=d^{ij}z_iz_j+Q_1^{(3)}(z_1,z_2),\]
and we obtain that $C_0=0$. This corresponds to the fact that even though the eigenvalues are resonant, the resonant monomial $\begin{pmatrix} 0 & y_1^2 \end{pmatrix}$ does not appear on the right hand side of (<ref>). We also obtain $b^{11}=-1.$ Thus:

\[y_1=a_1e^{2s}-a_1^2e^{4s}+O(e^{6s}),\ y_2=\big(a_2+a_1^2d^{11}\big)e^{4s}+O(e^{6s}).\]

Finally, in a neighborhood of the negatively stable critical point $(0,\pm\sqrt{3})$, the orbits of (<ref>) satisfy the expansion:

\[\alpha=a_1e^{2s}-a_1^2e^{4s}+O(e^{6s}),\]

\[\theta=\pm\sqrt{3}\mp\frac{a_1}{\sqrt{3}}e^{2s}+\bigg(\pm\frac{a_1^2}{\sqrt{3}}+a_1^2d^{11}+a_2\bigg)e^{4s}+O(e^{6s}),\]

for some $(a_1,a_2)\in\mathfrak{a}.$

§.§ Sectional curvature computation

In this section, we compute the sectional curvatures of $S^2$ and $\mathcal{Q}^{(1+1)}$ with respect to the metric:

\[\Tilde{g}=-\frac{1}{f}e^{2\nu}du^2-\frac{2}{f}e^{\nu+\lambda}dudr+\frac{1}{f}r^2d\sigma_2^2+f^2d\sigma_1^2.\]

We express these using the corresponding $4$-dimensional metric:

\[g=-e^{2\nu}du^2-2e^{\nu+\lambda}dudr+r^2d\sigma_2^2,\]

and the quantities defined in Section <ref>.

In general, we compute for $a,b,c,d\neq\psi:$

\[\tilde{R}_{abc}{ }^d=R_{abc}{ }^d-\delta^d{ }_{[a}\nabla_{b]}\nabla_c\log f+g_{c[a}\nabla_{b]}\nabla^d\log f+\frac{1}{2}\big(\nabla_{[a}\log f\big)\delta^d{ }_{b]}\nabla_c\log f-\]

\[-\frac{1}{2}\big(\nabla_{[a}\log f\big)g_{b]c}\nabla^d\log f-\frac{1}{2}g_{c[a}\delta^d{ }_{b]}g(\nabla\log f,\nabla\log f).\]

We apply this formula to compute:

\[\tilde{K}_{S^2}=\frac{f}{r^2}\tilde{R}_{\theta\varphi\theta}{ }^{\varphi}=\frac{f}{r^2}\bigg[R_{\theta\varphi\theta}{ }^{\varphi}+\nabla_{\theta}\nabla_{\theta}\log f-\frac{1}{4}r^2g(\nabla\log f,\nabla\log f)\bigg]=\]

\[=\frac{f}{r^2}\bigg[\frac{2m}{r}-k(1-\mu)(\theta+\zeta)-k^2r^2g(\nabla\phi,\nabla\phi)\bigg]=\frac{f}{r^2}\bigg[1-(1-\mu)(k\theta+1)(k\zeta+1)\bigg].\]

Similarly, we also compute:

\[\tilde{K}_{\mathcal{Q}^{(1+1)}}=\frac{f^2}{4}\tilde{R}(n,l,n,l)=f^2K_{\mathcal{Q}^{(1+1)}}+\frac{f^2}{2}\nabla_{l}\nabla_{n}\log f+\frac{f^2}{4}\nabla_{l}\log f\nabla_{n}\log f+\frac{f^2}{4}g(\nabla\log f,\nabla\log f)=\]

\[=\frac{f}{r^2}\bigg[1-(1-\mu)(k\theta+1)(k\zeta+1)\bigg]+(1+k^2)\frac{f^2}{r^2}(1-\mu)\theta\zeta.\]

§.§ Equations in double null coordinates

In this section, we write the Einstein-Scalar field equations in spherical symmetry for a self-similar solution in double null coordinates as a system of autonomous ODEs. We consider the metric:

\[g=-\Omega^2dudv+r^2d\sigma_2^2.\]

The Einstein-Scalar field equations are:

\begin{equation}\label{partial u}
\partial_u(\Omega^{-2}\partial_ur)=-r\Omega^{-2}(\partial_u\phi)^2
\end{equation}

\begin{equation}\label{partial v}
\partial_v(\Omega^{-2}\partial_vr)=-r\Omega^{-2}(\partial_v\phi)^2
\end{equation}

\begin{equation}\label{partial uv Omega}
\partial_u\partial_v\log\Omega=-\frac{1}{r}\partial_u\partial_v r-4\partial_u\phi\partial_v\phi
\end{equation}

\begin{equation}\label{partial uv r}
\partial_u\partial_vr^2=-\frac{1}{2}\Omega^2
\end{equation}

\begin{equation}\label{wave}
\end{equation}

We notice that equation (<ref>) is implied by the other equations. Moreover, we can rewrite the Raychaudhuri's equations (<ref>), (<ref>) as:

\begin{equation}\label{partial u m}
\partial_um=-2r^2\Omega^{-2}\partial_vr(\partial_u\phi)^2
\end{equation}

\begin{equation}\label{partial v m}
\partial_vm=-2r^2\Omega^{-2}\partial_ur(\partial_v\phi)^2
\end{equation}

where $m$ is the Hawking mass:

\begin{equation}\label{mass}
\end{equation}

Following Section <ref>, we assume we are in self-similar double null gauge, so $S=u\partial_u+v\partial_v$ and:

\[S\Omega=0,\ Sr=r,\ S\phi=-k.\]

As before, we consider the following quantities:

\[\theta=\frac{r\partial_v\phi}{\partial_vr},\ \zeta=\frac{r\partial_u\phi}{\partial_ur},\ \beta=\frac{v\partial_vr}{r}.\]

We introduce the self-similar coordinate:

\[y=-\frac{v}{u}.\]

The change of coordinates is given by:

\[\partial_v=-\frac{1}{u}\partial_y,\ \partial_u=\frac{v}{u^2}\partial_y+\partial_u.\]

Thus, in $(u,y)$ coordinates we have $S=u\partial_u.$ We recall the notation:

\[\chi=\phi+k\log(-u),\ R=-\frac{r}{u}.\]

The self-similar condition is equivalent to $\Omega=\Omega(y),\ R=R(y),\ \chi=\chi(y).$ We compute:

\[\beta=\frac{y\partial_yR}{R},\ \zeta=\frac{\theta\beta+k}{\beta-1},\ \Omega^2=-\frac{4}{1-\mu}\partial_yR\big(y\partial_yR-R\big).\]

We can rewrite the wave equation (<ref>) as:

\[y\partial_y\theta+\theta\bigg(1+\frac{\beta-1}{1-\mu}\bigg)+k=0.\]

Also, we can rewrite equation (<ref>) as:

\[y\partial_y\beta+\bigg(2-\frac{1}{1-\mu}\bigg)\beta(\beta-1)=0.\]

Using this in equations (<ref>),(<ref>) eventually gives:

\[\frac{1}{1-\mu}=\beta\theta^2-\frac{(\theta\beta+k)^2}{\beta-1}+1.\]

Finally, we conclude that the Einstein-Scalar field system is equivalent to:

\[\begin{dcases}
y\partial_y\theta=k\big(k\theta-1\big)+\beta\theta\big[(\theta+k)^2-(1+k^2)\big] \\
\end{dcases}\]

[1]
M. Anderson.
Remarks on evolution of spacetimes in 3 + 1 and 4 + 1 dimensions.
Class. Quant. Grav., 18:5199–5209, 2001.

[2]
M. Anderson.
Geometric aspects of the ads / cft correspondence.
IRMA Lect.Math.Theor.Phys., 8:1–31, 2005.

[3]
D. Christodoulou.
The formation of black holes and singularities in spherically
symmetric gravitational collapse.
Comm. Pure Appl. Math., 44:339–373, 1991.

[4]
D. Christodoulou.
Bounded variation solutions of the spherically symmetric
einstein-scalar field equations.
Comm. Pure Appl. Math., 46:1131–1220, 1993.

[5]
D. Christodoulou.
Examples of naked singularity formation in the gravitational collapse
of a scalar field.
Ann. Math., 140:607–653, 1994.

[6]
D. Christodoulou.
The instability of naked singularities in the gravitational collapse
of a scalar field.
Ann. Math., 149:183–217, 1999.

[7]
D. Christodoulou.
On the global initial value problem and the issue of singularities.
Class. Quantum Grav., 16(A23), 1999.

[8]
D. Christodoulou.
The Formation of Black Holes in General Relativity.
Monographs in Mathematics, European Math. Soc., 2009.

[9]
M. Dafermos.
On naked singularities and the collapse of self-gravitating higgs
Adv. Theor. Math. Phys., 9(4):575–591, 2005.

[10]
M. Dafermos.
Spherically symmetric spacetimes with a trapped surface.
Class. Quantum Grav., 22(11):2221–2232, 2005.

[11]
C. Fefferman and C.R. Graham.
Conformal invariants.
Astérisque, pages 95–116, 1985.

[12]
R. Geroch, E. Kronheimer, and R. Penrose.
Ideal points in spacetime.
Proc. R. Soc., A 327:545–567, 1972.

[13]
S. Hawking and G. Ellis.
The Large Scale Structure of Space-Time.
Cambridge University Press, 1973.

[14]
S. Hawking and D. Page.
Thermodynamics of black holes in anti-de sitter space.
Comm. Math. Phys., 87:577–588, 1983.

[15]
J. Kommemi.
The global structure of spherically symmetric charged scalar field
Comm. Math. Phys., 323(1):35–106, 2013.

[16]
J. Li and J. Liu.
Instability of spherical naked singularities of a scalar field under
gravitational perturbations.
J. Diff. Geom., 120:97–197, 2022.

[17]
J. Luk and I. Rodnianski.
Nonlinear interaction of impulsive gravitational waves for the vacuum
einstein equations.
Cam. J. Math., 5(4):435–570, 2017.

[18]
C. W. Misner.
Taub-nut space as a counterexample to almost anything.
Relativity Theory and Astrophysics I: Relativity and Cosmology,
8:160–169, 1967.

[19]
I. Rodnianski and Y. Shlapentokh-Rothman.
Naked singularities for the einstein vacuum equations: The exterior
solution, 2019.

[20]
Y. Shlapentokh-Rothman.
Naked singularities for the einstein vacuum equations: The interior
solution, 2022.

[21]
Y. Shlapentokh-Rothman.
Twisted self-similarity and the einstein vacuum equations, 2022.