Evolution of cancer stem cell lineage involving feedback regulation

Tumor emergence and progression is a complex phenomenon that assumes special molecular and cellular interactions. The hierarchical structuring and communication via feedback signaling of different cell types, which are categorized as the stem, progenitor, and differentiated cells in dependence of their maturity level, plays an important role. Under healthy conditions, these cells build a dynamical system that is responsible for facilitating the homeostatic regulation of the tissue. Generally, in this hierarchical setting, stem and progenitor cells are yet likely to undergo a mutation, when a cell divides into two daughter cells. This may lead to the development of abnormal characteristics, i.e. mutation in the cell, yielding an unrestrained number of cells. Therefore, the regulation of a stem cell’s proliferation and differentiation rate is crucial for maintaining the balance in the overall cell population. In this paper, a maturity based mathematical model with feedback regulation is formulated for healthy and mutated cell lineages. It is given in the form of coupled ordinary and partial differential equations. The focus is laid on the dynamical effects resulting from acquiring a mutation in the hierarchical structure of stem, progenitor and fully differentiated cells. Additionally, the effects of nonlinear feedback regulation from mature cells into both stem and progenitor cell populations have been inspected. The steady-state solutions of the model are derived analytically. Numerical simulations and results based on a finite volume scheme underpin various expected behavioral patterns of the homeostatic regulation and cancer evolution. For instance, it has been found that the mutated cells can experience significant growth even with a single somatic mutation, but under homeostatic regulation acquire a steady-state and thus, ensuing healthy cell population to either a steady-state or a lower cell concentration. Furthermore, the model behavior has been validated with different experimentally measured tumor values from the literature.


