One and two-phase cell cycle models

I. INTRODUCTION

The cell cycle is a series of events that take place in a cell leading to its replication. It is regulated by a complex network of protein interactions. For example, a relatively simple mathematical model of mammalian cell-cycle control consists of eighteen differential equations [32]. Usually the cell cycle is divided into four phases [2], [18], [31].

The first one is the growth phase G 1 with synthesis of various enzymes. The duration of the phase G 1 is highly variable even for cells from one species. The DNA synthesis takes place in the second phase S. In the third phase G 2 a significant protein synthesis occurs, which is required during the process of mitosis. The last phase M consists of nuclear division and cytoplasmic division. We consider also G 0 phase (quiescence). A cell can enter the G 0 phase from G 1 and may remain quiescent for a long period of time, possibly indefinitely, or after some period of time it can go back to the G 1 phase. The schematic model of the cell cycle is given in Fig. 1.

There are several mathematical models of the cell cycle. One can consider four-phase models [6], but the most popular are one or two-phase models. In one-phase models, we put together phases G 1 , S, G 2 , and M and neglect the phase G 0 . The second category are two-phase models. Biologists used to divide the whole cycle into interphase, which consists of G 1 , S and G 2 , and the mitotic phase M . From a mathematical point of view it is better to divide the cell cycle into the resting (or growth) phase A = G 1 with a random duration t A and the proliferating phase B which consists of the phases S, G 2 and M , and has an almost constant duration t B . In these models the phase G 0 is also neglected. There are also two-phase models with cells in the proliferating state (phases G 1 , S, G 2 and M ) and in the quiescent state [4], [41].

The mathematical models of the cell cycle are based on the concept of maturity, also called the physiological age. The maturity can be the size of a cell, its volume or DNA content. The core of the theory was formulated in the late sixties [24], [43]. A lot of new models appeared in the eighties and we can divide them into two groups. The first group contains discrete-time models (generational models) which describe the relation between the initial maturity of mother and daughter cells. In this group we have one-phase models [20] and two-phase models [57], [58].

The second group is formed by continuoustime models characterizing the time evolution of distribution of cell maturity. We consider two types of continuous-time one-phase models: with the division of a cell when the cell maturity has a given level [24], [42], [43], or with division at a random maturity (including size structured models) [5], [7], [11], [15], [29], [34], [52], [59]. The second type consists of two-phase models [1], [8], [27], [39], [56]. The above mentioned models describe the distribution of maturity in the whole population. We also investigate models given by piecewise deterministic Markov processes [54], which describe the evolution of maturity of consecutive descendants of a single cell [22], [29], [38]. Such models seem to be the most suitable for the description of the cell cycle because they do not include environmental components.

This paper provides an introduction to one and two-cycle models, with particular emphasis on models given by piecewise deterministic Markov processes. We also formulate some results concerning their long-time behaviour and compare asymptotic properties of discrete and continuous time models. The evolution of distribution of maturity in the discrete time models and in models given by piecewise deterministic Markov processes is described by stochastic (Markov) operators and semigroups [21], [47]. The results concerning asymptotic stability and sweeping of stochastic semigroups are based on papers [35], [36], [37]. We also present results concerning chaotic properties of some maturity structured models [49], [50].

II. ASYMPTOTIC PROPERTIES OF STOCHASTIC

OPERATORS AND SEMIGROUPS Now we introduce the notion of a stochastic operator and a stochastic semigroup. Then we present some results concerning asymptotic stability, sweeping and Foguel alternative for stochastic semigroups. These results will be applied to models of cell cycle presented in the next sections.

Let a triple (X, Σ, µ) be a σ-finite measure space. Denote by D the subset of the space

is called a stochastic semigroup if it is a strongly continuous semigroup and all operators P (t) are stochastic. Now, we introduce some notions which characterize the asymptotic behaviour of iterates of stochastic operators P n , n = 0, 1, 2, . . . , and stochastic semigroups {P (t)} t≥0 . The iterates of stochastic operators form a discrete-time semigroup and we can also use the notation P (t) = P t for their powers so that we formulate most of the definitions and results for both types of semigroups without distinguishing them. A stochastic semigroup {P (t)} t≥0 is asymptotically stable if there exists a density f * such that

