Numerical illustrations

In this section, we provide numerical simulations for the single-cell model (2.1), (2.2) and for the population model (2.4).

The single-cell model is a system of two ordinary differential equations (ODEs) which can be easily solved by the Runge–Kutta method, using ode45 in MATLAB. The population model (2.4) is es-sentially an integro-differential equation. The integrations in the Si

(

t

)

need to be carried out through quadrature rule (numerical integration); we shall use Simpson’s rule which has third order accu-racy. The solution of Eq. (2.4) is then obtained by using Lax–Friedrichs method [10,11]. Notice that

the asymptotic solution of the population model becomes singular for large time. In order to ob-tain highly accurate solution, refinement is definitely needed at the places where population density tends to grow, and the corresponding quadrature rule has to be redesigned; this we have done for the one-peak case, but not for the multi-peak cases: we stopped the numerical simulations after the asymptotic singular solutions are observed.

9.1. The single-cell model

In Fig. 7, we first demonstrate the single-cell model results. The parameters for (a)–(e) are chosen as those in [15], namely,

μ ₌

5 day⁻¹

, α

1

= α

2

=

^{5 day−1}

, σ

1

= σ

2

=

^{5 day−1}

,

(9.1) k1

=

^k2

=

¹

, ρ

1

= ρ

2

=

¹

, γ

1

=

¹

, γ

2

=

⁰

.

5

,

(9.2)

β

1

= β

2

=

⁰

.

05 day⁻¹

,

n

=

⁶

.

(9.3)

For these parameters, we take A1

=

^A2

=

²

.

01. Then

α

1n

˜

k1

·

¹

1

+

^A2

/ γ

2

< μ < α

1n

˜

k1

, α

2n

˜

k2

·

¹

1

+

^A1

/ γ

1

< μ < α

2n

˜

k2

,

so that the conditions (M1) and (M2) are not satisfied. Thus,

ˆ

_fi defined in (4.1) has a local minimum and a local maximum for i

=

¹

,

2.

In addition, conditions (B1), (B2) hold for the B_idefined in Section 4. This gives the flexibility for the system to be either monostable or bistable under different choices of S1and S2. For example, the system is

(a) monostable (MS) for S₁

=

⁰

.

05

,

S₂

=

⁰

.

025, with the choice of B₁

=

⁰

.

058 and B₂

=

⁰

.

035; in this case (B1), (B2) hold and

ˆ

_fi

(ˆ

p^M_i

) <

0 (i

=

¹

,

2);

(b) bistable (BS-ll,hl) for S₁

=

¹

.

2

,

S₂

=

⁰

.

025 with B₁

=

¹

.

181 and B₂

=

⁰

.

035; in this case (B1), (B2) hold, and

ˆ

_f1

(

p

ˆ

^m₁

) <

0

, ˇ

f1

(

p

ˇ

^M₁

) >

0,

ˆ

_f2

(

p

ˆ

^M₂

) <

0;

(c) bistable (BS-ll,lh) for S₁

=

⁰

.

05

,

S₂

=

¹

.

3 with B₁

=

⁰

.

058 and B₂

=

¹

.

493; in this case (B1), (B2) hold and

ˆ

_f2

(

p

ˆ

^m₂

) <

0

, ˇ

f2

(ˇ

p^M₂

) >

0,

ˆ

_f1

(

p

ˆ

^M₁

) <

0 hold.

Each of these cases is shown in Fig. 7, where we chose 36 different initial conditions

(

x1

(

0

),

x2

(

0

))

and depicted their evolution. The blue curve is the nullcline of f1while the red curve is the nullcline of f2. We can clearly see that the solutions converge to a single stable equilibrium in case (a) and to two stable equilibria in cases (b) and (c). The bistable-ll,hl (bistable-ll,lh) system with low x₁–low x₂ and high x1–low x2 states (low x1–low x2 and low x1–high x2 states), shown in case (b) ((c)) can become monostable with high x1–low x2 state (low x1–high x2), shown in (d) ((e)) by increasing the value of S₁

(

S₂

)

. Notice that (d) ((e)) satisfies conditions (B1),

ˇ

f₁

(ˇ

^pm₁

) >

0 and (B2),

ˆ

f₂

(

^p

ˆ

^M₂

) <

0 with B₁

=

¹

.

626 and B₂

=

⁰

.

035 ((B2),

ˇ

f₂

(

^p

ˇ

^m₂

) >

0, (B1),

ˆ

f₁

(

^p

ˆ

^M₁

) <

0 with B₁

=

⁰

.

058 and B₂

=

¹

.

626). It is also possible to switch from bistable-ll,hl (bistable-ll,lh) to monostable by decreasing S2(S1).

