A Simple Mathematical Model Based on the Cancer Stem Cell Hypothesis Suggests Kinetic Commonalities in Solid Tumor Growth

Background The Cancer Stem Cell (CSC) hypothesis has gained credibility within the cancer research community. According to this hypothesis, a small subpopulation of cells within cancerous tissues exhibits stem-cell-like characteristics and is responsible for the maintenance and proliferation of cancer. Methodologies/Principal Findings We present a simple compartmental pseudo-chemical mathematical model for tumor growth, based on the CSC hypothesis, and derived using a “chemical reaction” approach. We defined three cell subpopulations: CSCs, transit progenitor cells, and differentiated cells. Each event related to cell division, differentiation, or death is then modeled as a chemical reaction. The resulting set of ordinary differential equations was numerically integrated to describe the time evolution of each cell subpopulation and the overall tumor growth. The parameter space was explored to identify combinations of parameter values that produce biologically feasible and consistent scenarios. Conclusions/Significance Certain kinetic relationships apparently must be satisfied to sustain solid tumor growth and to maintain an approximate constant fraction of CSCs in the tumor lower than 0.01 (as experimentally observed): (a) the rate of symmetrical and asymmetrical CSC renewal must be in the same order of magnitude; (b) the intrinsic rate of renewal and differentiation of progenitor cells must be half an order of magnitude higher than the corresponding intrinsic rates for cancer stem cells; (c) the rates of apoptosis of the CSC, transit amplifying progenitor (P) cells, and terminally differentiated (D) cells must be progressively higher by approximately one order of magnitude. Simulation results were consistent with reports that have suggested that encouraging CSC differentiation could be an effective therapeutic strategy for fighting cancer in addition to selective killing or inhibition of symmetric division of CSCs.


Introduction
Fundamental and applied clinical research into cancer could greatly benefit from mathematical models that contribute to the basic understanding of this disease, to the planning of more efficient therapeutic strategies, or to the generation of accurate patient prognosis. This paper presents a general, simple, and flexible mathematical model, mechanistically based on the Cancer Stem Cell (CSC) hypothesis, that is capable of reproducing the dynamics observed during the exponential growth of a tumor.
Recently, the CSC hypothesis has gained credibility within the cancer research community [1][2][3][4][5]. In its simplest version, this hypothesis postulates that most tumors (if not all) arise by consecutive genetic changes in a small subpopulation of cells that have intrinsic characteristics similar to those of normal stem cells (SCs) [6][7][8][9]. A fast growing body of experimental evidence suggests that these so-called cancer stem cells (CSCs) are the drivers of cancer and are responsible for sustained tumor growth. Although no general consensus has yet been reached on several key aspects of the biology of CSCs, there is agreement in some of their distinctive features: (a) self-renewal capabilities, (b) potential for differentiation into the various cell subtypes of the original cancer, and (c) increased tumorigenesis [9][10][11][12][13][14].
Although tumor growth has been a subject of intensive mathematical modeling in the last two decades, the concept of existence of a CSC population within tumors has been only recently included as an element in describing tumor growth [30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45]. Among these examples, different modeling approaches have been used, ranging from stochastic [35,42,45] to deterministic modeling [37,41]. CSC-cancer modeling has frequently focused on the exploration of therapeutic strategies [36,37,41,43]. For instance, Dingli and Michor [36] used mathematical modeling to demonstrate the importance of selective targeting of CSCs to improve the efficiency of cancer therapies. Similarly, Ganguly and Puri [39] formulated a model to evaluate chemotherapeutic drug efficacy in arresting tumor growth, based on the cancer stem cell hypothesis. Their results suggested that the best response to chemotherapy occurs when a drug targets abnormal stem cells. CSC based mathematical models have also been used to forecast the effect of specific therapeutic agents (and combinatory therapies). Several contributions have explored different aspects of the treatment with imatinib [37,41,43]. Mathematical modeling has also been used to gain understanding of fundamental issues underlying CSC biology [31,32,40,42,44,45]. The biology of CSCs has not been fully elucidated and many questions still remain unresolved [16,45]. In particular, some of these uncertainties are related to the dynamics of tumor growth. As an illustration, little is known about the balance between the multiple and complex cellular events that occur during the early stages of tumor progression. One of the central objectives of this work is to identify if some commonalities (or universal features) may exist with respect to the kinetics of early tumor growth. Experimentally studying the balance between the different cellular events involved on the process of tumor growth is not a trivial matter. Mechanistically based mathematical modeling might be highly useful for simulating the dynamics of cancer initiation and progression, the response to different therapies, and the evolution of resistance to drugs [30], as well as for gaining further fundamental understanding on the underlying dynamics of tumor growth.
In the present manuscript, we present a simple mathematical model that is designed to study the role of CSCs in tumor growth, with the aim of understanding the kinetic relationships between the different processes leading to exponential growth in solid tumors and evaluating possible therapeutic strategies for cancer treatment. We attempted to capture the key features of the known biological behavior of CSCs in a pseudo-chemical model, where cell division and death of the three cell subtypes considered are represented as ''chemical reactions.'' The intrinsic rates at which these reactions (cellular events) occur are the parameters of the model (k j ) and are analogous to reaction rate kinetic constants. Based on an exploration of the parameter space of these kinetic constants, we derive conclusions related to their relative magnitudes. Some inferences regarding the fundamental biology of tumor growth and the effectiveness of some therapeutic strategies against cancer are discussed.