Introduction
A tissue structure is comprised of various cell types arranged in a hierarchy according to specific characteristics, properties and functionalities. Typically, stem cells have the inherent property of indefinite self-renewal and differentiation into specialized cells [1]. Self-renewal in a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 stem cells results in the production of the cells identical to the parent [2]. As sources of a lineage structure, stem cells produce progenitor cells via differentiation and their properties vary accordingly. For a given cell type, cell lineage has to undergo a number of maturity levels between the stem and differentiated cells. At the end of a cell line, the progenitor cells give rise to a mature cell population which does not possess the power to proliferate anymore, but can only experience apoptosis (the programmed cell death [3]). The specialised functions in the tissue are performed by mature cells, while the tissue homeostasis is preserved by regulating the ratio of stem cells' self-renewal rate to differentiation. According to tumor stem cell hypothesis [4], cancer invasion and maintenance is driven by a small number of cells possessing the properties of stem cells. It has been observed that cancer initiating cells are characterized by high proliferative potential, capability to differentiate into diverse phenotypes and strength to escape apoptosis [4,5]. In fact, these so-called "tumor-initiating cells" are stem cells that have acquired mutations [4], while the rest of the tumor cells are either mutated progenitor or differentiated cells. The latter can undergo apoptosis and are less likely responsible to invade and persist the tumor [6]. Therefore, it has been suggested to eradicate the cancer stem cells by treatment to completely eliminate the cancer [7]. This motivates particularly the study of stem cell dynamics and their role in the cancer evolution. In this sense, the present paper tends to develop a mathematical modeling framework, which is useful to predict the observed behavioral patterns of cancer evolution and, additionally, help in a purposeful impact by means of external inputs (e.g. radiation) which leads to mutation acquisition.
Tumor development results from acquiring mutations and escaping the enzyme-coded fixation process [8]. After acquiring a nonsense mutation, it can increase in number via cell division. Although not all mutations are harmful, certain mutations can contribute to malignant cell growth when acquired successively. While there exist various types of mutations, the ones which are crucial to cancer are characterized by enhanced proliferative potential, reduced apoptosis, genetic instability and reduced tumor suppression [5]. It has also been observed that typically one to ten mutations are required in a cell to revamp into a malignant one [5,9,10]. The mutated cells also possess a progeny, because these cells not only proliferate, but can also differentiate to successive cell types. In other words, there exists another hierarchical structure of mutated (i.e. cancer) cells besides the healthy one. Herein, the interesting aspect to study is the joint evolution of both progenies sharing the same environment.
The functionality of any multi-cellular organism as a whole depends greatly on the active feedback regulation process [11]. The loss of this homeostatic control escalates the growth of cells in the tissue which culminate in the advent of cancer. The precise nature of this feedback is not known [11]. In the literature, it is assumed that the mature cells secrete feedback signals which manipulate the stem cell's division strength in order to maintain the balance between its self-renewal and differentiation rate [11]. The escalating growth of the cell population may approach the steady-state due to the effects induced by the feedback [12]. Various cell lineage frameworks have been introduced in the literature to investigate the dynamics of tissue regulation via feedback loop [11][12][13]. For a structural inspection of the feedback in a system consisting of two different cell lines with distinct properties, it is necessary to consider a model of each sub-population. The latter is based on the assumption that in every lineage, there exist a chain of maturation stages, which is sequentially arranged [14,15].
Maturity represents a quantitative macroscopic measure which characterises cell differentiation. A variety of mathematical models have been formulated for explicit modeling of each cell subpopulation in a lineage using discrete [11,12,[16][17][18][19][20][21][22][23] cell maturity representations. The discrete modeling paradigm assumes that maturation occurs only during the division of a cell, yielding a sequence of maturation stages. However, it is becoming evident that differentiation inherently displays continuous transitions. For instance, in neurogenesis cell differentiation without cell divisions are reported [24] and in haematopoiesis, stem cell lineage commitment is described to be a continuous process [25]. Besides, in some tissues, such as the mammary gland, different differentiation stages are not well-identified [26]. In addition to these experimental evidences, a rationale for considering continuous maturity representation is provided by the fact that differentiation is controlled by intracellular biochemical processes, which are continuous in time, at least when averaged over a large number of cells. Finally, modeling continuous maturation is exceptionally relevant considering the heterogeneity observed in various cancers (e.g. in breast cancer [27]) and thus necessitates a continuous maturity structured modeling approach. In the literature, continuous maturity structured models have been proposed in [28][29][30], where maturity is defined by a continuous variable which stands for the remaining proliferative potential of the cell and its capability to perform cellular functions. In [29], a continuous maturity structured model of granulopoiesis has been developed using partial differential equations (PDEs) for bone marrow granulocyte precursors and ordinary differential equation (ODE) for the blood granulocytes. The population of stem cells is assumed to be constant. The proliferation and mobilization rates along with apoptosis were modeled as functions of cell maturity. The scaled maturity level lies between zero and one. While the authors focused on the identification of the fastest mutation sequences leading to emergence of the cancer, the feedback regulation from the mature cells was entirely neglected. Due to the lack of regulation, such structures produce unbounded growth of cell populations only, and in particular can not predict steady-state evolution. On the other hand, in [30], the authors have used a similar maturity based continuous model along with additional stem cell dynamics in the form of an ODE model. The model is rather general and supports hierarchical structures of cell lineage. As opposed to [30], we assign a separate sub-population to mature cells and introduce the feedback homeostatic regulation therefrom, which has been neglected in both [29,30]. Our model provides a generic framework to investigate the dynamics involved in the evolution of both normal and mutated cell populations under continuous maturation process and feedback regulation. The main motivation behind this model is to develop an insight into the process, while taking into account most relevant features of this multi-step process.
In the present paper, we consider the dynamical interaction of three different sub-populations: (i) the stem, (ii) progenitor and (iii) differentiated cells, while highlighting the effects of feedback regulation from the mature cell population. More specifically, we analyze the coupling of two progenies consisting of healthy and mutated cells, while our main interest lies in investigating the feedback regulation from the separately modeled dynamics of mature cell population into the stem and maturity structured progenitor cell populations. In our framework, the stem and differentiated cells are modeled using ODEs, assuming minimum and maximum maturity, respectively, while PDEs with continuous maturity distributions are used to predict the evolution of the progenitor cells. In particular, the differentiation rate is not assumed to be constant as in [30], it is rather considered to be a function of maturity. Although there exist several models with feedback regulation in the literature, to our best knowledge, this is a first attempt to cover the feedback regulation in a more generic framework of stem cell lineage with continuous maturity distribution along with the mutated cell lineage resulting from the mutation acquisition in healthy cells. Finally, it is also interesting to highlight that our mathematical model can predict the stem cell hypothesis, claiming that even a small number of mutated stem cells can invade the overall cell population.
The paper is organized as follows. In the sequel, we introduce our model of stem cell lineage for both healthy and mutated cell populations. Then, the steady-state solution are derived analytically in the following section. After that, we present the simulations of the developed model, its validation using the experimental data from the literature and the numerical scheme that has been implemented to find numerical solution. Finally, we conclude with the short outlook.

Mathematical model
The mathematical model of the stem cell line is rather complex as the cells vary continuously in course of maturation with time. In the present paper, instead of considering the evolution of net cell population, we split it into different sub-populations to account for their specific dynamics, as depicted in Fig 1. The very initial cell state, i.e., stem cells, has the potential to stay undifferentiated and not to divide frequently under the conditions of homeostatic regulation [1,31,32]. As a middle stage in the cell evolution from stem cells all the way to full differentiation, we discriminate the progenitor cells, which undergo proliferation at relatively high rates and give rise to the population of fully mature cells. In the process of maturation, the proliferative potential and mortality rate of progenitor cells vary drastically until the terminal differentiation. To capture these dynamical effects, it is necessary to consider the maturity distribution of progenitor cells, which is mathematically described by means of PDEs. The last transition stage in the cell line from progenitor cells refers to fully mature cells that are specialized to perform their functions in the respective tissue without further division, and undergo apoptosis after a short span of life.
The schematics of our model in Fig 1 depicts the possible interactions between subsequent cell types ensuing from symmetric/asymmetric self-renewal, mutation, differentiation and apoptosis. The two parallel cell lines refer to the healthy and mutated cells (perhaps cancer, if cells acquire a lethal mutation) with zero and one mutation, respectively. Notice that, we consider only one mutation to keep the model simple for investigation. The model can be scaled up easily to the acquisition of multiple mutations. The potential for self-renewal is labeled as the property only for the stem cells in both healthy and mutated states, while the differentiation of cells is undertaken by both stem cell and progenitor cell populations. However, apoptosis can occur at all transition states with a certain rate. Since the number of cells increases with each step of maturation [33], the evolution scheme of all cell states as described above may lead to abnormal growth tending towards an unbounded number of cells. To avoid such unrestrained behavior of cell growth, one has to introduce feedback regulation. The modeling scheme in Fig 1 enables investigation of the dynamical behavior (i.e., evolution and control) of the overall cell population with and without feedback homeostatic regulation. The discrete compartmental setting of the model facilitates the implementation and withdrawal of the feedback signals into and from the different transition states, respectively.
In the sequel, we explain the governing model equations for the cell lineage dynamics of healthy and mutated cells. Thereby, C 0 (t) and C 1 (t) refer to the number of stem cells with zero and one mutation, respectively. Similarly, P 0 (x, t) and P 1 (x, t) correspond to the progenitor cells with zero and one mutation, respectively, while M 0 (t) and M 1 (t) refer to the number of fully differentiated healthy and mutated cell populations. As the cells in the compartment of the stem and mature cells are assumed to behave alike, one can infer a modeling scheme in the form of ODEs. On the other hand, the healthy and mutated progenitor cells require a property space, with cell maturity as a property variable. Since all the intermediate transitioning cell states between stem and fully differentiated cells are modeled as progenitor cells, the cells therein continuously differentiate to higher maturity states. Thus, the progenitor cells in both healthy and mutated states require PDEs.
Stem cell population. Stem cells are assumed to possess zero maturity level and they define the boundary conditions for the progenitor cell population at minimum maturity. The primary characteristics of stem cells responsible for their evolution are self-renewal and differentiation. The self-renewal can occur in two different ways, symmetrically or asymmetrically. Either way, there is a probability of acquiring a mutation during the division process, this yields an influx into the mutated stem cell population. On the other hand, differentiation of stem cells without mutation acquisition results in healthy progenitor cells, and those with a mutation influence to mutated progenitor cells. The stem cell population increases by symmetric self-renewal only, whereas the other mechanisms, e.g., mutation acquisition and differentiation, cause a decrement. The dynamical behavior of the stem cells is then described by the following mathematical expressions The initial conditions of healthy and mutated stems cells are C 0 (0) = c 0 and C 1 (0) = c 1 , respectively. In the above equations, k 0 and k 1 are the proliferation rates of stem cell with zero and one mutation, respectively. The first term on the right-hand side in Eq (1) refers to symmetric self-renewal with probability a S 0 , which results either in a decrement in the stem cell population by one, if the stem cells acquire a mutation with rate m, i.e. À a S 0 ðsÞmk 0 C 0 , or increase the pool by one in case of no mutation, i.e. ð1 À mÞa S 0 ðsÞk 0 C 0 . The second term represents an asymmetric self-renewal of stem cells with probability a A 0 in which the stem cells decrease by one, while asymmetrically self-renewing and acquiring a mutation. The third term represents the differentiation of stem cells with probability a D 0 , which is always followed by a decrement in stem cell population by one. The resulting progeny from the differentiation of cells will influx into the progenitor cell population. Note that the feedback signal s has been introduced into the stem cell probability of self-renewal (symmetric/asymmetric) and differentiation to maintain tissue homeostasis. This feedback signal is determined by the mature cell population, as shown below in Eq (12). Finally, the fourth term describes the death of stem cells with a rate of d C 0 , which reduces the stem cell population by one. In Eq (2), the first three terms in a square bracket on the right-hand side describe the symmetric self-renewal, differentiation and death of mutated stem cells C 1 with the probability of a S 1 a D 1 , and d C 1 , respectively. The last two terms (in the second square bracket) correspond to the influx from healthy stem cell population C 0 as a consequence of mutations acquired during symmetric and asymmetric self-renewal, Eq (1).
Progenitor cell population. The maturity distribution of progenitor cells represented by P 0 (x, t) and P 1 (x, t) constitutes of all maturity stages between stem and mature cell populations with x as maturity variable. Obviously, R x 2 x 1 P i ðx; tÞdx, i = 0, 1, is equal to the number of cells between maturity x 1 and x 2 . The governing equations and the initial conditions for normal and mutated progenitor cells read: with initial conditions P 0 (x, 0) = p 0 (x), P 1 (x, 0) = p 1 (x) and boundary conditions for t > 0. The functions g 0 (x) and g 1 (x) stand for the differentiation rate of progenitor cells with zero and one mutation, respectively. On the right-hand side of Eq (3), the first and second terms in the square bracket represent the birth and loss of progenitor cells due to a mutation with the rate m 0 . The progenitor cells are assumed to acquire one mutation at a time. The proliferation rates β 0 (x, s) and β 1 (x, s) of healthy and mutated progenitor cells depend on the maturity level and tend to zero as the cells achieve the higher maturity level [34]. The third term describes the apoptosis of progenitor cells P 0 (x, t) with maturity dependent death rate μ 0 (x) and generally it gets higher as the cell matures. On the right-hand side of Eq (4), the first two terms in a square bracket represent the proliferation and death of the mutated progenitor cells with the rate β 1 (x, s) and μ 1 (x), respectively. The last term represents the influx from the healthy progenitor cells via mutation. Note that, the proliferation potential and rate of apoptosis for the progenitor cells is defined as function of maturity. The early progenitor cells have higher proliferation potential as compared to the late progenitor cells. On the other hand, as mentioned earlier the death rate of progenitor cells is meager for early progenitor cells and increases after differentiating to a certain maturity level [35,36], see Fig 2(a). Although this does not hold for all cell types but true for of haemopoietic cells. There are many choices which can be suitable for proliferation and death rates of progenitor cells. Here, we borrow from [30] the following functional forms for β i and μ i , where i = 0, 1, and b i and d i represent the maximum rate of proliferation and apoptosis, respectively. Furthermore, o b i represent the maturity level at which the progenitor cells proliferate at half of the maximum rate and r b i refers to the steepness of the proliferation switch. Similarly, the maturity at which progenitor cells die at half of the maximum rate is o m i and the steepness of the switch is r m i . Here, the feedback signal s is introduced in the proliferation rate of progenitor cells. The behavior of the functional forms of β i (for a fixed value of feedback, i.e., s = 1) and μ i are depicted in Fig 2(a). Next, we introduce the differentiation function g i (x) which describes the rate at which the cells mature. It is a strictly positive and continuously differentiable function. In the pool of progenitor cells, continuous differentiation takes place alongside cellular division (maturation process). In maturity-time representation, mitosis takes place at same maturity levels [37]. From the modeling viewpoint, it means that in an infinitesimal time interval (t, t + dt), a cell with maturity x either matures to level x + dx with probability g(x)dt or divides into two daughter cells with probability β(x, s)dt. The progenitor cell population is heterogeneous with respect to cell maturation velocities [37] and typically, the maturity rate decreases along with an increasing maturity level. In order to define g(x), we fix the maximum maturation velocity equal to one then we use a monotonically decreasing function of the following functional form g i (x) = exp(−h i x), i = 0, 1 in our model. Here, the differentiation functions g i (x) are bounded between 0 and 1, thus the parameters h i are to be selected in such a way that g i (x) should not get near zero within the specified range of maturity variable, i.e., x min � x � x max , see Fig 2(b). We assume lower differentiation potential for mutated cells because poor cellular differentiation is one of the important traits of cancer [38,39].
During a division process, progenitor cells can also undergo a mutation. In this model, the healthy progenitor cells P 0 (x, t) acquire a mutation with a mutation rate m 0 to either proliferate into a mutated progenitor cell population P 1 (x, t) or to differentiate into a mutated differentiated cell population M 1 (x, t).
Mature cell population. The mature cell population is constituted by all the cells that attain the maximum maturity level, i.e. x � . Here, maximum maturity is assumed to be same for both normal and mutated cells. Normally, cancers are graded as 'low-differentiated' and 'well-differentiated' cancers based on the level of maturation of the cells in an organ where the cancer arises. This implies that mutated cells are as mature as normal cells in 'well-differentiated' cancers. However, mutated cells may remain immature sometimes if mutations are taking place in early phases of maturation because mutated cells grow rapidly and divide before cells are fully mature. Nevertheless, it is plausible to assume that also the mutates cells can achieve the maximum maturation level. Mature cells do not possess any proliferating potential and only undergo apoptosis after a particular time. Therefore, mature cells cannot acquire any additional mutation and are specialized to perform their assigned functions in the tissue. The following equations describe the healthy and mutated density of mature cells represented by M 0 (t) and M 1 (t), respectively: In Eq (9), the first term on the right-hand side describes the inflow into healthy mature cells via terminal differentiation of progenitor cells with maturity rate g 0 (x � ), while the second term defines the apoptosis of the mature cells with a rate of d M 0 . The first term on the right-hand side of Eq (10) is the influx from fully differentiated mutated progenitor cells, and the second term represents the influx from healthy progenitor cells due to an acquired mutation. The last term involving the rate d M 1 refers to the death of mutated mature cell population. Feedback regulation. In the signaling mechanism among the cells, the growth response is modulated by cytokine proteins along with other proliferation regulating factors [11]. Cytokines bind to their specific membrane-associated receptors which results in the activation of signal transduction pathways [40]. Studies have shown that in order to maintain the number of cells in balance, these signals have to depend on the mature cell population [41,42]. The dynamics of cytokine signaling molecules S can be described by an ODE as: where υ is the maximum secretion rate of cytokine signals, δ S represents the natural decrement of the signals S, and γ is the rate by which the total mature cell population M = M 0 + M 1 (consisting of both healthy and mutated mature cells) regulate the cytokine signals. Substituting s = (δ S /υ)S and k s = γ/δ S , the above equation turns into Since the cytokine signals are typically secreted at a higher rate than that of proliferation and differentiation of the cells, these drift quickly to a steady state. Hence, using quasi-steady state approximation, the equilibrium state for the feedback signal intensity follows from Eq (11). This shows that in the absence of mature cell population, the signal intensity is maximal, i.e., s = 1, and it drops to a minimum with a significant increase in the mature cell population. The nature of this signal can be understood as a controlling parameter that identifies the need for proliferation based on how many cells are present in the vicinity. These signals get weaker and weaker with a larger number of cells, i.e., there are not enough resources critical for the division process. The probabilities a S i ðsÞ, a A i ðsÞ and a D i ðsÞ in Eqs (1) and (2) and the maximum birth rates b i (s) in the birth functions β i (x, s) of progenitor cells in Eqs (3) and (4) are assumed to change linearly in s, cf. Eq (7), given that slopes are positive which leads to the following linear forms a D i and � b i represent maximum symmetric self-renewal, maximum asymmetric self-renewal, maximum differentiation of stem cell and maximum birth rate of progenitor cells, respectively.

Steady-state solutions
In this section, we derive the steady-state solutions of our governing Eqs (1)- (10). For the sake of convenience, we hereby simplify the notation, reading for the populations of healthy and mutated cells, respectively, with γ 00 , γ 10 and γ 11 defined as In a similar manner, the PDEs that describe the progenitor cells read: where γ 0 (x, s) and γ 1 (x, s) are given by g 0 ðx; sÞ ≔ ð1 À 2m 0 Þb 0 ðx; sÞ À m 0 ðxÞ; g 1 ðx; sÞ ≔ b 1 ðx; sÞ À m 1 ðxÞ: Finally, the equations of mature cell population stay same as before in Eqs (9) and (10). In the sequel, we assume the following conditions: x � � � RÞ; m 0 ðxÞ; m 1 ðxÞ 2 L 1 ð½0; x � �Þ g 00 ðsÞ ≔ g 00 ðMÞ is a decreasing function; i:e:; g 00 ðþ1Þ < 0: To address the question of the existence of any steady-state under a homeostatic regulation, we need to solve the following system of equations, for the steady-state unknowns � C 0 , � C 1 , � P 0 , � P 1 , � M 0 and � M 1 : with � g 0 ðxÞ ≔ g 0 ðx; � MÞ, � g 1 ðxÞ ≔ g 1 ðx; � MÞ and � b 0 ðxÞ ¼ b 0 ðx; � MÞ and the boundary conditions given as The trivial steady-state, i.e., � (26). However, the system also admits a non-trivial steady-state under the assumption γ 00 (0) > 0. In this case, from Eq (21) we get immediately Now, using Eq (15) in the above relation, we obtain: where the probabilities a S 0 ð � MÞ, a A 0 ð � MÞ, and a D 0 ð � MÞ are defined as and � a S 0 , � a A 0 ; � a D 0 2 R >0 . By employing the above relations in Eq (30), we derive the relation for � M, which is: Next, to solve the ODEs (23) and (24) for progenitor cells, we compute the boundary conditions at the final maturity, i.e., x = x � from the Eqs (25) and (26) The steady-state values result by solving the ODEs (23) and (24) for the healthy and mutated progenitor cells � P 0 ðxÞ and � P 1 ðxÞ:

PLOS ONE
Further, we use the boundary conditions given in Eqs (27) and (28) to compute the steadystate values of healthy and mutated stem cells, respectively where Eventually, we derive the steady-state relation for mutated mature cells � M 1 from Eq (22): where l 2 ¼ 2m 0 ð1 À mÞe f 1 ð0Þ þ me f 0 ð0Þ À Þ e f 0 ð0Þ : Note that, the steady-state relation for healthy mature cells � M 0 can be easily determined utilizing Eqs (32) and (39). From the above derivation of steady-states, we can summarise the following observations. The steady-states of our coupled nonlinear model cannot be defined explicitly, but the sum of the steady-states of healthy and mutated mature cells, used to compute feedback can be represented by an explicit relation. Moreover, the steady-states of stem and progenitor cells highly depend on the steady-state of mature cells due to the feedback inclusion.

Model simulations
In this section, we present the model simulations for illustration purposes. The initial states and the used parameters are given in Table 1. The forthcoming results are computed by the numerical scheme given at the end of this section. The maturity variable x belongs to [0,5] with the value of maximum maturity x � set to be 5. The step sizes for time Δt and maturity Δx used in simulations are 0.01 and 0.05, respectively. The behavioral patterns of the model are investigated hereby with the objective to observe the evolution of all six sub-populations with the feedback regulation, which is determined from the total number of both healthy and mutated mature cells. In general, after acquiring a mutation, the mutated cell gains fitness and thus differ considerably from healthy cells [17]. Therefore, the probability of mutated stem cells to self-renew is greater in simulations relative to the healthy ones, while the death rate is reduced, see Table 1. The feedback influences the the probabilities of symmetric/asymmetric self-renewal and differentiation in stem cells whereas in progenitor cells, the feedback is influencing the maximum proliferation rate b i in the birth function β(x, s). The simulations are initialized with healthy stem cells as c 0 ≔ C 0 (0) = 18000mL −1 , while all the rest of the subpopulations have been set to zero. Initially, the feedback signal is maximum, i.e., equal to one, because no differentiated cells exist, while over the course of time, the increase in the healthy and mutated mature cell population leads to a reduction of the feedback signal, as shown in  Table 1. The exponential growth in healthy stem cell population results in the increase of all healthy and mutated cell types Fig 4(a). A steady-state is achieved in healthy stem cells and consequently in all other sub-populations, see Figs 4 and 5. The achieved steady-states coincide with the analytically calculated steady states of all cell states in Eqs (35)- (39). The feedback signal plays the central role in the stabilization of the model states. In the absence of this feedback signal, the exponential growth continues and thus results in an unnatural number of cells. Fig 6 depicts another behavior of the model in which all the parameters used are same as in Table 1 but the symmetric self-renewal rate of stem cells � a S 0 ¼ 0:175 and their death rate d C 0 ¼ 0:016 day −1 . Contrary to Fig 4, the healthy stem cells C 0 in Fig 6(a) start decreasing after a gradual increase for a while and eventually land to a very low number. Indeed, similar behavior has been shown by the healthy mature cells in Fig 6(c), whereas the mutated stem and mature cells still attain their respective steady-states, as shown in Fig 6(a) and 6(b). In accordance with the dynamics of healthy and mutated stem cell populations described in Eqs (1) and (2), respectively, the probabilities of symmetric/asymmetric self-renewal and differentiation rates are influenced by the feedback signal. The rapidly increasing healthy and mutated Table 1 (12). It can be easily observed from Eq (1) that, as the feedback signal drops, the probabilities of symmetric/asymmetric stem cell self-renewal a S 0 =a A 0 and differentiation a D 0 also decrease. Note that these probabilities vary linearly with feedback signal s

PLOS ONE
having a positive slope. As a consequence, with the temporal evolution of healthy stem cells, death rate dominates over self-renewal, and healthy stem cells start declining in number, cf. Eq (1). In Eq (2) for mutated stem cells, the decline in feedback signal reduces the probability of self-renewal and differentiation in mutated stem cells. Nevertheless, the mutated cells still manage to grow in a higher number due to the increased fitness as stated before and thus approaching a steady-state. This scenario, in which the healthy cell line deteriorates and only mutated cells prevail over time, could also be called 'pure cancerous steady-state'. The establishment of the steady-state via feedback regulation has already been suggested in the literature [11], where the authors have considered ODE settings for the discrete cell populations of the stem cells all the way to the differentiated cells. Moreover, it is mentioned that the whole dynamics of stem cell lineages can be controlled by a single negative feedback loop, i.e., cytokine signaling.  Table 1. It can be seen that with any number of initial healthy stem cells C 0 (0), the steady states are achieved at the same time in healthy stem C 0 (t) and mature M 0 (t) cells, Fig 7(a) and 7(c); however, the magnitudes of the steady states are different. On the other hand, in mutated cell populations, the effect of different initial conditions is reflected only in the rates at which the steady states are achieved, while the magnitude of steady-state is the same for all initial values, see Fig 7(b) and 7(d). It implies that no matter how many healthy stem cells are there at any particular age, the subsequent mutations can lead to a substantial amount of mutated cell populations.
Feedback signal as Hill function. Here we want to analyse the model behavior when we define the relation between feedback signal s and total number of mature cells M using Hill function as compared to the behavior produced by the relation in Eq (12). Since increase in the concentration of mature cells represses the feedback signal, we define their relation using the Hill function as follows: where K M is the mature cell concentration (2.6 × 10 6 mL −1 ) at which feedback signal is half a maximum and n is the Hill coefficient. Note that the Eq (40) coincides with Eq (12) Table 1, whereas in Fig 8(e)-8(h) the parameter values which have been varied are symmetric and asymmetric self renewal rate of stem cells as � a S0 ¼ 0:175 and � a A0 ¼ 0:6650, respectively. It is to be noted that the dynamics of the stem and progenitor cells are maintained in homeostasis by inducing feedback only in the self-renewal rates. One can also achieve homeostasis by inducing the feedback in the death rates [49]. However, death rates are kept constant in the current paper; see Table 1. Moreover, in our model, the division rates of stem cell populations {k i , i = 0, 1}, are not depending on feedback signal because it has been validated in [11,50] that an efficient control mechanism underlies the modulation of self-renewal and differentiation rates as compared to the maintenance of proliferation rates in stem cells. It is evident from the simulation results as well that even without feedback regulation in division rates, the mutated stem cell population does achieve a steady-state.

Model validation
In this section, we validate the behavior of our model with different experimental measurements taken from the literature. In Fig 9, we use the tumor volume measurements for three different cancers, namely prostate, breast and colon for validation purpose. The experimental data sets are taken from [51][52][53]. These data sets are obtained by establishing human tumor xenografts in mice. The only measurements available are tumor volumes from the day of implantation and during exponential growth. We fit our model's parameters for validation purposes since we do not have the experimentally measured parameter values from these experiments. To compare our model behavior with the tumor volume, we first compute the total number of mutated cells N(t), which is the sum of mutated stem, progenitor, and mature cells, as The tumor volume calculated from the cell count of the proposed model (blue lines) fits very well to the experimental data (black marks) in all three scenarios, see Fig 9. The grey shaded regions depict the predicted model behavior before and after the available experimental values. Our model attains a steady-state in all three simulations due to the feedback via cytokine signals. Note that, the steady-states may vary in reality for different cancer types and also individually but the proposed model is flexible enough to depict various steady-state scenarios. The healthy cell lineage is considered to be initially at a steady-state. The parameters used in Fig 9 are given in Table 2.  Table 1 are given in Table 3. The exponential growth and achievement of the steady-state requires enhanced proliferation and selfrenewal rates of stem cells. The sum of the probabilities for symmetric self-renewal, asymmetric self-renewal and differentiation is still kept equal to 1.

Implemented numerical scheme
In this section, the numerical method used to solve the governing nonlinear Eqs (1)-(10) is presented. We apply already developed finite volume method (FVM) with central upwind scheme for the flux approximation on hyperbolic partial differential equations in MATLAB. The domain of the problem has been discretized in both, space and time. The timeline is discretized into N t steps with equidistant interval Δt = t k + 1 − t k . The spatial stepsize is given by where N x is the maximum number of spatial nodes given by The discretized progenitor cell density associated with the j th spatial interval at time k reads The necessary Courant-Friedrichs-Lewy (CFL) condition for convergence of the solution requires max x2fx j g g i ðxÞ Dt Dx � 1: The PDEs (3)-(4) are hyperbolic in nature and with the discretization defined above, we can implement the following algorithm to solve the coupled differential equations.
First, the initial conditions are given as For each time step k, the feedback from mature cells is calculated as Then, the stem cell population at time t k+1 can be discretized as follows and the following relation result for healthy and mutated stem cells, respectively The boundary conditions for progenitor cells at j = 0 are given as P kþ1 The discretization of the PDEs concerning the density of the progenitor cell populations is given in accordance with the central upwinding scheme as follows Finally, the discretized ODEs for mature cells are given as following which also involves the influx from progenitor cells P k 0 and P k 1 at the maximum maturity x = x � . The mature cell populations M 0 and M 1 will manipulate the feedback in the next time step and consequently, feedback will alter the dynamics of stem, progenitor or both cell populations to stabilize the exponential growth.

Discussion
The model predicts the evolution of healthy and mutated stem cell lineages in various case studies. Both lineages evolve with time and achieve a steady state under homeostatic regulation. According to the stem cell hypothesis, the persistence of cancer is regulated by a small number of cancer cells which share the same biological properties as the stem cells [4]. The results of this model are in accordance with the stem cell hypothesis because the mutated stem cells are responsible for the evolution and persistence of the whole mutated cell lineage due to their elevated self-renewal and differentiation potential, see Figs 4 and 5. Moreover, it is comprehensible that an efficient feedback mechanism must exist with heavy cross-talks between the cells themselves and the extracellular environment to robustly regulate the system comprising of cell lineages with various cell types. In the proposed model, the feedback in both the stem and progenitor cell populations, allows mimicking the intercellular interactions among the cells. Besides, the model provides an insight that the self-renewal rate of stem cells is a very sensitive parameter for the persistence and maintenance of healthy cell lines. As shown in Fig  6, the feedback's influence led to the extinction of the healthy cell line because the self-renewal rate of stem cells was reduced. In a nutshell, the critical ratio of stem cells' self-renewal rate to the differentiation rate should be preserved to maintain homeostasis in the healthy cell population.
The proposed model has some drawbacks too. First, the model assumes a single mutation leading to cancer evolution, which might not be accurate in many cases, but the model structure is flexible to incorporate more mutations and to predict the evolution of cancer depending on their individual effects. Secondly, the model assumes only cytokines feedback signaling but of course, there exist several feedback mechanisms, for instance, chalone [11], mechanosensing [58] etc. Furthermore, the simulations are performed assuming linear dependence on the feedback signal, which might not be very realistic; however, the precise nature of this feedback is still unknown.

Conclusion
We propose a generic modeling framework to investigate the coupled dynamics of the healthy and mutated cell lineages, entailing homeostatic regulation. We show that the model predicts familiar behavior and evolutionary patterns of cancer. For instance, the small number of mutated stem cells is responsible for evolving the whole mutated cell lineage. Moreover, the healthy cell line significantly declines in number due to the sensitivity of the symmetric selfrenewal rate of stem cells. Thus, the symmetric and asymmetric self-renewal rates of stem cells are crucial for the persistence and maintenance of both cell lines. The model is also validated with different experimental measurements of the tumor available in the literature. Concerning future work, it is possible to extend our model to include additional phenomena, e.g., cell dedifferentiation and other feedback mechanisms. Its architecture also enables heterogeneous type mutations to be introduced, which can be of interest to gain additional insights into the development of cancer and, additionally, in the faster emergence of cancer. An appealing future step concerns the stability analysis of the dynamical behavior and sensitivity analysis for the process parameters.