From (1) it follows immediately that f * is invariant with respect to {P (t)} t≥0 , i.e. P (t)f * = f * for each t ≥ 0. A stochastic semigroup {P (t)} t≥0 is called sweeping with respect to a set B ∈ Σ if for every f ∈ D In order to formulate a result concerning asymptotic stability of a stochastic semigroup we need the following definition. A stochastic semigroup {P (t)} t≥0 is called partially integral if there exists a measurable function

for some t > 0 and

Theorem 2.1 ( [35]): Let {P (t)} t≥0 be a continuous-time partially integral stochastic semigroup. Assume that the semigroup {P (t)} t≥0 has a unique invariant density f * . If f * > 0 a.e., then the semigroup {P (t)} t≥0 is asymptotically stable.

The next result concerns the Foguel alternative [21], that is, when a stochastic semigroup {P (t)} t≥0 is asymptotically stable or sweeping from all compact sets. We assume additionally that (X, ρ) is a separable metric space, Σ = B(X) is the σ-algebra of Borel subsets of X and that the semigroup {P (t)} t≥0 is partially integral with the kernel q which satisfies the following condition: (K) for every x 0 ∈ X there exist ε > 0, t > 0 and a measurable function η ≥ 0 such that η(y) µ(dy) > 0 and

We define condition (K) for a stochastic operator P in the same way, remembering the notation P (t) = P t . Condition (K) is satisfied if, for example, for every point x ∈ X there exist t > 0 and y ∈ X such that the kernel q(t, •, •) is continuous in a neighbourhood of (x, y) and q(t, x, y) > 0. We need an auxiliary definition. We say that a stochastic semigroup {P (t)} t≥0 overlaps supports if for every f, g ∈ D there exists t > 0 such that

The support of any measurable function f is defined up to a set of measure zero by the formula

Theorem 2.2: Assume that {P (t)} t≥0 satisfies (K) and overlaps supports. Then {P (t)} t≥0 is sweeping or {P (t)} t≥0 has an invariant density f * with a support A and there exists a positive linear functional α defined on L 1 (X, Σ, µ) such that (i) for every f ∈ L 1 (X, Σ, µ) we have

(ii) if Y = X \A, then for every f ∈ L 1 (X, Σ, µ) and for every compact set F we have

In particular, if {P (t)} t≥0 has an invariant density f * with the support A and X \ A is a subset of a compact set, then {P (t)} t≥0 is asymptotically stable. The proof of Theorem 2.2 is based on theorems on asymptotic decomposition of stochastic operators [36,Theorem 1] and stochastic semigroups [36,Theorem 2] and it is given in [38]. Another consequence of [36,Theorem 2] it the following.

Corollary 1:

Assume that a continuous-time stochastic semigroup {P (t)} t≥0 satisfies condition (K) and has no invariant densities. Then {P (t)} t≥0 is sweeping from compact sets.

III. RUBINOW-TYPE MODELS

In all models in this section we consider a sequence of consecutive descendants of a single cell in a single line. One of the oldest models of cell cycle, introduced by Rubinow [43], is based on the concept of maturity and maturation velocity. However, experimental observations concerning cell populations, cultured under identical conditions for each member, revealed high variability of l in the population. It means that the population is heterogeneous with respect to cell maturation velocities and therefore mathematical models of the cell cycle should take into account maturation velocities. A model of this type was proposed by Lebowitz and Rubinow [24]. In their model the cell cycle is determined by its maturation velocity v, which is fixed at the birth of the cell and constant during the cell cycle. The relation between the maturation velocities of mother's v and daughter's cellsv is given by a probability transition function P(v, dv). Denote by t n the time, when a cell from the nth-generation splits. Then the maturity and the maturation velocity of a cell from the nth-generation are described by a stochastic process

The process ξ(t), t ≥ 0, has jumps at points t 0 , t 1 , t 2 , . . . and between jumps the pair (m(t), v(t)) satisfies the following system of differential equations

At jump points we have m(t n ) = 0 and

