Asymptotic phases in a cell differentiation model

J. Differential Equations 247 (2009) 736–769

Contents lists available atScienceDirect

Journal of Differential Equations

www.elsevier.com/locate/jde

Asymptotic phases in a cell differentiation model

Avner Friedman

a

, Chiu-Yen Kao

a

, Chih-Wen Shih

b,∗

a_{Mathematical Biosciences Institute, Department of Mathematics, The Ohio State University, OH 43210, United States} b_{Department of Applied Mathematics, National Chiao Tung University, Hsinchu, Taiwan 300, ROC}

a r t i c l e i n f o a b s t r a c t

Article history:

Received 15 November 2008 Available online 29 April 2009 Keywords: Cell differentiation Th1/Th2 cells Conservation law Multistationary Integro-differential equation Transcription factors

T cells of the immune system, upon maturation, differentiate into either Th1 or Th2 cells that have different functions. The decision to which cell type to differentiate depends on the concentrations of transcription factors T-bet (x1) and GATA-3 (x2). The population

density of the T cells, φ(t,x1,x2), satisfies a conservation law

∂φ/∂t_{+ (∂/∂}x1)(f1φ)+ (∂/∂x2)(f2φ)=gφ where fi depends on

(t,x1,x2)and, in a nonlinear nonlocal way, onφ. It is proved that,

as t→ ∞,φ(t,x1,x2) converges to a linear combination of 1, 2,

or 4 Dirac measures. Numerical simulations and their biological implications are discussed.

©2009 Elsevier Inc. All rights reserved.

1. Introduction

The development of a multicellular organism from a single fertilized egg cell to specialized cells depends on programs of gene expression. Following the initial stage of cell determination is a mat-uration process called differentiation by which cells acquire specific recognizable phenotypes and functions. In particular, the T lymphocytes of the immune system, upon maturation, differentiate into either Th1 or Th2 cells that have different functions. The decision to which of the cell type to differ-entiate depends on the concentration of transcription factors T-bet

(

x1)and GATA-3

(

x2). If x1is high (low) and x2 is low (high), the T cell will differentiate into Th1 (Th2).

A mathematical model by Yates et al. [15] describes the differentiation process in terms of two differential equations

dxi

dt

=

fi

(

t

,

x1,x2, φ) (i

=

1

,

2

),

(1.1)

*

Corresponding author. Fax: +886 3 5724679. E-mail address:cwshih@math.nctu.edu.tw(C.-W. Shih).

0022-0396/$ – see front matter ©2009 Elsevier Inc. All rights reserved.

where

φ (

t

,

x1,x2)is the population density of cells with concentration

(

x1,x2)at time t;

φ

satisfies the conservation of mass law

∂φ

∂

t

+

∂

∂

x1

(

f1φ)

+

∂

∂

x2

(

f2φ)

=

g

φ,

(1.2)

where g is the growth rate. Here fi

(

t

,

x1,x2, φ)is a nonlinear, nonlocal function of

φ (

t

,

x1,x2). In this paper we analyze the asymptotic behavior of

φ (

t

,

x1,x2

)

as t

→ ∞

. We prove that

φ (

t

,

x1,x2)

→



ω

j

δ

_(aj 1,a

j 2)

as t

→ ∞,

(1.3)

where the limit is a linear combination of Dirac measures at

(

a₁j

,

a₂j

)

, and the number of terms in the linear combination is 1, 2 or 4, depending on the parameters which occur in the definition of the fi. Conservation laws of the form (1.2), but with very different velocity terms

(

f1,f2), were considered in [6, Chapter 3], [7,8,16] and [9, Chapter 3], and some asymptotic estimates were derived in [6,7,9]. A theoretical study of bistable switches appeared in [3]. An analytic approach in studying multista-tionary dynamics for neural networks was reported in [2,12,14]. We finally note that mathematical models of differentiation of T cell and other cells appeared in [4,5] and [13], respectively; see also [1, Chapter 9].

2. The mathematical model

Lymphocytes are white blood cells that play important roles in the immune system. T cells and B cells are two major types of lymphocytes. B cells produce antibodies against pathogens while T cells are involved in autoimmunity. Th lymphocytes represent a subtype of T cells that are identified by the presence of surface antigens called CD4; they are referred to as CD4+T cells. Other subtypes of T cells include cytotoxic T cells (CD8+) and regulatory T cells. Th cells are the most numerous of the T cells in a healthy person. After an initial antigenic stimulation, Th lymphocytes differentiate into either one of two distinct types of cells called Th1 and Th2. Th1 cells make IFN

γ

that combat intracellular pathogens, and this immune response, if abnormal, is associated with inflammatory and autoimmune diseases. Th2 cells produce cytokines that activate B cells to produce antibodies against extracellular pathogens; this response, if abnormal, is associated with allergies such as asthma. Whether a precur-sor Th cell (henceforth to be denoted by Th0) becomes Th1 or Th2 depends on ‘polarizing’ signals.

The Yates et al. [15] model of Th differentiation is based on the interaction of two transcription factors, T-bet and GATA-3. High protein level of T-bet or GATA-3 corresponds to the Th1 phenotype or the Th2 phenotype. We shall denote by S1 and S2 the Th1 and Th2 polarizing cytokines, and by x1 and x2 the concentrations of T-bet and GATA-3, respectively, in a Th0 cell. Then the dynamics of x1 and x2 is described by dx1 dt

= −

μ

x1

+



α

1 xn₁ kn₁

+

xn₁

+

σ

1 S1

ρ

1

+

S1



·

1 1

+

x2/

γ

2

+ β

1, (2.1) dx2 dt

= −

μ

x2

+



α

2 xn 2 kn₂

+

xn₂

+

σ

2 S2

ρ

2

+

S2



·

1 1

+

x1/

γ

1

+ β

2. (2.2)

The first term on the right-hand side of each equation represents the rate of protein degradation. The last term

β

i is the constant basal rate of protein synthesis. The autoactivation rate of protein xi is represented by the term

α

i

xn i

738 A. Friedman et al. / J. Differential Equations 247 (2009) 736–769

where n is the Hill exponent that tunes the sharpness of the activation switch. The contribution of external signaling to the rate of growth in xi is given by the term

