Asymptotic limit in a cell differentiation model with consideration of transcription

Contents lists available atSciVerse ScienceDirect

Journal of Differential Equations

www.elsevier.com/locate/jde

Asymptotic limit in a cell differentiation model with

consideration of transcription

Avner Friedman

a

, Chiu-Yen Kao

a

,

b

, Chih-Wen Shih

c

,

∗

a_{Mathematical Biosciences Institute, Department of Mathematics, The Ohio State University, OH 43202, United States} b_{Department of Mathematics and Computer Science, Claremont McKenna College, Claremont, CA 91711, United States} c_{Department of Applied Mathematics, National Chiao Tung University, Hsinchu, 300, Taiwan}

a r t i c l e

i n f o

a b s t r a c t

Article history:

Received 25 August 2011 Revised 21 January 2012 Available online 20 February 2012

Keywords: Cell differentiation Th1/Th2 cells Conservation law Multistationary Integro-differential equation Transcription factors mRNA

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). These factors

are translated by the mRNA whose levels of expression, y1and y2,

depend, respectively, on x1 and x2 in a nonlinear nonlocal way.

The population density of T cells,φ(t,x1,x2,y1,y2), satisfies a

hy-perbolic conservation law with coefficients depending nonlinearly and nonlocally on(t,x1,x2,y1,y2), while the xi, yi satisfy a

sys-tem of ordinary differential equations. We study the long time behavior ofφand show, under some conditions on the parameters of the system of differential equations, that the gene expressions in the T-cell population aggregate at one, two or four points, which connect to various cell differentiation scenarios.

©2012 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. For example, the T lymphocytes of the immune system, upon maturation, differentiate into either Th1 or Th2 cells. These cells are different by the repertoire of chemokines they produce. Th1 cells secrete IFNγ needed to combat intracellular pathogens and, if abnormal, are associated with inflammatory and autoimmune diseases. Th2 cells secrete cytokines that activate B cells to produce

*

Corresponding author. Fax: +886 3 5724679.

E-mail address:cwshih@math.nctu.edu.tw(C.-W. Shih).

0022-0396/$ – see front matter ©2012 Elsevier Inc. All rights reserved. doi:10.1016/j.jde.2012.02.006

antibodies against extracellular pathogens and, if abnormal, are associated with asthma and other allergies.

The variables of primary interest in a quantitative description of gene expression are the number of mRNA copies of a given gene and the number of transcription factors (proteins). The mRNA are translated into proteins, and transcription factors promote the mRNA transcription by genes. Hence in order to determine quantitatively the cellular concentration of mRNA and protein, we need a math-ematical model that connects these two concentrations. In terms of the balance equations, these concentrations are governed by

d

(

mRNA

)

dt

=

vtranscription

−

vmRNA degradation

,

(1.1) d

(

protein

)

dt

=

vtranslation

−

vprotein degradation

,

(1.2) where the v’s are the rates of transcription, translation, and degradation as indicated; cf. [9].

In the case of T cell differentiation, the decision to which cell type to differentiate, Th1 or Th2, depends on proteins x1and x2, and their mRNA y1and y2, where x1is the concentration of

transcrip-tion factor T-bet and x2is the concentration of transcription factor GATA-3; yiis the concentration of the mRNA which translates into xi. By (1.2), we then have

dxi

dt

=

viyi

−

τ

ixi

=:

gi

,

(1.3)

where vi

,

τ

i are constants. On the other hand, the rate of change dyi

/

dt is far more complex, since vtranscription depends on intrinsic signals from all the T cells and on extrinsic signals by IL4 and IL12.

Yates et al. [10] introduced the following model for the rate of the transcription of xi:

vtranscription

=



α

i xn_i kn_i

+

xn_i

+

σ

i Si

ρ

i

+

Si



·

1 1

+

xj

/

γ

j

+ β

i

,

where

α

i, ki,

σ

i,

ρ

i,

γ

i,

β

i are constants and j

=

2 if i

=

1, j

=

1 if i

=

2. Here Si is the combined intrinsic/extrinsic signal, and xj inhibits xi ( j

=

i); the autocatalytic process, given by

α

ixn_i

/(

kn_i

+

xn_i

)

, is modeled by Hill’s dynamics with exponents n



2. The first balance equation (1.1) then becomes

dyi dt

= −

μ

iyi

+



α

i xn_i kn_i

+

xn_i

+

σ

i Si

ρ

i

+

Si



·

1 1

+

xj

/

γ

j

+ β

i

=:

fi

,

(1.4) for

(

i

,

j

)

= (

1

,

2

)

and

(

i

,

j

)

= (

2

,

1

)

.

Introducing the population density of cells with concentration

(

x1

,

x2

,

y1

,

y2

)

at time t,

φ (

t

,

x1

,

x2

,

y1

,

y2

)

, the mass conservation law then yields

∂φ

∂

t

+

2



i=1

∂

∂

xi

(

gi

φ)

+

2



i=1

∂

∂

yi

(

fi

φ)

=

g∗

φ,

(1.5)

where g∗is the growth factor.

For a healthy normal individual in homeostasis, the expressions of mRNA/T-bet and mRNA/GATA-3 are at stationary levels, and, at intermediate times, Th0 does not differentiate into Th1 or Th2. How-ever, when a strong signal Si is generated in response to pathogens, the Th cells differentiate into either Th1 or Th2, but usually not both. In the present model, a single cell with high (low) con-centration of T-bet

(

x1

)

and low (high) concentration of GATA-3

(

x2

)

