The overshoot and phenotypic equilibrium in characterizing cancer dynamics of reversible phenotypic plasticity
Abstract
The paradigm of phenotypic plasticity indicates reversible relations of different cancer cell phenotypes, which extends the cellular hierarchy proposed by the classical cancer stem cell (CSC) theory. Since it is still questionable if the phenotypic plasticity is a crucial improvement to the hierarchical model or just a minor extension to it, it is worthwhile to explore the dynamic behavior characterizing the reversible phenotypic plasticity. In this study we compare the hierarchical model and the reversible model in predicting the cellstate dynamics observed in biological experiments. Our results show that the hierarchical model shows significant disadvantages over the reversible model in describing both longterm stability (phenotypic equilibrium) and shortterm transient dynamics (overshoot) of cancer cells. In a very specific case in which the total growth of population due to each cell type is identical, the hierarchical model predicts neither phenotypic equilibrium nor overshoot, whereas the reversible model succeeds in predicting both of them. Even though the performance of the hierarchical model can be improved by relaxing the specific assumption, its prediction to the phenotypic equilibrium strongly depends on a precondition that may be unrealistic in biological experiments, and it also fails to capture the overshoot of CSCs. By comparison, it is more likely for the reversible model to correctly describe the stability of the phenotypic mixture and various types of overshoot behavior.
Xiufang Chen, Yue Wang, Tianquan Feng, Ming Yi, Xingan Zhang, Da Zhou

School of Mathematical Sciences, Xiamen University, Xiamen 361005, P.R. China
(Cocorresponding Author, zhouda@xmu.edu.cn) 
School of Mathematics and Statistics, Central China Normal University, Wuhan 430079, P. R. China
(Cocorresponding Author, zhangxinan@hotmail.com) 
Department of Applied Mathematics, University of Washington, Seattle, WA 98195, USA

School of Teachers’ Education, Nanjing Normal University, Nanjing 210023, China

Department of Physics, College of Science, Huazhong Agricultural University, Wuhan, Hubei 430070, China

Key Laboratory of Magnetic Resonance in Biological Systems, Wuhan Institute of Physics and Mathematics, Chinese Academy of Sciences, Wuhan, 430071, P. R. China

