Stability and Hopf Bifurcation in a Mathematical Model of Pluripotent Stem Cell Dynamics

# Stability and Hopf Bifurcation in a Mathematical Model of Pluripotent Stem Cell Dynamics

Mostafa Adimy,  Fabien Crauste  and  Shigui Ruan111Research was partially supported by the NSERC of Canada and the College of Arts and Sciences at the University of Miami. On leave from Dalhousie University, Halifax, Canada.
Year 2004
###### Abstract

We study a mathematical model describing the dynamics of a pluripotent stem cell population involved in the blood production process in the bone marrow. This model is a differential equation with a time delay. The delay describes the cell cycle duration and is uniformly distributed on an interval. We obtain stability conditions independent of the delay. We also show that the distributed delay can destabilize the entire system. In particularly, it is shown that Hopf bifurcations can occur.

Laboratoire de Mathématiques Appliquées, FRE 2570,

Université de Pau et des Pays de l’Adour,

Avenue de l’université, 64000 Pau, France.

E-mail: mostafa.adimy@univ-pau.fr, fabien.crauste@univ-pau.fr

Department of Mathematics, University of Miami,

P. O. Box 249085, Coral Gables, FL 33124-4250, USA.

E-mail: ruan@math.miami.edu

Keywords: Blood production system, stem cells, delay differential equations, stability, Hopf bifurcation.

## 1 Introduction

Blood production process, called hematopoiesis, is one of the major biological phenomena occurring in human body. It takes place in the bone marrow where pluripotent stem cells give birth to mature cells. After ejecting their nuclei, these cells enter the bloodstream and become blood cells.

According to the study of Burns and Tannock [4], the population of pluripotent stem cells can be divided into two distinct groups: quiescent cells and proliferating cells. Mathematical models describing the dynamics of this cell population have been studied since the end of the seventies, in particularly by Mackey [9, 10]. We refer to the review articles by Haurie et al. [8] and Mackey et al. [12] for further study and more references on this topic. More recently, Pujo-Menjouet et al. [14] and Pujo-Menjouet and Mackey [15] proved the existence of a Hopf bifurcation for the hematopoiesis model proposed in [9]. In all these works, the authors assumed that the proliferating phase duration is constant. Mathematically, this means that the delay in their models is a discrete delay. However, experimental data (see Bradford et al. [3]) indicate that cells do not spend the same time in the proliferating phase.

In this paper, taking into account this assumption, we assume that the delay (or proliferating phase duration) is uniformly distributed on an interval. The main objective is to investigate the effect of time delay on the dynamical solutions. It is shown that there exist some critical values of time delay such that a local Hopf bifurcation occurs at the non-trivial equilibrium.

The paper is organized as follows. In section 2, we present our model, which is given in equation (1). In section 3, we derive stability conditions for the two equilibria of equation (1) which do not depend on the delay. We show the existence of Hopf bifurcations at the non-trivial equilibrium in section 4. A brief discussion is given in section 5.

## 2 The Model

Pluripotent stem cells can be either in a resting phase, also known as -phase, or in a proliferating phase. In the resting phase, they can die at a constant rate , which also includes the cellular differentiation, or be introduced in the proliferating phase at a rate . According to the work of Sachs [16], is assumed to depend on the resting phase population.

In the proliferating phase, which is in fact the so-called cell cycle, pluripotent stem cells are committed to divide and give birth to two daughter cells at the end of this phase. The two daughter cells enter directly the resting phase and complete the cycle. We assume that proliferating cells divide according to a uniform law on an interval with . This assumption comes from the fact that, even if only a little is known about phenomena involved in hematopoiesis, there are strong evidences (see Bradford et al. [3]) indicating that cells do not divide at the same age. The function is then defined by

 f(r)=⎧⎨⎩1τ−τmin, if r∈[τmin,τ],0, otherwise.

Let denote the pluripotent stem cell population density (cells/kg) at time . It satisfies the nonlinear delay differential equation

 x′(t)=−(δ+β(x(t)))x(t)+2τ−τmin∫ττminβ(x(t−r))x(t−r)dr. (1)

The first term in the right-hand side of equation (1) accounts for the cellular loss due to mortality and cellular differentiation, , and introduction in the cell cycle, . The second term is for the division of proliferating cells into two daughter cells during mitosis. Proliferating cells are in fact resting cells introduced in the proliferating phase one generation earlier, so that the quantity appears with a time delay. The factor 2 is, of course, for the division of each proliferating cell into two daughter cells.

