Clonal pattern dynamics in tumor: the concept of cancer stem cells

We present a multiphase model for solid tumor initiation and progression focusing on the properties of cancer stem cells (CSC). CSCs are a small and singular cell sub-population having outstanding capacities: high proliferation rate, self-renewal and extreme therapy resistance. Our model takes all these factors into account under a recent perspective: the possibility of phenotype switching of differentiated cancer cells (DC) to the stem cell state, mediated by chemical activators. This plasticity of cancerous cells complicates the complete eradication of CSCs and the tumor suppression. The model in itself requires a sophisticated treatment of population dynamics driven by chemical factors. We analytically demonstrate that the rather important number of parameters, inherent to any biological complexity, is reduced to three pivotal quantities.Three fixed points guide the dynamics, and two of them may lead to an optimistic issue, predicting either a control of the cancerous cell population or a complete eradication. The space environment, critical for the tumor outcome, is introduced via a density formalism. Disordered patterns are obtained inside a stable growing contour driven by the CSC. Somewhat surprisingly, despite the patterning instability, the contour maintains its circular shape but ceases to grow for a typical size independently of segregation patterns or obstacles located inside.

Discrete 1-3 and continuous 4 models have been proposed recently to study solid tumor progression within the cancer stem cells (CSCs) hypothesis. According to this concept 5 , CSCs are a small subpopulation within the tumor that are the main cause for its progression. These cells have high proliferation rates, self-renewal and quasi-immortality capacities 6,7 . Consequently, they are very difficult to detect and to eradicate. The cancer stem cell hypothesis has been matter of strong debates during the last forty years 8 . Nevertheless it explains several biological evidences such as the relapse of tumor post-surgery and treatments 9 or the ability to generate tumors in xenotransplantation 10 . Nowadays, surface biomarkers, in charge of identification and detection of these cells, are more and more evaluated for different classes of cancer which allows new therapeutics in clinical 11 but also experimental measurements of their physical properties 12 .
The CSC model 13 enters into a more general class, usually referred as hierarchical models where cancer cells originate from a precursor which divides either symmetrically giving two identical CSCs or two DCs (differentiated cancer cells) or asymmetrically (one CSC and one DC). In these models, only the complete extinction of the CSC population eradicates the tumor. The stochastic model, instead, considers the possibility that every cell can be a tumor-initiating one and so the best eradication strategy becomes a more complex issue. Our work lies between both concepts: indeed, as stated by Plaks et al. 14 , this distinction is often a false dichotomy because of the phenotype plasticity of cancerous cells 15 which allows to recover the stem traits by dedifferentiation. The same feature also occurs with healthy differentiated cells which are known to gain a stem state under traumatic conditions, but not only, and this plasticity process seems much more common than previously believed 8 . Within these ideas, we propose a rather simple dynamical tumor initiation model in the continuous limit, for a mixed population, composed of 3 subpopulations: cancer stem cells (CSCs), differentiated cancer cells (DCs) and all other cells that a tumor harbors such as quiescent, dead or immune cells. Each cell population, represented by its density, enters a nonlinear dynamical system based on proliferation rates controlled by chemical pathways and physical interactions. We demonstrate that the dynamics is governed by the existence and stability of 3 fixed points showing the possibility of full tumor extinction, but also of over-shooting dynamics after a period of remission or relapse of the tumor post surgery.These features, observed in practice, are not commonly recovered in cancer growth modeling 16,17 when biological parameters are fixed. 1 Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Str. 38, D-01187, Dresden, Germany. 2 Laboratoire de Physique de l'Ecole normale supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, F-75005, Paris, France. 3 Institut Universitaire de Cancérologie, Faculté de médecine, Sorbonne Université, 91 Bd de l'Hôpital, 75013, Paris, France. *email: benamar@lps.ens.fr In the plastic region, d < q 0 , the value of the parameters are intimately related, m 2 being always close to m 0 indicating an easy phenotypic switching and so plasticity. Close to  2 , tumor cells are then in the highest plastic state. In the following, Table 2 classifies all the fixed points in each region of the phase diagram. At this stage, the spatial variation is ignored and the time unit t 0 has been chosen as the inverse of the mitotic rate of the CSCs. CSCs are considered as immortal, as opposed to DCs and to the chemicals a or m, so no death rate is added for S. According to the second equation of Eq. (1), the DC proliferation and death are included in the single parameter d. Excessive growth of D leads to the disappearance of CSCs, the so-called Allee effect, demonstrated by Konstorum et al. 4 and to an optimistic issue for tumor elimination. Unfortunately, as CSCs tend to disappear, the de-differentiation (or plasticity) occurs from DC to CSC, controlled by m as shown in Fig. (1). Then, the Allee effect may be compromised by the re-population of CSCs, and ultimately of DCs, which happens for m greater than a typical value m 0 . To ensure a sharp variation of q(m), the parameter σ which enters the last equation of Eq. (2) is a tiny quantity and a real feedback DCs versus CSCs is possible if γ/α > m 0 . In practice, in vivo or in vitro, the population of CSCs is always a tiny fraction of the entire tumor population (few per cents). When S is below this tiny value, the activator m starts to grow as q(m), leading to an increase of S. Notice that the time-scale for de-differentiation q 0 , Eq. (2), has to be compared to the symmetric mitotic rate but the issue of this dynamical system depends on the comparison of q 0 with d (the time death rate of DCs) and with α (time elimination rate of both chemicals a and m which are chosen equal for simplification). As in 4 , the dynamics of a results from the positive competition between the activator itself, the CSC's population S and the spontaneous degradation rate α. To conclude, in our model, the selected biological ingredients of the CSCs hypothesis leads to a dynamical system of population of fourth order with a minimal set of 8 parameters, two of them having small values as S 0 and σ. As demonstrated in the next section and suggested by the phase diagram in Fig.(1b), only 3 parameters really control the full dynamics and the tumor progression.