School of Computer Science and Information Engineering, Qilu Institute of Technology, Jinan, Shandong 250000, China
1 Introduction
The cancer stem cell theory has provided a hierarchical model of how diverse cancer cells being organized [1, 2, 3]. Similar to the stem cell theory in normal tissues, this hierarchical model assumes that a small number of stemlike cancer cells (termed cancer stem cells, CSCs) are capable of selfrenewal and differentiation into other more committed cancer cells (termed nonstem cancer cells, NSCCs) but not vice versa. That is, CSCs are thought to be at the apex of this cellular hierarchy. However, some recent researches may extend this unidirectional relation of cancer cells. It has been reported that cancer cells can convert from NSCC phenotype to CSC phenotype (e.g. breast cancer [4, 5], melanoma [6], colon cancer [7] and glioblastoma multiforme [8]). Furthermore, the interconversions among cancer cell phenotypes have also been found (breast cancer [9]). All these works indicate reversible relations of different cancer cells (termed phenotypic plasticity, which has long been an issue concerned in bacterial populations [10]).
Very recently special attention has already been paid to the reversible cancer models with the phenotypic plasticity by theoreticians. In particular, dos Santos and da Silva explained the variable frequencies of CSCs in tumors by establishing a model with stochastic cell plasticity [11, 12]. Leder et al investigated the model of reversible conversions between the stemlike resistant cells (SLRCs) and the differentiated sensitive cells (DSCs) in glioblastomas [13]. Wang et al showed how tumor heterogeneity arises in the model of cooperating CSC hierarchy with cell plasticity [14]. Zhou et al showed that the dedifferentiation from NSCC phenotype to CSC phenotype was essential for explaining the transient increase of the minority populations of CSCs observed in cancer cell lines [15, 16]. Chen et al studied stochastic models that capture the transitions between endocrine therapy responsive and resistant states of breast cancer cells [17]. Zhou et al investigated nonequilibrium dynamics with phenotype transitions of cancer cells [18]. Jilkine and Gutenkunst studied the effect of dedifferentiation on time to mutation acquisition in cancers [19].
However, it is still questionable if the phenotypic plasticity is a crucial improvement to the hierarchical model or just a minor extension to it [21, 22]. Thus a rigorous analysis on the characteristics owned by the model with the phenotypic plasticity is necessary for model validation. In this study, we try to investigate this issue by giving a comparative study of the reversible model and the hierarchical model. Note that Gupta et al studied SUM149 and SUM159 breast cancer cell lines [9], in which two interesting phenomena were observed: phenotypic equilibrium and overshoot. For the phenotypic equilibrium, they found that the breast cancer cell lines will tend to a stable phenotypic mixture over time regardless of initial states. Similar results on the phenotypic equilibrium were also reported in [23, 24, 7]. For the overshoot, they found that shortly after the cell sorting, the proportion of the minority subpopulation will increase transiently above the final equilibrium level, and then decrease to it (Fig. 1). Enlightened by these observations, we show that the reversible model is more capable of capturing the phenotypic equilibrium and overshoot than the hierarchical model. Under a very specific assumption that the total population growth due to each cell phenotype is identical, the hierarchical model can perform neither phenotypic equilibrium nor overshoot, whereas the reversible model succeeds in predicting both the phenotypic equilibrium and three types of overshoots (asynchronous, synchronous and oscillating overshoots). Even though the performance of the hierarchical model can be improved by relaxing the specific assumption, it is still not good enough to correctly capture the cellstate dynamics. On one hand, the phenotypic equilibrium predicted by the hierarchical model strongly depends on the condition that the selfcontributed growth rate by CSCs is faster than that of more committed cancer cells. However, it has been reported that this condition cannot be satisfied in some cancers [25, 26]. On the other hand, by fitting the two models to experimental data in [9], the hierarchical model can only fit the overshoot of NSCCs, but it cannot capture the overshoot of CSCs. By contrast, the reversible model can fit both of them. Therefore, our results generally imply that the reversible model shows distinct advantages over the hierarchical model in predicting both longterm and transient dynamics of cancer cells.
The paper is organized as follows. The mathematical model is presented in Section 2. Main results are shown in Section 3, where we give a comparative study of the reversible and hierarchical models. Conclusions are presented in Section 4.
2 Model
In this section we describe the assumptions of our model. We model cancer as population dynamics of cancer cells, each cancer cell can be assigned to one of the following cell phenotypes: , , , …, and . Let , ,…, be the cell numbers of , , …, at time respectively, then the model can be described as a cellular dynamical system of . In particular, the model will reduce to a twophenotypic model when . That is, besides CSCphenotype, all the other cancer cells are grouped into one whole NSCCphenotype. This CSCNSCC model has extensively been investigated in previous literature [21, 11, 12, 15, 13, 14]. Here we pay more attention to the multiphenotypic case (). As a starting point, the case with should be a natural and favorable candidate for theoretical analysis. Meanwhile, motivated by Gupta et al [9] where three phenotypes (stemlike, basal and luminal) were identified, we are concerned about  model in light of both theoretical and experimental reasons. Unlike the threephenotypic cell lineages investigated in [27, 28, 29] where cancer stem cells, progenitor cells, and terminally differentiated cells are hierarchically cascaded, here we are more interested in the cell population structure investigated in [9]. In their work, CSCs can differentiate into and , respectively. They were not assumed to be cascaded. Moreover, by the phenotypic plasticity, they can interconvert into each other and dedifferentiate into CSCs.
We now present the cellular processes included in our model. For , it can not only divide symmetrically into two identical daughters, but also divide asymmetrically into and ( or ) [30, 31, 32]. That is,