corresponds to the polarization

and mRNA/GATA-3 may aggregate at one or several points y1

/

x1 and y2

/

x2, respectively. For those

cells whose expressions aggregate at the point with low-x1 and low-x2, cell differentiation does not

occur, while the cells whose expressions aggregate at the point with high (low) x1 and low (high) x2,

differentiate into Th1 (Th2). The model parameters (although similar to those in Yates et al. [10]) are not experimentally known; hence our aim is to show that with a specific choice of parameters, the present model illustrates the main biological phenomena on cell differentiation. The fact that we end up with 1, 2, or 4 limit aggregations may not be biologically significant; the model with other param-eters may end up with different number of aggregation points. What is important is that although there may be a number of limit points, only points with significant contrast of protein concentrations, i.e., x2



x1 or x1



x2, indicate cell differentiation. In a recent paper, we studied the asymptotic

behavior of the reduced system (with yi

≡

xi) dxi dt

= −

μ

ixi

+



α

i xn_i kn_i

+

xn_i

+

σ

i Si

ρ

i

+

Si



·

1 1

+

xj

/

γ

j

+ β

i

=: ˜

fi

,

(1.6)

for

(

i

,

j

)

= (

1

,

2

)

and

(

i

,

j

)

= (

2

,

1

)

, with the conservation law

∂φ

∂

t

+

2



i=1

∂

∂

xi

( ˜

fi

φ)

=

g∗

φ,

(1.7)

where

φ

= φ(

t

,

x1

,

x2

)

, and proved under some conditions on the parameters of (1.6) that

φ (

t

,

x1

,

x2

)

converges to a linear combination of one, two, or four Dirac functions, as t

→ ∞

.

In the present paper, we consider the more general model (1.3), (1.4), (1.5) and establish similar asymptotic behaviors for the population density of T cells,

φ (

t

,

x1

,

x2

,

y1

,

y2

)

. The proof,

however, involves a far deeper analysis than the analysis we used in the reduced case of (1.6) and (1.7).

2. The mathematical model

Denote x1 and x2as the concentrations of transcription factors T-bet and GATA-3, respectively, and

by y1and y2 their respective mRNA concentrations. By combining the models of Yates et al. [10] and

Mariani et al. [9] (see also [1]), we obtain the following system:

⎧

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎨

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎩

dy1 dt

= −

μ

1y1

+



α

1 xn 1 kn₁

+

xn₁

+

σ

1 S1

ρ

1

+

S1



1 1

+

x2

/

γ

2

+ β

1

=:

f1

(

x1

,

x2

,

y1

,

S1

),

dy2 dt

= −

μ

2y2

+



α

2 xn₂ kn₂

+

xn₂

+

σ

2 S2

ρ

2

+

S2



1 1

+

x1

/

γ

1

+ β

2

=:

f2

(

x1

,

x2

,

y2

,

S2

),

dx1 dt

=

ν

1y1

−

τ

1x1

=:

g1

(

x1

,

y1

),

dx2 dt

=

ν

2y2

−

τ

2x2

=:

g2

(

x2

,

y2

).

(2.1)

The first term on the right-hand side of the yi-equation represents the rate of mRNA degradation, and

β

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

α

i xn_i kn_i

+

xn_i

where n



2 is the Hill exponent that tunes the sharpness of the activation switch. The contribution of combined cytokine signaling to the rate of growth in yiis given by the term

σ

i Si

ρ

i

+

Si

.

The cross-inhibition between y1 and y2 occurs at both the autoactivation level and the cytokine

(membrane) signaling level, and is represented by the factors

1 1

+

xj

/

γ

j

.

The parameter

γ

jis the value of xjat which the ratio of production of yi(i

=

j), due to the combined autoactivation and cytokine signaling, is halved.

We denote by

φ (

t

,

x1

,

x2

,

y1

,

y2

)

the population density of T cells with protein concentration

(

x1

,

x2

)

and mRNA concentration

(

y1

,

y2

)

at time t. Then the total levels of expression of T-bet and

GATA-3, at time t in the cell population are given by

xi

˜φ(

t

,

x1

,

x2

)

dx1dx2

,

for i

=

1 and i

=

2, respectively, where

˜φ(

t

,

x1

,

x2

)

=

φ (

t

,

x1

,

x2

,

y1

,

y2

)

dy1dy2. If we denote by

Ei

(

t

)

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

(

t

)

=

Ei

(

t

)

+

xi

˜φ(

t

,

x1

,

x2

)

dx1dx2

˜

φ(

t

,

x1

,

x2

)

dx1dx2

,

i

=

1

,

2

.

(2.2)

Here, a normalization by total cell numbers is adopted in order 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

(

g1

φ)

+

∂

∂

x2

(

g2

φ)

+

∂

∂

y1

(

f1

φ)

+

∂

∂

y2

(

f2

φ)

=

g∗

φ,

(2.3)

where g∗is a growth factor. Note that (2.3) is associated with the velocity field described by

dxi

(

t

)

dt

=

gi



xi

(

t

),

yi

(

t

)



,

(2.4) dyi

(

t

)

dt

=

fi



t

,

xi

(

t

),

yi

(

t

),

Si

(

t

)



,

(2.5)

where fiand giare defined in (2.1). We shall consider system (2.4)–(2.5) in the rectangular region

Ω

= {

0



x1



B1

,

0



x2



B2

,

0



y1



A1

,

0



y2



A2

}

where Bi

=

ν

i

τ

i Ai

,

i

=

1

,

2

,