. It is not difficult to check that ξ(t), t ≥ 0, is a homogeneous Markov process. Observe that there is an increasing sequence of random times (t n ), called jump times, such that the sample paths (trajectories) of ξ(t) are defined in a deterministic way in each interval (t n , t n+1 ).

A process which has such properties is called piecewise deterministic.

The Lebowitz and Rubinow model can be identified with a one-dimensional stochastic billiard on the interval [0, 1]. Namely, consider a particle moving in the interval [0, 1] with a constant velocity. We assume that when the particle hits the boundary points 0 and 1, it changes its velocity according to the probability measures P 0 (−v, B) = P(v, B) and P 1 (v, −B) = P(v, B), respectively, where v > 0 and B is a Borel subset of (0, ∞). Observe that the PDMP defined in the Lebowitz-Rubinow model is given by

where m(t) and v(t) represent position and velocity of the moving particle at time t. Asymptotic properties of the general onedimensional stochastic billiard were studied in [30]. Based on that paper we briefly present properties of the Lebowitz-Rubinow model: for

The infinitesimal generator A of the semigroup {P (t)} t≥0 is given by the formula

The functions f from the domain of the operator A satisfy the boundary condition:

Observe that if this semigroup has an invariant density f * , then Af

given by

Since the boundary condition contains a kernel operator, one can check that the semigroup {P (t)} t≥0 is partially integral. We assume that K is irreducible which is equivalent to say that

for each measurable set B ⊆ (0, 1] of the Lebesgue measure 0 < |B| < 1, see [55, p. 334]. From irreducibility it follows that if an invariant density f * exists then it is unique and f * (v) > 0 for va.e. The question of the existence of an invariant density is nontrivial. If for example we assume that there exist C > 0 and γ > 0 such that

then an invariant density exists [30]. It means that irreducibility of K and condition (8) imply asymptotic stability of the semigroup. Now we consider the case when the semigroup {P (t)} t≥0 has no invariant density. Assume that the kernel q is continuous and bounded. Then the semigroup satisfies condition (K) and from Corollary 1 it follows that the semigroup is sweeping from compact sets. It means that

for every density f and every ε > 0. The sweeping property in this case means that the length of the cell cycle tends to infinity in the sense of distribution, which is not so a rare phenomenon in tissue cells. For example if q ≡ 1, then the semigroup has no invariant density, and consequently is sweeping from compact sets. It is interesting that in this example we have

where c is some constant. The Lebowitz-Rubinow model is a special case of the Rotenberg model [42]. In the Rotenberg model the maturation velocity can also change during the cell cycle. A new born cell inherits the initial maturation velocity from its mother according to a transition probability P(v, dv), as in the Lebowitz-Rubinow model. During the cell cycle it can change its maturation velocity with intensity ϕ(m, v), i.e., a cell with parameters (m, v) can change the maturation velocity in a small time interval of length ∆t with probability ϕ(m, v)∆t + o(∆t). We suppose that if (m, v) is the state of the cell at the moment of changing of the maturation velocity, then a new maturation velocity is drawn from a distribution P(m, v, dv). The process ξ(t) = (m(t), v(t)) describing consecutive descendants of a single cell is a PDMP which has jumps when cells split and random jumps during their cell cycles. Between jumps the pair (m(t), v(t)) satisfies system (6). If a jump is at the moment of division, then it is given by the same formula as in the Lebowitz-Rubinow model. If a jump is during the cell cycle, then

for each Borel subset B of (0, ∞). Sample graphs of maturity in the Rubinow, Lebowitz-Rubinow and Rotenberg models are presented in Fig. 2.

IV. BELL-ANDERSON-TYPE MODELS

Now we consider one-phase models which are based on two main assumptions: the maturity m grows with velocity g(m) and a cell can splits at a rate ϕ(m), i.e. a cell with maturity m divides during a small time interval of length ∆t with We assume that lim m→∞ Q(m) = ∞, which guaranties that each cell splits with probability one. Let ξ be the maturity of the cell at the moment of division and let η be a positive random variable with density e −x . Since P(ξ > m) = e Q(m0)−Q(m) = P(η > Q(m)−Q(m 0 )),

