In this section, we provide numerical simulations for the single cell model (2.1) and the population model (2.3). The parameters used in some of the simulations do not satisfy the assumptions made in the previous theorems.
The single cell model is a system of four ODEs which can be easily solved by the Runge–Kutta method, using ode45 in MATLAB. The population model (2.3) is essentially an integro-differential equation. The integrations in the Si
(
t)
need to be carried out through quadrature rule (numerical integration); we shall use midpoint rule which has second-order accuracy. The solution of Eq. (2.3) is then obtained by using standard Lax–Friedrichs method [3,7]. Notice that the asymptotic solution of the population model becomes singular for large time.6.1. The single cell model
In Fig. 6, we first show the dynamics of ODE system which exhibit single limit point, with param-eters
n
=
6, α
1= α
2=
5,
k1=
k2=
5, ρ
1= ρ
2=
1, β
1= β
2=
0.
05,
σ
1= σ
2=
5, μ
1= μ
2=
5, γ
1= γ
2=
1, ν
1= ν
2=
10, τ
1= τ
2=
5.
(6.1) All the parameters are as in [10] except for ki(the level of T-bet or GATA-3 at which the rate of auto-stimulation is half-maximum) which is altered from 1 to 5 to make sure that the conditions in (3.18) are satisfied. We choose six different initial pointsx1
(
0),
x2(
0),
y1(
0),
y2(
0)
=
i5B1
, (
5−
i)
5 B2,
i5A1
, (
5−
i)
5 A2,
for integer i from 0 to 5. This means that there are six cells with different initial levels of transcription factors and mRNA. The polarizing cytokines are S1
=
B1/
2 and S2=
B2/
2. (The same choice is used in the ODEs’ simulations unless otherwise mentioned.) It is clear that all the points from different initial locations converge to a single limit point which has low concentration of T-bet (x1) and low concentration of GATA-3 (x2). Thus all six cells do not differentiate. The plot of(
x1,
x2,
y1,
y2)
versusFig. 7. ODE simulation (a) x1,x2,y1,y2versus t, (b) the trajectories in(x1,x2)-plane, (c) the trajectories in(y1,y2)-plane. (For interpretation of the references to color, the reader is referred to the web version of this article.)
time is shown in Fig. 6(a). The trajectories in
(
x1,
x2)-plane and(
y1,
y2)-plane are shown in Figs. 6(b) and 6(c) respectively. The trajectories may cross each other on these projection planes. However, they do not cross each other in the four-dimensional space.In Fig. 7, we illustrate the case of two limit points with parameters
n
=
6, α
1= α
2=
5,
k1=
0.
9,
k2=
0.
6, ρ
1= ρ
2=
1, β
1= β
2=
0.
05, σ
1= σ
2=
2, μ
1= μ
2=
5, γ
1= γ
2=
30, ν
1= ν
2=
5, τ
1= τ
2=
5.
The corresponding time evolution and projection on
(
x1,
x2)-plane and(
y1,
y2)-plane are given in Figs. 7(a), 7(b), and 7(c) for 16 different initial conditions. This case demonstrates that there are two stable limit points (blue dot and green dot). The trajectories which converge to blue (green) limit point are colored as blue (green). In this parameter setting, the green point has larger attracting basin. Thus there are more cells with low concentrations of T-bet (x1) and high concentrations of GATA-3 (x2) and they differentiate into Th2 cells.In order to generate quadstable phase, we choose the parameters as
n
=
6, α
1= α
2=
5,
k1=
k2=
0.
6, ρ
1= ρ
2=
1, β
1= β
2=
0.
05, σ
1= σ
2=
2, μ
1= μ
2=
5, γ
1= γ
2=
30, ν
1= ν
2=
5, τ
1= τ
2=
5.
We can see that the system is quadstable, as illustrated in Fig. 8. The corresponding time evolution and projection on
(
x1,
x2)-plane and(
y1,
y2)-plane are given in Figs. 8(a), 8(b), and 8(c). We start with 81 different initial conditions; the trajectories converge to one of the four stable limit points.For example, there are nine different trajectories emitting from the center point (0
.
705,
0.
705) on(
x1,
x2)-plane because their initial locations on(
y1,
y2)-plane are different. One, two, two, and four trajectories tend to blue, red, green, and cyan limit points, respectively, which indicate no differenti-ation (low-x1, low-x2), Th1 (high-x1, low-x2) differentiation, Th2 (low-x1, high-x2) differentiation and undetermined (abnormal). It is hard to distinguish two of the trajectories (colored as mixed cyan and blue) because they overlap each other.In Fig. 9, we use the same parameters as in [10]:
n
=
6, α
1= α
2=
5,
k1=
k2=
1, ρ
1= ρ
2=
1, β
1= β
2=
0.
05, σ
1= σ
2=
5, μ
1= μ
2=
5, γ
1=
1, γ
2=
0.
5, ν
1= ν
2=
10, τ
1= τ
2=
5.
The parameters do not satisfy the conditions (Ma1), (Ma2), (B1) and (B2). The number of limit points vary with respect to the polarizing cytokines S1 and S2. When S1
=
0.
05, S2=
0.
025, the trajectoriesFig. 8. ODE simulation (a) x1,x2,y1,y2versus t, (b) the trajectories in(x1,x2)-plane, (c) the trajectories in(y1,y2)-plane. (For interpretation of the references to color, the reader is referred to the web version of this article.)
approach to a low–low single limit point (no differentiation). When S1
=
0.
15, S2=
1.
0, the trajec-tories approach two limit points. If we further increase S2 to 1.7, it becomes one single limit point again. However, the limit point is at a low–high level which indicates the differentiation toward Th2.6.2. The population model
The main difference between single cell model and population model is the coupling dynamics generated from the total signals S1and S2defined in (2.2) for population model. Since S1and S2are not constants, their evolutions depend on both the initial population of cells and the external signals E1
(
t)
, E2(
t)
. In [2], we have demonstrated, for the model (1.6), an interesting behavior, namely, the system may switch from one-peak to two-peak profile at intermediate times. Here we will focus on the singular behaviors which demonstrate one-, two- and four-peak solutions by choosing specific parameters.In the subsequent numerical simulations we assume that g
=
0,ψ
0= φ
0=
0 on∂Ω
, and Ei(
t) =
0.We take A1, A2, B1, B2 as in (2.6), (2.7) and choose the initial condition
ψ
0(
x1,
x2,
y1,
y2) =
constant=
1 A1A2B1B2(6.2)
so that N0
=
1. Even though the discretization method is the same as the one used in [2], the com-putations for the present four-dimensional model are much more intensive.6.2.1. Asymptotic one-peak solution
In Fig. 10, we show numerical results under conditions (Ma1) and (Ma2) which guarantee a sin-gle attracting point. The choice of parameters is as (6.1). In Figs. 10(a), 10(b), and 10(c), we plot
ψ
dy1dy2, at times t=
0.
5,
1,
5 respectively, because it is hard to visualize the original four-dimensional density function. Since (Ma1) and (Ma2) are satisfied no matter what S1 and S2 are, there is only one stable equilibrium point. The normalized population density gets more and more concentrated at an attracting point with low-x1and low-x2 so there is no polarization toward differ-entiation into Th1 or Th2.6.2.2. Asymptotic two-peak solution
Fig. 11 displays the bistable case (two-peak solution). We choose parameters
n
=
6, α
1= α
2=
5,
k1=
0.
9,
k2=
0.
6, ρ
1= ρ
2=
1, β
1= β
2=
0.
05,
σ
1= σ
2=
2, μ
1= μ
2=
5, γ
1= γ
2=
30, ν
1= ν
2=
5, τ
1= τ
2=
5.
Fig. 9. The trajectories in(x1,x2)-plane and(y1,y2)-plane for (a) S1=0.05, S2=0.025,(b) S1=0.15, S2=1.0,and (c) S1= 0.15, S2=1.7.
Fig. 10. Monostable:
ψdy1dy2at (a) t=0.5, (b) t=1, (c) t=5.
Fig. 11. Bistable (BS-ll,lh):
ψdy1dy2at (a) t=0.5, (b) t=1, (c) t=5.
Fig. 12. Quadstable:
ψdy1dy2at (a) t=0.5, (b) t=1, (c) t=5.
We see that integration of the population density with respect to y1 and y2 starts to accumulate at two attracting points and the population density is higher in low-x1 high-x2 state as time evolves.
Thus a large portion of cells differentiates into Th2 while others do not differentiate. The weights w1 and w2 in the asymptotic solution depend on the initial population density and the parameters. If most of the initial population is distributed in the attraction basin of low-x1 low-x2 state, then the weight for the Dirac function with center at low-x1 low-x2 state would be higher (not shown here).
6.2.3. Asymptotic four-peak solution
In Fig. 12, the population density becomes highly concentrated at four attracting points as expected from Theorem 5.3. In this example, the parameters are chosen as
n
=
6, α
1= α
2=
5,
k1=
k2=
0.
6, ρ
1= ρ
2=
1, β
1= β
2=
0.
05, σ
1= σ
2=
2, μ
1= μ
2=
5, γ
1= γ
2=
30, ν
1= ν
2=
5, τ
1= τ
2=
5.
The weights wll
,
wul,
wlu and wuu depend on the parameters of the system as well as on the ini-tial population density. In this case, the population density is higher in high-x1 high-x2 state as time evolves. Even though the population densities at the other three points are small and hardly notice-able, there is indeed some population. The parameters chosen satisfy conditions (B1) and (B2). Note that the mutual inhibition is small (i.e.,γ
1 andγ
2 are large), in this case.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 and Alfred P. Sloan Fellowship. 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 and 2011.
References
[1] B.D. Aguda, A. Friedman, Models of Cellular Regulation, Oxford University Press, 2008.
[2] A. Friedman, C.-Y. Kao, C.-W. Shih, Asymptotic phases in a cell differentiation model, J. Differential Equations 247 (3) (2009) 736–769.
[3] B. Gustafsson, H.-O. Kreiss, J. Oliger, Time-Dependent Problems and Difference Methods, Wiley–Interscience, New York, 1995.
[4] J. Hale, Ordinary Differential Equations, Robert E. Krieger Publishing Company, 1980.
[5] J. Hale, Asymptotic Behavior of Dissipative Systems, Math. Surveys Monogr., vol. 25, Amer. Math. Soc., Providence, RI, 1988.
[6] J. Hale, H. Kocak, Dynamics and Bifurcations, Springer-Verlag, 1991.
[7] R.J. LeVeque, Numerical Methods for Conservation Laws, Birkhäuser-Verlag, Basel, 1990.
[8] P. Magal, X.-Q. Zhao, Global attractors and steady states for uniformly persistent dynamical systems, SIAM J. Math.
Anal. 37 (1) (2005) 251–275.
[9] 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.
[10] 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.