(2.6) Ai

=

α

i

+

σ

i

+ β

i

μ

i

,

i

=

1

,

2

,

(2.7)

and set

˜

Ω

= {

0



x1



B1

,

0



x2



B2

}.

Then

Ω

is a positively invariant and an attracting set for (2.4)–(2.5). Therefore, in order to solve (2.3) for

(

x1

,

x2

,

y1

,

y2

)

in

Ω

, we need to assign both initial and boundary conditions to

φ

:

φ (

0

,

x1

,

x2

,

y1

,

y2

)

= φ

0

(

x1

,

x2

,

y1

,

y2

)

in

Ω,

(2.8)

φ (

t

,

x1

,

x2

,

y1

,

y2

)

|

∂Ω

=

0 for all t

>

0

.

(2.9) Assuming, for simplicity, that g∗

=

g∗

(

t

)

, and setting

G

(

t

)

=

t

0 g∗

(

s

)

ds

,

N0

=

Ω

φ

0

(

x1

,

x2

,

y1

,

y2

),

(2.10)

ψ (

t

,

x1

,

x2

,

y1

,

y2

)

=

e−G(t)

φ (

t

,

x1

,

x2

,

y1

,

y2

),

(2.11)

˜ψ(

t

,

x1

,

x2

)

=

ψ (

t

,

x1

,

x2

,

y1

,

y2

)

dy1dy2

,

we can replace (2.3) by the simpler equation

∂ψ

∂

t

+

∂

∂

x1

(

g1

ψ )

+

∂

∂

x2

(

g2

ψ )

+

∂

∂

y1

(

f1

ψ )

+

∂

∂

y2

(

f2

ψ )

=

0

,

(2.12)

and rewrite Si

(

t

)

in the form

Si

(

t

)

=

Ei

(

t

)

e−G(t) N0

+

xi

˜ψ(

t

,

x1

,

x2

)

dx1dx2 N0

,

(2.13)

where N0is the initial total population, and the integral in (2.13) is taken over

Ω.

˜

Let

Φ(

t

,

x1

,

x2

,

y1

,

y2

)

denote the solution map (flow map) of (2.4)–(2.5) and set

Ω(

t

)

= Φ(

t

, Ω)

.

Integrating the transport equation (2.12) over

Ω(

t

)

, we find that d

dt

Ω(t)

ψ (

t

,

x1

,

x2

,

y1

,

y2

)

dx1dx2dy1dy2

=

0

.

Furthermore, if

Ω(

t

)

→ (¯

a1

,

a

¯

2

,

a

¯

3

,

a

¯

4

)

as t

→ ∞

, then for any continuous function h

(

x1

,

x2

,

y1

,

y2

)

,

Ω

h

(

x1

,

x2

,

y1

,

y2

)ψ (

t

,

x1

,

x2

,

y1

,

y2

)

dx1dx2dy1dy2

→

h

(

a

¯

1

,

a

¯

2

,

a

¯

3

,

a

¯

4

)

N0 as t

→ ∞,

i.e.,

ψ (

t

,

x1

,

x2

,

y1

,

y2

)

→

N0

δ

(a¯1,a¯2,a¯3,a¯4) in measure, as t

→ ∞.

In the subsequent sections, we study the asymptotic behavior of the solutions of (2.4)–(2.5) in con-junction with the behavior of

Ω(

t

)

. Similarly to [2], one can prove that the system (2.3), (2.8), (2.9)

has a unique solution for all t



0. Hence we shall focus here only on the asymptotic behavior of the solution. We shall prove that

Ω(

t

)

converges to one, two, or four points, as t

→ ∞

, depending on the parameters of the dynamical system (2.4)–(2.5). The asymptotic study of dynamical system (2.4)–(2.5) will require a far deeper analysis than that developed for Eqs. (1.6)–(1.7) in [2].

3. Upper and lower dynamics

The system (2.1) can be written as a system of two second-order equations, d2_x 1 dt2

+ (

τ

1

+

μ

1

)

dx1 dt

=

h1



x1

,

x2

,

S1

(

t

)



,

(3.1) d2x2 dt2

+ (

τ

2

+

μ

2

)

dx2 dt

=

h2



x1

,

x2

,

S2

(

t

)



,

(3.2) where h1



x1

,

x2

,

S1

(

t

)



= −

μ

1

τ

1x1

+

ν

1



α

1 xn 1 kn₁

+

xn₁

+

σ

1 S1

(

t

)

ρ

1

+

S1

(

t

)



1 1

+

x2

/

γ

2

+

ν

1

β

1

,

h2



x1

,

x2

,

S2

(

t

)



= −

μ

2

τ

2x2

+

ν

2



α

2 xn₂ kn₂

+

xn₂

+

σ

2 S2

(

t

)

ρ

2

+

S2

(

t

)



1 1

+

x1

/

γ

1

+

ν

2

β

2

.

We introduce the upper boundsh

ˆ

i for the functions hi:

ˆ

hi

(

xi

)

= −

μ

i

τ

ixi

+

ν

i



α

i xn_i kn_i

+

xn_i

+

σ

i

ˆ

Si

ρ

i

+ ˆ

Si



+

ν

i

β

i for 0



xi

<

∞,

i

=

1

,

2

,

(3.3) where S

ˆ

i

=

supt>0Si

(

t

)

, and lower boundsh

ˇ

ifor hi:

ˇ

h1

(

x1

)

= −

μ

1

τ

1x1

+

ν

1



α

1 xn₁ kn₁

+

xn₁