we have

and therefore the random variables Q −1 (Q(m 0 ) + η) and ξ have the same distribution. It is easy to check that the random variable ξ has the density

where λ(m) = h −1 (m). If we assume that the distribution of the initial maturity of mother cells has a density f , from the above formula we infer that the initial maturity of the daughter cells has the density

Then P is a stochastic operator on the space

In the continuous version of the above model we consider a sequence of consecutive descendants of a single cell. The maturity of cells can be described by the following homogeneous piecewise deterministic Markov process ξ(t). Let t n be the time when a cell from the nth-generation divides. If t n−1 ≤ t < t n , then the maturity satisfies the equation ξ (t) = g(ξ(t)). The process ξ(t) has a jump at the moment of the division of the cell:

where π(t, m 0 ) is the solution at time t of the equation m = g(m) with the initial condition m(0) = m 0 .

Maturity structured models have been investigated in many papers. Usually, we are interested in the behaviour of the density of the maturity u(t, m). We should underline that u(t, •) is not density in our probabilistic sense because the integral of u with respect to m may be different than one. Since in these models we consider the whole population, beside the rate of division ϕ there is also the rate of death µ. These models coincide with the model given by the process ξ(t), when ϕ = µ and the rate of division in the definition of ξ(t) equals 1 2 ϕ. We briefly present some examples of these models. The basic model was proposed by Bell and Anderson [7]. They assume that the size of the cell is a number m ∈ (m min , 1), 0 < m min < 1. In order not to cross 1, it is assumed that 1 mmin ϕ(m) = ∞. This model was studied and generalized, for example, in [11], [15], [59]. Versions of this model with unequal division were investigated in [3], [16], [19], [52]. The papers [5], [33] are devoted to a general model that includes also age structure. There were also considered versions with unbounded growth of cells [29], [53].

V. TWO-PHASE MODEL

We start with a short biological description of two phase-cell cycle models. The cell cycle is divided into the resting and proliferating phase. The duration of the resting phase is random variable t A which depends on the maturity of a cell. The duration t B of the proliferating phase is almost constant. Therefore, we assume that t B = τ , where τ is a positive constant. A cell can move from the resting phase to the proliferating phase with rate ϕ(m). We assume that cells age with unitary velocity and mature with a velocity g 1 (m) in the resting phase and with a velocity g 2 (m) in the proliferating phase. The age variable a in the proliferating phase is assumed to range from a = 0 at the moment of entering the proliferating phase to a = τ at the point of cytokinesis. The maturity of the daughter cell m is a function of the maturity of the mother cell m, i.e. m = h(m) (see Fig. 3).

We consider a version of the model studied in [38]. Now we collect the assumptions concerning the model: Denote by π i (t, m 0 ) the solution of the equation

with the initial condition m(0) = m 0 ≥ 0. Now, we introduce two auxiliary functions. Let

According to (M4) lim m→∞ Q(m) = ∞, which guaranties that each cell enters the proliferating phase with probability one. Under this notation the initial maturity of the daughter cells has density P f , if f is the analogous density of the mother cells, where P is the stochastic operator given by (10). Now we present some results concerning the operator P .

Theorem 5.1: The operator P satisfies the Foguel alternative, i.e. P is asymptotically stable or sweeping from compact sets. Theorem 5.1 was proved in [38]. The results from Theorem 5.2 were proved, respectively, (a) in [14], (b) in [25], and (c) in [45]. The sweeping property in this case means that the maturity of cells statistically tends to infinity, which also means that the length of the cell cycle tends to infinity. Now we consider a continuous version of the model. The cell cycle can be described as a piecewise deterministic Markov process. We consider a sequence of consecutive descendants of a single cell. Let s n be the time, when a cell from the nthgeneration enters a resting phase and t n = s n + τ be the time of its division. If t n−1 ≤ t < t n then the state ξ(t) = (a(t), m(t), i(t)) of the n-th cell is described by the age a(t), maturity m(t) and the index i(t), where i = 1 if the cell is in the resting phase and i = 2 if it is in the proliferating phase. Random moments t 0 , s 1 , t 1 , s 2 , t 2 , . . . are called jump times. Between jump times the parameters change according to the following system of equations:

The process ξ(t) changes at the jump points according to the following rules:

If m(t n−1 ) = m 0 , then the cumulative distribution function F of s n − t n−1 is given by F (t) = 1 − e Q(m0)−Q(π1(t,m0)) .

Then ξ(t) is a homogeneous Markov process. If the distribution of ξ(0) is given by a density Having a homogeneous Markov process ξ(t) with the property that if the random variable ξ(0) has a density f 0 , then ξ(t) has a density f t , we can define a stochastic semigroup {P (t)} t≥0 corresponding to ξ(t) by P (t)f 0 = f t . The proper choice of the space X of the values of the process ξ(t) plays an important role in investigations of the process and the semigroup {P (t)} t≥0 . We define X = X 1 ∪ X 2 , where

and µ is the product of the twodimensional Lebesgue measure and the counting measure on the set {1, 2} (see Fig. 4).

We need two additional assumptions: (14) and

for somem > m P . Condition (14) guarantees that with a positive probability each cell will have a descendant with a sufficiently small maturity in some generation and thanks to that property each two states from the interior of the set X communicate. Condition (15) seems to be technical but if

for all m ≥ m P , then all descendants of a single cell in the same generation have the same maturity at a given time t. It means that the cells have synchronous growth and we cannot expect the model to be asymptotically stable. In particular, if g 1 ≡ g 2 and h(m) = m/2, then (15) reduces to 2g 2 (m) = g 2 (2m) for some m > π 2 (τ, m P ).

A similar condition appears in many papers concerning size-structured models [5], [11], [15], [52], [54].

The following results are proved in [38]. and there is ε > 0 such that ϕ(m) ≥ ε for sufficiently large m, then the semigroup {P (t)} t≥0 is asymptotically stable. If Q(λ(m)) − Q(m) ≤ 1 for sufficiently large m and ϕ is bounded, then the semigroup {P (t)} t≥0 is sweeping. Remark 1: One can give an example of the operator P which is asymptotically stable but the semigroup {P (t)} t≥0 is sweeping. Such a case can happen when lim m→∞ ϕ(m) = 0. The explanation of this phenomenon is that in this example the rate of entering the proliferating phase is very small for large m. Then the mean length of the resting phase can be large and more and more mature cells dominate the population as t → ∞. In [48], [60] one can find the comparison of the discrete time model presented here with a two-phase model of maturity structured population considered in the paper [27] and briefly presented in the next section.

VI. TWO-PHASE POPULATION MODELS

Now we recall a two-phase maturity structured model of a cellular population from the paper [27]. The model is based on the same biological assumption as that of Section V, but we include also the mortality rates µ r (m) and µ p (m) in both phases. Denote by r(t, m, a) and p(t, m, a) the maturity-age distribution of resting and proliferating cells, respectively. We also assume that the rate of entering the proliferating phase ϕ depends on m and the total number of cells in the resting phaseR(t) = R(t, m) dm, where R(t, m) = r(t, m, a) da. The time evolution of p and r is described by the following system of equations:

with the boundary conditions

Additionally, we assume that µ p , µ r and ϕ do not depend on m. Integrating the above equations over the age variable a and using the boundary conditions we obtain ∂R ∂t

We recall that λ(m) = π 2 (−τ, h −1 (m)) is the maturity of the mother cell at the moment of entering proliferating phase, if the new born daughter cell has maturity m. Integrating both sides of (16) over m we obtain

The following result is proved in [27]. Theorem 6.1: Assume that equation 17has a constant solutionR 0 > 0 andR 0 is globally asymptotically stable. If

then there exists a stationary solution R 0 (m) of equation (16) and for every solution R(t, m) of it we have

