# Effect of Delay on Selection Dynamics in Long-Term Sphere Culture of Cancer Stem Cells

Authors: Peng Tang; Aijun Fan; Jianquan Li; Jun Jiang; Kaifa Wang
Journal: Discrete Dynamics in Nature and Society
 (2013)

---

## Abstract

To quantitatively study the effect of delay on selection dynamics in long-term sphere culture of cancer stem cells (CSCs), a selection dynamic model with time delay is proposed. Theoretical results show that the ubiquitous time delay in cell proliferation may be one of the important factors to induce fluctuation, and numerical simulations indicate that the proposed selection dynamical model with time delay can provide a better fitting effect for the experiment of a long-term sphere culture of CSCs. Thus, it is valuable to consider the delay effect in the future study on the dynamics of nongenetic heterogeneity of clonal cell populations.


---

## Body

## 1. Introduction

In the past years research on cancer stem cells (CSCs) has become a focus of cancer research, because CSCs have self-renewing and multidirectional differentiation capability and may result in tumors [1–6]. Recently, cell state dynamics due to non-genetic heterogeneity of clonal cell populations also has received more and more attention [7–10].In order to expand CSCs, sphere culture is performed by experimental cell biologists [11, 12]. However, whether long-term sphere culture can maintain a high ratio of CSCs is unclear. For this question, it is interesting that [13, 14] obtain a similar quantitative result through different mathematical model; that is, the ratio of CSCs will towards an apparent equilibrium state in a long-term sphere culture. Concretely, [13] proposed a kinetic model using ordinary differential equations that considered the symmetric and asymmetric division of CSCs, as well as the proliferation and transformation of differentiated cancer cells (DCCs). And [14] puts forward a Markov model in which cells transition stochastically between states. However, the time delay due to the maturation of individual cells has been ignored in [13, 14].In the present paper, based on the kinetic model in [13], we further explore the effect of time delay on selection dynamics in long-term sphere culture. The results show that the ubiquitous time delay in cell proliferation may be one of the important factors to induce fluctuation, and the proposed selection dynamical model with time delay can provide a better fitting effect for the experiment of a long-term sphere culture of CSCs in [13]. The organization of the paper is as follows. In Section 2, we formulate the selection dynamical model with time delay in long-term sphere culture of CSCs. Section 3 first gives the analytic analysis on our proposed model and then presents numerical simulations to compare the effect of time delay. Finally, some predictive conclusions with biological implications are given in Section 4.


## 2. Model Description

Letx(t) and y(t) denote the population sizes of CSCs and DCCs at time t in long-term sphere culture, respectively. Our previous work [13] proposed the following mathematical model to describe the interactive growth of the CSCs and DCCs in a long-term sphere culture:
(1)dx(t)dt=bx(1-βx)x(t)+byβyy(t)≜x(t)f1(x(t),y(t)),dy(t)dt=bxβxx(t)+by(1-βy)y(t)≜y(t)f2(x(t),y(t)).
Here the constants bx,by are called the net birth rate or intrinsic growth rate of population x,y, respectively. βx denotes the conversion rate from CSCs to DCCs in the process of CSCs proliferation, and βy denotes the conversion rate from DCCs to CSCs in the process of DCCs proliferation.Note that time delay may play an important role in many biological models. As shown in [15], the maturation of individual cells may need a period of time τ; that is, the number of these cells at time t may depend on the population at a previous time t-τ. Under the assumption of equal discrete retarded cell proliferation, model (1) can be modified to
(2)dx(t)dt=x(t-τ)f1(x(t-τ),y(t-τ)),dy(t)dt=y(t-τ)f2(x(t-τ),y(t-τ)).
Here τ is the time delay due to maturation time.For (2), inspired by [16], we give the following average fitness of the population:
(3)ϕ=bxx(t-τ)+byy(t-τ).
Thus, the selection dynamics in long-term sphere culture of CSCs can be written as
(4)dx(t)dt=x(t-τ)(f1(x(t-τ),y(t-τ))-ϕ),dy(t)dt=y(t-τ)(f2(x(t-τ),y(t-τ))-ϕ).
For (4), let N(t)=x(t)+y(t). We have
(5)dN(t)dt=(1-N(t-τ))ϕ.
Therefore, N(t)→1 as t→∞; that is, the total population size remains constant. Hence x(t) and y(t) in (4) can be understood as the frequency of CSCs and DCCs, respectively. Furthermore, since y(t) can be replaced by 1-x(t), system (4) describes only a single differential equation; that is,
(6)dx(t)dt=a1+a2x(t-τ)+a3x2(t-τ)≜f(x(t-τ)),
in which a1=byβy,a2=bx(1-βx)-by(1+βy), and a3=by-bx. Note that (6) is the final selection dynamical model with time delay in long-term sphere culture of CSCs.