+

σ

1

ˇ

S1

ρ

1

+ ˇ

S1



·

1 1

+

B2

/

γ

2

+

ν

1

β

1

,

(3.4)

ˇ

h2

(

x2

)

= −

μ

2

τ

2x2

+

ν

2



α

2 xn₂ kn 2

+

xn2

+

σ

2

ˇ

S2

ρ

2

+ ˇ

S2



·

1 1

+

B1

/

γ

1

+

ν

2

β

2

,

(3.5)

where S

ˇ

i

=

inft>0Si

(

t

)

. Clearly,

ˆ

hi

(

0

) >

0

,

h

ˆ

i

(

0

) <

0

,

h

ˆ

i

(

xi

) <

0

,

ˇ

hi

(

0

) >

0

,

h

ˇ

i

(

0

) <

0

,

h

ˇ

i

(

xi

) <

0

,

for Bi



xi

<

∞

. Also

ˆ

h_i

(

xi

)

= −

μ

i

τ

i

+

ν

i

α

i nkn_ixn_i−1

(

kn_i

+

xn_i

)

2

,

and, as easily verified, the maximum of the last term is attained at the point

˜ξ

i

=

ki



n

−

1 n

+

1



1/n

,

and h

ˆ

_i

( ˜ξ

i

)

= −

μ

i

τ

i

+

ν

i

α

in

˜

ki (3.6)

where

˜

n

= (

n

+

1

)

1+1/n

(

n

−

1

)

1−1/n

/

4n

.

The maximum ofh

ˇ

_i is also attained at the same point

˜ξ

iwith

ˇ

h_i

( ˜ξ

i

)

= −

μ

i

τ

i

+

ν

i

α

in

˜

ki

·

1 1

+

Bj

/

γ

j for

(

i

,

j

)

= (

1

,

2

)

and

(

i

,

j

)

= (

2

,

1

)

. Clearly,h

ˇ

_i

(ξ ) < ˆ

h_i

(ξ )

for all

ξ

.

The systems d2x

ˆ

i dt2

+ (

τ

i

+

μ

i

)

d

ˆ

xi dt

= ˆ

hi

(

x

ˆ

i

) (

i

=

1

,

2

),

(3.7) d2x

ˇ

i dt2

+ (

τ

i

+

μ

i

)

d

ˇ

xi dt

= ˇ

hi

(

x

ˇ

i

) (

i

=

1

,

2

),

(3.8)

will be used to provide the upper and lower bounds for the dynamics of (3.1)–(3.2). It will be convenient to use a change of variables

(

x1

,

y1

,

x2

,

y2

)

↔ (

x1

,

v1

,

x2

,

v2

)

where

v1

=

ν

1y1

−

τ

1x1

,

v2

=

ν

2y2

−

τ

2x2

so that the system (2.1) can be rewritten in the form dxi

dt

=

vi

,

(3.9)

dvi

dt

= −(

τ

i

+

μ

i

)

vi

+

hi

(

x1

,

x2

,

Si

),

(3.10)

i

=

1

,

2, in the transformed region

Ω

∗

= {

0



x1



B1

,

−

τ

1x1



v1



ν

1A1

−

τ

1x1

,

0



x2



B2

,

−

τ

2x2



v2



ν

2A2

−

τ

2x2

}.

Notice that

Ω

∗remains positively invariant under (3.9)–(3.10).

We need several lemmas to study the asymptotic behavior of (3.9)–(3.10). The first one deals with a system

du

dt

=

v

,

(3.11)

dv

dt

= −δ

v

+

q

(

u

),

(3.12)

where

δ

is a positive constant and q is a continuously differentiable function on

[

0

,

∞)

. We shall consider (3.11)–(3.12) on a region D

⊆ [

0

,

∞) × R

, which is positively invariant under the flow

Ψ

t generated from the system. Let B

(

u

,

v

)

⊂ R

2 be an open disc with center

(

u

,

v

)

and radius



, and K be a compact set in D.

Lemma 3.1. Assume that limu→∞q

(

u

)

= −∞

. Then the following holds.

(i) Every solution of (3.11)–(3.12) tends to the set

{(

u

,

0

)

∈

D: q

(

u

)

=

0

}

, as t

→ ∞

; if, in addition, the set of zeros for q is finite, then each solution of (3.11)–(3.12) tends to a single point in set

{(

u

,

0

)

∈

D: q

(

u

)

=

0

}

, as t

→ ∞

.

(ii) If q has a unique zero a with q

(

a

) <

0 and

(

a

,

0

)

∈

D, then

A := {(

a

,

0

)

}

is the global attractor; thus, for any small



>

0, there exists a T such that

Ψ

t

(

K

)

⊂

B

(

a

,

0

)

, for all t



T .

(iii) If q has exactly three zeros a

,

b

,

c with a

<

b

<

c, q

(

a

) <

0, q

(

b

) >

0, q

(

c

) <

0, and

(

a

,

0

), (

b

,

0

),

(

c

,

0

)

∈

D, then

A :=

Wu

(

b

,

0

)

∪ {(

a

,

0

)

} ∪ {(

c

,

0

)

}

is the global attractor, where Wu

(

b

,

0

)

is the unstable manifold of

(

b

,

0

)

. Moreover, for any small



>

0, there exists a T

>

0 such that

|Ψ

t

(

K

\

Ws



(

b

,

0

))

− {(

a

,

0

), (

c

,

0

)

}| <



, for all t



T , where Ws

(

b

,

0

)