Condition (18) has an interesting biological interpretation. It shows that the stability of the population depends on the dynamics of immature (small) cells. The term (h −1 ) (0) describes the relation between the maturity of the mother and daughter cells. If m is the maturity of a small mother cell at the moment of entering the proliferating phase, then the maturity of a new-born daughter cell is m/(h −1 ) (0). The term c = µ r + ϕ(R 0 ) is the rate of leaving of the resting phase (by being lost or by entering the proliferating phase). Since g 1 (0) is the rate at which small cells mature, condition (18) means that the maturity of a large part of small cells will increase in the next generation. Now we present a model from the paper [28]. In this model we assume that the rate of entering proliferating phase ϕ for cell with maturity m depends on the total number of cells with this maturity R(t, m). We consider a simplified version of this model with g 1 (m) = g 1 m, g 2 (m) = g 2 m, h(m) = hm, where g 1 , g 2 , h are positive constants. We have also λ(m) = λm, λ > 0.

Then equation 16becomes a special case of the following nonlinear equation:

where u = R, x = m. Equations of the form (20) were used in description of cellular models in the papers [9], [12], [40]. We assume that the fuctions g :

We consider the solution of (20) with the initial condition

Now we consider the following delay differential equation associated with (20):

The following theorem plays a central role in investigations of equation (20).

Theorem 6.2: Let u(t, x) be a solution of (20). Let z(t) be the solution of (22) with the initial condition z(t) = u(t, 0) for t ∈ [−τ, 0]. Then for every t 0 ≥ 0 and ε > 0 there exist t 1 > 0 and another solutionū(t, x) of (20)

The proof of this result is given in [28]. From Theorem 6.2 the entire strategy of studying of equation (20) becomes clear. Namely, if z 0 (t) is a globally asymptotically stable solution of (22) and u 0 (x, t) = z 0 (t) is a locally asymptotically stable solution of (20), then u 0 (x, t) is globally asymptotically stable solution of (20). Thus, rather surprisingly, the question of determining the global stability of a solution of (20) can be reduced to the problem of examining the global stability of the corresponding differential delay equation (22) and the local stability of (20). Therefore, in the general case it is sufficient to focus on the global stability of the associated differential delay equation 22, which is itself usually quite difficult, and the local stability of (20), which is often easier.

We now turn to considerations of the local stability of the full partial differential equation (20). We assume that the function f does not depend on t. Then equation (20) takes the form

where u τ = u(t − τ, λ(x)). Letū(x, t) be a given solution of (23) and let A be a subset of C([0, 1] × [−τ, 0]). We say that the solutionū of (23) is exponentially stable on the set A if there exists µ > 0 such that for every ψ ∈ A the solution of the problem (23), (21) satisfies the inequality

where C is a constant which depends only on ψ. Let

We say thatū is locally exponentially stable if there exist an ε > 0, µ and C such that condition (24) holds for every solution of the problem (23), (21) with ψ ∈ A ε .

Theorem 6.3: Let w be a constant such that f (w, w) = 0 and ∂f ∂u (w, w) < − ∂f ∂u τ (w, w) .

Then the solutionū(x, t) ≡ w of (23) is locally exponentially stable. Let us summarize the results. Consider the associated delay differential equation corresponding to (23):

Let ϕ ∈ C[−τ, 0] and denote by z ϕ the solution of (26) satisfying the initial condition z ϕ (t) = ϕ(t) for t ∈ [−τ, 0]. Let w ∈ R be a constant such that f (w, w) = 0. Then w is a stationary solution of (26) In paper [28] one can find applications of this theory to maturation structured models and to blood production systems.

VII. CHAOS

Thus far, we have restricted our mathematical results to study such asymptotic properties of the models as the asymptotic stability and sweeping. But some of our models can have more complicated behavior which can be studied using theoretical methods of dynamical systems. Now we present some results concerning chaotic and ergodic properties.

It is not widely known that solutions of simple linear partial differential equations may behave in a chaotic way. The following equation with the initial condition: (27) defines a semiflow on the space

given by

which is chaotic, practically in each sense of the meaning of this word. For example, for λ > 0 there exists a Gaussian measure with the support X invariant under {S t } t≥0 and the system is mixing [44] (see also review papers [46], [51]). This implies the topological chaos: the existence of dense trajectories (topological transitivity) and instability of trajectories. The invariant measure µ for the semiflow {S t } t≥0 can be given by the formula µ(A) = P(ξ x ∈ A), where ξ x = w x 2c and w x is the standard Wiener process. Since equation (27) can describe the evolution of the distribution of maturity in a cellular population, one can prefer to consider the semiflow {S t } t≥0 restricted to the set X + = {v ∈ X : v ≥ 0}.