## 3. Results

### 3.1. Dynamic Analysis

The objective of this subsection is to analyze the dynamical behavior of (6). In order to explore the effect of the delay, we split this into two cases.


#### 3.1.1. Case ofτ=0

In this case, we focus on the dynamic analysis if the delay is nonexistent; that is,τ=0 in (6). We start by studying the existence of nonnegative equilibria in the interval [0,1]. Let f(x)=a1+a2x+a3x2=0. Clearly, the discriminant of the quadratic equation is
(7)Δ=a22-4a1a3=(by-bx+bxβx-byβy)2+4bxbyβxβy>0.
Hence there are two different real roots for f(x)=0 if a3≠0. Furthermore, since
(8)f(1)=a1+a2+a3=-bxβx<0,
we know that there is a unique positive equilibrium x*∈(0,1) for (6) if a3≠0 (see Figures 1(a) and 1(b)). When a3=0 (see Figure 1(c)), it is clear that there is only one positive equilibrium
(9)x*=-a1a2=βxβx+βy<1.Illustrations of function image off(x)=0 under different cases ((a), (b), and (c)), and the phase diagram (d) of system (6). Here (a) a3>0, (b) a3<0, and (c) a3=0.
(a)
(b)
(c)
(d)The combination of the above results and the phase diagram (see Figure1(d)) of system (6) yields the following result.Proposition 1.
For system (6), when τ=0, a unique positive equilibrium x*∈(0,1) always exists, and it is globally asymptotically stable.


#### 3.1.2. Case ofτ>0

In this case, we focus on the dynamic analysis if the delay is existent; that is,τ>0 in (6). Clearly, the unique positive equilibrium x*∈(0,1) still remains for (6) in spite of the delay. To study the stability of the equilibrium x*, we first translate x* to the origin. Let
(10)x~=x-x*.
Then (6) becomes, after replacing x~ by x again,
(11)dx(t)dt=(a2+2a3x*)x(t-τ)+a3x2(t-τ).The variational system of (11) at the origin is given by
(12)dx(t)dt=(a2+2a3x*)x(t-τ).
The characteristic equation of linear system (12) is given by
(13)λ-(a2+2a3x*)e-τλ=0.
If we let λ=α+iβ, then (13) becomes
(14)α-(a2+2a3x*)e-ταcosβτ=0,β+(a2+2a3x*)e-ταsinβτ=0.
By setting α=0 in (14), we have
(15)cosβτ=0,sinβτ=-βa2+2a3x*.Solving the first algebraic equation in (15), we have
(16)βτ=π2+2nπ,n=0,1,2,….
Substituting (16) into the second equation of (15), we have
(17)β=-(a2+2a3x*).
Now, substituting (17) into (16), we have
(18)τn=-π2(a2+2a3x*)-2nπa2+2a3x*,n=0,1,2,….Next, we computeα′(τn). Differentiating (13) with respect to τ, we have
(19)dλ(τ)dτ=-(a2+2a3x*)λe-τλ1+τ(a2+2a3x*)e-τλ,
thus,
(20)α′(τn)=Re(dλ(τ)dτ|τ=τn)=β21+τn2β2>0.Therefore, similar to [17], according to the results in [18, 19] or [15, Theorem  2.2], we can obtain the following results on (6).Proposition 2.
Supposeτ>0. Then system (6) has a Hopf bifurcation at
(21)τ=τn=-π2(a2+2a3x*)-2nπa2+2a3x*,n=0,1,2,….Furthermore, according to the results in [15, 20, 21], we have the following.Proposition 3.
Supposeτ>0 in (6). Then the unique positive equilibrium x*∈(0,1) is stable if 0<τ<τ0 and unstable if τ>τ0.