In the following, the rate of reintroduction in the proliferating compartment is taken to be a monotone and decreasing Hill function, given by

 β(x)=β0θnθn+xn for x≥0.

The coefficient is the maximum rate of reintroduction, is the -phase population density for which the rate of re-entry attains its maximum rate of change with respect to the resting phase population, and describes the sensitivity of with changes in the population. This function was firstly used in hematopoiesis models by Mackey [9] in 1978.

In [9] and [11], Mackey gave values of the above parameters for a normal human body production. These values are

 δ=0.05 d−1,β0=1.77 d−1 % and n=3. (2)

The value of is usually However, since we shall study the qualitative behavior of the pluripotent stem cells population, the value of is not really important and could be normalized without loss of generality.

Now if we consider an initial continuous nonnegative function defined on , then the equation (1) has a unique continuous and nonnegative solution , defined for , such that

 xφ(s)=φ(s) for s∈[−τ,0].

This can be obtained by using the results in Hale and Verduyn Lunel [7].

Notice that equation (1) has at most two equilibria, the trivial equilibrium and a non-trivial positive equilibrium . The trivial equilibrium always exists and corresponds to the extinction of the population.

###### Proposition 2.1.

Equation (1) has a non-trivial positive equilibrium if and only if

 β0>δ>0. (3)

In this case, is explicitly given by

 x∗=θ(β0δ−1)1/n.
###### Proof.

Let be an equilibrium of equation (1). Then satisfies

 x∗(β(x∗)−δ)=0.

Consequently, equation (1) has a non-trivial equilibrium if and only if the equation

 β(x∗)=δ

has a non-trivial solution. Since the function is decreasing and positive with , then equation (1) has a non-trivial equilibrium if and only if condition (3) holds. ∎

In the next section, we shall study the stability of the two equilibria of equation (1).

## 3 Stability

Throughout this section, we are interested in the stability of the equilibria of equation (1), in particularly the stability of the non-trivial equilibrium . We start by giving a result on the global stability of the trivial equilibrium of (1).

###### Theorem 3.1.

The trivial equilibrium of equation (1) is globally stable if

 β0<δ.
###### Proof.

The proof uses a similar technique employed by Adimy and Crauste [1]. It is based on the construction of a Lyapunov functional.

Denote by the space of all continuous nonnegative functions on . Let be the function defined by

 B(x)=∫x0β(s)s ds for x≥0.

Consider the mapping defined, for , by

 J(φ)=B(φ(0))+1τ−τmin∫ττmin∫0−r(β(φ(a))φ(a))2dadr.

Then,

 ˙J(φ)=˙φ(0)β(φ(0))φ(0)+1τ−τmin∫ττmin(β(φ(0))φ(0))2−(β(φ(−r))φ(−r))2dr.

Since

 ˙φ(0)=−(δ+β(φ(0)))φ(0)+2τ−τmin∫ττminβ(φ(−r))φ(−r)dr,

we obtain that

 ˙J(φ)=−(δ+β(φ(0)))β(φ(0))φ2(0)+2τ−τmin∫ττmin(β(φ(0))φ(0))2dr−1τ−τmin∫ττmin(β(φ(0))φ(0)−β(φ(−r))φ(−r))2dr.

Hence,

 ˙J(φ)≤−(δ−β(φ(0)))β(φ(0))φ(0)2.

Let be the function defined, for , by

 α(x)=(δ−β(x))β(x)x2.

Assume that Since is a decreasing function, it follows that the function is positive for . Hence, is nonnegative on and if and only if . Consequently, the mapping is a Lyapunov functional when . We then deduce that the trivial equilibrium of (1) is globally stable. ∎

The result in Theorem 3.1 describes the fact that when is the only equilibrium of (1), the population is doomed to extinction except when .

Now we focus on the stability of the positive equilibrium of equation (1). To ensure the existence of the equilibrium , we assume that condition (3) holds; that is,

 β0>δ>0.

We do not expect to obtain conditions for the global stability of . However, local stability results can be obtained by linearizing equation (1) about . Set

 β∗:=ddx(β(x)x)∣∣x=x∗=δ(1−nβ0−δβ0). (4)

The linearization of equation (1) at is

 x′(t)=−(δ+β∗)x(t)+2β∗τ−τmin∫ττminx(t−r)dr.

The characteristic equation of (1) is given by

 Δ(λ):=λ+δ+β∗−2β∗τ−τmin∫ττmine−λrdr=0. (5)

We now state and prove our first result on the stability of .

###### Theorem 3.2.

Assume that

 nβ0−δβ0≤1.