.

.

^{1}^{1}1It should be pointed out that, the symmetric differentiation “” is also an important type of transition [33] that should be considered in the model. However, it can be shown that the main results we are interested in this work can still hold (with only minor changes in proofs) even if adding this transition into Eq. (1). Therefore, to keep the model as minimal as possible, the symmetric differentiation is not included in presented model..
For and , besides the phenotypic plasticity, they can also perform symmetric divisions and cell death. That is,

;

;

;

;

;

;

;

.
Based on the above cellular processes, the dynamics of can be captured by the following ordinary differential equations (ODEs) {linenomath*}
(1) 
By letting {linenomath*}
(2) 
Eq. (1) can be expressed as {linenomath*}
(3) 
For this model, we have the following remarks:
(1)
Each element of can be seen as the per
capita growth rate of phenotype contributed by
phenotype. Let us take the first row of as an example.
is the symmetric division rate of , thus
the growth rate of due to themselves can be expressed
as ; and are the
dedifferentiation rates
from and respectively,
so the growth rate of
due to and is .
In this way, it is easy to explain the biological meaning of
“”, i.e. the rate of change of
is the sum of the growth rates contributed by all the three cell phenotypes in the
population. Similarly, we can explain the other two equations in Eq. (1).
(2)
Note that , , and are the parameters associated with the phenotypic plasticity,
the model will reduce to the hierarchical model by letting them be zero. Accordingly,
will become to a lowertriangular matrix
{linenomath*}
(4) 
(3) Let be the total number of the population, then {linenomath*}
(5) 
Note that the coefficient of each on righthand side of Eq. (5) is just the the sum of the corresponding column in , i.e. the per capita growth rate of the whole population due to . In particular, when {linenomath*}
(6) 
the growth rate of the whole population due to each phenotype will be identical. This special case implies meaningful biological significance and has been reported in some cancer cell lines (breast cancer [9] and colon cancer [15]). In this case Eq. (5) will reduce to a very simple form {linenomath*}
that is, will grow exponentially with constant rate. We will show that, even in this special case the reversible model and the hierarchical models are quite different in predicting the phenotypic equilibrium and overshoot. For convenience, we denote {linenomath*}
(7) 
(8) 
(9) 
(4) Note that all the parameters are assumed to be constant. In other words, there is no feedback control included in presented model. This is of course a simplification of the biological facts. To model real biological systems, feedback mechanisms are indispensable [34, 35, 36, 28]. However, we will show that our model is not overly simplistic to characterize the phenotypic plasticity in comparison to the cellular hierarchy.
In the remaining of this section, we will convert the cell number equations Eq. (1) into the following cell proportion equations Eq. (10). This is because in reality, relative cell numbers, rather than absolute cell numbers, are usually measured by fluorescenceactivated cell sorting (FACS) experiments. Let , and , by variable substitutions in Eq. (1) we have {linenomath*}
(10) 
Note that , the freedom of Eq. (10) is actually two. By letting , Eq. (10) can be transformed into {linenomath*}
(11) 
where {linenomath*}
(12) 
(13) 
(14) 
(15) 
When , i.e. the phenotypic conversion rates are set to be zero, we have {linenomath*}
(16) 
where {linenomath*}
(17) 
(18) 
(19) 
(20) 
We term Eq. (11) the reversible model, and term Eq. (16) the hierarchical model. The comparison of these two models is the main task in the following section.
One trivial difference between the two models happens when there is no CSC in the initial population. In this case, since in the hierarchical model CSCs cannot be spontaneously created by other cell phenotypes, there is no CSC in the population all the time, whereas newborn CSCs can be generated from NSCCs in the reversible model. However, in reality it is quite impossible to completely purify cancer cell lines without any CSCs by cell sorting. Based on this fact, we only discuss the cases with positive initial states of CSCs afterwards. Both longterm stable behavior and shortterm transient dynamics are investigated when comparing the two models. Enlightened by [9], we are particularly concerned about their predictions to the phenotypic equilibrium and overshoot.
2.1 Special case
Before exploring the general models in Eqs. (11) and (16), we first discuss a very specific case with . Investigating this special case can provide some valuable insights to the general analysis.
For the hierarchical model in this case, {linenomath*}
(21) 
It is easy to give their explicit solutions {linenomath*}
(22) 
where and are the initial states. So we have