σ

i

Si

ρ

i

+

Si

.

The cross-inhibition between x1 and x2 occurs at both the autoactivation level and external (mem-brane) signaling level, and is represented by the cross-inhibition factors

1 1

+

xi

/

γ

i

.

The parameter

γ

i represents the value of xi at which the ratio of production of xj, i

=

j (due to the combined autoactivation and external signaling) is halved.

We denote by

φ (

t

,

x1,x2) the population density of CD4+ T cells with concentration

(

x1

,

x2) at time t. Then the total levels of expression of T-bet and GATA-3, at time t in the cell population are given, respectively, by



xi

φ (

t

,

x1,x2)dx1dx2, i

=

1

,

2

.

If we denote by Ci

(

t

)

the exogenous (non-T cell) signals that stimulate T-bet and GATA-3 expressions, then the total signal Siis given by

Si

(

t

)

=

Ci

(

t

)

+



xi

φ (

t

,

x1,x2)dx1dx2



φ (

t

,

x1,x2)dx1dx2

,

i

=

1

,

2

.

(2.3)

Here, a normalization by total cell numbers is adopted to impose the limitation of access to cytokines due to cell crowding. The evolution of the population density is then derived from the equation of continuity, or mass conservation law:

∂φ

∂

t

+

∂

∂

x1

(

f1φ)

+

∂

∂

x2

(

f2φ)

=

g

φ,

(2.4) where f1



x1,x2,S1(t

)



= −

μ

x1

+



α

1 xn₁ kn₁

+

xn₁

+

σ

1 S1(t

)

ρ

1

+

S1(t

)



·

1 1

+

x2/

γ

2

+ β

1, (2.5) f2



x1,x2,S2(t

)



= −

μ

x2

+



α

2 xn₂ kn₂

+

xn₂

+

σ

2 S2(t

)

ρ

2

+

S2(t

)



·

1 1

+

x1/

γ

1

+ β

2. (2.6)

In [15], the extrinsic and intrinsic cytokine interactions during the differentiation process were described in detail. Several numerical simulations have been made there to illustrate the changes of percentage of population under varying magnitudes of stimulus. Switches of population between Th0 to Th2 (high GATA-3) or from Th1 (high T-bet) to Th0, and then to Th2, under various levels of stimulus by extrinsic cytokines IL4 and IL12 were demonstrated.

The primary aim of the present paper is to analyze the behavior of the dynamical system (2.1)– (2.2) and the associated conservation law (2.4). We prove that when the parameters in (2.1)–(2.2) belong to a well-defined regime Pi, 1



i



6, the solution

φ (

t

,

x1,x2) will tend to 1-peak Dirac measure if i

=

1, 2-peak Dirac measures if i

=

2

,

3

,

4

,

5 and 4-peak Dirac measure if i

=

6. We use numerical simulation to examine the intermediate behavior of

φ (

t

,

x1,x2), and to draw biological implications.

Note that (2.4) is associated with the velocity field described by dx1(t

)

dt

=

f1



x1

(

t

),

x2(t

),

S1(t

)



,

(2.7) dx2(t

)

dt

=

f2



x1

(

t

),

x2(t

),

S2(t

)



.

(2.8)

We consider (2.4) on a (closed) domain

Ω

= [

0

,

A1

] × [

0

,

A2

]

which is an attracting set for (2.7)–(2.8); for convenience, we choose

Ai

=

α

i

+

σ

i

+ β

i

μ

,

i

=

1

,

2

.

(2.9) We assume that

φ (

0

,

x1,x2)

|

∂Ω

=

0

,

and

φ (

t

,

x1,x2)

|

∂Ω

=

0 for all t

>

0

.

Assuming that g

=

g

(

t

)

, and setting G

(

t

)

=



₀tg

(

s

)

ds,

ψ(

t

,

x1,x2)

=

e−G(t)

φ (

t

,

x1,x2), we can replace (2.4) by

∂ψ

∂

t

+

∂

∂

x1

(

f1ψ)

+

∂

∂

x2

(

f2ψ)

=

0

,

(2.10) with Si

(

t

)

=

Ci

(

t

)

e−G(t) N0

+



xi

ψ(

t

,

x1,x2)dx1dx2 N0

,

(2.11)

where N0is the initial total population and the integral is taken over

Ω.

Let

Φ(

t

,

x1,x2) be the solution map (flow map) of (2.7)–(2.8) and let

Ω(

t

)

= Φ(

t

, Ω)

. Then the transport equation (2.10) yields

d dt



Ω(t)

ψ(

t

,

x1,x2)dx1dx2

=

0

.

Furthermore, if

Ω(

t

)

→ (¯

a1,a

¯

2)as t

→ ∞

then for any continuous function h

(

x1,x2),



Ω

740 A. Friedman et al. / J. Differential Equations 247 (2009) 736–769

i.e.,

ψ(

t

,

x1,x2)

→

N0δ(a¯1,a¯2) in measure as t

→ ∞.

(2.12)

In the subsequent sections we study the behavior of the solution of (2.7), (2.8) in conjunction with the behavior of

Ω(

t

)

.

In Section 3 we prove existence and uniqueness for the initial value problem of Eq. (2.10). In Sections 4–8, we establish the assertion (1.3) under some assumptions on the parameters of (2.5)– (2.6). Numerical simulations illustrating the dynamics of the single-cell model and the formation of peak-solutions as t increases are given in Section 9. In the concluding Section 10, we give a biological interpretation of our results.

3. Existence and uniqueness

We shall prove the existence and uniqueness for Eq. (2.10) with initial values

ψ

|

t=0

= ψ

0(x1,x2) in

Ω,

(3.1)

where

ψ0

vanishes on

∂Ω,



_Ω(0)

ψ0

=

N0,

ψ0,

∇ψ

0 are continuous functions in

Ω,

G

(

t

)

and Ci

(

t

)

are continuous functions for t



0

.

(3.2) Set f

= (

f1,f2)and write

f

=

f

(

t

,

x

, ψ)

=

F

(

x

)

+

H



t

,

x

, ψ(

t

,

·)



.

(3.3) The characteristic curves of (2.10) are given by

