1 2 3 4 5 6 7 8 9 10 11 12 13 14 Q1 15 Q3 16 17 18 19 20 21 22 23 24 25 Q4 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 Q5 56 57 58 59 60 61 62 63 64 65 66

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

A multi-phenotypic cancer model with cell plasticity Da Zhou a, Yue Wang b, Bin Wu c a b c

School of Mathematical Sciences, Xiamen University, Xiamen 361005, PR China Department of Applied Mathematics, University of Washington, Seattle, WA 98195, USA Evolutionary Theory Group, Max-Planck-Institute for Evolutionary Biology, August-Thienemann-Straβe 2, 24306 Plö n, Germany

H I G H L I G H T S

We established a multi-phenotypic cancer model with cell plasticity. Prove the stability of the nonlinear model of multi-phenotypic proportions. Cell plasticity explains the transient increase of cancer stem cells proportion. Overshooting of CSCs proportion arises only in multi-phenotypic case.

art ic l e i nf o

a b s t r a c t

Article history: Received 28 November 2013 Received in revised form 27 March 2014 Accepted 30 April 2014

The conventional cancer stem cell (CSC) theory indicates a hierarchy of CSCs and non-stem cancer cells (NSCCs), that is, CSCs can differentiate into NSCCs but not vice versa. However, an alternative paradigm of CSC theory with reversible cell plasticity among cancer cells has received much attention very recently. Here we present a generalized multi-phenotypic cancer model by integrating cell plasticity with the conventional hierarchical structure of cancer cells. We prove that under very weak assumption, the nonlinear dynamics of multi-phenotypic proportions in our model has only one stable steady state and no stable limit cycle. This result theoretically explains the phenotypic equilibrium phenomena reported in various cancer cell lines. Furthermore, according to the transient analysis of our model, it is found that cancer cell plasticity plays an essential role in maintaining the phenotypic diversity in cancer especially during the transient dynamics. Two biological examples with experimental data show that the phenotypic conversions from NCSSs to CSCs greatly contribute to the transient growth of CSCs proportion shortly after the drastic reduction of it. In particular, an interesting overshooting phenomenon of CSCs proportion arises in three-phenotypic example. Our work may pave the way for modeling and analyzing the multi-phenotypic cell population dynamics with cell plasticity. & 2014 Published by Elsevier Ltd.

Keywords: Cancer stem cell theory Cell-state conversions Multi-phenotype model

1. Introduction The cancer stem cell (CSC) hypothesis (Reya et al., 2001; Jordan et al., 2006) states that tumors or hematological cancers arise from a small number of stem-like cancer cells with the abilities of selfrenewal and differentiation into other non-stem cancer cells (NSCCs). That is, the conventional cancer stem cell theory suggests a cellular hierarchy where CSCs are at the apex (Dalerba et al., 2007). Based on this paradigm, cancer stem cell models were widely investigated in previous literature on theoretical biology (Colijn and Mackey, 2005; Johnston et al., 2007; Boman et al., 2007; Dingli et al., 2007; Johnston et al., 2010; Antal and Krapivsky, 2011; Sottoriva et al., 2011; Werner et al., 2011; Stiehl and Marciniak-Czochra, 2012; Molina-Peña and Álvarez, 2012).

E-mail address: [email protected] (D. Zhou).

However, recent studies have highlighted the complexities and challenges in the evolving concept of CSC (Nguyen et al., 2012; Visvader and Lindeman, 2012). In particular, it was reported that reversible phenotypic changes can occur between stem-like cancer cells and more differentiated cancer cells. Meyer et al. showed that interconversion occurred both in vivo and in vitro between noninvasive, epithelial-like CD44 þ CD24 þ cells and invasive, mesenchymal CD44 þ CD24 cells in breast cancer (Meyer et al., 2009). Chaffer et al. (2011) showed both in vivo and in vitro that transformed (oncogenic) CD44lo-HMECs (human mammary epithelial cells) could spontaneously convert to CD44hi-CSCs. The results by Quintana et al. (2010) indicated that phenotypically diverse cancer cells in both primary and metastatic melanomas can undergo reversible phenotypic conversions in vivo. Scafﬁdi et al. showed that stem-like cancer cells can be generated in vitro from transformed (oncogenic) ﬁbroblasts during neoplastic transformation (Scafﬁdi and Misteli, 2011). The conversion from NSCCs to CSCs was also in situ visualized

http://dx.doi.org/10.1016/j.jtbi.2014.04.039 0022-5193/& 2014 Published by Elsevier Ltd.

Please cite this article as: Zhou, D., et al., A multi-phenotypic cancer model with cell plasticity. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.04.039i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87

2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

D. Zhou et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

by Yang et al. (2012) in SW620 colon cancer cell line. Besides, the phenotypic transitions between different NSCCs in breast cancer were reported both in vivo and in vitro (Gupta et al., 2011). These various types of phenotypic conversions among cancer cells, also known as cancer cell plasticity mechanisms (French and Clarkson, 2012), provide new thinking about the CSC hypothesis and related therapeutic strategy (Pisco et al., 2013). Special attention has recently been paid to the mathematical models concerning cancer cell plasticity. Gupta et al. (2011) introduced a simple Markov chain model of stochastic transitions among stem-like, basal and luminal cells in breast cancer cell lines. Zapperi and La Porta (2012) compared mathematical models for cancer cell proliferation that take phenotypic switching and imperfect biomarker into account. A series works by dos Santos and da Silva (2013a,b) developed a model with the effects of stochastic noise and cell plasticity for explaining the variable frequencies of CSCs in tumors. Zhou et al. (2013) compared the transient dynamics of the bidirectional and unidirectional models of CSCs and NSCCs. Wang et al. (2014) showed that tumor heterogeneity may exist in the model with both CSC hierarchy and cell plasticity. Leder et al.'s (2014) model described the reversible phenotypic interconversions between the stem-like resistant cells (SLRCs) and the differentiated sensitive cells (DSCs) in glioblastomas, which revealed optimized radiation dosing schedules. These works demonstrated that cell plasticity provides new insight into cancer cell population dynamics. To further explore how cell plasticity challenges the hierarchical cancer stem cell scenario, and in particular how cell plasticity inﬂuences tumor heterogeneity, one should incorporate cell plasticity into the development of cancer which is full of biological complexities. One particular and crucial complexity arises from highly diverse phenotypes in the population of cancer cells. The aforementioned mathematical models were mainly focused on the relation between CSCs and NSCCs (Zapperi and La Porta, 2012; dos Santos and da Silva, 2013a,b; Zhou et al., 2013; Wang et al., 2014; Leder et al., 2014), that is, their models simply classiﬁed the cancer cells into two “opposite” phenotypes. This simpliﬁcation is effective for studying the reversible conversions between NSCCs and CSCs, but covers up various phenotype switchings between different cancer cells that are worthy of studying. As an exceptional case, Gupta et al.'s (2011) model consists of three phenotypes (stem-like, basal and luminal cells), but rigorous mathematical analysis for general multi-phenotypic models with cell plasticity is still lacking. In this study, we try to provide a multi-phenotypic framework for integrating cell plasticity with conventional growth model of cancer cells. A generalized model comprising 1 þ m cellular phenotypes (one CSC phenotype and m different NSCC phenotypes) is put forward. Besides cell-state conversions from NSCCs to CSCs and phenotype switchings between different NSCC phenotypes, the cellular processes in classical cancer stem cell models (i.e., asymmetric cell division, symmetric cell division and cell death) are also included in our model. When the cell population size is not large enough and subject to stochastic ﬂuctuations, our model is formulated by a continuous-time high-dimensional Markov process. In the limit of large population size, the model can be governed by a system of linear ordinary differential equations (ODEs). Moreover, to investigate the dynamics of phenotypic proportions, the population model is converted into a nonlinear frequency one. It is shown that under very weak assumption, the nonlinear frequency model has only one stable steady state and no stable limit cycle. Not only does this result theoretically explain the phenotypic equilibrium phenomena reported in various cancer cell lines (Chaffer et al., 2011; Yang et al., 2012; Gupta et al., 2011; Iliopoulos et al., 2011), but it is also predicted that the phenotypic equilibrium should be universal in the population of multi-phenotypic cancer cells.