Methods
A pseudo-chemical model for tumor growth: underlying biological concepts Tumors are a heterogeneous mix of cells, some of which exhibit SC-like characteristics [14,16,17,34,[46][47][48]. It is probably more accurate to say that a tumor possesses a continuous spectrum of cell types, ranging from CSCs to more differentiated cells. In most of the previous modeling studies, the complexity of tumor tissue has been addressed by defining several cell subpopulations (typically from two to four), leading to compartment models [31,32,37,40]. In order to reduce the complexity of the resulting model, only three subtypes of cells are considered: CSCs, transit amplifying progenitor cells (P), and terminally differentiated cells (D) (Fig. 1). This assumption is consistent with several experimental reports that simplify the cell heterogeneity found in cancer, in which three main cell subtypes are indentified [23,48] with some variants in nomenclature; i.e., holoclones, meroclones, and paraclones [49][50][51].
In our model, events related to CSC self-renewal, to maturation of CSCs into P cells, to further differentiation to D cells, and to death of all cell subtypes, are represented as ''chemical reactions'' and are mediated by specific rate constants. These reactions occur in a system that has no nutrient limitations during the phase of exponential tumor growth. This assumption presumes that angiogenesis occurs at a rate that ensures the accessibility of nutrients sustain constant growth.
Framed in this way, the time evolution of all cellular subpopulations can be represented by a set of ordinary differential equations that have an analytical solution. In the following paragraphs, we introduce each of the cellular events considered for the construction of the model, and their representation in the form of ''chemical reactions.'' Expansion of SCs can be accomplished through symmetric division [52,53], whereby one CSC originates two CSCs: Alternatively, a CSC can undergo asymmetric division (whereby one CSC gives rise to another CSC and a more differentiated progenitor (P) cell). This P cell possesses intermediate properties between CSCs and differentiated (D) cells [2,51,[54][55][56]: Both symmetrical and asymmetrical cell divisions of CSCs have been experimentally documented by staining of nuclear Oct-4 (a stem cell marker) [57]. Regulation of the ratio between symmetric and asymmetric division might possibly be crucial for the development and progression of cancer [44,52]. CSCs may also differentiate to P cells by symmetric division [44,58]: P cells can either self-renew, with a decreased capacity compared to CSCs, or they can differentiate to D cells [21,31,59,60]: D cells do not have the capacity to proliferate [20,51,62], so their corresponding kinetic constant should be very small (here considered negligible). In addition, all cellular subtypes can undergo cell death: To establish the model, we followed a classical strategy used in chemical reaction engineering to describe a system of chemical reactions. For each cellular event, a ''reaction rate'' can be successively established: r 1 = k 1 CSC; r 2 = k 2 CSC; r 3 = k 3 CSC; r 4 = k 4 P; r 5 = k 5 P; r 6 = k 6 CSC; r 7 = k 7 P; r 8 = k 8 D.
Neglecting all terms related to transport of cells (to and from the tumor), the rate of accumulation of a cell subtype ''j'', dj/dt, is equivalent to the net generation of that cell subtype, given by the addition of all reaction rates where that cell subtype is involved (i.e. produced or consumed).
For example, the cell species CSC is involved in the cellular events R1, R2, R3, and R6. In reaction R1, CSC is produced at a rate equivalent to 2r 1 , and consumed at a rate r 1 . In the cellular event R2, CSC is produced at a rate r 2 , and consumed at the very same rate r 2 . In the cellular event R3, CSC is consumed at a rate r 3 due to differentiation into P. Similarly, in the cellular event R6, CSC is consumed at a rate r 6 due to cell death.
The accumulation of CSCs within the system will be given by Eq. 1: And, This simple system of ordinary differential equations can be analytically solved to obtain the populations of each cell type (CSC, P, D) at any time, provided that a set of initial cell populations is specified (CSCo, Po and D).
CSC~CSC 0 e at ð4Þ Where The total number of tumor cells (N) will be the sum of the three cellular subtypes, assuming that dead cells are reabsorbed: The tumor volume can be calculated, considering that the effective volume contribution of a spherically shaped cell in a spherical tumor is 4.18610 26 mm 3 /cell [62], as: This assumption implies that the tumor grows at constant cellular density, which might be reasonable during exponential tumor growth, when we have assumed that no nutrient transport limitations exist, and space constraints are not significant [63].

Significance of the parameters and constraints of the model
The model has one kinetic parameter per cellular event (or ''reaction''). A brief discussion of the physical significance of these kinetic constants is pertinent here. In a typical elementary chemical reaction, the rate of appearance of a chemical species is proportional to the concentration of the reagents through a proportionality constant, the specific rate of reaction. Analogously, for our cellular system, the ''rate of reaction'' of each cellular event depends on both the number of precursor cells for that event and the proportionality constant k j that will multiply that number. For example, the rate of disappearance of CSCs, due exclusively to the occurrence of the cellular event R1, is r 1 , and is mediated by the kinetic parameter k 1 . Therefore, k 1 is an intrinsic reaction rate constant that indicates the natural predisposition of a cell, in this case a CSC, to divide symmetrically to originate two CSCs. In our model definition, k j is not equal to the growth rate of a cell subtype, but rather to the intrinsic proliferation rate associated with the frequency at which that particular cell subtype generally divides. To calculate the rate of accumulation of a particular cellular species (or net growth rate of that species), all terms where this species appears or disappears should be considered. For example, for the particular case of CSC, the associated rate of accumulation involves r 1 , r 3 , and r 6 (see equation 1).

Kinetic parameter grouping
The proposed kinetic model has eight independent variables. We simplified the analysis of the parameter space by grouping the eight kinetic constants of the model into three groups. Group I includes k 1 , k 2 , and k 3 ; kinetic constants associated with cellular events that relate to CSC proliferation and differentiation. Group II includes k 4 and k 5 , since they mediate cellular events related to proliferation or differentiation of P cells. Finally, k 6 , k 7 and k 8 were included in Group III, as they describe cell death of each cell subtype.
In addition, some mathematical relationships among parameter groups and kinetic constants were defined. The value of k 1 , the rate constant for symmetric CSC renewal, was set equal to one. This is convenient, since all the rest of the k j values can now be defined (or scaled) relative to k 1 . In addition, it is helpful to define ratios between the kinetic parameters, namely W 2/1 , W 3/1 , W 4/1 , W 5/4 , W 6/1 , W 7/1 , and W 8/1 . For example, W 4/1 is the ratio between k 4 /k 1 ; biologically, it reflects the relative magnitude of the intrinsic kinetic constant governing symmetrical division of progenitor cells with respect to that related to symmetrical division of cancer stem cells. Similarly, W 5/4 (k 5 /k 4 ) indicates the relative magnitude between the rate constants associated to symmetrical division of progenitor cells to render two progenitor cells (R4) and symmetrical division of progenitor cells to produce two differentiated cells (R5). In this way, a vector of W i/j values will define a complete set of k j values, and therefore a biological scenario for tumor growth. Illustratively, once k 1 is set to the unit value (k 1 = 1), the vector Although, a priori, all k j values are plausible, some constraints based on biological knowledge can be applied. For instance, in the present paper, the intrinsic apoptosis rate was considered to increase as cells become progressively more differentiated. Consequently, the greatest death rate corresponds to the most differentiated phenotype (D). On the other hand, CSCs have an extremely low apoptotic index [21,53,[64][65][66]; therefore, the next constraints will be imposed for all simulations: In addition, experimentally, a minimum fraction of CSCs has been found to be maintained through the evolution of cancer [20,61]. Normally, this fraction is lower than 1% of the total cell number [20,25,67]. Therefore, CSC/N,0.01 and d[CSC/N]/dt<0 for all times. Finally, the fraction of P cells can be estimated from experiments in the literature to be approximately 0.2 [20,68]. It is also well known that D cells constitute the majority of tumor cells [7,20,37]. Accordingly, we set the following expressions as constraints: P=N&0:2 and D=N&0:8 ð15Þ The biology of CSCs has not been fully elucidated and many questions still remain unresolved [16,45]. In particular, some of these uncertainties are related to the dynamics of tumor growth.
As an illustration, little is known about the balance between the multiple and complex cellular events that occur during the early stages of tumor progression. One of the central objectives of this work is to identify if some commonalities (or universal features) may exist with respect to the kinetics of early tumor growth. Experimentally studying the balance between the different cellular events involved on the process of tumor growth is not a trivial matter.
The simple model that we proposed here allows for the study of the effect of variations in the relationships between the intrinsic kinetic values of each one of the cellular events defined (from R1 to R8). In our discussion, we place particular importance on the experimentally documented fact that the CSC fraction in a tumor is constant during tumor evolution, which is clear evidence of the crucial role of the CSC reservoir in tumor growth [20,22,25,51]. Therefore, the constraint d[CSC/N]/dt<0 becomes central to identifying biologically consistent and feasible solutions for the model.

Feasible model solutions
In principle, one would expect that a vast number of combinations would render dynamical behaviors that would be consistent with experimental observations of the evolution of solid tumors. Based on our experience testing the model, the imposed constraints (namely d[CSC/N]/dt<0; P/N<0.2; D/N<0.8; CSC/N,0.01) substantially limit the number of sets of parameter values that lead to feasible solutions.
As an illustration, let us consider the particular solution, obtained when the vector W i/j = [W 2/1 ,W 3/1 ,W 4/1 ,W 5/4 ,W 6/1 ,W 7/1 ,  Figure 2B). Indeed, while exploring feasible solutions, we found the CSC/N vs. time plot to be particularly useful (see Figure 2C). For example, Figure 2C shows the dependence of the dynamics of CSC/N with respect to the value of W 4/1 . In this illustrative exercise, the rest of the W i/j values remain constant, while W 4/1 was progressively increased within the range of 5.0 to 5.7 units. Only specific value of  Table 1 presents results from a series of simulation experiments where W i/j values were varied around those that produced the solution previously discussed (Exp 0 in Table 1).
Some solutions, although satisfy d[CSC/N]/dt = 0, differ importantly in terms of their resulting cellular fractions at the steady state. For example, the solution of Exp. 7 (Table 1) is conducive to a steady state in which the cellular fractions CSC/N, P/N, and D/N are 0.0459, 0.1703, and 0.7836, respectively. In this particular case, the modification consisted of increasing the value of W 5/4 (or k 5 /k 4 ) from 0.80 to 0.85. This implies an increase of only 6.25% in the value of the intrinsic reaction rate constant of differentiation versus self-renewal of the subpopulation of progenitor cells. Not intuitively, the fraction of cancer stem cells  Table 1). We found several solutions (see Exp. 15, 16, and 17 in Table 1), with steady state values in the vicinity of P/N<0.18, D/N<0.81, and CSC/N<0.0018, by modifying the value of two of the parameters W i/j with respect to the values of the reference set (Exp. 0 in Table 1). To do so, we selected displacements (DW i/j ) with opposite effects on the steady state P/N or CSC/N values (according to column No. 10 in Table 1). For example, in Exp. 15, an increase of 50% on the value of W 2/1 was compensated by a proportional decrease (50%) on the value W 3/1 . This was an expected result: an increase in the rate of asymmetrical differentiation (r 2 ) has to be balanced by a decrease in the rate of symmetrical differentiation (r 3 ). Similarly, the opposite statement should hold. A decrease of 20% in the intrinsic rate of asymmetrical differentiation was compensated by a 20% increase in the rate of symmetrical differentiation (Exp. 16 in Table 1). Less intuitively, in Exp. 17 (see Figure 3B), an increment in the W 2/1 value was balanced by a decrease on W 4/1 . In this case, a 50% increase in the rate of asymmetrical differentiation is equilibrated by a minor decrease (0.5%) in the intrinsic rate of both symmetrical proliferation (r 4 ) and differentiation of progenitor cells (r 5 ). Note that both the value of k 4 and k 5 are influenced by W 4/1 , since the value of W 5/4 was left unmodified.

Kinetic commonalities during exponential growth
While multiple sets of W i/j could produce results consistent with the proposed set of constraints, those sets must comply with some general characteristics. For example, as illustrated before, the parameters W 2/1 and W 3/1 (and consequently k 2 and k 3 ) are inversely related. To keep the overall fraction of cancer stem cells CSC/N constant over time, an increase on the value of k 2 must be balanced by a proportional decrease on k 3 , and vice versa.
This observation suggests a fine feedback biochemical control, and not necessarily a fixed ratio k 2 /k 3 . This result is relevant, since the experimental determination of the relative probability of occurrence of symmetric and asymmetric CSC division is difficult. The suggestion that regulation of the ratio between symmetric and asymmetric division may be crucial for the development and progression of cancer appears recurrently in the literature [33,46,52,54,55]. For example, Boman et al. [31], using a compartmental model, concluded that the only mechanism that can explain how CSC subpopulations can increase exponentially during colorectal cancer development involves an increase in symmetric SC cell division. This finding suggests that systemic therapies for effective treatment of cancers must act to control or eliminate symmetrical cancer SC division in tumors, while minimally affecting normal SC division in non-tumor tissues. In this respect, Turner et al. [45] have consistently concluded that symmetric division rates are the key in dictating brain tumor composition. Their results also suggested the importance of developing novel treatment strategies that specifically target the CSC subpopulation in brain tumors. Other kinetic relationships appear to be more stringent. Our results suggest that the ratio k 4 /k 1 has to be maintained in a very narrow band in order to sustain exponential tumor growth with constant CSC/N, P/N, and D/N fractions. For example, in our model, for the case where CSC/N<0.0018, P/N<0.18, and D/ N<0.81, the ratio k 4 /k 1 must be maintained at <5. 35. Variations of less than 2.5% around this value must be compensated, for example, by important modifications to the value of k 2 (or more generally stated W 2/1 ; see Figure 4A) or alternatively, k 3 (or more generally stated W 3/1 ; Figure 4B).
Indeed, we found a perfect linear correlation between the values of W 4/1 and W 2/1 that fulfill the conditions d[CSC/N]/dt = 0, CSC/N,0.0018, P/N<0.20, and D/N<0.80, namely W 2/1 = 219.997W 4/1 +107.98 ( Figure 4A). An analogous linear relationship exists between W 4/1 and W 3/1 , namely W 3/1 = 20.1947W 4/1 + 1.0515 ( Figure 4B). This suggests that the relative magnitude of the proliferation and differentiation rate constants for progenitor cells (k 4 +k 5 ) must be at least half an order of magnitude above the analogous parameters for stem cells (CSC) in order to maintain exponential tumor growth while keeping (d[CSC/N]/dt = 0). Our results also suggest that the ratio between the intrinsic constant rate for symmetric proliferation of progenitor cells and the analogous parameter for stem cells symmetric self-renewal (W 4/1 = k 4 /k 1 ) should be approximately half an order of magnitude (between 5.10 and 5.4), independently of the values of k 2 and k 3 . We also observe that, to maintain tumor growth, the constant for symmetric cancer stem cell renewal should be in the same order of magnitude as the sum of k 2 and k 3 .
The time evolution of the rates of proliferation of CSCs, progenitors, and terminally differentiated cells (dCSC/dt, dP/dt, and dD/dt respectively) can be calculated using equations (1), (2), and (3). Despite the fact that the k j values associated with CSC events are of the order of magnitude, and in some cases even higher, than k j values corresponding to P and D cells, the number  of stem-like cells is much smaller than numbers of progenitors or differentiated cells. Consequently, the product (k 1 CSC) and the global rate of proliferation (dCSC/dt) are one or two orders of magnitude smaller than those rates for more differentiated cells during all stages of the exponential growth of a tumor. In addition, early in the development of a tumor, the ratios between the global growth rates of the different cell subtypes achieve equilibrium. According to our simulation results, (dP/dt)/(dCSC/dt)<100 and (dD/dt)/(dP/dt)<4. Therefore, (dD/dt)/(dCSC/dt)<400. This equilibrium of growth rates between cell subtypes must be achieved if balanced tumor growth is to maintain constant ratios of cell subpopulations. Again, these ratios between growth rates of different cell subpopulations are not easy to calculate in vivo.

The role of cell death
The model is sensitive to the ratios between cancer stem cell self-renewal and cell death (W 6/1 , W 7/1 , W 8/1 ). The initial assumption of k 6 ,k 7 ,k 8 should be satisfied in order to observe constant cell fractions throughout exponential growth. This is consistent with several experimental reports. For instance, it has been observed that the apparent rate of death in prostate CSC is much lower than in other tumor cell subtypes [20,51]. We also found that the W 6/1 , W 7/1 ,

Model fitness to experimental data sets
The model is flexible enough to allow proper adjustment to a wide range of possible tumor growth scenarios. Here, to validate the model, we selected three different tumor growth experimental data sets available from previous literature [13,25,69]. In all cases, the same W i/j vector was used to describe the ratios between the kinetic parameters of the model, namely W i/j = [W 2/1 , W 3/1 , W 4/1 , W 5/4 , W 6/1 , W 7/1 , W 8/1 ] = [1.0, 0.01, 5.35, 0.8, 0.01, 0.1, 1.0]. After proper scaling of the y axes (by multiplying by a different scaling factor in each case), the model reproduces the three experimental data sets with an R-square greater than 0.97 ( Figure 6). This suggests that, although the values for the intrinsic kinetic rates of each tumor might be different, the relationships between them (all W i/j values) could be approximately common among different cancer types. For all simulations, all constraints were satisfied for each set of experimental data, particularly the condition of a constant CSC cell fraction lower than 1% (specifically 0.0018).
We should emphasize that our model is capable of reproducing the evolution of tumor growth only during its exponential phase. For the experimental sets that we had analyzed, this means up to a size of 1500 mm 3 ; clinically, a medium size solid tumor. Tumors of this size quite often contain necrotic tissue at their central core. Our model does not distinguish between living and dead tissue, and only provides an overall volume based on a very naive and simple approximation of equal spherical volume contribution of each cell (dead or alive) within the tumor. Although simple and unrealistic, this assumption has yielded a good agreement with experimental sets for our descriptive purposes.

Contrasting strategies to combat cancer
Most of the currently used chemotherapy and radiotherapy strategies against cancer are unable to distinguish between different tumor cell types, or even healthy and tumor cells, killing them all unselectively. Some reports also indicate that CSCs are particularly resistant to conventional therapeutic procedures [18,19,23,26,51]. Friel et al. [19] reported that CSCs isolated from human EnCa were particularly resistant to Paclitaxel, a widely used chemotherapeutic anticancer agent. Kang et al. [18] found that CSCs from GBM were radio-resistant when exposed to radiation dosages that killed other tumor cell subpopulations. Eyler and Rich [26] and Dylla et al. [23] reported additional evidence of CSC resistance to conventional therapies. More recently, Tan et al. demonstrated that holoclone forming cells from pancreatic tumors (with stem cell characteristics) exhibit much higher chemoresistance to gemcitabine and 5-FU than meroclones and paraclones [51].
This enhanced resistance of CSCs would explain the observed aggressive regeneration of tumors after treatment. The remaining CSC population would grow exponentially (having an abundant amount of nutrients) and would readily regenerate the other subpopulations that form the bulk of the tumor. Figure 7A shows the behavior of a solid tumor when treated with a non-specific therapy such that 90%, 99%, 99.9% or 99.99% of the solid tumor cells are eradicated at a particular time, let us say 120 days. For illustrative purposes, the base case we have chosen is the case illustrated in Figure 6C [13]. In all instances, after an incubation time, the tumor relapses to exhibit practically the same progression profile originally followed. Even in the case of the most effective treatment, which kills 99.99% of all tumor An intuitive therapy to eradicate tumors, widely suggested in recent literature [22,23,24,[26][27][28][29]31,46,[70][71][72], relates to the targeting of CSCs, killing them selectively through treatment (i.e., augmenting k 6 significantly over k 7 , and k 8 ) or, alternatively, preferentially killing stem cells at a particular time. Other authors have evaluated this strategy using mathematical models; for example, Dingli and Michor [36] designed a simple mathematical model to demonstrate the importance of eliminating tumor stem cells. The authors explored different therapeutic scenarios to illustrate the properties required in novel anti-cancer agents for successful tumor treatment. Their results indicated that successful therapy must eradicate tumor stem cells. Similarly, Ganguly and Puri [39] formulated a mathematical model to evaluate chemotherapeutic drug efficacy in arresting tumor growth based on the cancer stem cell hypothesis. Their results also suggest that the best response to chemotherapy occurs when a drug targets abnormal stem cells.
However, our simulations show that this therapeutic avenue, although superior to unselective treatments, is only effective if progenitor cells (P) are also targeted. This is, for the time window that our model is applicable, the exponential tumor growth phase, both cancer stem cells and progenitor cells have to be eradicated to stop tumor growth. This can be easily seen by examining equation 2 [dP/dt = (k 2 +2k 3 )CSC+(k 4 2k 5 2k 7 )P ]. In that equation, the rate of accumulation of the subpopulation of progenitor cells depends on two terms, one affected by the number of CSCs and the other dependent on the number of P cells. Therefore, even in the complete absence of CSCs, the term depending on P has to be negative for P to decrease. That term, namely (k 4 2k 5 2k 7 )P, can only be negative if k 4 ,k 5 +k 7 . For the set of W j/i values that we have chosen (based on the rational that we had explained before: d[CSC/N]/dt = 0 with CSC/N,0.01) this condition is not met. was multiplied by a different scaling factor (17.56W i/j ; 5.06W i/j ; and 7.06W i/j respectively). (A) A human prostate tumor transplanted into a mouse model [69]; it is assumed that the original percentage of CSCs was 1%; (B) 2000 CSCs from a breast primary tumor implanted in NOD/SCID mice [25]; and (C) 1610 5 colon CSCs isolated and implanted in NOD/SCID mice [13]. For all simulations, the vector W i/j = W i/j = [W 2/1 ,  The need to not only specifically target CSC, but also the progenitor subpopulation to effectively detain tumor progression has been suggested before in literature [73], also at the light of a mathematical argument.
Conceivably, in the long term, selective killing of all cancer stem cells will be sufficient to eradicate the tumor, since CSCs are the ultimate cell reservoir responsible for sustained growth. At the end of exponential growth phase, the relative weight of k 7 has to increase (or alternatively the ratio W 5/4 = k 5 /k 4 will change) due to causes such as oxygen mass transfer or nutrient limitations. We plan to include this functionality in a later version of our model. This last comment on the relative ratio W 5/4 = k 5 /k 4 is useful to introduce a less intuitive therapeutic strategy. Promotion of differentiation has been suggested for clinical practice [74,75]. Using the proposed model, we explored variants of this approach. Indeed, augmenting W 5/4 = k 5 /k 4 appears to be therapeutically promising. In biological terms, this means increasing the relative intrinsic rate of progenitor cell differentiation with respect to the rate of progenitor cell renewal. Interestingly, a very modest 10% increase in W 5/4 (from 0.8 to 0.9), while keeping the rest of the W i/j values fixed, retards tumor relapse more effectively than does the unselective tumor treatment (see Figure 7B). Similarly, increasing the ratio W 4/1 = k 4 /k 1 and promoting a tumor richer in P cells also is another effective therapeutic strategy. Indeed, the delay induced by an increase of 15% on W 4/1 = k 4 /k 1 is comparable to that caused by unselective eradication of 99.99% of the tumor cell mass ( Figure 7B).

Concluding remarks
In summary, in this study, we presented a first version of a conceptually simple model, with an analytical solution, that is capable of describing the basic kinetic features of tumor growth during the exponential phase and that is consistent with the Cancer Stem Cell hypothesis. Three cell subpopulations were considered: CSCs, progenitors (P), and terminally differentiated (D) cells. Each event related to cell division or death of each one of these subpopulations has been represented and modeled as a chemical reaction. This resulted in an analytically solvable set of ordinary differential equations that describes the time evolution of each cell subpopulation, as well as the overall tumor volume evolution during exponential growth.
Although, in principle, an infinite set of combinations of model parameter can be studied, we found that only a limited set of model solutions is feasible if some biologically sound constraints are imposed. For example, if we accept that the fraction of cancer stem cells during the exponential phase of the tumor is practically constant and no greater than 0.01 (namely d[CNC/N]/dt = 0 with CSC/N,0.01), then only a reduced set of solutions is feasible. The analysis of those sets suggests kinetic commonalities in solid tumor growth: (a) the rates of symmetrical and asymmetrical CSC renewal must be in the same order of magnitude; (b) the intrinsic rate of renewal and differentiation of progenitor cells should be half an order of magnitude higher than the corresponding intrinsic rates for cancer stem cells; (c) the rates of apoptosis of the CSC, P, and D cells are progressively higher by approximately one order of magnitude. The flexibility of the model was tested by fitting experimental data sets from three different tumor growth scenarios. After adequate scaling, a single set of kinetic parameters can be used for adequate reproduction of different tumor growth cases.
We do not claim that our model renders accurate information about the kinetics of solid tumor formation, but we observe that it does provide insight into several underlying kinetic behaviors of solid tumor growth that would be difficult to directly study experimentally. As illustrated in this work, the model can even be useful for anticipating the effect of different therapeutic strategies (available or potential) against cancer.