d

ξ

t,x d

τ

=

F



ξ

t,x

(

τ

)



+

H



t

, ξ

t,x

(

τ

), ψ(

τ

,

·)



,

0

<

τ

<

t

,

(3.4)

ξ

t,x

(

t

)

=

x

.

(3.5)

Note that if x

∈ Ω

then

ξ

t,x

(

τ

)

∈ Ω

for all 0



τ

<

t.

We introduce the space C1

(Ω)

of continuously differentiable functions

ψ(

x

)

with norm

ψ =

max

x∈Ω



ψ(

x

)

+ ∇

ψ(

x

)



and the space C1_T

(Ω)

of continuous functions

ψ(

t

,

x

)

in

Ω

T

= [

0

,

T

] × Ω

with continuous derivative

∇

x

ψ(

t

,

x

)

in

Ω

T, and with norm

ψ

T

=

max

x∈Ω,0tT



ψ(

t

,

x

)

+ ∇

x

ψ(

t

,

x

)



.

Theorem 3.1. Under the condition (3.2) there exists a unique solution of (2.10), (3.1), with fi, Sidefined by (2.5), (2.6), (2.11), for all t

>

0 such that

ψ

∈

C1

Proof. Take any constant M, M

>

ψ

0



, and introduce the set

XM

=

ψ

∈

C1_T

(Ω),

ψ

T



M



for T small to be determined. We define a mapping W from XM into itself and prove that it has a unique fixed point. Given any

ψ

∈

XM, set

¯ψ =

W

(ψ)

where

¯ψ

is the solution of

∂ ¯

ψ

∂

t

+

f

(

t

,

x

, ψ)

· ∇

x

¯ψ = −



∇

x

·

f

(

t

,

x

, ψ)

 ¯

ψ ,

x

∈ Ω,

0

<

t

<

T

,

(3.6)

¯ψ|

t=0

= ψ

0, x

∈ Ω.

(3.7)

Using the representation

¯ψ(

t

,

x

)

= ¯ψ



ξ

t,x

(

0

)



−

t



0



∇

x

·

f



τ

, ξ

t,x

(

τ

), ψ(

τ

,

·)

 ¯

ψ



τ

, ξ

t,x

(

τ

)



d

τ

,

(3.8) we get max x∈Ω,0tT

ψ(

t

,

x

)

|

ψ0

|

L∞(Ω)

+

C T

,

where C is a constant which is actually independent of M.

Differentiating (3.6) with respect to xi and applying the preceding argument, we obtain a similar bound on ∂ ¯_∂ψ_x

i, so that

 ¯ψ

T

 ψ

0

 +

C T

<

M

if T is small enough. Hence W maps XM into XM. We next claim that W is a contraction. Indeed, given two functions

ψ1, ψ2

in XM, denote by

ξ

t1,x,

ξ

t2,x, the corresponding characteristic curves, and

set

¯ψ

i

=

W

(ψ

i

)

,

ψ

= ψ

1

− ψ

2, ¯

ψ

= ¯ψ

1

− ¯ψ

2.By ODE theory and (3.3),

ξ

_t1_,x

(

τ

)

− ξ

_t2_,x

(

τ

)



C T



max

x∈Ω,0tT

ψ(

t

,

x

)

.

(3.9)

Using the representation (3.8) for each

¯ψ

i, we deduce that max x∈Ω,0tT

¯

ψ(

t

,

x

)



C T



max x∈Ω,0tT

ψ(

t

,

x

)

.

Similarly we obtain a bound on

∇ ¯ψ(

t

,

x

)

by differentiating (3.6) with respect to xi, applying the previous argument, and using (3.9). Hence

 ¯ψ

T



C T

ψ

T

,

so that W is a contraction if T is small enough, and thus existence and uniqueness for (2.10), (3.1) follows for 0



t



T .

We can extend the solution step-by-step to all t

>

0 provided we can derive an a priori bound, say

742 A. Friedman et al. / J. Differential Equations 247 (2009) 736–769

where C

,

α

are constants. From (3.8) with

¯ψ = ψ

and (3.3) we get, by Gronwall’s inequality, sup

x∈Ω

ψ(

t

,

x

)



C

+

C eαt

.

Similarly, by differentiating (3.6) with respect to xi

,

we derive

sup

x∈Ω

∇

ψ(

t

,

x

)



C

+

C eα

t

.

Hence (3.10) holds and the proof of Theorem 3.1 is complete.

2

4. Single cell

We consider the single-cell model (2.1)–(2.2) in which S1,S2 are regarded as nonnegative con-stants. As we shall see, under some regimes of the parameter space, the system admits monostable, bistable, and quadstable phases. In order to study the dynamics of a single-cell, we introduce upper bounds

ˆ

fifor the functions fiin (2.5), (2.6):

ˆ

fi

(

xi

)

= −

μ

xi

+



α

i xn_i kn i

+

xni

+

σ

i Si

ρ

i

+

Si



+ β

i

,

i

=

1

,

2

.

(4.1)

Then

ˆ

fihas the following properties:

ˆ

fi

(

0

) >

0

,

ˆ

fi

(

0

) <

0

,

ˆ

fi

(

xi

) <

0 for Ai



xi

<

∞.

(4.2)

Let Bi

∈ (

0

,

Ai

)

be greater than the largest zero of f

ˆ

i, i

=

1

,

2. We also introduce lower bounds

ˇ

fi for fi:

ˇ

f1(x1)

= −

μ

x1

+



α

1 xn₁ kn 1

+

xn1

+

σ

1 S1

ρ

1

+

S1



·

1 1

+

B2/

γ

2

+ β

1, (4.3)

ˇ

f2(x2)

= −

μ

x2

+



α

2 xn₂ kn 2

+

xn2

+

σ

2 S2

ρ

2

+

S2



·

1 1

+

B1/

γ

1

+ β

2. (4.4) Indeed,

ˇ

f1(x1)



f1(x1,x2), for

(

x1,x2)

∈ [

0

,

A1

] × [

0

,

B2

],

ˇ

f2(x2)



f2(x1,x2), for

(

x1,x2)

∈ [

0

,

B1

] × [

0

,

A2

].