Furthermore, it is also found that cancer cell plasticity greatly inﬂuences the transient proportions of cell phenotypes. In particular, two concrete examples with experimental data are presented, showing that the cell-state conversions play an essential role in the transient growth of CSCs proportion shortly after the drastic reduction of it. In particular, an interesting overshooting phenomenon of CSCs proportion arises in the S–B–L model with cell plasticity. Note that two-phenotypic models never perform overshooting (Zhou et al., 2013; Jia et al.), overshooting can be a result of interplay between cell plasticity and diversity of phenotype. Moreover, it has been investigated in ecology and population genetics that phenotypic variability can serve as an advantageous strategy for biological populations in ﬂuctuating environments (Wolf et al., 2005; Kussell and Leibler, 2005; Lu et al., 2007; Acar et al., 2008; Kaneko, 2012), our ﬁndings thus enrich this idea that cell plasticity as a surviving strategy might be more essential in maintaining the phenotype diversity (heterogeneity) of cancer especially during transient dynamics. The paper is organized as follows. The model framework is formulated in Section 2. In Section 3, we investigate the frequency model and phenotypic equilibrium. The roles of cell-state conversions in transient dynamics are discussed in Section 4. Conclusions are presented in Section 5. 2. Model description 2.1. Assumptions This section describes the assumptions of the model investigated in this study. Consider a population of cancer cells comprising 1 þ m phenotypes: CSC represents cancer stem cell, and NSCC1 , NSCC2 ; …; NSCCm represent m different phenotypes of non-stem cancer cells. In this model, cell plasticity is integrated with the growth model of cancer cells. According to conventional cancer stem cell scenario, not only can CSC divide asymmetrically into two unequal daughter cells (one CSC and one NSCC) (Reya et al., 2001), but it can also divide symmetrically into two daughter CSCs (Todaro et al., 2010).1 So for CSC we assume that α Symmetric division: CSC⟶ CSC þ CSC. α Asymmetric division: CSC⟶CSC þ NSCCj ð1 r j r mÞ. α Cell death: CSC⟶∅. 00

0j

0

For NSCCs, besides the symmetric division, two types of cell plasticity mechanisms that have been reported in previous biological literature are included in our model, i.e., the cell-state conversions from NSCCs to CSCs (termed as de-differentiation) (Yang et al., 2012) and phenotype switchings between different NCSSs (Gupta et al., 2011). In this way, for NSCCi ð1 r i rmÞ we assume that

αii

Symmetric division: NSCCi ⟶NSCCi þ NSCCi . αi0 De-differentiation: NSCCi ⟶CSC. αij Phenotype switching: NSCCi ⟶NSCCj ði a jÞ. αi Cell death: NSCCi ⟶∅.

We list the elements of the model in Table 1 (see Fig. 1 for the example of three phenotypes). If we denote X 0t , X 1t ; …; X m t as the

1 It should be noted that, another type of symmetric division that CSC divides into two daughter NSCCs (termed as symmetric differentiation, Morrison and Kimble, 2006) is not accounted for in our model, however it is shown in Appendix F that the main results achieved in this study are still valid for the model including symmetric differentiation, implying the kinetic equivalence between the two models.

Please cite this article as: Zhou, D., et al., A multi-phenotypic cancer model with cell plasticity. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.04.039i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

D. Zhou et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

3

Table 1 Parameters used in the model. Symbol

Parameter

Cellular process

α00

Symmetric division rate by CSC

CSC⟶CSC þ CSC

α0j

Asymmetric division rate by CSC

CSC⟶CSC þ NSCCj ð1r jr mÞ

α0

Death rate by CSC

CSC⟶∅

αii

Symmetric division rate by NSCCi

NSCCi ⟶NSCCi þ NSCCi ð1 r ir mÞ

αi

α00 α0j α0

αii

αi0

Phenotypic conversion rate from NSCCi to CSC

NSCCi ⟶CSC ð1 r i r mÞ

αij

Phenotypic conversion rate from NSCCi to NSCCj

NSCCi ⟶NSCCj ði a jÞ

αi

Death rate by NSCCi

NSCCi ⟶∅ ð1 r i r mÞ

0

αij αi

Fig. 1. The model illustrated by the three-phenotype example (CSC, NSCC1 , NSCC2 ). For CSC (ﬁrst row, left to right): symmetric division, asymmetric division (differentiation to NSCC1 ), asymmetric division (differentiation to NSCC2 ), and cell death. For NSCC1 (second row, left to right): symmetric division, phenotypic conversion to CSC, phenotypic conversion to NSCC2 , and cell death. For NSCC2 (third row, left to right): symmetric division, phenotypic conversion to CSC, phenotypic conversion to NSCC1 , and cell death.

cell numbers of CSC, NSCC1 ; …; NSCCm at time t respectively, then ! T X t ¼ ðX 0t ; X 1t ; …; X m is an (m þ1)-dimensional cellular system t Þ consisting of ðm þ 1Þ ðm þ 2Þ cellular processes. On the basis of the modeling principles in biochemical reactions (see Chapter 7 in Van Kampen, 1992), our cellular system can be formulated by both stochastic and deterministic models as follows:

from ðx0 ; …; xi ; …; xm Þ to all the other possible states, formally, d Prfx0 ; …; xi ; …; xm ; tg ¼ Fðx0 ; …; xi ; …; xm ; tÞ; dt

ð1Þ

where Fðx0 ; …; xi ; …; xm ; tÞ m

2.2. Stochastic model The randomness of cellular systems is essentially rooted in the stochastic nature of gene expression in individual cells (Cai et al., 2006). The emergence of determinism from randomness is a result of the Law of Large Numbers. In other words, when the population is not ! large enough and subject to stochastic ﬂuctuations, X ¼ ðX 0t ; X 1t ; …; m T X t Þ should be stochastic. According to the theory of Chemical Master Equation (CME) (see Chapter 11 in Beard and Qian, 2008), the stochastic model can be formulated by a continuous-time Markov process on Nm þ 1 . If we deﬁne Prfx0 ; …; xi ; …; xm ; tg as the probability ! of X ¼ ðx0 ; …; xi ; …; xm ÞT at time t, then the rate of change in Prfx0 ; …; xi ; …; xm ; tg is equal to the rate of transition from all the other possible states to ðx0 ; …; xi ; …; xm Þ minus the rate of transition

¼ ∑ ðxi 1Þαii Prfx0 ; …; xi 1; …; xm ; t g i¼0 m

þ ∑ ðxi þ 1Þαi Prfx0 ; …; xi þ 1; …; xm ; t g i¼0 m

þ ∑ ðx0 α0j Pr x0 ; …; xj 1; …; xm ; t Þ j¼1

þ

∑

ði;jÞ A f1;2;…;mgf0;1;…;mg with i a j

∑

ðxi þ 1Þαij Pr x0 ; …; xi þ 1; …; xj 1…; xm ; t

ði;jÞ A f0;1;…;mgf0;1;…;mg

ðxi αij ÞPrfx0 ; …; xi ; …; xm ; tg

m

∑ ðxi αi ÞPrfx0 ; …; xi ; …; xm ; tg: i¼0

The ﬁrst four items of Fðx0 ; …; xi ; …; xm ; tÞ describe the transitions from other states to ðx0 ; …; xi ; …; xm Þ. For example, the ﬁrst item

Please cite this article as: Zhou, D., et al., A multi-phenotypic cancer model with cell plasticity. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.04.039i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

D. Zhou et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

4

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

corresponds to the sum of the rates of transition from ðx0 ; …; xi 1; …; xm Þ to ðx0 ; …; xi ; …; xm Þ by symmetric division of ith cell phenotype, 0 r i r m. Similarly, the other three items can be explained by other types of cellular processes accordingly. The last two items of Fðx0 ; …; xi ; …; xm ; tÞ describe the transitions from ðx0 ; …; xi ; …; xm Þ to other states. Since it is difﬁcult to calculate the analytic solution of Prfx0 ; …; xi ; …; xm ; tg, stochastic simulation (called Gillespie algorithm, Gillespie, 2007) can often be used as an efﬁcient numerical approach to investigate the CME. When the population size is large enough, it is more convenient to use deterministic equations to capture the model: 2.3. Deterministic model ! ! T Let 〈X t 〉 ¼ ð〈X 0t 〉; 〈X 1t 〉; …; 〈X m t 〉Þ be the expectation of X t , where 〈X it 〉≔

∑

x0 ;…;xi ;…;xm

xi Prfx0 ; …; xi ; …; xm ; tg

ð0 r i r mÞ:

We now derive the equation governing the dynamics of 〈X it 〉 (our method is based on Section 5.8 in Van Kampen, 1992). We multiply xi on the both sides of (1), and then calculate the summation over all the indexes d Prfx0 ; …; xi ; …; xm ; tg ¼ ∑ ∑ xi ðxi Fðx0 ; …; xi ; …; xm ; tÞÞ; dt x0 ;…;xi ;…;xm x0 ;…;xi ;…;xm that is, d〈X it 〉 ¼ ∑ ðxi Fðx0 ; …; xi ; …; xm ; tÞÞ: dt x0 ;…;xi ;…;xm

(1) when i ¼0 m d〈X 0t 〉 ¼ ðα00 α0 Þ〈X 0t 〉 þ ∑ ðαi0 〈X it 〉ÞÞ: dt i¼1