In this case the invariant measure on X + can be induced by the process ξ x = |w x 2c | and we still have very strong ergodic and chaotic properties of this semiflow. Similar results can be obtained for semiflows generated by equations of the form

In structured population models we mainly investigate the long time behavior of densities. Although in some models we expect chaotic behavior on the set of densities, it is rather difficult to prove such a result. One of the reasons is that there are no good mathematical tools to investigate such problems. Analytic methods of studying chaos in infinite dimensional spaces, based on the paper [10], do not work in this case because these methods are strictly connected with semiflows acting on the whole linear space. Now we present one model from the paper [49], which is well biologically motivated, where, by using methods of ergodic theory, we are able to prove chaotic behaviour of a semiflow acting on the set of densities.

We consider a population of stem cells. These cells live in the bone marrow and they are precursors of any blood cells. They are subjects of two biological processes: maturation and division. Stem cells can be at different levels of morphological development called maturity. The maturity of a cell differs from its age, because we assume that a newly born cell is in the same morphological state as its mother at the point of division. We assume that maturity is a real number x ∈ [0, 1]. The function u(t, x) describes the density distribution function of cells with respect to their maturity. The maturity grows according to the equation x = g(x). When one cell reaches the maturity 1 it leaves the bone marrow, then one of cells from the bone marrow splits. This cell is chosen randomly according to the distribution given by the density u(t, x). It follows from the assumptions that a newly born cell has the same maturity as its mother cell and each cell can divide with the same probability (see Fig. 5).

Although our model is rather simple in comparison with other models for erythroid production (e.g. [26], [23]), it is based on the same continuous maturation-proliferation scheme. We neglect here the fact that with the growth of maturity cells pass through consecutive morphological compartments from pluripotential stem cell to erythrocytes. In our model we do not have any exterior regulatory system in which the production of erythrocytes is stimulated by the hormone erythropoietin and the system tries to keep the number of erythrocytes on a constant level. One can say that chaos appears if the exterior regulatory system does not work (a pathological case).

The model is described by a nonlinear semiflow induced by the equation ∂u ∂t + ∂ ∂x (g(x)u) = g(1)u(t, 1)u(t, x)

with the initial condition

The semiflow is defined on the space of densities. In the paper [49] it was shown that the semiflow generated by the initial problem (29)- (30) posses an invariant measure which is mixing and supported on the whole set of all densities. From this result there follows instability of all trajectories and topological transitivity. The main idea of the proof is to show that the semiflow generated by (29)-(30) is isomorphic to a semiflow generated by (28). Then we construct an invariant measure for the second semiflow and transfer it to the initial semiflow. We skip the precise proof here. Fig. 6 presents the spatial temporal plot of the solution of (29) with the initial condition u 0 (x) = x sin 2 1

x .

If we change a little the initial condition replacing u 0 byū 0 (x) = x sin 2 1 x + 0.05 we obtain the solution with the plot shown in Fig. 7. The second example is the Bell and Anderson model of size structured cellular population given by the equation ∂u ∂t + ∂ ∂x (g(x)u) =−(µ+b)u(t,x)+4bu(t,2x), (31) where x ∈ [0, 1] and we put u(t, 2x) = 0 if 2x > 1. It is well known that if g(2x) = 2g(x) at least for one x ∈ [0, 1], then the solutions of (31) have asynchronous exponential growth, i.e., there exist λ ∈ R, positive functions f * , and a constant c which depends on the initial condition u(0, x) such that e −λt u(t,

Having in mind this result, it is difficult to imagine that some versions of this model can be chaotic. It is interesting that if g(x) = ax, then the semiflow generated by the equation 31is chaotic. The chaotic behaviour of the semiflow generated by this equation was studied using analytic methods by Howard [17] and El Mourchid et al. [13]. In [50] it was shown that for this semiflow there exists a mixing invariant measure supported on the whole space. From this property we can deduce chaotic properties of the semiflow.