Dynamics: Fixed Points, Stability, Basin of Attraction
Existence of fixed points. The strategy of population dynamics consists in the search of fixed points, achieved by cancelling all time derivatives in Eq. (1). Followed by stability analysis, this part does not present any difficulty and details can be found in the Appendix (A). Here we only focus on the main results. Fixed points represent static and homogeneous solutions of the physical system, which is assumed to be spatially infinite. Defined by the 4 concentrations: i  = {S i , D i , a i , m i }, they share in common 2 relationships:  (1). In other words, this second fixed point only exists under two conditions: first, there is no self-proliferation of the DCs and second, cancerous cell plasticity occurs. The first condition then implies that the D cells are in an homeostatic state. In practice, both S and D must remain smaller than 1 just like their sum, but more importantly S and D must remain positive leading to a new inequality m 2 < γ/α. Then, the existence of the second fixed point imposes the following inequalities of parameters: 0 < d < q 0 and m 2 < γ/α.
Considering now a i ≠ 0, the third equation of Eq. (1) gives S 3 = α(1 + a 3 )/(βa 3 ), so γ = − m ae / S S 3 / 3 0 and finally D 3 = S 3 /d, see Eqs (16,25). However, the first equation must be verified leading to (2p(D 3 ,a 3 ) − 1)d + q = 0 so p ~ 1/2, for small S 0 and σ. An algebraic equation of second order for a 3 gives a positive solution if ψ < d(β/α) (see Appendix C). Then the existence of this third fixed point requires an efficient but not too strong retro-action originating from the DCs. Table 2. Fixed points of the four regions of the phase diagram of Fig. (1b). Every fixed point is classified with respect to the eigenvalues of the corresponding Jacobian matrix. Every eigenvalue is fully determined but for 3  , the pair of complex conjugate eigenvalues may have either a positive or negative real part (in the table). This will not change the global stability, since an eigenvalue with positive real part (λ 2 in the Table) always exists as demonstrated in the Appendix C. In the blue region (d > q 0 , d < ψα/β), only one fixed point exists with real and negative eigenvalues. Every initial concentration of cells and chemicals will then decay exponentially to 1  , which indicates a complete tumor regression. In the orange region (d < q 0 , d < ψα/β), the Jacobian near 1  has two positive and real eigenvalues, leading to a long time exponential relaxation for arbitrary initial concentration to 2  . In the green region (d > q 0 , d > ψα/β), every initial concentration in the basin of attraction of  1 , will exponentially relax to full regression, whilst increasing if it lays outside the basin. A complete discussion about this scenario has been given in the main text (Section (Dynamics)). In the red region of Fig. 1, (d < q 0 , d > ψα/β), the scenario is very rich, due to the existence of all fixed points. The Jacobian of 1  has two real and positive eigenvalues, making this point unstable. The tumor will never undergo regression. Instead,  2 is globally stable; an initial concentration of chemicals and cells in the vicinity of 2  will relax to a low population tumor. Differently from the orange region, the existence of  3 restricts the basin of attraction of 2  , making the tumor less likely to relax to  2 .
A summary of these conclusions can be found in the phase-diagram in Fig. (1b) which shows that 2 fixed points always exist,  1 and 3  , contrary to 2  which only exists if 0 < d < q 0 . Independently of the parameter restrictions, these fixed points are potentially relevant for tumor growth if they are stable with a finite basin of attraction. This is the reason why we consider now the dynamics of small perturbations, described by the linearization of the right-hand-side of Eq. (1).