ð2Þ

(2) when 1 r ir m

! d〈X it 〉 j ¼ ∑ αji 〈X t 〉 þ αii αi ∑ αij 〈X it 〉: dt jai jai

ð3Þ

! Thus the dynamics of 〈X t 〉 can be formulated by the following (m þ1)-dimensional linear ordinary differential equations (ODEs): ! ! d〈X t 〉 ¼ A〈X t 〉; ð4Þ dt where

α00 α0 B α01 B B ⋮ A¼B B Bα 0ðm 1Þ @ α0m

α10

α11 α1 ∑j a 1 α1j

α20 α21

⋮

⋮

⋯

⋯

⋯

⋯

ð6Þ In comparison of A and An, it is interesting that there seems to be a trade-off between the diagonal and off-diagonal elements, that is, A has larger off-diagonal elements but smaller diagonal elements than An. Note that the diagonal elements correspond to the growth contributions from the cells with the same phenotypes themselves while the off-diagonal elements are the contributions from other phenotypes, cancer cell plasticity can be interpreted as a altruism behavior among cancer cells in the sense that it makes cancer cells of one phenotype help other phenotypes to survive by sacriﬁcing their own phenotype (Axelrod et al., 2006). How this altruism behavior affects the population structure of cancer cells is one of the major tasks in the study of cell plasticity models. In the following two sections, we will analyze both the asymptotic and transient dynamics of the model, then apply the theoretical results to explain the problems arising from concrete biological experiments.

3. Phenotypic equilibrium and frequency model

Then we have

0

to the one governing the hierarchical cancer stem cell model: 0 1 α00 α0 0 0 ⋯ 0 B α C α11 α1 0 ⋯ 0 01 B C B C n C: ⋮ ⋮ ⋮ ⋮ ⋮ A ¼B B C Bα C ⋯ ⋯ αðm 1Þðm 1Þ αðm 1Þ 0 @ 0ðm 1Þ A α0m ⋯ ⋯ 0 αmm αm

It has been reported in various cancer cell lines (Chaffer et al., 2011; Yang et al., 2012; Gupta et al., 2011; Iliopoulos et al., 2011) that the cancer cell populations starting from different initial proportions can return towards equilibrium proportions over time. The arise of this phenotypic equilibrium has performed as a profound indication in support of the intrinsic homeostasis in cancer. To explain this phenomenon, we investigate the dynamics of phenotypic proportions in this section. The model in previous section describes the population dynamics of the absolute cell numbers of different phenotypes. However in reality, relative numbers, i.e. the proportions of cell phenotypes, are usually measured by ﬂuorescence-activated cell sorting (FACS) experiments. Thus one often converts the population model to a frequency model. Let Nt be the total number of the population, N t ¼ X 0t þ X 1t þ i i i ⋯ þ Xm t , the proportion of X t is deﬁned as xt ¼ X t =N t . Then ! 0 1 m T xt ¼ ðxt ; xt ; …; xt Þ is the vector describing the phenotypic pro! ! portions. Let 〈xt 〉 be the expectation of xt , from Eq. (4) we can obtain the nonlinear (second order) ODEs governing the dynamics

⋯ ⋯

αm0 αm1

⋮

⋮

αðm 1Þðm 1Þ αðm 1Þ ∑j a ðm 1Þ αðm 1Þj αmðm 1Þ αðm 1Þm αmm αm ∑j a m αmj

1 C C C C: C C A

ð5Þ

! of 〈xt 〉 as follows (see mathematical details in Appendix A): It is noteworthy that Eq. (4) can also be obtained based on the Law of Mass Action (LMA) (Koudrjavcev et al., 2001), indicating the equivalent kinetic foundation between the stochastic theory of CME and its deterministic counterpart of LMA. If we choose to ignore cancer cell plasticity in the model, i.e. let αij ¼ 0 ð1 r i rm; 0 rj r m; and i ajÞ, the form of A will reduce

! d〈xt 〉 ! ! ! ¼ A〈xt 〉 〈xt 〉eT A〈xt 〉; dt

ð7Þ

where e ¼ ð1; 1; …; 1ÞT . Note that x0t þ x1t þ …þ xm t ¼ 1, if replacing 0 1 m1 xm Þ, the dimension of Eq. (7) should t by 1 ðxt þxt þ ⋯ þxt reduce from m þ 1 to m (see Appendix A).

Please cite this article as: Zhou, D., et al., A multi-phenotypic cancer model with cell plasticity. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.04.039i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

D. Zhou et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

To show the stability of the frequency model Eq. (7), we have Theorem 1: Theorem 1. If the Perron–Frobenius eigenvalue λ1 of the matrix A in Eq. (5) is simple (i.e. the algebraic multiplicity of λ1 is equal to one), then Eq. (7) has one and only one stable ﬁxed point and no stable limit cycle. (Technically, in very few cases it is necessary to add a small perturbation to the initial state for completing the ﬁnal proof, see proof in Appendix B). The Perron–Frobenius eigenvalue λ1 is the largest real eigenvalue satisfying that Re λ o λ1 for every other λ (Appendix B). On account of the inevitability of perturbations in real world, it can be proved that almost surely the algebraic multiplicity of λ1 is equal to one (Appendix C). Especially when validating the model with biological experiments, the condition of Theorem 1 can easily be satisﬁed. Theorem 1 thus theoretically explains the phenotypic equilibrium phenomena observed in cancer cell lines, that is, starting from different initial proportions, the cancer cell population can return towards stable equilibrium proportions as time passes. Furthermore, note that Theorem 1 holds for the general cases with any ﬁnite phenotypes, it can be predicted that the phenotypic equilibrium should be universal in multi-phenotypic population of cancer cells. In next section we will turn our attention from long-term asymptotic behavior to transient dynamics of the model.

4. Transient analysis and cell plasticity In this section, we investigate how the cell plasticity inﬂuences the transient proportions of cell phenotypes and maintains the heterogeneity of cancer. To explore the proportion of any phenotype i (0 r ir m), from Eq. (7) we have m m m d〈xit 〉 ¼ ∑ ðaij 〈xjt 〉Þ 〈xit 〉 ∑ ∑ ðanj 〈xjt 〉Þ: dt n¼0j¼0 j¼0

ð8Þ

We are interested in how 〈xit 〉 changes shortly after the drastic reduction of it. This problem is of particular interest to cancer biologists because of two major lines: one is that in cell culture experiments, they care about the transient changes of cancer cells shortly after the cell sorting. The other is that they are concerned about how the residual cancer cells relapse after the treatments. Assume that the proportion of phenotype i at time t0 becomes very small, i.e. 〈xit0 〉 0. Then Eq. (8) can be approximately simpliﬁed as follows: m d〈xit 〉 ∑ ðaij 〈xjt0 〉Þ: ð9Þ dt t ¼ t 0 jai Eq. (9) indicates that the transient growth of 〈xit 〉 at time t0 is almost contributed by the cell-state conversions from other phenotypes provided that current proportion of phenotype i is very small. That is, as soon as the cancer cells of one phenotype dramatically decreases, cancer cell plasticity helps increase their phenotypic proportion, which serves as a mutual-helping mechanism in transient dynamics. In particular, suppose the great majority of CSCs have already been eliminated (e.g. by CSCs-targeted drugs or cell sorting), i.e. 〈x0t 0 〉 0. It is easy to distinguish the conventional CSC model from the one with cell plasticity by comparing their transient dynamics: in the conventional CSC model without cancer cell plasticity, the growth of CSCs relies only on the self-renewal of CSCs themselves, the initial growth of CSCs proportion should be very limited (constrained by the rate limitation of cell division cycle). In contrast, the de-differentiations from other NSCCs can effectively speed up the transient grow rate of CSCs proportion. In this way, there should be a disparity of the

5

transient increase between these two models shortly after the initial time t0. To further illustrate the impact of cell plasticity on transient change of CSCs proportion, we will discuss two concrete examples in next two subsections: CSC–NSCC model and S–B–L model. 4.1. CSC–NSCC model The CSC–NSCC model is the most widely studied cell plasticity model. In this model, besides CSCs, all the other non-CSCs are grouped into one whole phenotype. We now treat the CSC–NSCC model as a two-phenotypic example of the multi-phenotypic framework, which contains six cellular processes as follows: (1) (2) (3) (4) (5) (6)

α00

CSC⟶CSC þCSC; α01 CSC⟶CSC þNSCC; α0 CSC⟶∅; α10 NSCC⟶CSC; α11 NSCC⟶NSCC þ NSCC; α1 NSCC⟶∅.

Let st be CSCs proportion, and nt ¼ 1 st be NSCCs proportion, the equation governing the change of CSCs proportion is then given by d〈st 〉 ¼ C〈st 〉2 þ D〈st 〉þ α10 ; dt