However, the system with parameters (9.1)–(9.3) cannot be quadstable due to the strong mutual inhibition (i.e., small

γ

1,

γ

2). If we decrease the mutual inhibition by taking parameters

γ

1

= γ

2

=

^30,

σ

1

= σ

2

=

^{2, k}1

=

^k2

=

⁰

.

6, but keep all the other parameters the same, then conditions (B1), (B2),

ˆ

fi

(ˆ

p^m_i

) <

0

, ˇ

fi

(

p

ˇ

^M_i

) >

0, i

=

¹

,

2, are satisfied, and by Theorem 4.1(iv), the system is quadstable, as illustrated in Fig. 8.

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

Fig. 7. Single-cell model. (a) Monostable: the stable equilibrium is(x1,x2)≈ (⁰.055,0.033).(b) Bistable-ll,hl: the stable equilib-ria are(x1,x2)≈ (0.556,0.026)and(1.368,0.020); the unstable equilibrium is(x1,x2)≈ (0.976,0.022).(c) Bistable-ll,lh: the stable equilibria are(0.032,0.602)and(0.022,1.444); the unstable equilibrium is(0.027,0.904). (d) Monostable: the stable equilibrium is(1.549,0.020).(e) Monostable: the stable equilibrium is(0.042,1.544). Notice that the dark blue curve is the nullcline of f1while the red curve is the nullcline of f2. The black curves converge to ll state, the green curves converge to hl state, and the light blue curves converge to lh state. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Fig. 8. Single-cell model. Quadstable: S1=⁰.04, S2=¹.6, the stable equilibria are (x1,x2)≈ (⁰.211,0.211), (0.204,1.186), (1.186,0.204),(1.146,1.146)and the unstable equilibria(x1,x2)≈ (⁰.208,0.531),(0.531,0.208),(0.535,0.535),(0.543,1.171), (1.171,0.543).

9.2. The population model

Since in the population model (2.4) S₁ and S₂ are not constant and their evolution depends on both the initial population of cells and the external signals C1

(

t

)

, C2

(

t

)

, one may expect interesting be-havior; for example, the system may switch from one-peak to two-peak profile at intermediate times.

In the subsequent numerical simulations we adapt the normalized population density

ψ(

t

,

x1

,

x2

)

, take A₁, A₂as in Section 9.1, and choose the initial condition

ψ

0

(

x1

,

x2

) =

^const

=

¹

A₁A₂ (9.4)

so that N0

=

1. Although in (3.1) we assumed that

ψ

0

=

^{0 on}

∂Ω

, the results of Sections 6–8 do not actually use this assumption. Furthermore, the simulations given below do not significantly change if we modify (9.4) near the boundary

∂Ω

so as to make

ψ

0 vanish there.

In Sections 9.2.1–9.2.3 we take C_i

(

t

) =

0, i.e., there is no external stimulus. In Section 9.2.4 we examine the effect of the stimulus C_i

(

t

)

.

We first demonstrate one-, two-, and four-peak solutions by choosing specific parameters in the regimes we discussed in Theorem 4.1.

9.2.1. Asymptotic one-peak solution

In Fig. 9 we show numerical results under conditions (M1) and (M2) which guarantee a single attracting point. Notice that we choose k1

=

^k2

=

2 instead of k1

=

^k2

=

1 in [15] in order to satisfy conditions (M1) and (M2). In Fig. 9(a), 9(b), 9(c) we have plotted

ψ

and the corresponding vector field

(

f₁

,

f₂

)

at times t

=

⁰

.

05

,

0

.

2

,

5. Since (M1) and (M2) are satisfied no matter what S₁and S₂are, there is only one stable equilibrium point (although the sufficient condition (6.6) in Theorem 6.2 is not satisfied). The vectors

(

f1

,

f2

)

all point toward the attracting point. The normalized population density gets more and more concentrated at an attracting point, and

(

S1

(

t

),

S2

(

t

))

converges to

(¯

a1

,

a

¯

2

) ≈ (

0

.

054

,

0

.

081

)

.

9.2.2. Asymptotic two-peak solution

Fig. 10 displays bistable case (bistable-ll,lh) with two-peak solution. We choose parameters

σ

1

= σ

2

=

^2,

γ

1

=

^30,

γ

2

=

^{1, k}1

=

⁵

,

k2

=

⁰

.

6. We see that the population density starts to ac-cumulate at two attracting points and the population density is higher in low x1–high x2 state as proved in Section 7. The weights w₁ and w₂ in the asymptotic solution depend on the initial

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

Fig. 9. Monostable: k1=^k2=2 but all other parameters are as in (9.1)–(9.3). (a) t=⁰.05, (b) t=⁰.2, (c) t=^5.

lation density. If most of the population density is initially in the attraction basin of low x1–low x2 state, then the weight for the Dirac function with center at low x₁–low x₂ state would be higher (not shown here).