Since and are both monotonic functions of , they cannot perform overshoot.

Note that {linenomath*}
CSCs proportion will eventually be zero, and the equilibrium proportions of depends on the initial states. Note that the phenotypic equilibrium corresponds to the stabilization of the phenotypic mixture that is independent of the initial states [9], the hierarchical model fails to predict the phenotypic equilibrium.
We turn our attention to the reversible model {linenomath*}
(23) 
For the phenotypic equilibrium, we have the following theorem:
Theorem 1.
There exists unique stable fixed point in Eq. (23), where
The proof is put in A. Theorem 1 indicates that the phenotypic equilibrium naturally holds in the special reversible model. We now consider whether this model can perform overshoot. The main results are listed in Table 1 (see B for more details). In particular, Fig. 2 illustrates the predictions by case 1 ():

Nonovershoot (left panel). Both and are monotonic functions.

Asynchronous overshoot (middle panel). Either or can perform overshoot (The figure shows an example of overshoot), but they cannot perform overshoot simultaneously.

Synchronous overshoot (right panel). and can perform overshoot simultaneously.
Case 2 ( and is diagonalizable) is the star case (see B), overshoot can never happen. For case 3 ( and is not diagonalizable), it belongs to the nodal case (see B), both asynchronous and synchronous overshoots can happen. Besides, in case 4 (), the model performs damped oscillatory dynamics, corresponding to an interesting oscillating overshoot (Fig. 3).
According to the above comparisons between Eqs. (21) and (23), the reversible model shows richer dynamics than the hierarchical model in performing both the phenotypic equilibrium and overshoot.
Parameters  Behavior  