ð10Þ

where C ¼ ðα00 þ α01 α0 Þ ðα11 α1 Þ and D ¼ ðα00 α0 Þ ðα11 þ α10 α1 Þ. Note that α10 is the transition rate from NSCCs to CSCs, the model will reduce to the conventional CSC model without cell plasticity by letting α10 ¼ 0. Fig. 2 shows the comparison of the models with and without cell plasticity by validating them to the published data on SW620 colon cancer cell line (Yang et al., 2012). The method for parameter ﬁtting we performed is presented in Appendix D. It is found that the major difference between the two models lies in puriﬁed NSCCs case (Fig. 2A): starting from very small CSCs proportion (0.6% purity), the initial growth of CSCs proportion predicted by the model without cell plasticity is very slow, the growth rate will gradually increase until CSCs proportion tends to its equilibrium level, which shows a typical sigmoidal growth pattern (Rodriguez-Brenes et al.). However, the model with cell plasticity predicted a transient increase of CSCs proportion shortly after the initiation, which is in line with the experimental data. Even though both models predict the ﬁnal phenotypic equilibrium, only can the model with cell plasticity ﬁt the transient change on the data. This indicates that cancer cell plasticity could play more essential role in transient dynamics than in long-term behavior of CSCs proportion. A similar result is obtained if the stochastic version of the model is simulated using Gillespie algorithm (see Appendix E). 4.2. S–B–L model We now consider a three-phenotypic example of our model. Enlightened by Gupta et al.'s (2011) work, the three-phenotypic model comprises stem-like (S), basal (B) and luminal cells (L). Based on our framework, the S–B–L model contains 12 reactions as follows:

(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)

α00

S⟶S þ S; α01 S⟶S þ B; α02 S⟶S þ L; α0 S⟶∅; α11 B⟶B þ B; α10 B⟶S; α12 B⟶L; α1 B⟶∅; α22 L⟶L þ L; α20 L⟶S;

Please cite this article as: Zhou, D., et al., A multi-phenotypic cancer model with cell plasticity. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.04.039i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

D. Zhou et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

6

1

Proportion of CSCs

Proportion of CSCs

1 0.8 0.6 0.4 0.2 0

0

2

4

6

8

10

0.8 0.6 0.4 0.2 0

12

0

2

4

Time

8

10

12

Proportion of CSCs

1

0.8 0.6 0.4 0.2 0

6

Time

1

Proportion of CSCs

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 Q8 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

0

2

4

6

8

10

12

0.8 0.6 0.4

Experimental data Model with cell plasticity Model without cell plasticity

0.2 0

0

2

Time

4

6

8

10

12

Time

Fig. 2. The comparison of CSC–NSCC models (10) with and without cell plasticity by ﬁtting to experimental data on SW620 colon cancer cell line: the black dots are the experimental data in Yang et al. (2012), the blue dashed lines are the predictions of the model without cell plasticity (α10 ¼ 0), and the red dashed lines are the predictions of the model with cell plasticity (α10 40). There are four groups of data in the experiment: (A) Puriﬁed NSCCs case (0.6% CSCs þ 99.4% NCSSs); (B) 70% CSCs þ 30% NCSSs; (C) unsorted case (65.4% CSCs þ34.6% NCSSs); (D) puriﬁed CSCs case (99.4% CSCs þ 0.6% NCSSs). By least squares method, we predicted the values of four different cases together with the best ﬁtting parameters (see methods in Appendix D). It is shown that in (B), (C) and (D), the predicted trajectories by both models are in good accord with the data. In (A), however, compared to the sigmoidal growth pattern predicted by the model without cell plasticity, the prediction by the model with cell plasticity shows a transient increase in the ﬁrst few days, which is in line with the experimental data. (For interpretation of the references to color in this ﬁgure caption, the reader is referred to the web version of this paper.)

α21

(11) L⟶B. α2 (12) L⟶∅.

population genetics (Kussell and Leibler, 2005; Acar et al., 2008) regarding phenotypic variability as a surviving strategy for biological populations in ﬂuctuating environments.

Let ðst ; bt ; lt ÞT be the vector describing the proportions of (S, B, L), then we have 8 d〈st 〉 > > > ¼ ðα00 α0 Þ〈st 〉 þ α10 〈bt 〉 þ α20 〈lt 〉 st ½A1 〈st 〉 þ A2 〈bt 〉 þ A3 〈lt 〉 > > dt > > < d〈b 〉 t ¼ α01 〈st 〉 þ ðα11 α10 α12 α1 Þ〈bt 〉 þ α21 〈lt 〉 〈bt 〉½A1 〈st 〉 þ A2 〈bt 〉 þ A3 〈lt 〉 > dt > > > > d〈l 〉 > t ¼ α 〈s 〉 þ α 〈b 〉 þ ðα α α α Þ〈l 〉 〈l 〉½A 〈s 〉 þ A 〈b 〉 þ A 〈l 〉 > : t 02 t 12 t 22 20 21 2 t 1 t 2 t 3 t dt

5. Conclusions

ð11Þ where A1 ¼ α00 þ α01 þ α02 α0 , A2 ¼ α11 α1 , and A3 ¼ α22 α2 . Fig. 3 shows the predictions of S–B–L model by ﬁtting to cell-state dynamics on data of SUM159 breast cancer cell line (Gupta et al., 2011) (see methods in Appendix D). Starting from very small 〈s0 〉, in contrast to the gradual increase of 〈st 〉 predicted by the model with α20 ¼ 0, there is an overshooting of 〈st 〉 predicted by the model with α20 4 0, which in line with the cell-state dynamics on data. That is, the proportion of stem-like cells is rapidly elevated by de-differentiation to the value above the equilibrium level, and then returns towards the ﬁnal equilibrium as time passes. It has been reported that the CSC– NSCC model can never perform overshooting (Zhou et al., 2013; Jia et al.), which implies that the overshooting phenomenon is rooted in the diversity of phenotype in the multi-phenotypic model. Moreover, note that overshooting behavior is a result of biological response to environmental stimulus through evolution (Ma et al., 2009), cancer cell plasticity may have meaningful implications for the evolution of cancer. One possible explanation is that the reversible phenotype switching between cells provides an advantageous cooperative strategy for cancer cell populations during their competitions with various anticancer mechanisms (Axelrod et al., 2006), thus increasing their collective ﬁtness. This idea is in line with the literature in ecology and

Instead of the hierarchical structure of cancer cells proposed by cancer stem cell theory, the updated CSC paradigm with cell plasticity suggests a reversible nature between different cancer cell phenotypes. In this study we proposed a generalized multi-phenotypic cancer model for providing a mathematical framework to investigate the interplay between cell plasticity and diversity of phenotypes in cancer. Based on our model, some interesting and insightful results have been achieved. On one hand, it was shown (in Theorem 1) that the phenotypic equilibrium phenomena can be predicted by the stability of our frequency model. On the other hand, by comparing the models with and without cell plasticity, it was found that cell plasticity may play more essential role in transient dynamics. It has been shown that only can the models with cell plasticity predict the transient increase and overshooting phenomena in CSCs proportion. In particular, the overshooting phenomenon performed by the S–B–L model was shown to be one of the salient features in the multi-phenotypic models with cell plasticity. With the further study of the multi-phenotypic models with cell plasticity, more interesting and insightful phenomena would be found in future. Cancer cell plasticity has enriched the theory of cancer stem cell and comes along with new fundamental problems in cancer biology. Firstly, its molecular mechanism is poorly understood. It was reported that TGF-β might have important roles in the process of cell-state conversions through activating epithelial mesenchymal transition (EMT) (Yang et al., 2012). Further studies on related molecular models should be important tasks. Accordingly, besides

Please cite this article as: Zhou, D., et al., A multi-phenotypic cancer model with cell plasticity. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.04.039i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

D. Zhou et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Predictions by S−B−L model with α =0

Predictions by S−B−L model with α >0

20

20

1

1

0.9

0.9 0.8

0.8

S* B* L* S B L

Overshooting of CSCs proportion

Proportions of cells

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

7

0.1

0.7

0.7

0.09

S*

0.08

0.6

0.6

0.07

0.5

0.5

0.06 0.05

0.4

0.4

0.04 0.03

0.3

0.3

0.02 0

2

4

6

8

10

0.2

0.2

0.1

0.1 0

0 0

1

2

3

4

5

6

7

8

9

10

11

0

2

4

6

Time

8

10

Time