Then the non-trivial equilibrium of equation (1) is locally asymptotically stable.

###### Proof.

We show, in fact, that is stable when . By the definition of given by (4), it follows that

 β∗≥0 if and only if nβ0−δβ0≤1.

So we assume that .

We first assume that , given by (5), is a real function. Then, is continuously differentiable and its first derivative is given by

 dΔdλ(λ)=1+2β∗τ−τmin∫ττminre−λrdr. (6)

One can see that is positive for as soon as . Moreover,

 limλ→−∞Δ(λ)=−∞ and limλ→+∞Δ(λ)=+∞.

Consequently, has a unique real root . Since

 Δ(0)=δ−β∗=nδ(1−δβ0)>0,

we deduce that .

Now, we show that if is a characteristic root of equation (5), then . By contradiction, we assume that there exists a characteristic root of equation (5) such that . By considering the real part of , we obtain that

 μ+δ+β∗−2β∗τ−τmin∫ττmine−μrcos(ωr)dr=0.

Consequently,

This yields a contradiction. We conclude that every characteristic root of (5) is such that . Hence, all characteristic roots of (5) have negative real parts and the equilibrium is locally asymptotically stable. ∎

When , that is, when

 1

the stability cannot occur for all values of and . In particularly, we shall show that a Hopf bifurcation can occur (see Theorem 4.1). However, we can still have the stability of the non-trivial equilibrium for values of , and if is not too large. This will be considered in the next theorem.

To present the results, without loss of generality we assume that

 τmin=0.

We want to point out that the results we are going to show remain true when , but the proof is more complicated.

Define a function , for , by

 K(x)=sin(x)x (7)

and let be the unique solution of the equation

 x1=tan(x1),x1∈(π,3π2).

Set

 u0:=cos(x1)∈(−1,0).

Then and

 u0=K(x1)=minx≥0K(x). (8)

We have the following local stability theorem.

###### Theorem 3.3.

Assume that

 1

Then the non-trivial equilibrium of equation (1) is locally asymptotically stable.

###### Proof.

Let us assume that (9) holds. Then , and

 δ+β∗2β∗

By contradiction, assume that there exists a characteristic root of (5) with . Then,

 μ=−(δ+β∗)+2β∗ωτ∫ωτ0e−μωrcos(r)dr

and

 ω=−2β∗ωτ∫ωτ0e−μωrsin(r)dr.

Integrating by parts, we obtain that

 2μ=−(δ+β∗)+2β∗e−μτK(ωτ).

Consequently,

 μ<−(δ+β∗)+2β∗e−μτK(ωτ).

If is such that

 sin(ωτ)≥0,

then

 2β∗K(ωτ)<0

and

 μ<−(δ+β∗)≤0.

So we obtain a contradiction.

Similarly, if is such that

 sin(ωτ)<0,

then, from (8) and (10), we deduce that

 δ+β∗2β∗≤K(ωτ).

It implies that

 δ+β∗≥2β∗K(ωτ)>2β∗e−μτK(ωτ).

Therefore,

 μ<−(δ+β∗)+2β∗K(ωτ)≤0.

Again we obtain a contradiction. Hence, all characteristic roots of (5) are such that .

Now, we assume that (5) has a purely imaginary characteristic root . Then and satisfy

 K(ωτ)=δ+β∗2β∗. (11)

Using (8) and (10), we obtain a contradiction. Consequently, (11) has no solution and equation (5) does not have purely imaginary roots. We conclude that all characteristic roots of (5) have negative real parts and is locally asymptotically stable. ∎

From Theorems 3.2 and 3.3, it follows that the non-trivial equilibrium of equation (1) is locally asymptotically stable when

 0≤nβ0−δβ0<2(1−u0)1−2u0. (12)

We are going to show that as soon as condition (12) does not hold, then the equilibrium can be destabilized. In the next section, we shall show that if condition (10) does not hold, then a Hopf bifurcation indeed occurs at .

## 4 Hopf Bifurcations

In this section we are going to show that the non-trivial equilibrium of equation (1) can be destabilized via Hopf bifurcations. The time delay will be used as a bifurcation parameter. This result is obtained in Theorem 4.1.

Recall that the non-trivial equilibrium of equation (1) exists if and only if In the following, without loss of generality we assume that

 τmin=0.

Again the results still hold when , but the proof is easier to understand when .