Note that

ˇ

fi

(

0

) >

0

,

ˇ

fi

(

0

) <

0

,

ˇ

fi

(

Bi

) <

0 for i

=

1

,

2

.

(4.5)

The functions

ˆ

fi

, ˇ

fi, extended to xi

∈ (

Ai

,

∞)

by the right-hand sides of (4.1), (4.3), (4.4), have a unique inflection point

˜ξ

i, given by

˜ξ

i

=

ki



n

−

1 n

+

1



1/n

,

where the slopes of

ˆ

fiand of

ˇ

fiare maximal. Therefore, if

ˇ

f_i

( ˜ξ

i

) <

0, then

ˇ

f_i

(

xi

)

cannot take positive values. Set

˜

n

= (

n

+

1

)

1+1/n

(

n

−

1

)

1−1/n

/

4n

.

We consider the following parameter regimes:

Condition (M1):

μ

>

α

1n

˜

k1

,

Condition (M2):

μ

>

α

2n

˜

k2

,

Condition (B1):

μ

<

α

1n

˜

k1

·

1 1

+

B2/

γ

2

,

Condition (B2):

μ

<

α

2n

˜

k2

·

1 1

+

B1/

γ

1

.

Condition (Mi) is equivalent to the inequality

ˆ

f_i

( ˜ξ

i

) <

0, i

=

1

,

2. Under this condition both

ˆ

fi and

ˇ

fiare strictly decreasing functions and have a unique zero.

Condition (Bi) is equivalent to

ˇ

f_i

( ˜ξ

i

) >

0 and, in that case, if

˜ξ

i

<

Ai then each of

ˆ

fi,

ˇ

fi has two critical points. Let p

ˆ

m

i

,

p

ˆ

Mi (respectively p

ˇ

mi

,

p

ˇ

Mi

)

be the local minimum and maximum of

ˆ

fi (respec-tively

ˇ

fi). Then,p

ˇ

m_i

<

p

ˇ

_iM,p

ˆ

m_i

<

p

ˆ

M_i , and

ˇ

fi



ˇ

pm_i



< ˆ

fi



ˆ

pm_i



,

ˇ

fi



ˇ

pM_i



< ˆ

fi



ˆ

pM_i



.

We shall consider only the following cases as illustrated in Fig. 1. (Note that if

˜ξ

i

>

Ai for i

=

1 or

i

=

2, then only case (Mi) can occur for this i.) (a) (Mi) holds for i

=

1

,

2;

(b) (Bi) holds and

ˆ

fi

(

p

ˆ

Mi

) <

0 for i

=

1

,

2; (c) (Bi) holds and

ˇ

fi

(

p

ˇ

mi

) >

0 for i

=

1

,

2;

(d) (Bi) holds and

ˆ

fi

(

p

ˆ

m_i

) <

0

, ˇ

fi

(

p

ˇ

M_i

) >

0 for i

=

1

,

2.

In cases (a), (b), and (c),

ˆ

fiand

ˇ

fihave a unique zero denoted bya

ˆ

ianda

ˇ

i, respectively. In case (d),

ˆ

fiand

ˇ

fi have three zeros, denoted by (a

ˆ

i

, ˆ

bi

,

c

ˆ

i) and (a

ˇ

i

, ˇ

bi

,

c

ˇ

i), respectively. We shall establish the following dynamical phases for (2.1)–(2.2): Monostable (MS): low x1–low x2; low x1–high x2; high x1–low x2;

high x1–high x2 states.

Bistable (BS-ll,lh): low x1–low x2state and low x1–high x2state; (BS-ll,hl): low x1–low x2 state and high x1–low x2 state; (BS-hl,hh): high x1–low x2 state and high x1–high x2 state; (BS-lh,hh): low x1–high x2 state and high x1–high x2 state. Quadstable (QS): low x1–low x2state, high x1–low x2 state,

low x1–high x2state, and high x1–high x2state.

These notions of ‘low’ and ‘high’ express only relative magnitude relations between x1 and x2. It will be shown that there exist six parameter regimes so that (2.1)–(2.2), with parameters in each of these regimes admit, respectively, a unique stable equilibrium; two stable equilibria and one unstable equilibrium; and four stable equilibria and five unstable equilibria. Moreover, every solution which is

744 A. Friedman et al. / J. Differential Equations 247 (2009) 736–769

Fig. 1. ˆf1and ˇf1have one zero in cases (a), (b), (c), and three zeros in case (d).

initially not an unstable equilibrium point converges to one of the stable equilibria as time tends to infinity.

In order to guarantee the convergence to equilibrium, we impose the following condition:

(

α

1

+

σ

1)

γ

2

·

(

α

2

+

σ

2)

γ

1

<

μ

−

α

1n

˜

k1

·

μ

−

α

2n

˜

k2

.

(4.6)

Theorem 4.1. Assume that condition (4.6) holds. Then

(i) phase (MS) takes place under conditions (M1) and (M2), or conditions (B1), (B2) with either

ˆ

f1(p

ˆ

M 1

) <

0,

ˆ