Fig. 3. Predictions of S–B–L model: we ﬁtted Eq. (11) to the cell-state dynamics on the data of SUM159 breast cancer cell line (Fig. 3 in Gupta et al., 2011). S, B and L are the predictions by the model, while Sn, Bn and Ln are the cell-state dynamics on data. There are three groups of data: (A) puriﬁed stem-like cells case; (B) puriﬁed basal cells case; (C) puriﬁed luminal cells case. We predicted the time points of three different cases together with the best ﬁtting parameters. This ﬁgure shows the predictions in puriﬁed luminal case, see Appendix D for the other two cases. The left panel is the predictions by the model without phenotypic conversion from luminal cells to stem-like cells (α20 ¼ 0), while the right panel is the predictions by the model with phenotypic conversion from luminal cells to stem-like cells (α20 4 0). It is shown that only can the model with cell plasticity predict the overshooting of the proportion of stem-like cells. The small window in left panel speciﬁcally focuses on the overshooting phenomenon of stem-like cells.

the population-level models studied in this work, molecular-level models should also be accounted for in future. Secondly, even though cell plasticity is found to be intimately related to the heterogeneity of cancer, cancer complexity is a result of multifactorial mechanisms (e. g. clonal evolution, tumor microenvironment, and reversible cell plasticity, Meacham and Morrison, 2013), and it is not clear in which cancer and to what extent the disease progression arises from the cell plasticity. That is, how to distinguish cell plasticity from other sources of cancer heterogeneity should be a central problem in future, for both theoretical and experimental researchers.

m

m

m

¼ ∑ ðaij 〈xjt 〉Þ 〈xit 〉 ∑ ∑ ðanj 〈xjt 〉Þ: n¼0j¼0

j¼0

Then we have Eq. (7). 0 1 m1 By letting xm Þ, t ¼ 1 ðxt þ xt þ ⋯ þ xt 0 r i rm 1)

m1

m

m1

n¼0

j¼0

i ¼ ∑ ðaij 〈xjt 〉Þ þ aim 〈xm t 〉 〈xt 〉 ∑ j¼0

m1

m

m1

j¼0

n¼0

j¼0

Acknowledgments

m

Appendix A. Derivation of Eq. (7) For Eq. (4) d〈X it 〉 ¼ ai0 〈X 0t 〉 þai1 〈X 1t 〉 þ ⋯ þ aim 〈X m t 〉: dt Note that xit ¼

X it 0 1 X t þ X t þ⋯ þ X m t

¼

X it ; Nt

d〈X it 〉 dð〈xit 〉〈N t 〉Þ d d〈xi 〉 ¼ ¼ 〈xit 〉 〈N t 〉 þ 〈N t 〉 t ; dt dt dt dt

have

m m m d〈xit 〉 ¼ ∑ ðaij 〈xjt 〉Þ 〈xit 〉 ∑ ∑ ðanj 〈xjt 〉Þ dt n¼0j¼0 j¼0

¼ ∑ ðaij aim Þ〈xjt 〉 〈xit 〉 ∑

We thank the referees for constructive comments. We also thank Prof. Minping Qian, Drs. Hao Ge, Chen Jia, Linyuan Liu, Hong Qian, Xin Tong, Benjamin Werner, Zhen Xie, Gen Yang, Michael Q Zhang, and Ruixiang Zhang, for helpful discussions. B. W. greatly acknowledges the generous sponsorship from Max-Planck Society. This work is supported by the Fundamental Research Funds for the Central Universities in China.

we

¼ 〈xit 〉 ∑

(for

!

∑ ðanj 〈xjt 〉Þ þ anm 〈xm t 〉 !

∑ ðanj anm Þ〈xjt 〉 þ anm þ aim

m1

m1

m

∑ ðanj anm Þ〈xjt 〉 〈xit 〉 ∑ anm þ ∑ ðaij aim Þ〈xjt 〉

þaim : |ﬄﬄﬄ{zﬄﬄ ﬄ} j¼0 |ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ}|ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ}constant term n¼0 j¼0

n¼0

nonlinear term

linear term

Appendix B. Proof of Theorem 1 Let us ﬁrst consider the properties of Eq. (4). We will show that there exists a largest real eigenvalue λ1 of A satisfying Re λ o λ1 for every other λ. Note that the off-diagonal elements of the matrix A are non-negative, A þ κ I becomes a non-negative matrix (I is the identity matrix) when κ is large enough. By Perron–Frobenius theory (see 0 Chapter 1 in Seneta, 1981), Aþ κ I has an real eigenvalue λ1 satisfying 0 0 0 Re λ o λ1 for any other eigenvalue λ of A þ κ I. It is easy to show that λ01 κ is just the largest real eigenvalue (called Perron–Frobenius eigenvalue λ1) satisfying that Re λ o λ1 for every other λ of A.2 Suppose λ1 is simple, the solution of Eq. (4) can be expressed as nj nj ! n ! ! 〈X t 〉 ¼ c1;1 u eλ1 t þ ∑ ∑ cj;l ∑ r jl;i t i 1 eλj t ; j¼2l¼1

i¼1

where λ1 ; λ2 ; …λn are the different eigenvalues of A, nj is the

then d〈xit 〉 1 d〈X it 〉 〈xit 〉 d〈Nt 〉 ¼ dt 〈N t 〉 dt 〈N t 〉 dt

0

0

2 If λ is an eigenvalue of A þ κ I, then detjA þ κ I λ Ij ¼ 0. For matrix A, we have 0 0 detjA ðλ κ ÞIj ¼ 0, namely λ κ is an eigenvalue of A. Then A has an real 0 0 0 0 eigenvalue λ1 κ satisfying Reðλ κ Þ o λ1 κ for any other eigenvalue λ κ of A.

Please cite this article as: Zhou, D., et al., A multi-phenotypic cancer model with cell plasticity. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.04.039i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

D. Zhou et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

8

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

! algebraic multiplicity of λj, u is the ! normalized (i.e. u0 þ u2 þ ⋯ þ um ¼ 1) right eigenvector of λ1, r jl;i is the corresponding eigenvector of λj, cj;l is a constant determined by initial values. Note that Re λi o λ1 ði a 1Þ, when c1;1 a 0, ! nj n cj;l nj ! 〈X t 〉 ! ! ¼ uþ ∑ ∑ ∑ r jl;i t i 1 eðλj λ1 Þt - u : λ t c1;1 e 1 j ¼ 2 l ¼ 1 c1;1 i ¼ 1 Thus ! ! ! 〈X t 〉 X t =c1;1 eλ1 t u ! ¼ 〈xt 〉 ¼ 〈Nt 〉 〈N t 〉=c1;1 eλ1 t u0 þ u2 þ ⋯ þum ! ! ! Since u has been normalized, 〈xt 〉 will tend to u when t- þ 1. In the following we will handle the case of c1;1 ¼ 0. Let t¼ 0, we get the linear equations of cj;l : nj ! n ! ! c1;1 u þ ∑ ∑ cj;l r jl;1 ¼ 〈 X 0 〉T j¼2l¼1

! ! ! ! Denote the coefﬁcient matrix by B ¼ ½ u r 21;1 r 22;1 ⋯r nnn ;1 , and let Bn ! be B with its ﬁrst column replaced by 〈 X 0 〉T . Notice that columns of B are linear independent. Then from Cramer's rule c1;1 ¼

detjBn j : detjBj

When c1;1 ¼ 0, namely detjBn j ¼ 0, we can ﬁnd an n-dimen! sional vector v which is linearly independent with the last n 1 ! ! columns of B. Actually, with a small perturbation ε v to 〈 X 0 〉T , all n the columns of B are linearly independent, and c1;1 a 0. Note that ﬂuctuations are inevitable in real world, c1;1 a 0 always holds in reality.

Appendix C. The algebraic multiplicity of Perron–Frobenius eigenvalue It was shown above that the Perron–Frobenius eigenvalue λ1 being simple is essential for proving Theorem 1. Here we will show that this condition is easily satisﬁed in real world. In fact, a more general result can be obtained as follows: Proposition 1. Denote Rnn as the space of n n real square matrices, equipped with Lebesgue measure. Then its subset of the matrices with repeated eigenvalues has measure zero. To prove this proposition, we need following lemma ﬁrst. Lemma 1. For any non-zero polynomial f ðx1 ; x2 ; …; xm Þ ðm Z1Þ, its zero set3 has measure 0 (w.r.t. Lebesgue measure on Rm ). Proof. We use inductive method. For m ¼ 1, the polynomial has only ﬁnite roots, whose zero set has measure 0. Assume that we have proved the lemma for m ¼ k 1. For m ¼k, let us take f ðx1 ; x2 ; …; xk Þ as the polynomial of ðx2 ; …; xk Þ, then the coefﬁcients can be regarded as polynomials of x1. It is easy to see that at most ﬁnite x1 A R can make all coefﬁcients equal to 0. Denote such set of x1 by S, then fðx1 ; x2 ; …; xk j x1 A Sg has measure 0 in Rk . For any ﬁxed x1 2 = S, f ðx1 ; x2 ; …; xk Þ is a non-zero polynomial of ðx2 ; ⋯; xk Þ. From inductive assumption we know its zero set has measure 0 in Rk 1 . By Fubini–Tonelli theorem, we know that fðx1 ; x2 ; …; xk Þj x1 2 = S; f ðx1 ; x2 ; ⋯; xk Þ ¼ 0g has measure 0 in Rk . Thus fðx1 ; x2 ; …; xk Þj f ðx1 ; x2 ; …; xk Þ ¼ 0g has measure 0 in Rk . □ Proof of Proposition 1. It is easy to show that for any square matrix in Rnn , the coefﬁcients of its characteristic polynomial are polynomials of matrix entries. Thus the discriminant of the 3 The zero set Rm jf ðx1 ; x2 ; ⋯; xm Þ ¼ 0g.

