Effect of feedback regulation on stem cell fractions in tissues and tumors: understanding chemo-resistance in bladder cancer

While resistance mutations are often implicated in the failure of cancer therapy, lack of response also occurs without such mutants. In bladder cancer mouse xenografts, repeated chemotherapy cycles have resulted in cancer stem cell (CSC) enrichment, and consequent loss of therapy response due to the reduced susceptibility of CSCs to drugs. A particular feedback loop present in the xenografts has been shown to promote CSC enrichment in this system. Yet, many other regulatory loops might also be operational and might promote CSC enrichment. Their identification is central to improving therapy response. Here, we perform a comprehensive mathematical analysis to define what types of regulatory feedback loops can and cannot contribute to CSC enrichment, providing guidance to the experimental identification of feedback molecules. We derive a formula that reveals whether or not the cell population experiences CSC enrichment over time, based on the properties of the feedback. We find that negative feedback on the CSC division rate or positive feedback on differentiated cell death rate can lead to CSC enrichment. Further, the feedback mediators that achieve CSC enrichment can be secreted by either CSCs or by more differentiated cells. The extent of enrichment is determined by the CSC death rate, the CSC self-renewal probability, and by feedback strength. Defining these general characteristics of feedback loops can guide the experimental screening for and identification of feedback mediators that can promote CSC enrichment in bladder cancer and potentially other tumors. This can help understand and overcome the phenomenon of CSC-based therapy resistance.