Fig. 10. Bistable (BS-ll,lh): σ1=σ2=^2, γ1=^30, γ2=^{1, k1}=^{5, k2}=⁰.6 and all other parameters are as in (9.1)–(9.3).

(a) t=0.05, (b) t=0.2, (c) t=1.

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

9.2.3. Asymptotic four-peak solution

In Fig. 11, the population density becomes highly concentrated at four attracting points as we expect from Theorem 8.1. The weights w11

,

wu1

,

w1u and wuu depend on the parameters of the system as well as on initial population density. The parameters chosen satisfy the condition (B1) and (B2) (and (6.6) is also satisfied). Note that the mutual inhibition is small (i.e.,

γ

1 and

γ

2are large).

9.2.4. Effect of the stimulus

In the previous subsections we have assumed that Ci

(

t

) ≡

0 (no external stimulus). We now want to examine the effect of these stimuli. We take the parameters as in (9.1)–(9.3): Fig. 12 shows how with no stimuli (i.e., with C1

(

t

) ≡

^C2

(

t

) ≡

0) the uniform populations begin to evolve and move into low x1–low x2 peak; this is interpreted biologically as no cell differentiation. In Fig. 13 we choose C₁

(

t

)

exp^−G(t)

=

⁰

.

5 and C₂

(

t

)

exp^−G(t)

=

¹

.

5 for all t

>

0. We see that the solution develops a two-peak solution. Due to the larger stimulus of x2 (i.e., C2

(

t

) >

C1

(

t

)

), as well as the stronger inhibition of x1 by x2, the low x1–high x2 peak appears instead of high x1–low x2 peak.

In Fig. 14 we use the same stimuli as in Fig. 13, but have taken

ψ

0 to be constant for x₁

<

A₁

/

5 and zero elsewhere. Thus we give GATA-3 initial density advantage as well as stimulus advantage. We see that the population density moves again toward two-peak solution, low x1–low x2 and low x1– high x2, but the population density at the low x1–high x2 is larger than in Fig. 13. In both Figs. 13 and 14, the low x1–high x2can be interpreted biologically as a population of differentiated Th2 cells.

10. Conclusions

In this paper, we considered a conservative law of the form

∂φ

∂

t

+ ∂

∂

x₁

(

f1

φ) + ∂

∂

x₂

(

f2

φ) =

^g

φ, φ = φ(

^t

,

x1

,

x2

),

(10.1) where the velocity vector f

= (

^f1

,

f₂

)

is a nonlinear nonlocal function of

φ

. This equation arises as a model of T cell differentiation where x

= (

^x1

,

x2

)

, and x1

,

x2are the concentrations of transcription factors T-bet and GATA-3, respectively. A precursor T cell growing at rate g with x1 large (small) and x₂ small (large) will differentiate into Th1 (Th2) T cell. Th1 and Th2 have different functions: Th1 T cells combat intracellular pathogens while Th2 T cells induce the activation of B cells to combat extracellular pathogens. A ‘good’ balance between these two populations of cells is maintained in homeostasis. The function

φ (

t

,

x1

,

x2

)

represents the population density of T cells with concentrations

(

x1

,

x2

)

at time t. Within an individual cell the concentrations of x1 and x2 vary according to the equations

dx_i dt

=

^fi



t

,

x1

,

x2

, φ (

t

, ·) 

,

i

=

¹

,

2

,

(10.2)

and the dependence on t and

φ (

t

, ·)

arises from stimuli Si

(

t

)

consisting of a stimulus which arises from within the entire population of the T cells and of an external stimulus Ci.

It is natural to ask what is the behavior of

φ

at intermediate and large times and how this depends on Ci

(

t

)

and on the initial condition. In this paper, we have depicted six regions from the space of parameters that are introduced in the definition of the f_i. We proved that for the first regime the function

φ (

t

,

x1

,

x2

)

converges to a 1-peak solution as t

→ ∞

; for regimes 2, 3, 4, and 5,

φ

converges to a 2-peak solution, and for regime 6,

φ

converges to a 4-peak solution; this was illustrated in Figs. 9–11 when C_i

(

t

) ≡

^0.

Numerical simulations given in Figs. 12–14, show how the location of these peaks depends on the external signals Ci

(

t

)

and the initial conditions. We interpret a peak centered at

(

x⁰₁

,

x⁰₂

)

with x⁰₁

,

x⁰₂ small as a population of T cell that do not differentiate. A peak with x⁰₁ small, x⁰₂ large represents a population of T cells that differentiate into Th2. In a similar way, we interpret the case of x⁰₁ large, x⁰₂ small. Finally, a situation where both x⁰₁ and x⁰₂ are large in viewed as abnormal: Since the con-centrations of both T-bet and GATA-3 are large, the cell receive conflicting instructions to differentiate