1  no overshoot, asynchronous and synchronous overshoots  
2  and is diagonalizable  no overshoot 
3  and is not diagonalizable  asynchronous and synchronous overshoots 
4  oscillating overshoot 
2.2 General case
We now remove the constraint of and further consider the general case in this section. For the general hierarchical model Eq. (16), the following theorem states the condition under which the phenotypic mixture will finally stabilize:
Theorem 2.
There exist three singular points in Eq. (16): , and , where {linenomath*}
Assume that and , is the unique stable fixed point.
The proof is put in C. Note that , and are the diagonal elements of the matrix in Eq. (4), representing the per capita growth rates due to the corresponding cell phenotypes themselves. Theorem 2 thus indicates that the general hierarchical model is capable of predicting the phenotypic equilibrium only provided that the selfcontributed growth rate by CSCs is larger than that of more committed cancer cells. However, this condition may not be satisfied in reality. It has been reported that more committed cancer cells can have faster cycling time than stemlike cancer cells in some cancers [25, 26, 14]. In contrast, the following theorem shows the phenotypic equilibrium of the general reversible model Eq. (10):
Theorem 3.
The proof is put in D. Theorem 2 shows that the result in the special case (Theorem 1) can be extended to the general case. By comparison of Theorems 2 and 2, we can see that it is more likely for the reversible model to correctly capture the phenotypic equilibrium.
In the remaining of this section, we discuss the overshoot performances of the two models by fitting them to the cellstate dynamics in [9]. We fitted Eqs. (16) and (10) to the SUM159 data shown in the Fig. 3 of [9]. From Fig. 4, we can see that the reversible model fits the data better than the hierarchical model. In particular, note that both the overshoots of luminal cells and stemlike cells were presented in the SUM159 data (see the left panel and right panel in the Fig. 3 of [9] respectively). Our result shows that, the reversible model can capture both of them, whereas the hierarchical model can only fit the overshoot of luminal cells but it fails to capture the overshoot of stemlike cells. In the left panel of Fig. 4, both the blue solid line (hierarchical model) and blue discrete trianglemarker line (reversible model) fit the overshoot of luminal cells. Meanwhile, in the right panel of Fig. 4, the red solid line (hierarchical model) shows a monotonic way up to the final equilibrium (nonovershoot), but the red discrete trianglemarker line (reversible model) predict the overshoot of stemlike cells.
Actually, it has already been reported in previous literature that the overshoot of NSCCs proportion can happen in classical CSC model [37], but the overshoot of CSCs proportion we think has added significance. Note that in the hierarchical model the increase of CSCs relies only on their selfrenewals. When the initial fraction of CSCs is very limited, the transient increase of CSCs proportion should be very slow (constraint by the limitation of cell division cycle), so intuitively it is quite unlikely for the hierarchical model to capture the overshoot of CSCs proportion. In contrast, the dedifferentiation from NSCCs can effectively speed up the accumulation of CSCs, it is more likely for the reversible model to capture the overshoot of CSCs proportion. In other words, the overshoot of CSCs should be a better candidate than that of NSCCs to characterize the phenotypic plasticity. Our numerical result just supports this idea.
3 Conclusions
In this paper, we have tried to characterize the phenotypic plasticity by giving a comparative study between the reversible model and hierarchical model. According to our results, the reversible model shows richer dynamics and better predictions to experimental phenomena than the hierarchical model.
It should be noted that, to focus our attention to the reversible phenotypic plasticity, the presented model does not include the biologically complex mechanisms such as feedback controls and time delays. It is undeniable that, by including the feedback and delay mechanisms, the hierarchical model can also show complex dynamics that are beyond the reach of the presented model [38, 39, 40, 34, 28]. For example, Bernard et al showed that oscillations can occur in their hierarchical model with feedbacks and delays [38]. Stiehl et al showed that the hierarchical model of hematopoiesis with nonlinear feedback can reproduce the overshoot and damped oscillations observed in clinical data [39]. We also think that the interplay between the feedback (or delay) mechanisms and the phenotypic plasticity should be a very interesting issue. As a staring point, it is suggested by our presented model that, the phenotypic plasticity facilitates the heterogeneity of cancer in two ways: it helps to sustain the coexistence of multiple phenotypes in cancer, and it also accelerates the recovery process of cancer stem cells in the population. These findings may shed some lights on the researches of the models incorporating the biological complexities with the phenotypic plasticity.
Moreover, the presented model is deterministic, whereby stochastic effects are not accounted for. Therefore, the branching model [41] concerning both the stochastic phenotypic conversions and proliferations of cancer cells is another interesting research direction. Finally, it is worth noting that, in addition to the longterm stabilization, the shortterm transient dynamics also deserve special attention in the study of biological population dynamics [42]. In this work we have seen the significance of the overshoot in characterizing the phenotypic plasticity. It is interesting to explore more types of transient dynamics and their biological significance in future researches.
Acknowledgements
X. C. and X. Z. acknowledge the support by NSFC (No. 11371161). M. Y. acknowledges NSFC (Nos. 11275259 and 91330113). D. Z. acknowledges the generous sponsorship from the National Natural Sciences Foundation of China (No. 11401499), the Natural Science Foundation of Fujian Province of China (No. 2015J05016),and the Fundamental Research Funds for the Central Universities in China (No. 20720140524).
Appendix A Proof of Theorem 1
Proof.
Let {linenomath*}
(24) 
where , , , . Note that {linenomath*}
it is easy to see that is the only fixed point. To proof the stability of , let , , then Eq. (23) is transformed to {linenomath*}
(25) 
Let and {linenomath*}
(26) 
Eq. (25) can be expressed as {linenomath*}
becomes the only fixed point and its characteristic equation is {linenomath*}
(27) 
Since {linenomath*}
the real parts of and are both negative. Hence is stable in Eq. (25), and then is the unique stable fixed point of Eq. (23). The proof is completed. ∎
Appendix B The overshoot predicted by the special reversible model
Let
(28) 
be the discriminant of (27). Since it has been shown that the real parts of and are negative, the special reversible model Eq. (23) can be classified into the following four cases (For the terminologies used below, please refer to Chapter 5 in [43]):