### 3.2. Numerical Simulations

For model (1), we designed a long-term sphere culture of human breast cancer MCF-7 stem cells [13]. Based on the experimental data, using an adaptive Metropolis-Hastings (M-H) algorithm to carry out an extensive Markov-chain Monte-Carlo (MCMC) simulation, we obtained the estimated parameter values as follows:
(22)bx=2.7506×10-1,by=3.2635×10-1,βx=1.4407×10-2,βy=2.3288×10-3.When retarded cell proliferation was considered, based on the induced selection dynamic model (6) and the experimental data in [13], using extensive MCMC simulation again, we can obtain the estimated delay τ=5.1401 (Figure 2).MCMC analysis of parameterτ based on (6), (22), and the experimental data in [13]. (a) is the random series, and (b) is its histogram. The algorithm ran for 104 iterations with a burn-in of 3000 iterations. The initial conditions were τ=0.9.
(a)
(b)Using the estimated values, we can plot the best-fit solution by fitting model (1) and (6) to the experimental data, respectively, (Figure 3). From Figure 3, we find that there is a better simulation effect in model (6) than that in model (1). In fact, the sum of squares of the deviations (SSD) in (1) is SSD(1)=7.2943, whereas SSD(6)=4.6387×10-2 in (6). Note that SSD(6) is far less than SSD(1). These results quantitatively confirmed that the induced selection dynamic model (6) with time delay can provide a better fitting effect in long-term sphere culture of CSCs.Figure 3
Simulations of the dynamical behaviors of the frequency of CSCs. Experimental data are represented by open circles. The blue dashed line denotes the best fit of model (1), and the red solid line denotes the best fit of model (6).


## 3.1. Dynamic Analysis

The objective of this subsection is to analyze the dynamical behavior of (6). In order to explore the effect of the delay, we split this into two cases.


### 3.1.1. Case ofτ=0

In this case, we focus on the dynamic analysis if the delay is nonexistent; that is,τ=0 in (6). We start by studying the existence of nonnegative equilibria in the interval [0,1]. Let f(x)=a1+a2x+a3x2=0. Clearly, the discriminant of the quadratic equation is
(7)Δ=a22-4a1a3=(by-bx+bxβx-byβy)2+4bxbyβxβy>0.
Hence there are two different real roots for f(x)=0 if a3≠0. Furthermore, since
(8)f(1)=a1+a2+a3=-bxβx<0,
we know that there is a unique positive equilibrium x*∈(0,1) for (6) if a3≠0 (see Figures 1(a) and 1(b)). When a3=0 (see Figure 1(c)), it is clear that there is only one positive equilibrium
(9)x*=-a1a2=βxβx+βy<1.Illustrations of function image off(x)=0 under different cases ((a), (b), and (c)), and the phase diagram (d) of system (6). Here (a) a3>0, (b) a3<0, and (c) a3=0.
(a)
(b)
(c)
(d)The combination of the above results and the phase diagram (see Figure1(d)) of system (6) yields the following result.Proposition 1.
For system (6), when τ=0, a unique positive equilibrium x*∈(0,1) always exists, and it is globally asymptotically stable.


### 3.1.2. Case ofτ>0