Fig. 11. Quadstable:σ1=σ2=2,γ1=γ2=30, k1=k2=0.6 and all other parameters are as in (9.1)–(9.3). (a) t=0.05, (b) t=0.2, (c) t=1.

simultaneously to Th1 and Th2. This situation arises in Fig. 11 where the mutual inhibition is weak (namely,

γ

1

= γ

2

=

30). Hence one of the conclusions of our simulations is that, in homeostasis, the mutual inhibition cannot be too weak.

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

Fig. 12. Monostable result: the parameters are as in (9.1)–(9.3). The population density moves toward low x1–low x2 state at (a) t=⁰.05, (b) t=¹.0 and (c) t=^5.

The results of the paper are obtained by approximating the full dynamical system (10.2) from above and below by a sequence of dynamical systems where in each step of approximation the total signaling is constant but is ‘sharper’ than in the previous step. This method is quite general and could be applied to more general functions f

(

t

,

x

, φ ( ·))

and in any number of dimensions for the x variable.

Fig. 13. The parameters are as in (9.1)–(9.3). A uniform density evolves toward two stable points under external stimulus C1(t)e^−G(t)=⁰.5, C1(t)e^−G(t)=¹.5. (a) t=⁰.01, (b) t=⁰.1 and (c) at t=^1.

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

Fig. 14. The parameters are as in (9.1)–(9.3). A uniform density at x1<A1/5 evolves toward two stable points under external stimulus C1(t)e^−G(t)=⁰.5, C1(t)e^−G(t)=¹.5. (a) t=⁰.01, (b) t=⁰.1 and (c) at t=^1.

Acknowledgments

A. Friedman was partially supported by the National Science Foundation under Agreement No. 0112050. C.-Y. Kao was partially supported by the National Science Foundation grant DMS-0811003. C.-W. Shih was partially supported by the NSC and NCTS of Taiwan; he is grateful to the MBI at OSU for the hospitality and support during his visit in 2007–2008; he also thanks Y.-H. Wan for helpful discussion.

References

[1] B.D. Aguda, A. Friedman, Models of Cellular Regulation, Oxford Univ. Press, 2008.

[2] C.Y. Cheng, K.H. Lin, C.W. Shih, Multistability in recurrent neural networks, SIAM J. Appl. Math. 66 (4) (2006) 1301–1320.

[3] J. Cherry, F. Adler, How to make a biological switch, J. Theoret. Biol. 203 (2) (2000) 117–133.

[4] O. Cinquin, J. Demongeot, High-dimensional switches and the modeling of cellular differentiation, J. Theoret. Biol. 233 (2005) 391–411.

[5] O. Cinquin, K. Page, Generalized switch-like competitive heterodimerization networks, Bull. Math. Biol. 69 (2007) 483–494.

[6] A. Friedman, Mathematics in Industrial Problems, IMA Vol. Math. Appl., vol. 16, Springer-Verlag, New York, 1988.

[7] A. Friedman, B. Ou, A model of crystal precipitation, J. Math. Anal. Appl. 137 (1989) 550–575.

[8] A. Friedman, B. Ou, D.S. Ross, Crystal precipitation with discrete initial data, J. Math. Anal. Appl. 137 (1989) 576–590.

[9] A. Friedman, D.S. Ross, Mathematical Models in Photographic Science, Springer-Verlag, Berlin, 2002.

[10] B. Gustafsson, H.-O. Kreiss, J. Oliger, Time-Dependent Problems and Difference Methods, Wiley–Interscience, New York, 1995.

[11] R.J. LeVeque, Numerical Methods for Conservation Laws, Birkhäuser, Basel, 1990.

[12] K.H. Lin, C.W. Shih, Multiple almost periodic solutions in nonautonomous delayed neural network, Neural Comput. 19 (12) (2007) 3392–3420.

[13] L. Mariani, M. Lohning, A. Radbruch, T. Hofer, Transcriptional control networks of cell differentiation: Insights from helper T lymphocytes, Biophys. Mol. Biol. 86 (2004) 45–76.

[14] C.W. Shih, J.P. Tseng, Convergent dynamics for multistable delayed neural networks, Nonlinearity 21 (2008) 2361–2389.

[15] A. Yates, R. Callard, J. Stark, Combining cytokine signaling with T-bet and GATA-3 regulation in Th1 and Th2 differentiation:

A model for cellular decision-making, J. Theoret. Biol. 231 (2004) 181–196.

[16] J. Zhang, A nonlinear nonlocal multidimensional conservation law, J. Math. Anal. Appl. 204 (1996) 353–388.

在文檔中 Asymptotic phases in a cell differentiation model (頁 23-34)