Stability of the fixed points.
At linear order, rewritten in matrix form, it reads:  www.nature.com/scientificreports www.nature.com/scientificreports/ parameters inside the basin of attraction of 2  , the streamlines and the dynamics of both S(t) and D(t) are presented. The concentrations exhibit an over-shooting behavior before reaching the equilibrium values (see Fig. (2b)) and that the streamlines are disordered in Fig. (2c). (c) The dynamics of 3  is more complex and requires appropriate approximations. Assuming ψ < dβ/α, it is at least a saddle point, possibly an unstable point, see Fig. (2d) where it is a saddle fixed point.
To conclude, among all parameters of the model, only d and q 0 play a pivotal role. It exists always a stable fixed point with low cancerous cell concentration: either  1 or 2  . When  2 appears, for d < q 0 , 1  becomes unstable. When 3  exists, it is never stable. For a clinical perspective, to reach either  1 or 2  takes a great importance. However when the feedback and plasticity of DCs is too weak compared to their death rate, the emergence of 3  may compromise these optimistic solutions. One can wonder if the scenario depicted here is maintained when the feedback parameter q 0 becomes a time-dependent variable, q 0 = q 0 (t) or if some chemical values as a or m become stochastic. This may occur for multiple reasons, the main coming from the genetic of cancerous cells, mutation after mutation. Another reason may be environmental factors, changing in the chemicals activity as drug for example. In the Appendix, D, we perform numerical simulations of the dynamical model, Eq. (1), with a time-dependent q 0 (t) and we give preliminary results on stochastic effects.
Strategy for tumor regression. The previous analysis shows the dramatic role of two quantities: d and q(m) representing respectively the DCs proliferation rate (d < 0) or the DCs death rate (d > 0) and q(m) the DCs dedifferentiation rate. A parameter d < 0 (by absorption of nutrients, for example) induces an exponential growth in both populations of CSCs and DCs. Even for d > 0, (by senescence or death via drug treatments), it may not be www.nature.com/scientificreports www.nature.com/scientificreports/ sufficient even if the death rate d is stronger than the plasticity represented by q(m). In other words, if d > q(m) holds, it will not imply automatically a regression in all cases due to the vicinity of the third repulsive fixed point  3 .
The ideal case would be to remain in the blue region of the diagram, (Fig.(1b)), where there exists a unique stable fixed point 1  , insisting that, in this case, q(m 1 ) < d < ψα/β. The most obvious strategy will consist in reducing q(m) as much as possible so that the first inequality becomes easier to reach. Since m 1 = γ/α (see Eq. (3)) if m 1 < m 0 , then q(m 1 ) ~ 0, the first inequality is checked and there is hope to reach complete extinction of the tumor by eradicating both CSCs and DCs. Notice that, it requires also the second inequality to be verified. In this case, one should control the degradation rate α of the activators which plays a crucial role. Increasing α will be the best way to lead the tumor to extinction. Moreover, controlling the aggressiveness of the self-renewal activator via the parameter β as well as the inhibitory feedback of the DCs via ψ (Eq. (2) makes this situation more likely to occur (see Fig. (6a,b) in the Appendix D.1). The inset of Fig. (2c) representing the (S, a) subspace at fixed value d = 0 shows the small size of the basin of attraction, meaning that even at low levels of S and a, the tumor cannot undergo full tumor regression and differentiated cells may be produced again. Considering now m 0s and m 1 fixed, varying σ will not affect too much q(m 1 ) since σ only controls the smoothness of the function q(m).
Suppose now that m 1 > m 0 , then q(m 1 ) ~ q 0 . The worst scenario occurs when d < q 0 because the stable attractive fixed point will be the second one  2 . What we would really like to do is to decrease as much as possible the value of S 2 and D 2 , concentrations of CSCs and DCs respectively. But, S 2 is of order S 0 , D 2 is of order S 0 /d. It will be difficult to decrease these cancerous populations since they are dependent on biological parameters that we cannot control. Estimation of the time scale for proliferation of CSCs depends on the cancer origin. In 27 , the rate of cell division is about 1.8 month −1 , and estimation of the de-differentiation factor is of order 0.6 per month. In pancreatic xenograft tumors, the rate of CSCs division is about 10 month −1 and the level of activator, normalized to 1 for healthy tissue, is about 4 for unsorted cancer cells, becoming 46 for selected and sorted CSCs 28 . Introducing these to Eq. (1), we derive α/β ~ 0.2 and d ~ 0.3. In Table 1, experimental results have been collected. Even if the biological literature on the subject is very rich these last years, it is unfortunate that not so much quantitative physical or chemical informations are available.Finally in 29 , the inhibitor i has been identified explaining why patients with high level in lung tumor have a better prognosis.
Thus, the best strategy seems to target the chemical activators rather than the cancerous cells. Indeed, the presence of the cancer stem cells effectively drives the tumor, even in the optimistic case where we hopefully kill all the cancer cells, still CSCs will differentiate again and start a new cycle making the treatment very unlikely to succeed when the chemical activators are not controlled, (compare the two figures Fig. (6a,b) in the Appendix D.1 which shows how the basin of attraction of a stable point can be reduced). In medical centers nowadays, this strategy is commonly employed post-surgery to avoid the relapse of a new tumor at the same place of the primary one which has been withdrawn. Long-time drug administration over years aims at suppressing these activators.
The analysis of this model is very rich. However, it assumes an infinite geometry and the results must be confronted to the spatial structure of the host tissue. It is the reason why we now focus on the spatial evolution of a tumor seed, with initial conditions in the vicinity of the fixed point  2 .

Spatio-Temporal Dynamics of Tumor Growth
Population dynamics modeling gives only a first overview of tumoro-genesis, far from the real complexity driven by cell-cell interactions, friction with other cells or with the stroma and physical perturbations due to obstacles, that are always present in the host organ. In addition, in many cancers, the primary tumor, once detected, is not unique and may be accompanied by smaller satellites. So the spatial component cannot be ignored as well. As shown in 30 , for inhomogeneous in vivo tumor, the density of each constituent varies locally, and the spatio-temporal behavior is completely different from the one predicted by a dynamical system. It is even worse when a spinodal decomposition 19 occurs leading to a a phase separation of Cahn-Hilliard type inside the tumor 31,32 . Furthemore, the possible existence of niches, area of localization of S cells, introduces a new complexity that we will not consider here. They do not seem to be present in every tumor and are believed to be dependent on the patient's age, as reported by Ferraro and Lo Clesio 33 . Besides, most of the physical characteristics of the cancerous niche such as their localization and density in the host tissue have not been reported even in recent review papers as 6,8,26 which focus mostly on biological pathways.
Given the diversity of situations for human beings, a unique model cannot exist and the same strategy is adopted here for CSCs and DCs. Our dynamical system must be completed by adding the standard advective For the chemicals, a or m, we restrict on diffusion and The average velocity of cells → v i will be derived from the Rayleighian variation, counter-part of the Lagrangian for highly dissipative systems 34,35 , taking into account the constraints of the Galilean Invariance: in addition to S + D + C = 1 (see the Appendix E). Then, cell velocities result from proliferation and mechanical forces. CSCs being much less numerous and more individualist with a propensity to escape and make metastasis, we neglect their energy of cohesion. Thus, the free energy of the cell mixing E is dominated by the DCs 36 , and we only consider the mixing energy of DCs:   www.nature.com/scientificreports www.nature.com/scientificreports/ cell-cell interaction of CSCs with DCs, but the mixing energy is scaled by the square of cell densities giving a net advantage to the DCs, responsible for the tumor cohesion. The cellular elasticity is discarded: recent AFM 38,39 or nano-indentation experiments 40 indicate that cancerous cells are less stiff than their healthy counterparts by 70% or 80%, due to an alteration and regression of their cytoskeleton. Since inertia is negligible, viscous dissipation controls the dynamics which is determined by a variational process acting on the Rayleighian (for technical details see Appendix E). Taking into account the principle of extremum of this quantity with the incompressibility constraint gives the following velocities for → v D and → v S : From Eq. (8), we can define the cell-cell adhesion in units of τ which gives a length unit l 0 as the square root of the cell-cell adhesion multiplied by t 0 , the time unit, so hereafter all quantities will be dimensionless. For simplicity, we do not change the notations. Different options can be made for χ D but must respect physical objectives  www.nature.com/scientificreports www.nature.com/scientificreports/ such as balance of the stresses at the tumor border, itself defined by a border density D b although this condition assumes no density gradient at the border. In addition, at low densities, χ is close to zero, being negative and strongly increases for growing density above D b . A standard choice is the Cahn-Hilliard potential,(see Appendix E) but here, by reason of numerical convergence the potential introduced by Byrne and Preziosi 41,42 has been preferred. Its derivative reads: Heterogeneity under tumor growth. Adopting the same strategy as before, we extend the stability analysis of each dynamically stable fixed point,  1 and  2 , to spatially periodic perturbations such that each cell or chemical density becomes φ(x, t) = φ homo + (δφ)e λt cos(kx). Each Jacobian (Eqs (5) and (41)) is then dependent on k 2 and the characteristic polynomial λ λ = ∑ = P Ck ( ) ( ) j j j 0 3 2 which fixes the eigenvalues can be found in the Appendix F. Static and periodic patterns are found for λ = dλ/dk = 0 which imposes that C 0 = dC 0 (k)/dk = 0 as well.
At linear order, advection does not play any role for  1 since S 1 = D 1 = 0, and the diffusion of chemicals being stabilizing, this fixed point remains stable leading to a homogeneous tumor. As predicted by the linear stability analysis explained in Appendix F and shown in Fig. (3a,b), there exist patterns for 2  . The selected wavenumber of the spatial pattern is found in the neighborhood of , so for D 2 ~ 0.44 according to our choice given by Eq. (9). One can observe that low values of D 2 (i.e. below 0.44) make segregation more likely to occur. Patterns in a square are shown in Fig. (3c,d) and in Fig. (4) in a growing tumor. Notice that in the vicinity of the selected D 2 values, the numerical patterns exhibit either unstable labyrinths or dots. Dots may represent the nests (clinical vocabulary) observed in many tumors. Assuming the soft tissue to be homogeneous, the tumoral seed develops with a phase segregation inside, both for CSCs and DCs, (see Fig.(4a,b)). Surprisingly, despite the interior patterning, the contour remains rather well defined with no contour instability except when obstacles are introduced. With a small density of random obstacles, (see Fig. (4c,d)) the interface is tortuous but in (e) it surprisingly recovers regularity at a high density. CSCs are more located at border except the case of low level of obstacles. This CSCs concentration at the border of tumors is also observed in mice after injection of human CSCS 43 . The fact that CSC have a higher density at the border also indicates that they are good candidate for metastasis. Finally, in Fig. (5), 3 seeds are introduced to mimic more closely the clinical reality. Details on the numerics are given in Section (Methods), however one must mention that modelling of growing tumors is numerically much more challenging when not achieved in a square with periodic boundary conditions.

Discussion and Conclusion
In this paper, we presented, analyzed and simulated a simple multiphase model for tumor initiation driven by CSC cells. We focused both on the differentiation process of cancer stem cells and on the de-differentiation of differentiated cancer cells within the tumor. We analyzed the effect of chemical activators responsible for these processes and the possible therapeutic strategies for tumor extinction. Our conclusion is that, in this situation, the most effective strategies are the ones that target chemical activators rather than cells themselves, due to a vicious feedback loop between CSCs and DCs. New therapeutics emerge in this direction 29,44,45 . With numerical simulations, we understood how modifications of the environment, in which the tumor expands, can abruptly change its long time behavior. In the future, a clear understanding of the mechanical and physical properties of CSCs can improve the model, especially when considering adhesion energy, elastic deformation and cell-cell interaction 41,46 . In this spirit, a more complicated multiphase model can be constructed, taking into account immune cells, the host cells and their www.nature.com/scientificreports www.nature.com/scientificreports/ interaction with tumor cells. Because of a lack of information, the existence of niches is dogged here: they may be at the origin of localized sources of D cells inside the tumor and would require separation between mobile and static DC populations. Our model can be extended to other cellular cell types which can be reprogrammed in cancer initiated cells 47 or in cancer stem cells after responses to therapeutics 48 . From the theoretical viewpoint, we have treated here an ecological model where the competition between species depends on a singular tiny sub-population. It will apply probably to other communities than the cancerous cells, opening future directions in nonlinear physics.

Methods
The system corresponding to Fig. (6) was numerically solved with the XMDS2 software package 49 . We used an adaptive Runge-Kutta of fourth-fifth order with a tolerance 10 −5 . All spatial operators were evaluated via spectral methods. The simulations were performed on a (400 × 400) grid with a space discretization of dx = 10 −1 and with periodic boundary conditions for Fig. (3) which shows the emergence of different patterns with respect to the fixed point values  2 . The most relevant parameters d and α were varied, keeping fixed in the following the ensemble {η, ψ, q 0 , m 0 , σ, β, γ} = {1, 1, 1, 0.5, 0.05, 1, 1} and {ε D , ε S , ε C , Λ, D b , α D , D m , D a } = {0.1, 0, 0, 2, 0.6, 2, 1, 1}. In Fig. (3c) different choices for chemical degradation rates and cell death result in labyrinth-like patterns for both S and D. In Fig. (3d), upon increasing the value 2  , a honeycomb-like pattern emerges. The origin of the different patterns may be due to the highest value of D 2 which increases the effective strength of the cell-cell DCs attraction given by Eq. (9). In Fig. (4) an initially localized concentration of CSCs and DCs grows in space and time. Numerical simulations were performed for different boxes and initial conditions. In Fig. (4) from a to f, S, D are initialized as a small and localized population with random perturbations around  2 . Similar initial conditions apply for the chemicals, but distributed on the entire grid. Concentration fields for both S,D are evaluated up to t = 30 -the time scale being in unit of differentiation of the CSCs -and plotted against each other. In Fig.  (4a,b), CSCs and DCs expand in a homogeneous tissue. In particular, Fig. (4a) highlights the relevance of CSCs in this process. CSCs form a ring of higher concentration with respect to bulk values which effectively drives the tumor. In Fig. (4) from c to f, obstacles are introduced inside the box. Obstacles are tumor regions not accessible to tumor cells. In order to model them, simulations were performed on a grid with randomly distributed holes. Numerically, at each step of the iterative process, we cancel the cell density inside the holes corresponding to obstacles and then we integrate the following discretized equation . O i is a (1,0) binary variable, extracted from a Bernoulli distribution with a different mean: μ = 0.98 in Fig. (4c,d) and μ = 0.95 in Fig.  (4e,f). In Fig. (4c,d) the tumor grows in a tortuous fashion, breaking the cylindrical symmetry. This effect is two-sided: first, the addition of obstacles breaks the circular ring and second, the values of  2 are the same as the honeycomb pattern, giving rise to a strong effective attraction between cells. In Fig. (4e), despite a higher number of holes in the grid, the only relevant effect to the tumor growth is the breakdown of the circular ring of CSCs. The size and form of the tumor are qualitatively unaffected.
Our dynamical model presents 3 fixed points, 2 of them are of special interest since they concern a low concentration of cancerous cells. Their stability and basin of attraction have a crucial importance for treatments: they can induce the complete disappearance of the cancer cells. Each fixed point F i is defined by 4 concentrations values i  = (S i , D i , a i , m i ) which satisfy the dynamical system Eq. (1), without time dependence. The notation φ j is used for either S, D, a or m. Following a very classical strategy of dynamical systems 50,51 , each concentration φ is perturbed according to where δ j is a small quantity and λ j the growth rate of the perturbation. The stability is entirely governed by the Jacobians of our system, Eq. (10), having coefficient matrix given by ∂Γ k /∂φ j where Γ k has been introduced to represent dφ k /dt. Being local, the stability analysis must be achieved for each fixed point, independently of the others.

Appendix
A. Existence and stability of fixed points. The Jacobian matrix, calculated with the densities of each fixed point, is then:  = (0, 0, 0, m 1 ). This fixed point corresponds to complete eradication of the tumor. The Jacobian reads: For this setting, the eigenvalues are:  One eigenvalue obviously being λ 1 = −α, the three others are determined by solving the third order characteristic polynomial: . Since d < q 0 , as a condition for 2  to exist, and m 2 < γ/α, B is negative and all coefficients in the characteristic polynomial are then positive. In the case of a polynomial written as λ 3 + pλ 2 + qλ + r, with all coefficients p, q, r being positive, the condition for negative real parts of the 3 roots is simply pq > r according to the Routh-Hurwitz result 52 . This condition is always satisfied in our case. The second fixed point is then always stable and is an attractor. With the Sturm Theorem, it is possible to prove that there exists only one negative real eigenvalue and a pair of two complex conjugates eigenvalues with negative real part. Let us proceed now to the analysis of the last fixed point. 3  . The existence of the third fixed point corresponds to q(m 3 ) ~ 0 and p ~ 1/2, which allows to get an approximate value of a 3 by solving p(a 3 , D 3 ) = 1/2. It follows an algebraic second order equation for a 3 :

C. Existence and stability analysis of the fixed point
where ψ 0 = ψα/(dβ). A unique positive root can be found: , we can deduce also that:  The treatment is more difficult for this fixed point if we do not make the appropriate approximations which lead to its condition of existence. Indeed, without approximation, from Eq. (10), the characteristic polynomial is expected to be of 4th order, then requiring a numerical treatment. Assuming that ψ 0 is smaller than 1 but not too close to 1, we can consider that −ẽ 0 S S / 3 0 , since S 0 is a tiny quantity compared to = α β + S a a 3 (1 ) 3 3 . In addition, q(m 3 ) ~ 0 and p ~ 1/2. This simplifies the Jacobian which reads: where τ 0 = (1 + a 3 )(1 + ηa 3 )/a 3 2 and q′ = q′(m 3 ). We can conclude that, there is at least a negative eigenvalue −α. The characteristic polynomial for the eigenvalues of the third fixed point, independently of λ 0 = −α reads: which also indicates a saddle point. However, due to the number of independent parameters, the study of the roots of the characteristic polynomial of third order remains cumbersome. Indeed, the product of the 3 roots being positive indicates that at least one real eigenvalue is positive, suggesting that the third fixed point is always a saddle point, once it exists (ψ 0 < 1), independently of all parameters.

D. Numerical investigation of the dynamical system.
D.1 First case: d > q 0 q(m 1 ) ~ q 0 . We select numerical values of parameters such that both  1 and  3 exist, the first one being stable but in the neighborhood of  3 which is unstable. In this case and as shown in Fig. (1b),  1 is the only attractive fixed point which gives the optimistic issue of the full tumor regression. However, 3  can possibly prevent the cancer regression, even in such a favorable situation. We numerically focus on the set {α, β, d, ψ, η} = {1, 3, 2, 1, 1}, giving ψ 0 = αψ/(dβ) = 1/6 which is smaller than 1 as required for the existence of www.nature.com/scientificreports www.nature.com/scientificreports/ The value of σ is chosen to mimic conveniently a Heaviside function for q(m). The computed eigenvalues of the Jacobian from Eq. (10) corresponds perfectly to Eq. (23) which validates our approximated analysis of the stability of the third fixed point. It also proves that, for the determination of the eigenvalues which fixes the dynamics, the only pertinent parameters are {α, β, d, ψ, η} and that the precise values of {S 0 , σ, γ} have much less importance. For this particular choice, we have 3 directions of attraction, 2 being spirals and one being a repulsive direction, see Fig. (2d). As shown analytically, this is a saddle point. In Fig. (2a,b), the population of DC and CSC cells are given as a function of time in a log-scale where a time-scale relaxation can be estimated and compared to the smallest (in absolute value) growth rate predicted by Eq. (13), λ 1 (4) ~ −0.536, see Fig. (7a).
D.2 Second case: d < q 0 . Here,  1 will be repulsive, while 2  will be attractive and  3 is a saddle fixed point (c.f Fig.(1b) 24). The situation looks similar to the previous case but the main and important difference relies on the existence of two repulsive fixed points,  1 and 3  (see Fig. (1b)). As shown in Fig. (2b), one striking feature is the overshoot of the dynamical variables induced by the feedback activator m. The tumor will not go under regression most of the time, although  2 is an attractive fixed point, Fig. (2c). Moreover even when the system reaches a steady state, it can maintain a high level of differentiated cells if d is small. The tumor will not proliferate but may stabilize at dangerous levels, which is not an optimistic issue in practice.
Time dependent perturbation of the dynamical system: Stochasticity of the chemicals. We proceed by investigating the dynamics of the system numerically. We perform numerical simulations of Eq. (1) when q 0 (t) is an explicit time-dependent function. Of particular interest is the case where q 0 (t) oscillates around d, in the area of the bisector. In this regime, we would expect the system to exhibit an oscillatory behavior, due to the switching in the stability of the fixed points:  1 <=> 2  . In particular,  1 will change from a stable fixed point (tumor regression) to an unstable fixed point. The concentration of cells will then be pushed away from its decreasing value when at a particular time t 1 , q 0 (t 1 ) > d and pushed back toward regression when q 0 (t 2 ) < d. In Fig. (7), we show numerical results for q 0 = d + f(t), f(t) being a periodic function with period T. Despite the symmetry around d, the system spending equal time in the different regions of the phase space, the tumor tends to go under regression at long times. Intuitively, the system relaxes to the different attractors in each regions ( 1  and 2  ) with a time scale governed by the slowest decaying mode. We compare the slowest decaying modes in the two extreme cases, q 0 (max) = d + f max and q 0 (min) = d + f min . In the first case, being q 0 (max) > d, the stable fixed point is  2 , with slowest decaying mode: λ = maxλ i ~ −0.2 and associated time-scale τ = λ | |~5 1 . In the second case, q 0 (min) < d, the stable fixed point is  1 , with slowest decaying mode: λ = maxλ i ~ −0.04 and associated time-scale τ = λ | |~2 5 1 . We would expect a faster relaxation toward  2 than  1 , but this is not observed in Fig. (7). In particular, when the Jacobian matrices are not normal the instantaneous response -which determines the dynamic for a time varying q 0 (t)-is given by the eigenvalues of its hermitian part: H J = (J + J T )/2 53 . The hermitian matrix, H J , has positive eigenvalues in both scenarios, making an analytical understanding of Fig. (7) challenging. Considering a and m as a stochastic variables, our preliminary study indicates that the noise stabilizes  1 which seems to be selected rather than 2  even in area where 2  is stable without noise. This is especially verified above the bisector. Below the bisector, the answer is more complex but either  1 or 2  are selected. A more extensive study of stochasticity can be found in 54 . The results obtained for sinusoidal perturbations in the asymptotic may be due to the choice of initial data. e. Spatio-dynamical equations deduced from the Rayleighian principle. The Rayleighian  represents the energy variation per unit time for a system whose dynamics results only from dissipation, as for Stokes flow. When inertia is discarded, a variational principle of extremum of dissipation applies and the constraints acting on the system are included via Lagrange parameters Q i . It reads:  is the interaction energy between phases, the viscous dissipation  includes both interphase friction and viscosity and Ω is the tumor volume at time t. Here, only the constraint of global incompressibility is applied (via the function Q). Let us define more precisely each term in our case. The energy density  is the sum of the energy of each phase S, D and C, the energy of interaction between phases and a penalty imposed to strong density gradient variations: The energy variation is then: where τ is a mobility coefficient and μ the viscosity. The time variation of the free-energy is decomposed in two terms for a growing domain Ω(t): a surface (resp. a volume) integral calculated on Ω and a line (resp. a surface) integral along ∂Ω(t).
Summing the three equations of Eq. (31) gives ∇Q: Finally, assuming μ → 0 we derive: Classical choices for the free energy are inspired by the classical Cahn-Hilliard potential: which allows to determine order of magnitudes of velocities in the realistic limit: D ≫ S. Indeed, ∂χ D /∂D ~ α D D while ∂χ S /∂S ~ α S S and ∂χ DS /∂S ~ α DS D and χ DS /∂D ~ α S . It shows that we can keep only χ D and it reads: In this case, both velocities are opposite in direction with quite the same order of magnitude. If α DS ⋙ α d and α s , then both velocities have the same direction but the stem cells velocity dominates.

E.2 Boundary conditions at the tumor border.
First of all, it must be emphasized that our treatment is rather general and can be applied to any cell-cell interaction potential. In our simulations, no boundary conditions are required since our simulations correspond to an initial value problem and we do not take into account the exterior of the tumor assumed like a gel or a liquid. However, for completeness and future works, we briefly remind how they are determined from the variational process. They are helpful for selecting the densities D b and S b which enter the free energies Eq. (37). In few cases, as for melanoma, the cancerous cell densities at border may be measured and the experimental values are chosen for D b and S b . At the frontier, the interface velocity is the averaged velocity of the cancerous cells: if there is no phase separation. Indeed, it is possible that the mixing is fully inhomogeneous and that the phases split up from each other. Considering the terms which have been left by integration by part in the variation with respect to → v D , we are faced with a contour integral: where the first term comes from the contour integral of Eq. (29) while the last two terms come from the first integral of the same relation, once integrated by parts using Eqs (26) and (29). A similar equation is verified for S cells. Then we have 2 relations at the border: one gives the Lagrange parameter Q and the second, one (among 2) density value.
Assuming that the border is the same for D and S cells (full mixing, no phase segregation) imposes = . So boundary conditions cannot be checked without the help of boundary layers of different thicknesses for both the S and D populations. A naive relationship will be found for the thicknesses of both boundary layers: Having now the full spatio-dynamical equations, we can now investigate the stability analysis.

F. Stability analysis in the vicinity of dynamically stable fixed points.
A spatial perturbation is added to any constituent of our system: φ(x, t) = φ homo + δφ(t)cos(kx). The stability of the first fixed point will not change, for any wave-vector. For the second fixed point 2  , the Jacobian will be: J 2 k = J 2 + k 2 /τδJ, J 2 has yet been defined in Eq. (14) and δJ reads:  The characteristic polynomial of third order is then transformed into: with all coefficients being dependent of k. This equation accepts 3 roots λ and a spatial static instability is obtained for the double condition C 0 = −λ 1 λ 2 λ 3 = 0 and dC 0 /dk = 0 which gives us the value of the threshold for instability and the wavelength at threshold. Due to the number of independent parameters, it is more convenient to fix the 3 parameters: α, B, d 0 which partially define the initial fixed point  2 up to the value of D 2 , so all the curves C 0 (k) will start at the same value and will differ because of D 2 . A representation of a spinodal decomposition can be found in the Section (Spatio-Temporal), see Fig. (3c,d)) with the Byrne-Preziosi potential 41,42 .