In this case, we focus on the dynamic analysis if the delay is existent; that is,τ>0 in (6). Clearly, the unique positive equilibrium x*∈(0,1) still remains for (6) in spite of the delay. To study the stability of the equilibrium x*, we first translate x* to the origin. Let
(10)x~=x-x*.
Then (6) becomes, after replacing x~ by x again,
(11)dx(t)dt=(a2+2a3x*)x(t-τ)+a3x2(t-τ).The variational system of (11) at the origin is given by
(12)dx(t)dt=(a2+2a3x*)x(t-τ).
The characteristic equation of linear system (12) is given by
(13)λ-(a2+2a3x*)e-τλ=0.
If we let λ=α+iβ, then (13) becomes
(14)α-(a2+2a3x*)e-ταcosβτ=0,β+(a2+2a3x*)e-ταsinβτ=0.
By setting α=0 in (14), we have
(15)cosβτ=0,sinβτ=-βa2+2a3x*.Solving the first algebraic equation in (15), we have
(16)βτ=π2+2nπ,n=0,1,2,….
Substituting (16) into the second equation of (15), we have
(17)β=-(a2+2a3x*).
Now, substituting (17) into (16), we have
(18)τn=-π2(a2+2a3x*)-2nπa2+2a3x*,n=0,1,2,….Next, we computeα′(τn). Differentiating (13) with respect to τ, we have
(19)dλ(τ)dτ=-(a2+2a3x*)λe-τλ1+τ(a2+2a3x*)e-τλ,
thus,
(20)α′(τn)=Re(dλ(τ)dτ|τ=τn)=β21+τn2β2>0.Therefore, similar to [17], according to the results in [18, 19] or [15, Theorem  2.2], we can obtain the following results on (6).Proposition 2.
Supposeτ>0. Then system (6) has a Hopf bifurcation at
(21)τ=τn=-π2(a2+2a3x*)-2nπa2+2a3x*,n=0,1,2,….Furthermore, according to the results in [15, 20, 21], we have the following.Proposition 3.
Supposeτ>0 in (6). Then the unique positive equilibrium x*∈(0,1) is stable if 0<τ<τ0 and unstable if τ>τ0.


## 3.1.1. Case ofτ=0

In this case, we focus on the dynamic analysis if the delay is nonexistent; that is,τ=0 in (6). We start by studying the existence of nonnegative equilibria in the interval [0,1]. Let f(x)=a1+a2x+a3x2=0. Clearly, the discriminant of the quadratic equation is
(7)Δ=a22-4a1a3=(by-bx+bxβx-byβy)2+4bxbyβxβy>0.
Hence there are two different real roots for f(x)=0 if a3≠0. Furthermore, since
(8)f(1)=a1+a2+a3=-bxβx<0,
we know that there is a unique positive equilibrium x*∈(0,1) for (6) if a3≠0 (see Figures 1(a) and 1(b)). When a3=0 (see Figure 1(c)), it is clear that there is only one positive equilibrium
(9)x*=-a1a2=βxβx+βy<1.Illustrations of function image off(x)=0 under different cases ((a), (b), and (c)), and the phase diagram (d) of system (6). Here (a) a3>0, (b) a3<0, and (c) a3=0.
(a)
(b)
(c)
(d)The combination of the above results and the phase diagram (see Figure1(d)) of system (6) yields the following result.Proposition 1.
For system (6), when τ=0, a unique positive equilibrium x*∈(0,1) always exists, and it is globally asymptotically stable.


## 3.1.2. Case ofτ>0

In this case, we focus on the dynamic analysis if the delay is existent; that is,τ>0 in (6). Clearly, the unique positive equilibrium x*∈(0,1) still remains for (6) in spite of the delay. To study the stability of the equilibrium x*, we first translate x* to the origin. Let
(10)x~=x-x*.
Then (6) becomes, after replacing x~ by x again,
(11)dx(t)dt=(a2+2a3x*)x(t-τ)+a3x2(t-τ).The variational system of (11) at the origin is given by
(12)dx(t)dt=(a2+2a3x*)x(t-τ).
The characteristic equation of linear system (12) is given by
(13)λ-(a2+2a3x*)e-τλ=0.
If we let λ=α+iβ, then (13) becomes
(14)α-(a2+2a3x*)e-ταcosβτ=0,β+(a2+2a3x*)e-ταsinβτ=0.
By setting α=0 in (14), we have
(15)cosβτ=0,sinβτ=-βa2+2a3x*.Solving the first algebraic equation in (15), we have
(16)βτ=π2+2nπ,n=0,1,2,….
Substituting (16) into the second equation of (15), we have
(17)β=-(a2+2a3x*).
Now, substituting (17) into (16), we have
(18)τn=-π2(a2+2a3x*)-2nπa2+2a3x*,n=0,1,2,….Next, we computeα′(τn). Differentiating (13) with respect to τ, we have
(19)dλ(τ)dτ=-(a2+2a3x*)λe-τλ1+τ(a2+2a3x*)e-τλ,
thus,
(20)α′(τn)=Re(dλ(τ)dτ|τ=τn)=β21+τn2β2>0.Therefore, similar to [17], according to the results in [18, 19] or [15, Theorem  2.2], we can obtain the following results on (6).Proposition 2.
Supposeτ>0. Then system (6) has a Hopf bifurcation at
(21)τ=τn=-π2(a2+2a3x*)-2nπa2+2a3x*,n=0,1,2,….Furthermore, according to the results in [15, 20, 21], we have the following.Proposition 3.
Supposeτ>0 in (6). Then the unique positive equilibrium x*∈(0,1) is stable if 0<τ<τ0 and unstable if τ>τ0.