. In this case the matrix in Eq. (26) has two unequal negative eigenvalues. Thus is the nodal point of Eq. (25), correspondingly is the nodal point of Eq. (23). The solution can be expressed as {linenomath*}
(29) where
(30) is a matrix determined by the coefficients of Eq. (25), both and are determined by the initial states. Different types of transient dynamics can happen for different parameter ranges. By numerical simulations, it is shown in Fig. 2 that the model can perform nonovershoot, asynchronous overshoot and synchronous overshoot respectively.

and is diagonalizable. Since is the star point in this case (Fig. B.5). There is no overshoot of and .

and is not diagonalizable. In this case is also called nodal point. From Fig. B.6, we can see that the model can perform both asynchronous and synchronous overshoots. If , the model can perform asynchronous overshoot (Fig. B.6(a)). If , the model can perform synchronous overshoot (Fig. B.6(b)).

. has two conjugate complex eigenvalues with negative real part. is the focal point, resulting in oscillating overshoot (Fig. 3).
Appendix C Proof of Theorem 2
Theorem (page 2, Theorem 2).
There exist three singular points in Eq. (16): , and , where {linenomath*}
Assume that and , is the unique stable fixed point.
Proof.
To obtain the fixed points of Eq. (16), we consider {linenomath*}
(31) 
For the first equation
,
we have either or
(i) When , based on the second equation we have or
that is, we have two fixed points
{linenomath*}
(ii) When the fixed point is {linenomath*}
where {linenomath*}
Since and , we have and Let{linenomath*}
then
The Jacobian matrix of Eq. (16) is given by {linenomath*}
Let us first consider the stability of , {linenomath*}
its characteristic equation is {linenomath*}
Note that
then , implying that is unstable.
Similarly, for
, is also unstable.
For , {linenomath*}
where {linenomath*}
Then we have {linenomath*}
Then {linenomath*}
It is shown that is a stable nodal point, and then it is the unique stable fixed point in Eq. (16). The proof is completed. ∎
Appendix D Proof of Theorem 2
Proof.
Our proof is valid for general multiphenotypic cases, so in this section the dimension of the matrix in Eq (2) is not restricted to three any more, but in general. Without loss of generality, we only prove the theorem when is positive (all the elements of are positive) ^{3}^{3}3Note that the offdiagonal elements of are positive, ( is the identity matrix) can be positive provided is large enough. It is easy to show that when is an eigenvalue of , is an eigenvalue of accordingly. That is, the eigenvalue structure of is the same as that of . So we just need to consider instead of .. According to the wellknown PerronFrobenius theory (see chapter 1 in [44]), has a real eigenvalue satisfying for any other eigenvalue of (called the PerronFrobenius eigenvalue) and is simple (a simple root of the characteristic equation of ).
Since is simple, the solution of Eq. (1) can be expressed as {linenomath*}
(32) 
where are the different eigenvalues of , is the algebraic multiplicity of , is the normalized () right eigenvector of , is the corresponding eigenvector of , is determined by initial states. Suppose , since Re, {linenomath*}
Thus {linenomath*}
Noticing that PerronFrobenius theory ensures that all the components of are positive, i.e. .
Before we complete the proof, we need to discuss . In this case, the above method does not work. However, since fluctuations are inevitable in real world, will hardly happen in reality. To show this, let in Eq. (32) {linenomath*}
This is a linear equation of . For we have {linenomath*}
where , is just with its first column replaced by . It is easy to add a small perturbation to , so that all the columns of are linear independent, hence