:= {(

u

,

v

)

:

|(

u

,

v

)

−

Ws

(

b

,

0

)

| <



}

is the

-neighborhood of W

s

₍

_b

_,

₀

₎

_.

Proof. (i) Consider the Lyapunov function

V

(

u

,

v

)

=

1 2v 2

₋

u

0 q

(

s

)

ds

.

Then

˙

V

(

u

,

v

)

=

v

·

−δ

v

+

q

(

u

)



−

q

(

u

)

·

v

= −δ

v2



0

,

(3.13)

andV

˙

=

0 if and only if v

=

0. All solutions are bounded in forward time due to limu→∞q

(

u

)

= −∞

. By LaSalle’s invariance principle [4,5], every solution of (3.11)–(3.12) tends to the maximal invariant set in



(

u

,

v

)

: V

˙

(

u

,

v

)

=

0



=



(

u

,

0

)



which is the set



(

u

,

0

)

: q

(

u

)

=

0



,

as t

→ ∞.

Since the

ω-limit set of an orbit is connected, if q has a finite number of zeros, then the

ω-limit set for an orbit of (3.11)–(3.12) is a single point

(

u

,

0

)

, where u is a zero of q.

(ii) If q has a unique zero a with q

(

a

) <

0, then

(

a

,

0

)

is a sink. From (3.13), it follows that

{(

a

,

0

)

}

is the global attractor for (3.11)–(3.12). The assertion about

Ψ

t

(

K

)

⊂

B

(

a

,

0

)

for t



T follows from [5,8].

(iii) If q has exactly three zeros a

,

b

,

c with a

<

b

<

c and q

(

a

) <

0, q

(

b

) >

0, q

(

c

) <

0, then

(

a

,

0

), (

c

,

0

)

are both sinks, and

(

b

,

0

)

is a saddle, for system (3.11)–(3.12). By (3.13), the level curve analysis, and Poincare–Bendixson Theorem, the unstable manifold Wu

(

b

,

0

)

for

(

b

,

0

)

consists of het-eroclinic orbits connecting

(

b

,

0

)

with

(

a

,

0

)

and with

(

c

,

0

)

respectively; cf. Fig. 1. It follows that

A :=

Wu

(

b

,

0

)

∪ {(

a

,

0

)

} ∪ {(

c

,

0

)

}

is the global attractor for (3.11)–(3.12); cf. [6, p. 395]. Therefore, for any



>

0 there exists a T

>

0, so that

Ψ

t

(

K

)

falls within a distance



>

0 from

A

, for all t



T . Moreover, for every point

(

u

,

v

)

in compact set K

\

Ws



(

b

,

0

)

,

Φ

t

(

u

,

v

)

approaches

(

a

,

0

)

or

(

c

,

0

)

, as t tends to infinity. By the continuity with respect to initial condition and the compactness of K , there exists a T

>

0 such that

Ψ

t

(

u

,

v

)

∈

Bε

(

a

,

0

)

or Bε

(

c

,

0

)

, for all

(

u

,

v

)

∈

K

\

Ws

(

b

,

0

)

, for all t



T .

2

Fig. 1. Stable manifold Ws

(b,0)(indicated as blue lines) and unstable manifold Wu

(b,0)(indicated as red lines) of the saddle point(b,0)of (3.11)–(3.12). The green dots allocate three equilibria(a,0), (b,0), (c,0). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Lemma 3.1 applies to system (3.7)–(3.8) on the domain Di, i

=

1

,

2, respectively, where Di

:=



(

xi

,

vi

)

: 0



xi



Bi

,

−

τ

ixi



vi



ν

iAi



⊂ R

2

_.

In particular, every solution x

ˆ

i

(

t

)

(resp., x

ˇ

i

(

t

)

) of (3.7) (resp., (3.8)) tends to the set of zeros of h

ˆ

i (resp.,h

ˇ

i), as t

→ ∞

, i

=

1

,

2.

Lemma 3.2. Consider the non-autonomous equation

d2_z

dt2

+ (

τ

+

μ

)

dz

dt

+

a

(

t

)

z

=

f

(

t

),

0

<

t

<

∞,

(3.14)

where a

(

t

)



α

,

f

(

t

)



0 for 0



t

<

∞

, and f

(

0

) >

0, 0

<

4α

< (

τ

+

μ

)

2. If z

(

0

)

= (

dz

/

dt

)(

0

)

=

0, then z

(

t

)



0 for all t



0.

Proof. We rewrite (3.14) in the form

d2_z dt2

+ (

τ

+

μ

)

dz dt

+

α

z

=

f

(

t

)

+



α

−

a

(

t

)



z

(

t

).

(3.15)

Eq. (3.15) has two linearly independent homogeneous solutions eλ1t_{, e}λ2t _where

λ

1,2

=

−(

τ

+

μ

)

±



(

τ

+

μ

)

2

₋

₄

_α

2

and, by assumption,

λ

1

< λ

2

<

0. Since z

(

0

)

= (

dz

/

dt

)(

0

)

=

0, we can represent z, by the variation of

constant formula, in the form

z

(

t

)

=

t

0 eλ2(t−s)

−

_eλ1(t−s)

(λ

2

− λ

1

)

f

(

s

)

+



α

−

a

(

s

)



z

(

s

)



ds

;

(3.16)

indeed observe that the right-hand side vanishes at t

=

0 together with its first derivative.

Since f

(

0

)

− (

α

−

a

(

0

))

z

(

0

)

=

f

(

0

) >

0, z

(

t

) >

0 for small t. We claim that z