We look for purely imaginary roots of . Of course, we assume that , otherwise is locally asymptotically stable. Let , with , be a purely imaginary characteristic root of equation (5). Then, and satisfy the following system

 {δ+β∗(1−2C(τ,ω))=0,ω+2β∗S(τ,ω)=0, (13)

where

 C(τ,ω)=1τ∫τ0cos(ωr)dr,S(τ,ω)=1τ∫τ0sin(ωr)dr.

First, one can see that cannot be a solution of (13). Otherwise

 δ=β∗<0.

Moreover, if is a solution of system (13), then is also a solution of (13). Hence, we only look for positive solutions .

One can check that and are given, for and , by

 C(τ,ω)=sin(ωτ)ωτ=K(ωτ),S(τ,ω)=1−cos(ωτ)ωτ,

where the function is defined by (7). Consequently, system (13) can be rewritten as

 K(ωτ) = δ+β∗2β∗, (14) cos(ωτ)−1(ωτ)2 = 12β∗τ. (15)

Consider the sequence

 {xk}k∈N:={x≥0 ; x=tan(x)}, (16)

with

 0=x0

In fact, one can check that

 {xk}k∈N={x≥0 ; K′(x)=0}.

Moreover, for all ,

 xk∈(kπ,kπ+π2).

Define two sequences and , for , by

 uk:=cos(x2k+1)<0,vk:=cos(x2k)>0.

Using the definition of , one can see that

 uk=K(x2k+1) and vk=K(x2k).

Thus, the sequence is increasing with and the sequence is decreasing with and for (see Figure 1).

Moreover,

 limk→+∞uk=limk→+∞vk=0.

Furthermore, one can check that, as soon as ,

 δ+β∗2β∗<1=v0.

Finally, define a function , for , by

 h(x)=2(1−x)1−2x

and set

 h(v0)=+∞.

We have the following results about the properties of the function

###### Lemma 4.1.

Suppose that

 h(u0)≤nβ0−δβ0 and δ+β∗≠0.

(i) If , then there exists such that

 h(uk)≤nβ0−δβ0

(ii) If , then there exists such that

 h(vk+1)≤nβ0−δβ0
###### Proof.

Since the function is increasing on the interval , we can see that

 h(uk)≤nβ0−δβ0

is equivalent to

 uk≤δ+β∗2β∗

and

 h(vk+1)≤nβ0−δβ0

is equivalent to

 vk+1≤δ+β∗2β∗

The lemma now follows. ∎

###### Proposition 4.1.

(i) If

 h(uk)

then system (14)-(15) has exactly solutions , …, and , …, with

 {ωl,1τl,1∈((2l−1)π,x2l−1), for l=1,…,k+1,ωl,2τl,2∈(x2l−1,2lπ), for l=1,…,k+1

and

 0<τ1,1<⋯<τk+1,1<τk+1,2<⋯<τ1,2.

(ii) If

 nβ0−δβ0=h(uk),k∈N,

then system (14)-(15) has exactly solutions , …, and , …, with

 ⎧⎪⎨⎪⎩ωl,1τl,1∈((2l−1)π,x2l−1), for l=1,…,k,ωl,2τl,2∈(x2l−1,2lπ), for l=1,…,k,ωk+1,1τk+1,1=x2k+1

and

 0<τ1,1<⋯<τk+1,1<τk,2<⋯<τ1,2.

(iii) If

 h(vk+1)

then system (14)-(15) has exactly solutions , …, and , …, with

 ⎧⎪⎨⎪⎩ω1,1τ1,1∈(π/2,π),ωl,1τl,1∈(xl+1,(l+2)π),for l=2,…,k+1,ωl,2τl,2∈((l+1)π,xl+1),for l=1,…,k

and

 0<τ1,1<⋯<τk+1,1<τk,2<⋯<τ1,2.

(iv) If

 nβ0−δβ0=h(vk),k∈N∗,

then system (14)-(15) has exactly solutions , …, and , …, , with

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩ω1,1τ1,1∈(π/2,π),ωl,1τl,1∈(xl,(l+1)π),for l=2,…,k,ωl,2τl,2∈((l+1)π,xl+1),for l=1,…,k−1,ωk,2τk,2=x2k

and

 0<τ1,1<⋯<τk,1<τk,2<⋯<τ1,2.

(v) If

 h(v1)

then system (14)-(15) has a unique solution such that and

 ω1τ1∈(0,π).
###### Proof.

We only prove when . The other cases can be deduced similarly. Assume that

 h(u0)

This is equivalent to

 u0<δ+β∗2β∗

The function is strictly negative and decreasing on with (see Figure 1). So the equation

 K(y)=δ+β∗2β∗

has a unique solution on the interval . Set

 τ1,1=(y1)22β∗(cos(y1)−1)

and . Then, is a unique solution of system (14)-(15) satisfying .

Moreover, the function is strictly negative and increasing on with , so the equation has a unique solution on the interval . Set

 τ1,2=(y2)22β∗(cos(y2)−1)

and . Then, is a unique solution of system (14)-(15) which satisfies .

Furthermore, the function is nonnegative on and

 u1=K(x3)=minx≥2πK(x).

Therefore, system (14)-(15) has two solutions, and .

Finally, using the fact that

 cos(y1)≤cos(y2),

we obtain that

 τ1,1<τ1,2.

This completes the proof. ∎

Lemma 4.1 and Proposition 4.1 give conditions for the existence of pairs of purely imaginary roots of equation (5). In the next proposition, we study the properties of the purely imaginary roots of (5).

###### Proposition 4.2.

Assume that there exists a such that equation (5) has a pair of purely imaginary roots for with . If

 ωcτc≠xk for all k∈N,

where the sequence is defined by (16), then are simple roots such that

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩dRe(λ)dτ∣∣τ=τc>0, if ωcτc∈(x2k,x2k+1),dRe(λ)dτ∣∣τ=τc<0, if ωcτc∈(x2k+1,x2k+2),k∈N.
###### Proof.

Assume that there exists a such that equation (5) has a pair of purely imaginary roots for with . Then, satisfies system (14)-(15).

Assume that

 ωcτc≠xk for all k∈N.

Let us show that are simple characteristic roots of (5). Using (6), one can see that are simple roots of (5) if

 1+2β∗∂S∂ω(τc,ωc)≠0 or ∂C∂ω(τc,ωc)≠0.

We will show that

 ∂C∂ω(τc,ωc)≠0.

A simple computation shows that

 ∂C∂ω(τc,ωc)=g(ωcτc)ω2cτc,

where the function is defined by

 g(x)=xcos(x)−sin(x)for x≥0.

One can check that if and only if there exists a such that . Moreover,

 g(x)>0 if and only if x∈(x2k+1,x2k+2), k∈N.

This yields that

 ∂C∂ω(τc,ωc)<0 if ωcτc∈(x2k,x2k+1)

and

 ∂C∂ω(τc,ωc)>0 if ωcτc∈(x2k+1,x2k+2).

Hence, are simple characteristic roots of (5).

Let be a characteristic root of (5) such that . By separating the real and imaginary parts, we obtain that

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩μ(τ)+δ+β∗−2β∗τ∫τ0e−μ(τ)rcos(ω(τ)r)dr=0,ω(τ)+2β∗τ∫τ0e−μ(τ)rsin(ω(τ)r)dr=0.

We denote by (respectively ) the first derivative of (respectively ) with respect to . For , we obtain that

 μ′(τc)[1+2β∗∂S∂ω(τc,ωc)]=2β∗∂C∂ω(τc,ωc)ω′(τc)\vspace1ex+2β∗τc(cos(ωcτc)−C(τc,ωc)) (17)

and

 ω′(τc)[1+2β∗∂S∂ω(τc,ωc)]=−2β∗∂C∂ω(τc,ωc)μ′(τc)\vspace1ex+2β∗τc(S(τc,ωc)−sin(ωcτc)). (18)

We consider two cases. First, assume that

 1+2β∗∂S∂ω(τc,ωc)=0. (19)

One can verify that

 1+2β∗∂S∂ω(τc,ωc)=2+(δ+β∗)τc.

Consequently, (19) is equivalent to

 τc=−2δ+β∗. (20)

Then, it follows from equation (18) that

 ∂C∂ω(τc,ωc)μ′(τc)=S(τc,ωc)−sin(ωcτc)τc.

Moreover, by using (14) and (15), we have

 S(τc,ωc)−sin(ωcτc)τc=1−(cos(ωcτc)+ωcτcsin(ωcτc))ωcτ2c=−δ+β∗4β∗ωc.

Hence, (20) implies that

 ∂C∂ω(τc,ωc)μ′(τc)=ωc2β∗τc<0.

Since we have

 μ′(τc)≠0.

Furthermore, the sign of is the same as the sign of

We now assume that

 1+2β∗∂S∂ω(τc,ωc)≠0.

Then, by using (17) and (18), we obtain that satisfies

 μ′(τc)[(1+2β∗∂S∂ω(τc,ωc))2+(2β∗∂C∂ω(τc,ωc))2]=2β∗τc[2β∗∂C∂ω(τc,ωc)(S(τc,ωc)−sin(ωcτ