## 3.2. Numerical Simulations

For model (1), we designed a long-term sphere culture of human breast cancer MCF-7 stem cells [13]. Based on the experimental data, using an adaptive Metropolis-Hastings (M-H) algorithm to carry out an extensive Markov-chain Monte-Carlo (MCMC) simulation, we obtained the estimated parameter values as follows:
(22)bx=2.7506×10-1,by=3.2635×10-1,βx=1.4407×10-2,βy=2.3288×10-3.When retarded cell proliferation was considered, based on the induced selection dynamic model (6) and the experimental data in [13], using extensive MCMC simulation again, we can obtain the estimated delay τ=5.1401 (Figure 2).MCMC analysis of parameterτ based on (6), (22), and the experimental data in [13]. (a) is the random series, and (b) is its histogram. The algorithm ran for 104 iterations with a burn-in of 3000 iterations. The initial conditions were τ=0.9.
(a)
(b)Using the estimated values, we can plot the best-fit solution by fitting model (1) and (6) to the experimental data, respectively, (Figure 3). From Figure 3, we find that there is a better simulation effect in model (6) than that in model (1). In fact, the sum of squares of the deviations (SSD) in (1) is SSD(1)=7.2943, whereas SSD(6)=4.6387×10-2 in (6). Note that SSD(6) is far less than SSD(1). These results quantitatively confirmed that the induced selection dynamic model (6) with time delay can provide a better fitting effect in long-term sphere culture of CSCs.Figure 3
Simulations of the dynamical behaviors of the frequency of CSCs. Experimental data are represented by open circles. The blue dashed line denotes the best fit of model (1), and the red solid line denotes the best fit of model (6).


## 4. Conclusions

In order to demonstrate the interesting facts about the structural heterogeneity of cancer (the stable ratio between CSCs and DCCs), many studies have been reported because it is helpful for the cancer community to elucidate the controversy about the CSC hypothesis and the clone evolution theory of cancer [10, 13, 14] and references cited therein. In the present paper, a selection dynamic model with time delay is proposed, and its dynamical behavior is studied. Based on the theoretical analysis and numerical simulations, we can conclude the following predictive conclusions.(i) The maturation of individual cells may produce a significant effect on the dynamic behavior of the selection dynamics. When the delay is nonexistent, the frequency of CSCs will tend to a stable size because the unique positive equilibrium is globally asymptotically stable (Proposition1). Conversely, if the delay is existent, the unique positive equilibrium may not always maintain its stability and a Hopf bifurcation may be induced (Propositions 2 and 3); that is, an oscillated phenomenon may be induced by the maturation of individual cells.(ii) Since the induced selection dynamic model (6) with time delay can provide a better fitting effect in long-term sphere culture of CSCs (Figure 3), it is reasonable to consider the delay effect in the future study on the dynamics of non-genetic heterogeneity of clonal cell populations.Since mathematical models can be at best approximate the behavior of real biological process, the results presented here may extend those studies on the structural heterogeneity of cancer. Note that distributed delay may be more tractable and realistic than discrete delay in the applications of biology. Hence it is a worthwhile study in future work to better understand these topics based on the idea of [15].



---

*Source: 606250-2013-03-20.xml*