(

t

) >

0 for all t

>

0. Indeed, otherwise there exists a smallest time t

=

t0 such that z

(

t

) >

0 if 0

<

t

<

t0 and z

(

t0

)

=

0.

But since

α

−

a

(

t

)



0, we have f

(

t

)

+ (

α

−

a

(

t

))

z

(

t

) >

0 for 0

<

t

<

t0 and from (3.16) we obtain

Remark 3.1. By approximation, the lemma remains true if f

(

0

)

=

0 and 4α

 (

τ

+

μ

)

2.

Lemma 3.2 can be used to compare solutions

(

x1

(

t

),

x2

(

t

))

of (3.1)–(3.2) with solutionsx

ˆ

1

(

t

),

x

ˆ

2

(

t

)

of (3.7), andx

ˇ

1

(

t

),

ˇ

x2

(

t

)

of (3.8), provided they have the same initial conditions.

Lemma 3.3. Let

(

x1

(

t

),

x2

(

t

))

be a solution of (3.1), (3.2). Suppose

min

ˇ

h_i

(

η

)

:

η

∈ [

0

,

Bi

]





μ

i

τ

i

,

i

=

1 or i

=

2

.

(i) If a solution

ˆ

xi

(

t

)

of (3.7) satisfies

ˆ

xi

(

0

)

=

xi

(

0

),

(

dx

ˆ

i

/

dt

)(

0

)

= (

dxi

/

dt

)(

0

)

then

ˆ

xi

(

t

)



xi

(

t

)

for all t

>

0

.

(ii) If a solution

ˇ

xi

(

t

)

of (3.8) satisfies

ˇ

xi

(

0

)

=

xi

(

0

),

(

dx

ˇ

i

/

dt

)(

0

)

= (

dxi

/

dt

)(

0

)

then

ˇ

xi

(

t

)



xi

(

t

)

for all t

>

0

.

(iii) If solutions

ˆ

xi

(

t

)

,x

ˇ

i

(

t

)

of (3.7)–(3.8) satisfy

ˆ

xi

(

0

)

= ˇ

xi

(

0

),

(

d

ˆ

xi

/

dt

)(

0

)

= (

dx

ˇ

i

/

dt

)(

0

)

then

ˆ

xi

(

t

)

 ˇ

xi

(

t

)

for all t

>

0

.

Proof. From (3.3), (3.4), and (3.5), it follows that

ˆ

h_i

(

η

)

 −

μ

i

τ

i

,

h

ˇ

i

(

η

)

 −

μ

i

τ

i

.

(3.17) Consider case (i). The function X

= ˆ

xi

−

xisatisfies

d2_X

dt2

+ (

τ

i

+

μ

i

)

d X

dt

= ˆ

hi

(

ˆ

xi

)

−

hi

(

x1

,

x2

)

and the right-hand side is equal to

ˆ

hi

(

x

ˆ

i

)

− ˆ

hi

(

xi

)

+ ˆ

hi

(

xi

)

−

hi

(

x1

,

x2

)

= ˆ

hi

(

η

i

)

X

+ ˆ

hi

(

xi

)

−

hi

(

x1

,

x2

)

where

η

i

=

η

i

(

t

)

lies between xiandx

ˆ

i, by the mean value theorem. Hence d2X

dt2

+ (

τ

i

+

μ

i

)

d X

where

a

(

t

)

= −ˆ

h_i

(

η

i

)



μ

i

τ

i

by (3.17) andh

ˆ

i

(

xi

)

−

hi

(

xi

,

xj

)



0,

(

i

,

j

)

= (

1

,

2

)

,

(

i

,

j

)

= (

2

,

1

)

. Applying Lemma 3.2 and Remark 3.1, we conclude that X

(

t

)



0 for all t

>

0. Hence

ˆ

xi

(

t

)



xi

(

t

)

for all t

>

0

.

The proofs of cases (ii) and (iii) are similar.

2

3.1. Single equilibrium

In this section, let us discuss the conditions under which h

ˆ

i (resp., h

ˇ

i) has a single zero and, consequently, by Lemma 3.1, all solutions to (3.7) (resp., (3.8)) converge to a single point

(

a

ˆ

i

,

0

)

(resp.,

(

a

ˇ

i

,

0

)

), as t

→ ∞

. According to (3.6), if

μ

i

τ

i

>

ν

i

α

in

˜

ki (3.18)

thenh

ˆ

_i

( ˜ξ

i

) <

0 and, consequently,

−

μ

i

τ

i

 ˆ

hi

(

xi

) <

0 for all 0



xi



Bi

;

then also

ˇ

h_i

(

xi

) <

0 and

∂

hi

(

x1

,

x2

)

∂

xi

<

0 for 0



xi



Bi

.

Note that

∂

hi

(

x1

,

x2

)/∂

xi(with x2fixed if i

=

1 and x2 fixed if i

=

2) attains its maximum at the same

point xi

= ˜ξ

i whereh

ˆ

i

(

xi

)

attains its maximum.

In addition to condition (3.18), we consider other situations which are more of biological interest. Analogously to [2], we assume that, for a given i (i

=

1 or i

=

2),

μ

i

τ

i

<

ν

i

α

in

˜

ki

·

1 1

+

Bj

/

γ

j

,

j

=

i

.

(3.19)

These conditions are equivalent to h

ˇ

_i

( ˜ξ

i

) >

0 and, in that case, if

˜ξ

i

<

Bi then each of h

ˆ

i, h

ˇ

i has two critical points. Let p