of

f ðx1 ; x2 ; ⋯; xm Þ

is

deﬁned

as

fðx1 ; x2 ; ⋯; xm Þ A

characteristic polynomial4 is again a polynomial of matrix entries. Applying Lemma 3 to the discriminant of characteristic polynomial, we have that the zero set of the discriminant has measure zero. And then according to the fact that a polynomial has repeated root if and only if its discriminant is 0 (Dickenstein and Emiris, 2005), the set of the matrices with repeated eigenvalues has measure zero. The proof is completed. □ Proposition 1 implies that almost surely the Perron–Frobenius eigenvalue λ1 should be simple. In particular, perturbations are inevitable in the real work, so λ1 can easily be simple. To further illustrate this result, we consider a 3 3 example: 0 1 3 2 3 B C A ¼ @ 0 3 5 A; ðC:1Þ 0 0 2 its eigenvalues are λ1 ¼ λ2 ¼ 3, λ3 ¼ 2. In this case the Perron– Frobenius eigenvalue λ1 is repeated. However, if adding perturbations to A, e.g. 0 1 3:23 2:19 3:22 B 3:78 5:11 C A¼@ 0 ðC:2Þ A; 0 0 2:31 its eigenvalues become λ1 ¼ 3:78, λ2 ¼ 3:23, λ3 ¼ 2:31. In this way the Perron–Frobenius eigenvalue λ1 become simple. Since the perturbations from experimental measurements and natural ﬂuctuations are inevitable in reality, the Perron–Frobenius eigenvalue λ1 in A is seldom to be repeated. For example, consider the S–B–L model in Section 4.2, let α00 ¼ 0:34, α01 ¼ 0:29, α02 ¼ 0:34, α0 ¼ 0:11, α11 ¼ 0:43, α10 ¼ 0:07, α12 ¼ 0:13, α1 ¼ 0:09, α22 ¼ 0:36, α20 ¼ 0:06, α21 ¼ 0:17, α2 ¼ 0:10, then the matrix A turns to 0 1 0:23 0:09 0:10 B C A ¼ @ 0:29 0:14 0:06 A: ðC:3Þ 0:34 0:07 0:03 The eigenvalues of A are λ1 ¼ 0:4391, λ2 ¼ 0:0844, λ3 ¼ 0:0453. The Perron–Frobenius eigenvalue λ1 is simple. Appendix D. Methods for parameter ﬁtting We used least squares method to estimate the parameters. Let ! !! y ¼ ðy0 ; …; yn ÞT be the data on n time points, x ð θ Þ ¼ ðx0 ; …; xn ÞT ! be the discrete trajectory of the model with parameter vector θ . The best-ﬁtting parameters are obtained by solving the following least squares problem: n

min ∑ ðxi yi Þ2 : ! i¼0

θ

Speciﬁcally, we ﬁrst discuss the parameter ﬁtting of the CSC–NSCC model to the data on SW620 cell line (Yang et al., 2012). There are four groups of data in their experiment, each group contains 13 time points (from day 0 to day 24, the data points were recorded every two days). That is, there are a total of 13 4 ¼ 52 time points. We estimated the parameters in Eq. (10) by ﬁtting to these 52 data points of four groups all together. It should be noted that, based on the population-level data measured by FACS experiment, only can we estimate the values of C, D and

α10 in Eq. (10), we cannot estimate all the αij or αi

4 The discriminant of m-order polynomial am xm þ am 1 xm 1 ⋯þ a1 x þ a0 is a polynomial function (see Dickenstein and Emiris, 2005, pp. 16–19 for concrete form) of its coefﬁcients am ; am 1 ; ⋯a1 ; a0 . A polynomial has repeated root if and only if its discriminant is 0 (Dickenstein and Emiris, 2005). For example, the 2 discriminant of quadratic polynomial ax2 þ bx þ c is b 4ac. ax2 þ bx þ c has 2 repeated root if and only if b 4ac ¼ 0.

Please cite this article as: Zhou, D., et al., A multi-phenotypic cancer model with cell plasticity. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.04.039i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

D. Zhou et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

individually. However, it is interesting that α10 is presented as an independent coefﬁcient, so we can still compare the models with and without cell plasticity based on the population-level analysis. We list the main points for parameter ﬁtting as follows: 1. Fit four groups data together: Let yji be the ith time point in jth groups, then we need to solve the following least squares problem: 4

min

12

∑ ∑ ðxji yji Þ2 :

ðC;D;α10 Þ j ¼ 1 i ¼ 0

9

case and puriﬁed luminal cells case). In each case, 12 time points were recorded for each cell phenotype. That is, a total of 3 12 3 ¼ 108 data points were recorded. Note that st þ bt þ lt ¼ 1, the effective number of the data is actually reduced to 72. By letting lt ¼ 1 st bt , we have 8 d〈st 〉 > > > > dt ¼ A4 〈st 〉 þ α10 〈bt 〉 þ α20 〈lt 〉 st ½A1 〈st 〉 þA2 〈bt 〉 þ A3 〈lt 〉 > < d〈bt 〉 ¼ α01 〈st 〉þ A5 〈bt 〉þ α21 〈lt 〉 〈bt 〉½A1 〈st 〉 þ A2 〈bt 〉 þ A3 〈lt 〉 > > > dt > > : 〈lt 〉 ¼ 1 st bt ðD:2Þ

2. Fixed initial states: We set that the initial states of Eq. (10) are the same as those on data, that is, let xj0 ¼ yj0 ð1 rj r 4Þ. 3. Parameter range: Based on the rate limitation of human cell proliferation cycle (Cowan et al., 2004), we assume that each cell requires at least one day to ﬁnish a cell cycle. To rescale the discrete time point to continuous time and note that the data was recorded every two days, we have 0 r αij r 2 log 2. Meanwhile, for the non-extinction of the cell population, we assume that cell death rate is smaller than symmetric division rate, i.e. 0 r αi o αii r 2 log 2. We performed fmincon algorithm in Matlab to calculate optimal values of ðC n ; Dn ; αn10 Þ. To test the uniqueness of the solutions, 100 different starting values of ðC 0 ; D0 ; α010 Þ were uniformly selected on the basis of the above parameter range. It was shown that the values of ðC; D; α10 Þ converge to C n ð10 6 ; 10 5 Þ, Dn ¼ 0:4002 7 0:0001, and αn10 ¼ 0:2630 7 0:0002. For the goodness of ﬁt, the coefﬁcient of determination is used to measure how well the observations on data are predicted by the models, that is, SSres 0:0588 ¼ 0:9407; ¼ 1 R ¼ 1 0:9921 SStot

where A4 ¼ α00 α0 , A5 ¼ α11 α10 α12 α1 . Similar to the CSC– NSCC model, on the basis of population level data, we cannot estimate all the 12 parameters αij and αi in the S–B–L model individually, since the freedom of the parameters in Eq. (D.2) is only nine. Note that the parameters for de-differentiation (α20 and α10) are presented as independent coefﬁcients, it is feasible to investigate how dedifferentiation inﬂuences the proportion of stem-like cells in transient dynamics. By applying similar numerical scheme we used in the CSC– NSCC model,6 we ﬁtted Eq. (D.2) to the cell-state dynamics on data. Fig. D1 shows the predictions by the model with α20 ¼ 0, and Fig. D2 shows the predictions by the model with α20 4 0. Note that the cellstate dynamic data in Gupta et al. (2011) was produced by the Markov chain, it is not surprising that the coefﬁcients of determination of the two models are both close to 1. But it is note worthy that the main difference between the two cases lies in the puriﬁed luminal case, where only can the model with α20 4 0 predict the overshooting of stem-like cells proportion. Appendix E. Stochastic simulations

2

where SSres is the sum of squared residuals and SStot is the total sum of squares (proportional to the sample variance) of data. By letting α10 ¼ 0, with the same numerical scheme we have C n ¼ 1:3865 70:0002, Dn ¼ 0:9550 7 0:0001, and the coefﬁcient of determination is R2 ¼ 1