Introduction
Tumor stem cells are thought to be important for the initiation and maintenance of cancers. In addition, it is becoming clear that tumor stem cells can also contribute to a reduced response to treatments, especially if the fraction of the stem cells in the tumor is high [1]. Cancer stem cells (CSCs) are intrinsically less responsive to drug treatment [2][3][4] due to protective properties they share with normal stem cells, including higher expression of drug-efflux pumps [1], better DNA-repair capacity [5], and enhanced protection against reactive oxygen species [6]. High CSC fractions can therefore give rise to poor a response to therapy even in the absence of resistance-inducing mutations.
CSC-based resistance falls in the more general category of non-genetic drug resistance. There is strong evidence that tumors contain non-genetic, heritable variation, and that this contributes to clonal evolutionary processes [7]. The importance of the interplay between genetic and non-genetic heterogeneity for clonal evolutionary processes has been demonstrated with mathematical models, especially in the context of therapeutic interventions [8,9].
The importance of stem cell-based resistance in vivo has been established in patient-derived bladder cancer mouse xenografts [10]. In clinically realistic chemotherapy regimes, it has been demonstrated that stem cell fractions increased during successive treatment cycles, and that this increase in the stem cell fractions correlated with a reduced response to chemotherapy during the next therapy cycle. Experiments identified a positive feedback loop that operates within the tumor as an important component of stem cell enrichment during the treatment phase [10]. Chemotherapy-induced death of more differentiated cells resulted in a wound healing type response that is mediated by PGE2 and that resulted in the activation and proliferation of CSCs during treatment. A mathematical model of these treatment dynamics confirmed that this wound-healing response could contribute to the amplification of the CSC population during chemotherapy [11]. At the same time, however, the model showed that this CSC enrichment was not maintained following cessation of therapy. Instead, CSC fractions were predicted to return to pre-treatment levels after treatment cessation. In other words, the positive feedback wound-healing response alone could not account for a pronounced reduction in the response to chemotherapy during the next cycle following an off-therapy phase. This situation changed when we added negative feedback from differentiated cells onto the division rate of stem cells in the mathematical model [11]. We found that in this case, the enriched stem cell fractions were predicted to be maintained in the long term during the off-treatment phase. The model with the negative feedback loop could thus reproduce the experimentally observed result that CSC enrichment during a round of chemotherapy could reduce the response to a subsequent chemotherapy cycle.
The notion that regulatory feedback loops originating from healthy tissue remain partially active in tumors is not well-established, nor is the notion that they determine the disease course or the response to therapy. The PGE 2 -based wound healing response in bladder cancer xeografts is the best experimental support for these ideas. Negative feedback loops have so far not been implicated in a similar way, but there is clear evidence that negative feedback regulation likely plays and important role in healthy tissue dynamics, e.g., in the olfactory epithelium, where GDF11 and Activin βB negatively regulate self-renewal rates in progenitor and stem cells [12,13]. Some of this negative feedback regulation could thus very well continue to be operational to a certain extent in tumors, and indeed the analysis of tumor growth data has found patterns that are difficult to account for in the absence of negative feedback regulation [14,15].
While our previous modeling showed that one particular negative feedback loop could lead to sustained CSC enrichment following chemotherapy of bladder cancer xenografts, many other regulatory feedback loops can potentially be present within the tumor. Some of these might also be able to contribute to sustained CSC enrichment and to a loss of therapy response, while others might not contribute to this phenomenon. A targeted experimental search for feedback factors that promote loss of treatment responses in bladder cancer requires us to know which types of feedback mechanisms intrinsic to the cell population can in principle promote sustained CSC enrichment. The aim of this paper is to identify such candidate feedback loops, based on the analysis of mathematical models, thus providing a theoretical basis for the experimental identification of relevant, specific feedback molecules.
We focus on two basic scenarios. First, we assume that partial breakage of feedback regulation results in temporary cell growth towards a new equilibrium, characterized by an overall larger number of cells. This could correspond to a single step in step-wise tumor progression. We investigate the conditions required for the stem cell fraction to be larger at the new compared to the old equilibrium. In particular we study how remaining feedback loops determine the stem cell fraction. Second, we consider unbounded tumor growth and investigate how different feedback mechanisms that remain in a growing tumor cell population can determine whether or not CSC enrichment occurs during growth. Much of this work is done using ordinary differential equations. In the context of unbounded tumor growth, the effect of spatial growth patterns on stem cell enrichment is explored using a stochastic agent-based model and an analytical approximation.
The present study provides a solid theoretical basis for implicating the presence of feedback regulatory loops as a determinant of responses to cancer therapy. This adds to the mathematical literature quantifying the role of feedback regulation for tissue and tumor dynamics [11,13,[15][16][17][18][19][20][21][22][23][24], and builds upon the wider mathematical literature concerned with the dynamics of hierarchically structured cell populations, e.g. [25][26][27][28][29][30][31] and stem cell fractions [32]. Two major ap-proaches can be mentioned as most relevant in the present context. One approach uses spatial, agent-based or hybrid models to investigate SC dynamics, including questions of SC enrichment. Reference [32] provides an excellent review of the relevant literature. In particular, a versatile model of CSC and non-SCs has been developed [33,34], where the cellular expansion and competition dynamics could be explored. It was found that an increased proliferation capacity of non-stem cells results in encapsulation of SCs and tumor dormancy, while an increase in migration may lead to ``liberation" and expansion of SCs. In [34] it was demonstrated that the CSC fraction of a tumor population can vary by multiple orders of magnitude as a function of the generational life span of the nonstem cancer cells. In [27], the time-dependent SC fraction was studied in the context of a similar model, and it was further found that spontaneous cell death yields a higher SC fraction; Reference [35] studied dynamics and localization of SCs by using a hybrid model.
The second approach is non-spatial, where the dynamics of SCs and their differentiation are represented by using ODEs with or without feedback. For example, in [24], a systematic analysis of feedback is performed by using delay differential equations. Reference [8] considers a two-compartment, multi-clone model with a specific type of negative feedback on SC self-renewal, and shows that clonal evolution selects for fast reproducing and highly self-renewing cells at primary diagnosis, while relapse following therapy-induced remission is associated with highly self-renewing but slowly proliferating cells. A multicompartment model of differentiation was introduced axiomatically in [25,26,36] and the steady states were studied, both for general functional forms and for a particular case of regulatory feedback functions. In [37], a similar model was used to compare cellular properties of leukemic stem cells to those of their benign counterparts, by deriving conditions for expansion of malignant cell clones. Reference [23] investigated the impact of feedback loops and their breakage on cancer progression.
While questions of SC enrichment were not directly addressed in the above studies, reference [38] presented a hierarchical model of cell growth and differentiation, and showed that in exponential growth (and in the absence of any control loops) the fraction of SCs reaches a steady state level, which is a decreasing function of the number of compartments and an increasing function of SC proliferation rate.
While a wealth of results of SC fraction dynamics have been obtained in the spatial (agent based and hybrid) models, analytical understanding or generalizations are more difficult in the context of these models; also feedback factors were not usually considered explicitly. On the other hand, a high level of analytical understanding has been reached in ODE models, but SC enrichment dynamics has not been specifically studied, except in the simpler cases in the absence of feedback loops. In the present study, we develop an axiomatic model of SC dynamics where feedback loops are included explicitly, and the functional form of the controls (which remains unknown) is kept as general as possible. We investigate conditions that lead to SC enrichment under different types of growth dy-namics, and define what types of regulatory feedback loops can and cannot contribute to CSC enrichment, providing guidance to further experimental inquiries into specific signaling mechanisms.

The basic mathematical modeling approach
An ordinary differential equation model has been used to describe tissue hierarchy dynamics in a healthy tissue [13,16], and the models presented here build on these approaches. While cell lineages consist of stem cells, transit amplifying cells, and terminally differentiated cells, our models make a simplification and take into account only stem cells (which encompass all the proliferating cells) and It has been shown that introduction of negative feedback loops can result in more realistic behavior, where a stable equilibrium is attained for P>0.5 [13].
This was shown in the context of two specific feedback loops, and subsequently generalized to comprehensively list all possible (positive and negative) feedback loops compatible with stability [19,39,40]. Here, we also use a general model to assume different kinds of feedback on the rate of cell division, L, the rate of cell death, D, and the probability of self-renewal, P. We also add the possibility that stem cells die with a rate δ (which can also be subject to feedback). In the context of our model, feedback is equivalent to a dependence of rates and probabili-ties on the population sizes, x and/or y. Hence, the model is given by the following ODEs: The division rate L, the death rates, D and δ, and the probability of self-renewal, P, are now functions of either the number of stem cells, x, or the number of differentiated cells, y, or both.
Evolution can result in the generation of mutant cell populations that are characterized by a higher self renewal probability, given by P 2 , which can be viewed as an early step towards carcinogenesis. Hence, we now have two stem and differentiated cell populations denoted by subscripts 1 and 2 for wild type and mutant types, respectively. The equations are thus given by where x=x 1 +x 2 and y=y 1 +y 2 . The two cell populations are in competition with each other, mediated by the feedback factors that are shared between the two populations. Table 1 summarizes all the variables used in this paper (both in this section and in the later sections). These equations are structurally similar to previously published models that were explored in different contexts [37].
Note that we only consider mutants with a higher self-renewal probability, because these are the only mutant types that can grow from low numbers and invade in this type of model. Mutants in other parameters fail to grow and will die out, as described in a different context in reference [21]. In principle, it is possible that a mutation causes a simultaneous change in more than one parameter, which, however, does not lead to a change in conclusion and is not pursued further.

Cell growth towards a new equilibrium.
The first important scenario happens when the mutant cell population gains a selective advantage, outcompetes the original, healthy cell population, and grows towards a new and higher equilibrium level. This is achieved by assuming that P 2 >P, along with other conditions on the rate functions that are specified in Section 1 of SI We examine the conditions under which the stem cell fraction at the new equilibrium is increased compared to that at the original equilibrium. In terms of the model notation, we are interested in the quantity v(x,y) = x / y, i.e. the ratio of stem to differentiated cells. We investigate how different combinations of feedback loops that remain in the mutant cell population impact stem cell enrichment.
In the following we assume that the division and death rates in the above models are monotonic functions of the number of stem and/or differentiated cells. We further assume that in the absence of mutants, the system is at equilibrium, characterized by the pair (x , y) , which satisfies:

At this equilibrium, the fraction of stem to differentiated cells is given by
We assume that the mutant cell population can invade from low numbers and displace the original cell population. This occurs if P 2 (x,y) > P(x,y) (this is a sufficient condition). Assuming that a new equilibrium is reached, it is characterized by and the new fraction of stem cells is given by We examine under what conditions ν * >ν , i.e. when the stem cell ratio at the newly obtained mutant equilibrium exceeds that at the original equilibrium.
We observe two qualitatively distinct outcomes, see We will illustrate these points by using some specific examples. The first example is of SC enrichment. Consider a system where δ=0, the self renewal probabilities P and P 2 are given by decreasing functions of y (see solid lines in Fig 1(a)), and the division rate L is also a function of y, solid line in Fig 1(b). The cell dynamics and control loops corresponding to this system are schematically shown in Fig   2(a). In this and other such diagrams, blue arrows correspond to cellular processes (characterized by kinetic rates, such as L and D), and cell fate decisions, which are probabilities (P or P 2 ). The red negative arrows originating in the DC circle represent a negative dependence of both functions L and P on y (the number of DCs). In this first example, SC enrichment is predicted to occur. Fig 2(b) shows the cell dynamics, once a mutant is introduced at 100 time units. While the wild type population goes extinct, the mutants rise to a new equilibrium characterized by a significantly higher ratio x/y compared to the original equilibrium.
The second example is SC depletion. It is given by a system with rate functions  Fig 2(a). First, we consider the death rate of the stem cells, δ. The amount of enrichment is smaller for nonzero stem cell death rates, compared to the case of δ=0. In Fig 3(b) we can see that in the presence of SC death, the increase in the enrichment parameter, x/y, is more modest than that of Fig 2(b).
The second modification is a smaller self-renewal probability, P 2 , of the mutant cell population, Fig. 3(c). The self-renewal probability P 2 of the mutant cells is that given by the dashed line in Fig 1(a). Again, this results in a more modest increase in the stem cell fraction compared to Fig 2(b).
The third modification is a less pronounced feedback on the division rate. The result is that the SC enrichment becomes smaller. In Fig 3(

Stem cell enrichment in non-equilibrium situations
Next, we study the scenarios where the mutant population grows from low numbers and does not reach a new equilibrium, but instead, continues to grow indefinitely. This would correspond to tumor growth. The definition of enrichment in this context and the relevant methodology will be somewhat different in this case. We will consider the mutant population alone, and study the growth of x 2 and y 2 , in order to find the dynamics of the quantity ν=x 2 /y 2 . We will say that stem cell enrichment occurs if the quantity ν increases during population growth, either infinitely, or temporarily. For simplicity we will assume that the cancer stem cell (CSC) population does not die (δ=0). We further assume that the probability of self-renewal of mutants, P 2 , is a monotonic function of the population size (stem cells or differentiated cells) that satisfies P 2 > 1/2, and that in the limit of large populations, it approaches a limiting value, P > 1/2. When we consider the mutant dynamics, we will drop the subscript 2, and simply study the equations: If we assume that both D and L depend on a single population, that is, L=L(x), D=D(x), or L=L(y),D=D(y), the following approximations can be derived (see Section 4 of SI). As the cell population expands, the DC population behaves as and ratio of stem to differentiated cells, ν, in the limit of large times is given by: Equation (4)  ). In this case, a temporary reduction in the ratio ν can be observed during tumor growth, until ν converges to a constant. In other words, we observe a reduction in the CSC fraction over time, until it reaches a limiting value. Note that the SC fraction can never decrease to zero, the minimum fraction is given by To summarize, the following three types of SC fraction dynamics are predicted: (1) If L is subject to negative feedback and decays to zero and/or D is subject to positive feedback and increases without bound, then we have unlimited SC enrichment, such that the content of SCs grows to 100% in the log run. However, if the feedback on D is (effectively) greater such that L/D is decreasing and bounded above 0, then this would correspond to type (2) and result in limited SC enrichment, where SC fraction grows and saturates below 100%.

SC enrichment in spatially structured populations
According to the above results, stem cell enrichment during growth requires certain feedback mechanism to be present in the tumor cell population, such that the ratio L/D is reduced as the population size is increased. Typically such feedback can occur through signaling factors that are secreted from stem or differentiated cells. Another way to achieve a similar result can be spatially restricted reproduction of cells. In such scenarios, cells experience range expansion in two dimensions, or grow as expanding sphere-like structures in 3D. Inside the expanding population, divisions must be balanced with deaths because free space is limited.
In a way this works similarly to control loops affecting division and/or death rates, which were discussed earlier in the paper; in the case of spatially restricted growth, control is essentially competition for space, which leads to slower divisions/ higher death as the density increases.
To explore this, we first considered a two-dimensional stochastic agent- is randomly picked as a target for one of the daughter cells. If the chosen spot is already filled, the division event is aborted, otherwise it proceeds. If division proceeds, both daughter cells will be stem cells with a probability P 0 (self-renewal).
With probability 1-P 0 , both daughter cells will be differentiated cells. If the sampled spot contains a differentiated cell, death occurs with a probability D 0 . No explicit feedback processes were included in the model. The simulation was started with 9 stem cells (and no differentiated cells). Assuming that stem cells do not die, the resulting average growth curve is shown in Figure 7(a). While initially, the differentiated cells grow to be more abundant than the stem cells, the stem cell population enriches over time and eventually becomes dominant as the cell population grows. Figure 7(b) shows the same kind of simulation, but assuming that stem cells die with a rate that is smaller than the death rate of differentiated cells.
Consistent with the results obtained for explicit feedback mechanisms, we find that the degree of stem cell enrichment is reduced in the presence of stem cell death (higher rates of stem cell death lead to less enrichment).
A simple mean-field model that takes account of the space limitations inside an expanding population can explain these results (see Section 5 of SI). Indeed, the system in the interim of the expanding globe reaches a dynamic equilibrium state where the density of SCs and DCs is dictated by the balance of division and death rates. Using the mean-field model, we have also calculated the degree of SC enrichment. The higher the SC death rate, the lower the overall density and the higher the fraction of DCs in the population. As time increases, the fraction of SCs reaches the value given by In particular, in the absence of SC death (δ 0 =0), the equilibrium value of v is infinity, which means that SCs exclude DCs everywhere in the core of the expanding range, such that DCs only concentrate in the exterior part of the colony. This corresponds to unlimited SC enrichment.
Agent-based simulations were repeated assuming a 3-dimensional space, where offspring cells could be placed into one of the 27 nearest neighboring spots (Section 5 of SI). Because there are more neighbors in three dimensions, the system is more mixed, and the extent of stem cell enrichment is generally lower. (As mentioned in the beginning of this paper, no stem cell enrichment occurs in a perfectly mixed system in the absence of explicit feedback mechanisms).

Discussion and Conclusions
The experimental observation that a wound-healing type mechanism in bladder cancer can modulate the responsiveness of the tumor to chemotherapy [10] is to our knowledge the first clear case where a specific feedback loop has been implicated in determining therapy outcome. At the same time, however, this interesting observation has brought up a number of questions that can best be answered by using mathematical models to interpret the data. Mathematical modeling [11] suggested that the wound-healing response alone cannot result in CSC enrichment that is sustained beyond the treatment phase, and hence cannot ex- It is further important to note that the models discussed here are based on a set of core assumptions that we consider critical when addressing the particular questions that were the focus of this investigation. Further complexities can be included in future work when studying more advanced questions. An important example is additional tumor heterogeneity that is certainly present in bladder caners as well as in most other tumors. Variation is likely to occur in division rates, death rates, and the ability to secrete or respond to feedback signals; it will be important to extend the models with this in mind. Furthermore, while we con-     At t=100, mutant stem cells are introduced at a low level, resulting in the extinction of wild type cells, and convergence to a new equilibrium, (0,0,x * ,y * ). In the insets, the SC fraction, x/y, is depicted as a function of time. (a,b) Negative con-trol by DCs, resulting in SC enrichment: the functions L, P, and P 2 are given by the solid lines in Fig 1(a,b). (c,d) Positive control by SCs, resulting in SC depletion: the functions L, P, and P 2 are specified in Fig 1(c,d).   Fig 1(a). (d) Shallower L: same parameters as in Fig 2(b), except the division rate L is given by the dashed line in Fig 1(b). We observe that the SC enrichment is (b-d) is less pronounced compared to Fig 2(b). Positive control by SCs; (d) Positive control by SCs. Unlike in Fig 2(a,c), the mu-tant self-renewal probability is assumed to be constant (and thus the control loop is completely severed in mutants), leading to unlimited growth. The 2 top panels correspond to Fig 5(a,b); the 2 bottom panels correspond to Fig 6(a,b). Notations are as in Fig 2(b,d). The purple line in each panel plots the approximation for y (formula (3)), and in the insets the approximation for x/y (formula (4)).
(a) Negative control by DCs: parameters are as in Fig 2( Negative control by SCs: similar to (a), except the rate functions depend on SCs: , and P 2 = P =0.6. In both panes, unlimited SC  Notations are as in Fig 2(b,d). The purple line in each panel plots the approximation for y (formula (3)), and in the insets the approximation for x/y (formula (4)).