ˆ

m_i

,

p

ˆ

M_i (resp., p

ˇ

m_i

,

p

ˇ

M_i

)

denote the points where h

ˆ

i (resp.,h

ˇ

i) achieves its local minimum and maximum. Each of functionsh

ˆ

i

, ˇ

hi may have one or three zeros as illustrated in Fig. 2.

We consider the following cases for i

=

1 or i

=

2:

(

Mai

) ˆ

hi

(

p

ˆ

Mi

) <

0;

(

Mbi

) ˇ

hi

(

p

ˇ

m_i

) >

0;

(Bi) h

ˆ

i

(

p

ˆ

m_i

) <

0,h

ˇ

i

(

p

ˇ

M_i

) >

0.

Fig. 2.h1ˆ andh1ˇ have one zero in cases (a), (b), (c), and three zeros in case (d).

Proposition 3.4. Suppose one of the conditions (3.18) or (3.19) with

(

Mai

)

, or (3.19) with

(

Mbi

)

holds for i

=

1 or i

=

2. Then every solution of (3.7) (resp., (3.8)) converges to a single equilibrium

(

a

ˆ

i

,

0

)

(resp.,

(

a

ˇ

i

,

0

)

). 3.2. Multiple equilibria

In this section we assume that (3.19) and (Bi) hold where i

=

1 or i

=

2. Then the dynamics (3.7) (resp., (3.8)) has three equilibrium points:

(

a

ˆ

i

,

0

), (ˆ

bi

,

0

), (

ˆ

ci

,

0

)

(resp.,

(

a

ˇ

i

,

0

), (ˇ

bi

,

0

), (

c

ˇ

i

,

0

)

) wherea

ˆ

i

<

ˆ

bi

<

c

ˆ

i(resp.,a

ˇ

i

< ˇ

bi

<

c

ˇ

i) anda

ˇ

i

<

a

ˆ

i,c

ˇ

i

<

c

ˆ

i, but

ˇ

bi

> ˆ

bi

.

(3.20)

Proposition 3.5. Under the conditions (3.19) and

(

Bi

)

, every solution of (3.7) (resp., (3.8)) converges to one of the equilibrium points

(

a

ˆ

i

,

0

)

,

(ˆ

bi

,

0

)

,

(

ˆ

ci

,

0

)

(resp.,

(

a

ˇ

i

,

0

)

,

(ˇ

bi

,

0

)

,

(

c

ˇ

i

,

0

)

).

It can easily be computed that the equilibrium points

(

a

ˆ

i

,

0

), (

c

ˆ

i

,

0

)

(resp.,

(

a

ˇ

i

,

0

), (

c

ˇ

i

,

0

)

) of (3.7) (resp., (3.8)) are both sinks, whereas the equilibrium

(ˆ

bi

,

0

)

(resp.,

(ˇ

bi

,

0

)

) is a saddle. In addition, one branch of the unstable manifold for

(ˆ

bi

,

0

)

(resp.,

(ˇ

bi

,

0

)

) converges to

(

a

ˆ

i

,

0

)

(resp.,

(

a

ˇ

i

,

0

))

, and the other branch converges to

(

c

ˆ

i

,

0

)

(resp.,

(

ˇ

ci

,

0

)

), as t

→ ∞

, cf. Fig. 1. We denote by Ws

(ˆ

bi

)

(resp., Ws

_(ˇ

_b

i

)

) the (one-dimensional) stable manifold for

(ˆ

bi

,

0

)

(resp.,

(ˇ

bi

,

0

)

) and set

ˆ

yi

=

dx

ˆ

i

/

dt,

ˇ

yi

=

d

ˇ

xi

/

dt. We partition the phase plane for (3.7) and (3.8) respectively



(

ˆ

xi

,

ˆ

yi

)

: 0

 ˆ

xi



Bi

,

y

ˆ

i

∈ R



=

Ws

(ˆ

bi

)

∪

U

(

a

ˆ

i

)

∪

U

(

c

ˆ

i

),



(

ˇ

xi

,

ˇ

yi

)

: 0

 ˇ

xi



Bi

,

y

ˇ

i

∈ R



=

Ws

(ˇ

bi

)

∪

U

(

a

ˇ

i

)

∪

U

(

c

ˇ

i

),

where U

(

p

)

is the basin of attraction for sink

(

p

,

0

)

= (ˆ

ai

,

0

)

,

(

ˆ

ci

,

0

)

,

(

a

ˇ

i

,

0

)

,

(

c

ˇ

i

,

0

)

. Notice that Ws

(ˆ

bi

)

and Ws

(ˇ

bi

)

do not intersect. Indeed, if they intersect at one point

(

u0

,

v0

)

, then we can apply

Lemma 3.3(iii) with initial point

(

u0

,

v0

)

and deduce that b

ˆ

i

< ˇ

bi, a contradiction to (3.20). In ad-dition, Ws

_(ˆ

_b

i

)

lies on the left-hand side of Ws

(ˇ

bi

)

, again by Lemma 3.3(iii). Moreover, Ws

(ˆ

bi

)

(resp., Ws

(ˇ

bi

)

) is tangent at

(ˆ

bi

,

0

)

(resp.,

(ˇ

bi

,

0

)

) to the stable subspace Eswhich is given, respectively, by

Es

(ˆ

bi

,

0

)

=

span



1

,

−(

τ

i

+

μ

i

)

−



(

τ

i

+

μ

i

)

2

+

4h

ˆ

_i

(ˆ

bi

)

2



,