f2(p

ˆ

M

(ii) phase (BS-ll,lh) takes place under conditions (B2),

ˆ

f2(p

ˆ

m₂

) <

0

, ˇ

f2

(

p

ˇ

M₂

) >

0, and condition (M1), or (B1)

and

ˆ

f1(p

ˆ

M₁

) <

0;

(iii) phase (BS-ll,hl) takes place under condition (B1),

ˆ

f1(p

ˆ

m

1

) <

0

, ˇ

f1(p

ˇ

M1

) >

0, and condition (M2), or (B2),

ˆ

f2(p

ˆ

M 2

) <

0;

(iv) phase (BS-hl,hh) takes place under conditions (B2),

ˆ

f2(p

ˆ

m

2

) <

0

, ˇ

f2(p

ˇ

M2

) >

0, and condition (M1), or (B1) and

ˇ

f1(p

ˇ

m₁

) >

0;

(v) phase (BS-lh,hh) takes place under condition (B1),

ˆ

f1(p

ˆ

m₁

) <

0

, ˇ

f1(p

ˇ

M₁

) >

0, and condition (M2), or (B2),

ˇ

f2(p

ˇ

m₂

) >

0;

(vi) phase (QS) takes place under conditions (B1), (B2),

ˆ

fi

(

p

ˆ

m_i

) <

0

, ˇ

fi

(

p

ˇ

M_i

) >

0, for i

=

1

,

2.

The proof of Theorem 4.1 follows from an iteration scheme which is similar to that introduced in Sections 5–8; in order to avoid repetition, the proof is omitted.

Remark 4.1. Note that

Condition (B1) :

μ

<

α

1n

˜

k1

·

1 1

+

A2/

γ

2

,

Condition (B2) :

μ

<

α

2n

˜

k2

·

1 1

+

A1/

γ

1

imply, respectively (B1) and (B2). Moreover, with Ai defined in (2.9), if conditions (B1) and (B2) are satisfied then (4.6) holds. However, these conditions are more restrictive than conditions (B1), (B2), and are not involved with the cytokine rates

σ

1,

σ

2.

Remark 4.2. The conditions expressed by the signs of

ˆ

fi

(

p

ˆ

mi

), ˇ

fi

(

p

ˇ

Mi

)

depend on the levels of cy-tokines S1,S2. There exist parameters so that phase (QS) takes place if both S1 and S2are sufficiently large. With the same parameters, the dynamics reduces to phase (BS-ll,lh) (respectively (BS-ll,hl)) if S2 (respectively S1) is sufficiently small and reduces to phase (MS) if both S1 and S2 are sufficiently small. We shall illustrate this situation numerically in Section 9.

5. The population model

In the subsequent sections we shall consider the asymptotic behavior of

ψ(

t

,

x1,x2) and of the corresponding dynamical system (2.7)–(2.8) in case Si

=

Si

(

t

)

is defined by (2.11). Typically g

(

t

)

=

2 day−1 for some time t

<

t0 and g

(

t

)

=

0 if t

>

t0, but Ci

(

t

)

may not vanish for large t. Throughout this paper we assume that

Ci

(

t

)

→

Ci

(

∞) 

0

,

G

(

t

)

→

G

(

∞) >

0 as t

→ ∞.

(5.1) The derivation of the asymptotic behavior will be based on a sequence of approximations by means of upper bounds

ˆ

f_i(k) and lower bounds

ˇ

f_i(k) of fi

(

x1,x2,Si

(

t

))

. In this section we construct these functions for the case k

=

0. As in the discussion in Section 4, we introduce an upper bound for

fi

(

x1,x2,Si

(

t

))

:

ˆ

fi

(

xi

)

= −

μ

xi

+



α

i xn i kn_i

+

xn_i

+

σ

i

ˆ

Ci

+

Ai

ρ

i

+ ˆ

Ci

+

Ai



+ β

i

,

whereC

ˆ

i

=

sup

{

Ci

(

t

)

e−G(t)

/

N0: t

∈ [

0

,

∞)}

;

ˆ

fi clearly satisfies (4.2). Let Bibe the largest zero of

ˆ

fi. Thus,

[

0

,

B1

] × [

0

,

B2

]

is an attracting set for (2.7)–(2.8).

746 A. Friedman et al. / J. Differential Equations 247 (2009) 736–769

Fig. 2. ˆf1and ˇf1have two zeros.

Next we define a lower bound for f1 on

R × [

0

,

B2

]

and a lower bound for f2 on

[

0

,

B1

] × R

, respectively:

ˇ

f1(x1)

= −

μ

x1

+



α

1 xn 1 kn₁

+

xn₁

+

σ

1

ˇ

C1

ρ

1

+ ˇ

C1



·

1 1

+

B2/

γ

2

+ β

1,

ˇ

f2(x2)

= −

μ

x2

+



α

2 xn₂ kn 2

+

xn2

+

σ

2

ˇ

C2

ρ

2

+ ˇ

C2



·

1 1

+

B1/

γ

1

+ β

2,

whereC

ˇ

i

=

inf

{

Ci

(

t

)

e−G(t)

/

N0: t

∈ [

0

,

∞)}

, i

=

1

,

2;

ˇ

fi clearly satisfies (4.5). The functions

ˆ

fi

, ˇ

fi share other properties with those defined in Section 4. Indeed, under conditions (Mi), (Bi) with

ˆ

fi

(

p

ˆ

Mi

) <

0, or (Bi) with

ˇ

fi

(

p

ˇ

mi

) >

0, both

ˆ

fi and

ˇ

fi have a unique zero, denoted respectively bya

ˆ

i

,

a

ˇ

i; under conditions (Bi), each of

ˆ

fiand

ˇ

fihas a local minimum and a local maximum, denoted byp

ˆ

m_i

,

p

ˆ

M_i , and

ˇ

pm

i

,

p

ˇ

Mi , respectively, and it can be computed that

ˇ

fi

(

p

ˇ

mi

) < ˆ

fi

(

p

ˆ

mi

)

and

ˇ

fi

(

p

ˇ

iM

) < ˆ

fi

(

ˆ

pMi

)

. Furthermore, under conditions (Bi), and

ˆ

fi

(

p

ˆ

m_i

) <

0

, ˇ

fi

(

p

ˇ

M_i

) >

0, both

ˆ

fi and

ˇ

fi have three zeros, denoted by

(

a

ˆ

i

, ˆ

bi

,

c

ˆ

i

)

,

(

a

ˇ

i

, ˇ

bi

,

ˇ

ci

)

, respectively; cf. Fig. 2. Set Smin_i

(

t

)

=

inf

Si

(

s

)

: s

∈ [

t

,

∞)



,

Smax_i

(

t

)

=

sup

Si

(

s

)

: s

∈ [

t

,

∞)



for i

=

1

,

2 and t



0. Then Smin_i

(

t

)

 ˇ

Ci, Smax_i

(

t

)

 ˆ

Ci

+

Ai, and Smin_i

(

t

)



Si

(

t

)



Smax_i

(

t

)

. Note that

Smin_i

(

t

)

is nondecreasing, Smax_i

(

t

)

is nonincreasing, and

Smin i

(

t

)

ρ

i

+

Smini

(

t

)



Si

(

t

)

ρ

i

+

Si

(

t

)



Smax i

(

t

)

ρ

i

+

Smax_i

(

t

)

for i

=

1

,

2 and t



0

.

We formulate the first step for the iteration scheme via the functions

ˆ

f_i(0)

(

xi

)

= −

μ

xi

+



α

i xn_i kn i

+

xni

+

σ

i Smax_i

(

0

)

ρ

i

+

Smaxi

(

0

)



+ β

i for i

=

1

,

2

,

ˇ

f₁(0)

(

x1)

= −

μ

x1

+



α

1 xn₁ kn 1

+

xn1

+

σ

1 Smin₁

(

0

)

ρ

1

+

S1min

(

0

)



·

1 1

+

B2/

γ

2

+ β

1

,

ˇ

f₂(0)

(

x2)

= −

μ

x2

+



α

2 xn₂ kn 2

+

xn2

+

σ

2 Smin 2

(

0

)

ρ

2

+

S₂min

(

0

)



·

1 1

+

B1/

γ

1

+ β

2

.

Then

ˆ

f_i(0)

, ˇ

f_i(0)admit the same properties as in (4.2) and (4.5). Moreover,

ˇ

fi

(

xi

)

 ˇ

fi(0)

(

xi

) < ˆ

fi(0)

(

xi

)

 ˆ

fi

(

xi

),

i

=

1

,

2

.

Therefore,

ˆ

fi

(

p

ˆ

mi

) <

0 implies

ˆ

f (0) i

(

p

ˆ

mi

) <

0, whereas

ˇ

fi

(

ˇ

pMi

) >

0 implies

ˇ

f (0) i

(

p

ˇ

Mi

) >

0. In addition,

ˇ

f_i(0)

(

xi

) < ˆ

f(0)

i

(

xi

)

for all xi

∈ [

0

,

∞)

, and both of

ˆ

fi(0) and

ˇ

f

(0)

i have their inflection points at

˜ξ

i

=

ki

(

n_n+−1₁

)

1/nwhere they attain their largest slopes. Observe that

ˇ

f_i(0)

(

xi

)



fi



x1,x2,Si

(

t

)



 ˆ

f_i(0)

(

xi

)

(5.2)

for i

=

1

,

2 and

(

x1,x2)

∈ [

0

,

B1

] × [

0

,

B2

]

, t



0. In addition, for all t



0,

f1



x1,x2,S1(t

)



 ˆ

f₁(0)

(

x1) if

(

x1,x2)

∈ [

0

,

A1

] × [

B2,A2

],

(5.3) f2



x1,x2,S2(t

)



 ˆ

f₂(0)

(

x2) if

(

x1,x2)

∈ [

B1,A1

] × [

0

,

A2

].

(5.4) In the sequel, x

(

t

,

x0)denotes the solution of (2.7)–(2.8) starting from point x0 at t

=

0.

6. Asymptotic one-peak solution

Similarly to the case of Theorem 4.1(i) we assume that one of the following conditions holds:

(

M1

)

and

(

M2

)

;

(6.1)

(

B1

)

and

(

B2

)

with

ˆ

f1



ˆ

pM₁



<

0

, ˆ

f2



ˆ

pM₂



<

0

;

(6.2)

(

B1

)

and

(

B2

)

with

ˇ

f1



ˇ

pm₁



>

0

, ˇ

f2



ˇ

pm₂



>

0

.

(6.3)

Then each

ˆ

f_i(0)and

ˇ

f_i(0)has a unique zero which is denoted bya

ˆ

_i(0)anda

ˇ

_i(0), respectively. Let

ε

0

>

0 be small so that

ˆ

f_i(0)

(

xi

)

 ˆ

fi(0)



ˆ

a(_i0)

+

ε

0



<

0 for all xi

 ˆ

a(i0)

+

ε

0,

ˇ

f_i(0)

(

xi

)

 ˇ

fi(0)



ˇ

a(_i0)

−

ε

0



>

0 for all xi

 ˇ

a(i0)

−

ε

0,

for i

=

1

,

2; cf. Fig. 3. Combining these with inequalities (5.2)–(5.4), we deduce that there exists a T0

>

0 such that any solution x

(

t

,

x0) starting from a point x0

∈ [

0

,

A1

] × [

0

,

A2

]

falls into the rectangle

Ω

(0)

:=



a

ˇ

(₁0)

−

ε

0,a

ˆ

(₁0)

+

ε

0

×



a

ˇ

(₂0)

−

ε

0,a

ˆ

(₂0)

+

ε

0

⊂ [

0

,

B1

] × [

0

,

B2

]

748 A. Friedman et al. / J. Differential Equations 247 (2009) 736–769

Fig. 3. The configuration ofˆf₁(0), ˇf₁(0), ˆf₁(1), ˇf₁(1)and their zeros, under conditions (B1) andˆf1(pˆM1) <0.

for t



T0. Define

ˆ

f₁(1)

(

x1)

= −

μ

x1

+



α

1 xn₁ kn₁

+

xn₁

+

σ

1 Smax₁

(

T0)

ρ

1

+

Smax₁

(

T0)



·

1 1

+ (ˇ

a(₂0)

−

ε

0)/

γ

2

+ β

1,

ˇ

f₁(1)

(

x1)

= −

μ

x1

+



α

1 xn₁ kn₁

+

xn₁

+

σ

1 Smin₁

(

T0)

ρ

1

+

Smin1

(

T0)



·

1 1

+ (ˆ

a(₂0)

+

ε

0)/

γ

2

+ β

1,

ˆ

f₂(1)

(

x2)

= −

μ

x2

+



α

2 xn 2 kn 2

+

xn2

+

σ

2 Smax 2

(

T0)

ρ

2

+

Smax2

(

T0)



·

1 1

+ (ˇ

a(₁0)

−

ε

0)/

γ

1

+ β

2,

ˇ

f₂(1)

(

x2)

= −

μ

x2

+



α

2 xn₂ kn 2

+

xn2

+

σ

2 Smin₂

(

T0)

ρ

2

+

Smin2

(

T0)



·

1 1

+ (ˆ

a(₁0)

+

ε

0)/

γ

1

+ β

2.

Then

ˇ

f_i(0)

(

xi

) < ˇ

f_i(1)

(

xi

) < ˆ

f_i(1)

(

xi

) < ˆ

f_i(0)

(

xi

)

for xi

∈ [

0

,

Ai

]

, i

=

1

,

2. Let a

ˆ

(_i1) and a

ˇ

(_i1) denote the unique zeros of

ˆ

f_i(1)and

ˇ

f_i(1), respectively. Thena

ˆ

(_i1)

<

a

ˆ

(_i0)anda

ˇ

_i(1)

>

a

ˇ

(_i0). Furthermore,

ˇ

f_i(1)

(

xi

)



fi



x1,x2,Si

(

t

)



 ˆ

f_i(1)

(

xi

)

(6.4) for all

(

x1,x2)

∈ Ω

(0)_{, t}

_

_T

0, i

=

1

,

2, and

ˇ

f_i(1)

(

xi

) >

0 for xi

<

a

ˇ

_i(1),

ˆ

f_i(1)

(

xi

) <

0 for xi

>

a

ˆ

(_i1). Hence for any small

ε

1

>

0 there exist a T1

>

T0 such that any solution x

(

t

,

x0) starting from a point

x0

∈ [

0

,

A1

] × [

0

,

A2

]

falls into the region

Ω

(1)

:=



a

ˇ

(₁1)

−

ε

1,a

ˆ

(₁1)

+

ε

1

×



a

ˇ

(₂1)

−

ε

1,a

ˆ

(₂1)

+

ε

1

Fig. 4.Ω(0)_andΩ(1)_{, for one-peak case.}

for t



T1; cf. Fig. 4. We can proceed in a similar manner to define successively

ˆ

fi(k) and

ˇ

f

(k) i , k



2, by

ˆ

f₁(k)

(

x1)

= −

μ

x1

+



α

1 xn 1 kn₁

+

xn₁

+

σ

1 Smax 1

(

Tk−1)

ρ

1

+

Smax1

(

Tk−1)



·

1 1

+ (ˇ

a(₂k−1)

−

ε

k−1)/

γ

2

+ β

1,

ˇ

f₁(k)

(

x1)

= −

μ

x1

+



α

1 xn 1 kn 1

+

xn1

+

σ

1 Smin 1

(

Tk−1)

ρ

1

+

Smin₁

(

Tk−1

)



·

1 1

+ (ˆ

a(₂k−1)

+

ε

k−1)/

γ

2

+ β

1,

ˆ

f₂(k)

(

x2)

= −

μ

x2

+



α

2 xn₂ kn 2

+

xn2

+

σ

2 Smax₂

(

Tk−1

)

ρ

2

+

Smax2

(

Tk−1)



·

1 1

+ (ˇ

a(₁k−1)

−

ε

k−1)/

γ

1

+ β

2,

ˇ

f₂(k)

(

x2)

= −

μ

x2

+



α

2 xn₂ kn 2

+

xn2

+

σ

2 Smin₂

(

Tk−1

)

ρ

2

+

Smin2

(

Tk−1)



·

1 1

+ (ˆ

a(₁k−1)

+

ε

k−1)/

γ

1

+ β

2 and their zerosa

ˆ

(_ik)

,

a

ˇ

(_ik), i.e.,

ˆ

f_i(k+1)



a

ˆ

(_ik)



=

0

,

ˇ

f_i(k+1)



a

ˇ

_i(k)



=

0

.

(6.5) We may clearly assume that

ε

k

→

0 and Tk

→ ∞

as k

→ ∞

.

We can then prove that for any small

ε

k

>

0 there exists a Tk such that any solution x

(

t

,

x0) starting from a point x0

∈ [

0

,

A1

] × [

0

,

A2

]

falls into the region

Ω

(k)

:= [ˇ

a1(k)

−

ε

k

,

a

ˆ

(1k)

+

ε

k

] × [ˇ

a(2k)

−

ε

k

,

a

ˆ

(2k)

+

ε

k

] ⊂ Ω

(k−1)for t



Tk.

We shall need the following conditions:

(

α

2

+

σ

2)

γ

1

<

μ

−

α

1n

˜

k1

−

σ

1

ρ

1

,

(

α

1

+

σ

1)

γ

2

<

μ

−

α

2n

˜

k2

−

σ

2

ρ

2

.

(6.6)

Lemma 6.1. Under the conditions (6.6) and either (6.1), (6.2) or (6.3), the intersection



∞_k=1

Ω

(k)_{consists of}

750 A. Friedman et al. / J. Differential Equations 247 (2009) 736–769

Proof. Note that for each i

=

1

,

2,

{ˇ

a(_ik)

−

ε

k

}

is an increasing sequence,

{ˆ

a(_ik)

+

ε

k

}

is a decreasing sequence,a

ˇ

(_ik)

−

ε

k

<

a

ˆ

_i(k)

+

ε

kfor each k, and

ε

k

→

0 as k

→ ∞

. Hence

ˇ

a∗_i

=

lim k→∞a

ˇ

(k) i and a

ˆ

∗i

=

_klim_→∞a

ˆ

(k)

i exist

,

and a

ˇ

∗i

 ˆ

a∗i for i

=

1

,

2

.

Assuming that



∞_k=1

Ω

(k) _{is not a single point so thata}

_ˆ

∗

i

>

a

ˇ

∗i for either i

=

1 or i

=

2 (or both), we shall derive a contradiction.

By passing to the limit in (6.5) we get

−

μ

a

ˇ

∗₁

+



α

1

(

a

ˇ

∗₁

)

n kn₁

+ (ˇ

a∗₁

)

n

+

σ

1

ˇ

S1

ρ

1

+ ˇ

S1



·

1 1

+ ˆ

a∗₂

/

γ

2

+ β

1

=

0

,

(6.7)

−

μ

a

ˆ

∗₂

+



α

2

(

a

ˆ

∗₂

)

n kn₂

+ (ˆ

a∗₂

)

n

+

σ

2

ˆ

S2

ρ

2

+ ˆ

S2



·

1 1

+ ˇ

a∗₁

/

γ

1

+ β

2

=

0

,

(6.8)

−

μ

a

ˆ

∗₁

+



α

1

(

a

ˆ

∗₁

)

n kn 1

+ (ˆ

a∗1

)

n

+

σ

1

ˆ

S1

ρ

1

+ ˆ

S1



·

1 1

+ ˇ

a∗₂

/

γ

2

+ β

1

=

0

,

(6.9)

−

μ

a

ˇ

∗₂

+



α

2

(

a

ˇ

∗₂

)

n kn₂

+ (ˇ

a∗₂

)

n

+

σ

2

ˇ

S2

ρ

2

+ ˇ

S2



·

1 1

+ ˆ

a∗₁

/

γ

1

+ β

2

=

0

,

(6.10) where

ˆ

Si

=

lim t→∞S max i

(

t

),

S

ˇ

i

=

lim t→∞S min i

(

t

),

and

ˆ

S1

 ˆ

a∗1

+ ¯

C1, S

ˇ

1

 ˇ

a∗1

+ ¯

C1, (6.11)

ˆ

S2

 ˆ

a∗2

+ ¯

C2, S

ˇ

2

 ˇ

a∗2

+ ¯

C2, (6.12) with

¯

Ci

=

lim t→∞Ci

(

t

)

e −G(t)

_/

_N0.

Taking the difference of (6.7), (6.9) we obtain

μ



a

ˆ

∗₁

− ˇ

a∗₁



−

α

1



(

a

ˆ

∗₁

)

n kn₁

+ (ˆ

a∗₁

)

n

−

(

a

ˇ

∗₁

)

n kn₁

+ (ˇ

a∗₁

)

n



·

1 1

+ ˇ

a∗₂

/

γ

2

=



α

1

(

a

ˇ

∗₁

)

n kn₁

+ (ˇ

a∗₁

)

n

+

σ

1

ˇ

S1

ρ

1

+ ˇ

S1



·



1 1

+ ˇ

a∗₂

/

γ

2

−

1 1

+ ˆ

a∗₂

/

γ

2



+

σ

1



_ˆ

S1

ρ

1

+ ˆ

S1

−

S

ˇ

1

ρ

1

+ ˇ

S1



·

1 1

+ ˇ

a∗₂

/

γ

2

.

(6.13)

Thus, by the mean value theorem and the estimates (6.11) for

ˆ

S1,S

ˇ

1,

ˆ

a∗₁

− ˇ

a∗₁

·

μ

−

α

1n

˜

k1



(

α

1

+

σ

1)

γ

2

ˇ

a∗₂

− ˆ

a∗₂

+

σ

1

ρ

1

ˆ

a∗₁

− ˇ

a∗₁

,

or

ˆ

a∗₁

− ˇ

a∗₁

·



μ

−

α

1n

˜

k1

−

σ

1

ρ

1





(

α

1

+

σ

1)

γ

2

ˇ

a∗₂

− ˆ

a∗₂

.

(6.14) Similarly, from (6.8), (6.10), (6.12) we obtain

ˇ

a∗₂

− ˆ

a∗₂

·



μ

−

α

2n

˜

k2

−

σ

2

ρ

2





(

α

2

+

σ

2)

γ

1

ˆ

a∗₁

− ˇ

a∗₁

.

(6.15) Assuming that the LHS of (6.14) and (6.15) are positive, these two inequalities yield



μ

−

α

1n

˜

k1

−

σ

1

ρ

1



·



μ

−

α

2n

˜

k2

−

σ

2

ρ

2



<

(

α

2

+

σ

2)

γ

1

·

(

α

1

+

σ

1)

γ

2

,

(6.16)

which is a contradiction to (6.6). We thus conclude that a

ˇ

∗_i

= ˆ

a∗_i for i

=

1

,

2, which proves the

lemma.

2

From Lemma 6.1 it follows that the limit

(

a

¯

1,a

¯

2)of

Ω

(k)(as k

→ ∞

) satisfies the equations

−

μ

a1

+



α

1 an 1 kn₁

+

an₁

+

σ

1 a1

+ ¯

C1

ρ

1

+

a1



·

1 1

+

a2

/

γ

2

+ β

1

=

0

,

(6.17)

−

μ

a2

+



α

2 an 2 kn₂

+

an₂

+

σ

2 a2

+ ¯

C2

ρ

2

+

a2



·

1 1

+

a1/

γ

1

+ β

2

=

0

,

(6.18)

and the solution is unique. We have thus proved:

Theorem 6.2. If

(

6

.

6

)

and one of the conditions (6.1), (6.2), or (6.3) hold, then the solution

ψ

of (2.10), (3.1), with fi, Sidefined by (2.5), (2.6), (2.11), satisfies:

lim

t→∞

ψ(

t

,

x1,x2)

=

N0δ(a1,a2)

,

(6.19)

where

δ

(a1,a2)is the Dirac measure at point

(

a1,a2)which is uniquely determined from (6.17)–(6.18), and the

convergence in (6.19) is in the sense of convergence in measure as defined in (2.12).

7. Asymptotic two-peak solutions

Analogously to the case of Theorem 4.1(iii) we assume that condition (B1) holds

,

ˆ

f1



ˆ

pm₁



<

0

,

and

ˇ

f1



ˇ

pM₁



>

0

,

(7.1)

either condition (M2) holds

,

or (B2) and

ˆ

f2



ˆ

pM₂



<

0 hold

.

(7.2) Leta

ˆ

(₁0)

, ˆ

b(₁0)

,

c

ˆ

₁(0)(respectively a

ˇ

₁(0)

, ˇ

b(₁0)

,

c

ˇ

(₁0)) be the zeros of

ˆ

f₁(0) (respectively

ˇ

f₁(0)), anda

ˆ

(₂0)

,

a

ˇ

(₂0)be the zeros of

ˆ

f₂(0)

, ˇ

f₂(0), respectively; cf. Fig. 5.

Then, by (7.1), (7.2) and (5.2)–(5.4), for any small

ε

0

>

0 there exists a T0

>

0 such that any solution

x

(

t

,

x0)starting from a point x0

∈ [

0

,

A1

] × [

0

,

A2

] \

K(0)falls into the region

Ω

(0)

= Ω

_l(0)

∪ Ω

u(0)