α00 α0 B α01 þ 2β01 B B ⋮ Q ¼B B B α0ðm 1Þ þ 2β 0ðm 1Þ @ α0m þ 2β0m

Appendix F. A remark for the model By adding symmetric differentiation of CSC into our model,

SSres 0:5932 ¼ 0:4021: ¼ 1 0:9921 SStot

β0j

Fig. 2A shows that the disparity between the predictions by the two models mainly lies in the transient dynamics in puriﬁed NSCCs case, where only can the model with α10 4 0 predict the transient increase of CSCs proportion. Now we discuss the parameter ﬁtting of the S–B–L model. The time-series data we used is the cell-state dynamics on SUM159 breast cancer cell line (Fig. 3 in Gupta et al., 2011,5).There were three cases in their experiments (puriﬁed stem-like cells case, puriﬁed basal cells 0

The stochastic simulation of the CSC–NSCC model is illustrated in Fig. E1. The simulation scheme we used is based on Monte Carlo simulation for the CME (see Section 11.4.4 in Beard and Qian, 2008).

α10

α11 α1 ∑j a 1 α1j

α20 α21

⋮

⋮

⋯

⋯

⋯

⋯

CSC⟶NSCCj þ NSCCj ð1 r j rmÞ ! the dynamics of 〈X t 〉 can still be formed by the linear ODEs as follows: ! ! d〈X t 〉 ¼ Q 〈X t 〉; ðF:1Þ dt where

⋯

αm0 αm1

⋮

⋮

⋯

αðm 1Þðm 1Þ αðm 1Þ ∑j a ðm 1Þ αðm 1Þj αmðm 1Þ αðm 1Þm αmm αm ∑j a m αmj

5 It should be noted that this data was produced by discrete-time Markov chain with transition matrix: 0 1 0:58 0:35 0:07 B 0 C P ¼ @ 0:01 0:99 ðD:1Þ A: 0:04 0:49 0:47

1 C C C C: C C A

ðF:2Þ

(footnote continued) According to the methods in Gupta et al. (2011) this cell-state dynamic data was in good line with experimental results on SUM159 breast cancer cell line. 6 Since the time point in the S–B–L model was recorded every day, the parameter constraints are changed to 0r αij r log 2.

Please cite this article as: Zhou, D., et al., A multi-phenotypic cancer model with cell plasticity. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.04.039i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

D. Zhou et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1

1

1

0.9

0.9

0.9

0.8

0.8

0.8

0.7

0.7

0.7

0.6

0.6

0.6

0.5

0.5

0.5

0.4

0.4

0.4

0.3

0.3

0.3

0.2

0.2

0.2

0.1

0.1

0.1

0

0

1

2

3

4

5

6

7

8

9

10 11

S* B* L* S B L

0

0 0

1

2

3

4

5

6

7

8

9

0

10 11

1

2

3

4

Time

Time

5

6

7

8

9

10 11

Time

Fig. D1. Fit Eq. (D.2) to cell-state dynamics in Gupta et al. (2011) by letting α20 ¼ 0. The best ﬁtting parameters are A4 ¼ 0:032 7 0:002, α10 ¼ 0:016 7 0:003, α20 ¼ 0, A1 ¼ 0:578 7 0:002, A2 ¼ 0:581 7 0:002, A3 ¼ 0:550 7 0:001, α01 ¼ 0:416 7 0:001, A5 ¼ 0:567 7 0:002 and α21 ¼ 0:693 7 0:001. The sum of squared residuals SSres is about 0.06 and the coefﬁcient of determination R2 ¼ 0:9851. The unﬁtting part lies in the transient dynamics of the puriﬁed luminal cells case. In particular, the overshooting of the stem-like cells proportion cannot be predicted by the model. (a) Puriﬁed stem-like cells case, (b) Puriﬁed basal cells case and (c) Puriﬁed luminal cells case.

Proportions of cells

1

1

1

0.9

0.9

0.9

0.8

0.8

0.8

0.7

0.7

0.7

0.6

0.6

0.6

0.5

0.5

0.5

0.4

0.4

0.4

0.3

0.3

0.3

0.2

0.2

0.2

0.1

0.1

0.1

0

0 0

1

2

3

4

5

6

7

8

9

10 11

S* B* L* S B L

0 0

1

2

3

Time

4

5

6

7

8

9

10 11

0

1

2

3

4

Time

5

6

7

8

9

10 11

Time

Fig. D2. Fit Eq. (D.2) to cell-state dynamics in Gupta et al. (2011) by letting α20 40. The best ﬁtting parameters are A4 ¼ 0:0277 0:002, α10 ¼ 0:0137 0:001, α20 ¼ 0:0717 0:001, A1 ¼ 0:580 7 0:002, A2 ¼ 0:578 7 0:002, A3 ¼ 0:5797 0:001, α01 ¼ 0:4187 0:001, A5 ¼ 0:567 7 0:003, and α21 ¼ 0:688 7 0:001. The sum of squared residuals is about ð10 7 ; 10 8 Þ, so the coefﬁcient of determination R2 is very close to 1. (a) Puriﬁed stem-like cells case, (b) Puriﬁed basal cells case and (c) Puriﬁed luminal cells case.

It is easy to ﬁnd that the only difference between Q and A in Eq. (5) lies in the ﬁrst column, but this does not change the nonnegativity of the off-diagonal elements. Note that the main ingredient of the proof of Theorem 1 is the non-negativity of the off-diagonal elements in A (see Appendix B), Theorem 1 still holds for Q. Moreover, the transient properties discussed in our model are still valid for the model with Q here, by only replacing aij in Eq. (8) with qij here. Let us take the CSC–NSCC model as an example. The equation governing the change rate of CSCs proportion is given by

Stochastic simulations

0.45 0.4 0.35

Proportion of CSCs

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

Proportion of cells

10

0.3 with cell plasticity without cell plasticity

0.25 0.2 0.15

d〈st 〉 ¼ C n 〈st 〉2 þDn 〈st 〉 þ α10 ; dt

0.1 0.05 0 0

5

10

15

20

25

Time points

Fig. E1. Stochastic simulation of CSC–NSCC model: it is shown that starting from very small CSCs proportion, the stochastic model without cell plasticity shows sigmoidal pattern growth, while the model with cell plasticity predicts the transient increase of CSCs proportion. This is in accord with the prediction by Eq. (10) in Fig. 2A. The small window shows the strategy of our stochastic simulation. For given parameters, 20 stochastic samples (colorful thin lines) were produced by Monte carlo simulation. The thick black line is the representative trajectory by averaging these 20 stochastic samples. (For interpretation of the references to color in this ﬁgure caption, the reader is referred to the web version of this paper.)

ðF:3Þ

where C n ¼ ðα00 þ α01 þ β 01 α0 Þ ðα11 α1 Þ and Dn ¼ ðα00 α0 Þ ðα11 þ α10 α1 Þ. The property of transient increase is still dependent on α10, which shares the same nature as Eq. (10). References Acar, M., Mettetal, J.T., van Oudenaarden, A., 2008. Stochastic switching as a survival strategy in ﬂuctuating environments. Nat. Genet. 40 (4), 471–475. Antal, T., Krapivsky, P.L., 2011. Exact solution of a two-type branching process: models of tumor progression. J. Stat. Mech. 2011, P08018. Axelrod, R., Axelrod, D., Pienta, K., 2006. Evolution of cooperation among tumor cells. Proc. Natl. Acad. Sci. USA 103 (36), 13474–13479.

Please cite this article as: Zhou, D., et al., A multi-phenotypic cancer model with cell plasticity. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.04.039i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

D. Zhou et al. / Journal of Theoretical Biology ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1 2 Q7 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 Q6 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

Beard, D.A., Qian, H., 2008. Chemical Biophysics: Quantitative Analysis of Cellular Systems. Cambridge University Press. Boman, B., Wicha, M., Fields, J., Runquist, O., 2007. Symmetric division of cancer stem cells—a key mechanism in tumor growth that should be targeted in future therapeutic approaches. Clin. Pharmacol. Ther. 81 (6), 893–898. Cai, L., Friedman, N., Xie, X., 2006. Stochastic protein expression in individual cells at the single molecule level. Nature 440 (7082), 358–362. Chaffer, C., Brueckmann, I., Scheel, C., Kaestli, A., Wiggins, P., Rodrigues, L., Brooks, M., Reinhardt, F., Su, Y., Polyak, K., et al., 2011. Normal and neoplastic nonstem cells can spontaneously convert to a stem-like state. Proc. Natl. Acad. Sci. USA 108 (19), 7950–7955. Colijn, C., Mackey, M.C., 2005. A mathematical model of hematopoiesisłi. periodic chronic myelogenous leukemia. J. Theor. Biol. 237 (2), 117–132. Cowan, C., Klimanskaya, I., McMahon, J., Atienza, J., Witmyer, J., Zucker, J., Wang, S., Morton, C., McMahon, A., Powers, D., et al., 2004. Derivation of embryonic stem-cell lines from human blastocysts. N. Engl. J. Med. 350 (13), 1353–1356. Dalerba, P., Cho, R., Clarke, M., 2007. Cancer stem cells: models and concepts. Annu. Rev. Med. 58, 267–284. Dickenstein, A., Emiris, I.Z., 2005. Solving Polynomial Equations: Foundations Algorithms and Applications, vol. 14. Springer. Dingli, D., Traulsen, A., Michor, F., 2007. (a) symmetric stem cell replication and cancer. PLoS Comput. Biol. 3 (3), e53. dos Santos, R.V., da Silva, L.M., 2013a. A possible explanation for the variable frequencies of cancer stem cells in tumors. PloS One 8 (8), e69131. dos Santos, R.V., da Silva, L.M., 2013b. The noise and the kiss in the cancer stem cells niche. J. Theor. Biol. 335 (21), 79–87. French, R., Clarkson, R., 2012. The complex nature of breast cancer stem-like cells: heterogeneity and plasticity. J. Stem Cell Res. Ther. 7, 009 (Special). Gillespie, D.T., 2007. Stochastic simulation of chemical kinetics. Annu. Rev. Phys. Chem. 58, 35–55. Gupta, P., Fillmore, C., Jiang, G., Shapira, S., Tao, K., Kuperwasser, C., Lander, E., 2011. Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells. Cell 146 (4), 633–644. Iliopoulos, D., Hirsch, H., Wang, G., Struhl, K., 2011. Inducible formation of breast cancer stem cells and their dynamic equilibrium with non-stem cancer cells via il6 secretion. Proc. Natl. Acad. Sci. USA 108 (4), 1397–1402. Jia, C., Qian, M., Jiang, D., Overshoot in Biological Systems Modeled by Markov Chains: a Nonequilibrium Dynamic Phenomenon, arXiv preprint arXiv:1311.2260. Johnston, M.D., Edwards, C.M., Bodmer, W.F., Maini, P.K., Chapman, S.J., 2007. Mathematical modeling of cell population dynamics in the colonic crypt and in colorectal cancer. Proc. Natl. Acad. Sci. USA 104, 4008–4013. Johnston, M.D., Maini, P.K., Jonathan Chapman, S., Edwards, C.M., Bodmer, W.F., 2010. On the proportion of cancer stem cells in a tumour. J. Theor. Biol. 266 (4), 708–711. Jordan, C., Guzman, M., Noble, M., 2006. Cancer stem cells. N. Engl. J. Med. 355 (12), 1253–1261. Kaneko, K., 2012. Evolution of robustness and plasticity under environmental ﬂuctuation: formulation in terms of phenotypic variances. J. Stat. Phys. 148 (4), 686–704. Koudrjavcev, A.B., Jameson, R.F., Linert, W., 2001. The Law of Mass Action. Springer. Kussell, E., Leibler, S., 2005. Phenotypic diversity, population growth, and information in ﬂuctuating environments. Science 309 (5743), 2075–2078. Leder, K., Pitter, K., LaPlant, Q., Hambardzumyan, D., Ross, B.D., Chan, T.A., Holland, E.C., Michor, F., 2014. Mathematical modeling of pdgf-driven glioblastoma reveals optimized radiation dosing schedules. Cell 156 (3), 603–616. Lu, T., Shen, T., Bennett, M.R., Wolynes, P.G., Hasty, J., 2007. Phenotypic variability of growing cellular populations. Proc. Natl. Acad. Sci. USA 104 (48), 18982–18987.

11

47 Ma, W., Trusina, A., El-Samad, H., Lim, W., Tang, C., 2009. Deﬁning network topologies that can achieve biochemical adaptation. Cell 138 (4), 760–773. 48 Meacham, C.E., Morrison, S.J., 2013. Tumour heterogeneity and cancer cell plasticity. 49 Nature 501 (7467), 328–337. 50 Meyer, M., Fleming, J., Ali, M., Pesesky, M., Ginsburg, E., Vonderhaar, B., 2009. Dynamic regulation of cd24 and the invasive, cd44poscd24neg phenotype in 51 breast cancer cell lines. Breast Cancer Res. 11 (6), R82. 52 Molina-Peña, R., Alvarez, M., 2012. A simple mathematical model based on the 53 cancer stem cell hypothesis suggests kinetic commonalities in solid tumor growth. PloS One 7 (2), e26233. 54 Morrison, S.J., Kimble, J., 2006. Asymmetric and symmetric stem-cell divisions in 55 development and cancer. Nature 441 (7097), 1068–1074. 56 Nguyen, L., Vanner, R., Dirks, P., Eaves, C., 2002. Cancer stem cells: an evolving concept. Nat. Rev. Cancer 12 (2), 133–143. 57 Pisco, A.O., Brock, A., Zhou, J., Moor, A., Mojtahedi, M., Jackson, D., Huang, S., 2013. 58 Non-Darwinian dynamics in therapy-induced cancer drug resistance. Nat. 59 Commun. 4, 2467. Quintana, E., Shackleton, M., Foster, H.R., Fullen, D.R., Sabel, M.S., Johnson, T.M., 60 Morrison, S.J., 2010. Phenotypic heterogeneity among tumorigenic melanoma 61 cells from patients that is reversible and not hierarchically organized. Cancer 62 Cell 18 (5), 510–523. 63 Reya, T., Morrison, S., Clarke, M., Weissman, I., 2001. Stem cells, cancer, and cancer stem cells. Nature 414, 105–111. 64 Rodriguez-Brenes, I.A., Komarova, N.L., Wodarz, D., Tumor growth dynamics: insights into evolutionary processes, Trends. Ecol. Evol. Q2 65 66 Scafﬁdi, P., Misteli, T., 2011. in vitro generation of human cells with cancer stem cell properties. Nat. Cell Biol. 13 (9), 1051–1061. 67 Seneta, E., 1981. Non-negative Matrices and Markov Chains, second edition 68 Springer, New York. 69 Sottoriva, A., Vermeulen, L., Tavaré, S., 2011. Modeling evolutionary dynamics of epigenetic mutations in hierarchically organized tumors. PLoS Comput. Biol. 7 70 (5), e1001132. 71 Stiehl, T., Marciniak-Czochra, A., 2012. Mathematical modeling of leukemogenesis 72 and cancer stem cell dynamics. Math. Model. Nat. Phenom. 7 (01), 66–202. Todaro, M., Francipane, M.G., Medema, J.P., Stassi, G., 2010. Colon cancer stem cells: 73 promise of targeted therapy. Gastroenterology 138 (6), 2151–2162. 74 Van Kampen, N.G., 1992. Stochastic Processes in Physics and Chemistry, vol. 1. 75 Elsevier Science, Amsterdam. Visvader, J., Lindeman, G., 2012. Cancer stem cells: current status and evolving 76 complexities. Cell Stem Cell 10 (6), 717–728. 77 Wang, W., Quan, Y., Fu, Q., Liu, Y., Liang, Y., Wu, J., Yang, G., Luo, C., Ouyang, Q., 78 Wang, Y., 2014. Dynamics between cancer cell subpopulations reveals a model 79 coordinating with both hierarchical and stochastic concepts. PloS One 9 (1), e84654. 80 Werner, B., Dingli, D., Lenaerts, T., Pacheco, J., Traulsen, A., 2011. Dynamics of 81 mutant cells in hierarchical organized tissues. PLoS Comput. Biol. 7 (12), 82 e1002290. Wolf, M.D., Vazirani, V.V., Arkin, A.P., 2005. Diversity in times of adversity: 83 probabilistic strategies in microbial survival games. J. Theor. Biol. 234 (2), 84 227–253. 85 Yang, G., Quan, Y., Wang, W., Fu, Q., Wu, J., Mei, T., Li, J., Tang, Y., Luo, C., Ouyang, Q., et al., 2012. Dynamic equilibrium between cancer stem cells and non-stem 86 cancer cells in human sw620 and mcf-7 cancer cell populations. Br. J. Cancer 87 106 (9), 1512–1519. 88 Zapperi, S., La Porta, C., 2012. Do cancer cells undergo phenotypic switching? The case for imperfect cancer stem cell markers. Sci. Rep. 2, 441. 89 Zhou, D., Wu, D., Li, Z., Qian, M., Zhang, M.Q., 2013. Population dynamics of cancer 90 cells with cell state conversions. Quant. Biol. 1 (3), 201–208.

Please cite this article as: Zhou, D., et al., A multi-phenotypic cancer model with cell plasticity. J. Theor. Biol. (2014), http://dx.doi.org/ 10.1016/j.jtbi.2014.04.039i

91