Es

(ˇ

bi

,

0

)

=

span



1

,

−(

τ

i

+

μ

i

)

−



(

τ

i

+

μ

i

)

2

+

4h

ˇ

_i

(ˇ

bi

)

2



.

Let

(

x

ˆ

i

(

t

;

u0

,

v0

),

y

ˆ

i

(

t

;

u0

,

v0

))

(resp.,

(

x

ˇ

i

(

t

;

u0

,

v0

),

y

ˇ

i

(

t

;

u0

,

v0

))

) be the solution to (3.7) (resp., (3.8)),

starting from point

(

u0

,

v0

)

at t

=

0, i

=

1

,

2. Clearly, if

(

u0

,

v0

)

∈

U

(

a

ˆ

i

)

∩

U

(

a

ˇ

i

)

, then as t

→ ∞

,



ˆ

xi

(

t

;

u0

,

v0

),

y

ˆ

i

(

t

;

u0

,

v0

)



→ (ˆ

ai

,

0

),



ˇ

xi

(

t

;

u0

,

v0

),

y

ˇ

i

(

t

;

u0

,

v0

)



→ (ˇ

ai

,

0

)

;

if

(

u0

,

v0

)

∈

U

(

c

ˆ

i

)

∩

U

(

ˇ

ci

)

, then



ˆ

xi

(

t

;

u0

,

v0

),

ˆ

yi

(

t

;

u0

,

v0

)



→ (ˆ

ci

,

0

),



ˇ

xi

(

t

;

u0

,

v0

),

y

ˇ

i

(

t

;

u0

,

v0

)



→ (ˇ

ci

,

0

)

;

if

(

u0

,

v0

)

∈ [

U

(

c

ˆ

i

)

∩

U

(

a

ˇ

i

)

]

, then



ˆ

xi

(

t

;

u0

,

v0

),

ˆ

yi

(

t

;

u0

,

v0

)



→ (ˆ

ci

,

0

),



ˇ

xi

(

t

;

u0

,

v0

),

y

ˇ

i

(

t

;

u0

,

v0

)



→ (ˇ

ai

,

0

)

as t

→ ∞.

In addition, U

(

a

ˆ

i

)

∩

U

(

c

ˇ

i

)

= ∅

, according to Lemma 3.3. As seen in Fig. 3, Ws

(ˇ

bi

)

lies to the right of Ws

(ˆ

bi

)

. Orbits of (3.9)–(3.10) cannot enter the region bounded by Ws

(ˇ

bi

)

and Ws

(ˆ

bi

)

, but an orbit initially from this region may exit it.

Fig. 3. Stable manifolds Ws

(ˆb1), Ws

(ˇb1)(indicated as blue lines) and unstable manifolds Wu

(ˆb1), Wu

(ˇb1)(indicated as red lines) of(ˆb1,0)and(ˇb1,0). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

4. Asymptotic behavior: single limit point

As in [2], we shall introduce an iterative scheme to prove convergence to a single point; the last step in the convergence proof will require the following condition:

ν

1

(

α

1

+

σ

1

)

γ

2

<





μ

1

τ

1

−

ν

1

α

1n

˜

k1



 −

ν

1

σ

1

ρ

1

,

ν

2

(

α

2

+

σ

2

)

γ

1

<





μ

2

τ

2

−

ν

2

α

2n

˜

k2



 −

ν

2

σ

2

ρ

2

.

(4.1)

We also assume that functions G

(

t

)

(in (2.10)) and Ei

(

t

)

(in (2.13)) satisfy the following conditions:

lim

t→∞G

(

t

)

and limt→∞Ei

(

t

)

exist. (4.2)

Theorem 4.1. Assume that (3.18) holds for i

=

1 and i

=

2, and that (4.1) and (4.2) hold. Then every solution of (3.9)–(3.10) converges to a single point

(

a

¯

1

,

0

,

a

¯

2

,

0

)

, as t

→ ∞

.

Corollary 4.2. The solution

ψ

of (2.8)–(2.13) has the following asymptotic behavior:

ψ(

t

,

x1

,

x2

,

y1

,

y2

)

→

N0

δ

(a¯1,a¯2,τ1a¯1/ν1,τ2¯a2/ν2)in measure, as t

→ ∞

.

Proof of Theorem 4.1. Set, for t



0,

Smin_i

(

t

)

=

inf



Si

(

s

)

: s

∈ [

t

,

∞)



,

Smax_i

(

t

)

=

sup



Si

(

s

)

: s

∈ [

t

,

∞)



.

Then Smin

i

(

t

)



Si

(

t

)



Smaxi

(

t

)

. Note that Smini

(

t

)

is nondecreasing, Smaxi

(

t

)

is nonincreasing, and Smin_i

(

t

)

ρ

i

+

Smini

(

t

)



Si

(

t

)

ρ

i

+

Si

(

t

)



Smax_i

(

t

)

ρ

i

+

Smax_i

(

t

)

for t



0

.

Under the condition (3.18),h

ˆ

i(resp.,h

ˇ

i) is a strictly decreasing function, and has a single zero, denoted by a

ˆ

i (resp.,a

ˇ

i). Let

(

x1

(

t

),

v1

(

t

),

x2

(

t

),

v2

(

t

))

be the solution to (3.9)–(3.10), starting from arbitrary

initial point

(

x1

(

0

),

v1

(

0

),

x2

(

0

),

v2

(

0

))

∈ Ω

∗. By Lemmas 3.1(ii), 3.3 and Proposition 3.4, for any small