Modeling the stem cell hypothesis: Investigating the effects of cancer stem cells and TGF−β on tumor growth

1.   Introduction

It's been almost 20 years since the renewed hypothesis and subsequent confirmation of cancer stem cells has fundamentally changed the way we view the development and treatment of tumors ^[1,2]. Soon after this discovery, cancer stem cells were confirmed in many tumor models (including both solid tumors and leukemias), cultured cell lines, and primary human tumors ^{[3,4,5,6,7,8,9,10,11]}. Cancer stem cell theory describes tumors as comprised of a heterogeneous population of cells, in which there is a subpopulation of tumor cells that behave similarly to the haematopoetic stem cells in normal tissues (reviewed in ^[12,13]). Tumor cells with this "stemness" property have low proliferation rates, but when they do divide they can undergo symmetric (creation of two identical daughter stem cells) or asymmetric cell division ^[14,15]. Asymmetric division creates one daughter cell that retains its "stemness" to become any type of cell and seed new tumors, while the other daughter cell becomes slightly more differentiated (yet still capable of great plasticity) to regenerate the tumor tissue ^[13,16]. Whether these cells arise from mutations of haematopoetic stem cells within adult tissue, de-differentiated tumor cells due to environmental pressures ^[16,17,18], former immune cells ^[19], or another mechanism, it is clear that-cells with cancer stem cell phenotypes contribute to drug and radiation resistance of tumors ^[20,21] as well as increased metastatic potential ^[22,23]. Recent approaches indicate that targeting cancer stem cells as part of a combined therapeutic regime may increase patient survival outcomes ^[24,25,26].

Recent advances have also elucidated the critical importance of transforming growth factor $\beta$ (TGF-$\beta$) in the progression of cancer, as both a tumor repressor and later as a driver of metastasis and immune evasion (reviewed in ^[27,28]). The TGF-$\beta$ pathway is typically responsible for signaling cellular senescence and apoptosis, making it a key anti-tumor agent and powerful immune suppressor. Paradoxically, TGF-$\beta$ also promotes angiogenesis and metastasis, thus promoting tumor invasiveness. TGF-$\beta$ is often upregulated within the tumor microenvironment, and can impact not only survival of tumor cells but suppression of the immune response to tumor cells. For example, TGF-$\beta$ secreted by tumor cells specifically prevents CD8$^+$ t-cell release of perforin and granzymes ^[29], TGF-$\beta$ secreted by fibroblasts block CD8$^+$ t-cells from entering the tumor microenvironment ^[30], and TGF-$\beta$ signaling may be a mechanism that keeps cancer stem cells dormant and resistant to immune recognition ^[31].

Targeting resident CD8$^+$ T-cells, as well as repurposing anti-viral CD8$^+$ T-cells may be important in developing new cancer therapies ^[32,33]. However, the immune suppression of CD8$^+$ T-cells by TGF-$\beta$, lymphocyte exhaustion, and the repression of antigen-specific targeting of both cancer cells and cancer stem cells must be overcome for these therapies to be effective ^{[34,35,36,37]}. To explore this, we have created a mathematical model that encompasses both the contributions of cancer stem cells and TGF-$\beta$ in tumor growth and tumor progression, adding to the work done in ^[38,39]. We calibrate the model parameters using experimental data from ^[40], which present data for chemically induced skin papillomas, a type of HPV in mice where one randomized group of mice treated with 5-Fluorouracil (5-FU) and the control group were not treated.

Mathematical models have sought to help explain observed behaviors of tumor cell populations. Early mathematical models to include cancer stem cells were published in 2006 ^[41,42,43], were modified as more information was gathered ^{[44,45,46,47]} and remain a key component of tumor modeling today ^{[38,39,48,49,50]}. Similarly, mathematical modeling of the role of TGF-$\beta$ in tumor progression, which began in 1991 ^[51], continues to be modified as new information emerges ^[52,53].

2.   The mathematical model

In order to understand the combined contributions of cancer stem cells and TGF-$\beta$ on the progression of the cancer, we propose a coupled model that incorporates cancer stem cells ($S$), effector CD$8^+$ T-cells ($E$), regulatory T-cells ($R$), tumor cells ($T$), and transforming growth factor TGF-$\beta$ ($B$), which is produced by both regulatory T-cells (also called Tregs) and tumor cells in the model. The relationships between the five state variables are illustrated in Figure 1. We begin by describing the no-treatment model.

2.1.   The No-treatment model

In order to model the dynamics between immune cells, tumor cells and cancer stem cells, several simplifying assumptions must be made. To effectively fight tumor cells, the immune system must balance the necessity of being non-reactive to self-tissues (or risk auto-immune reactions) with the essential task of recognizing subtle differences in tumor antigens to remove aberrant cancerous cells. The immune cells responsible for this tumor recognition are primarily lymphocytes: natural killer cells (NK), natural killer T-cells (NKT), and CD8$^+$ T-cells, with help from other immune mediators (reviewed in ^[54,55,56]). However, chronic stimulation, such as those found in anti-viral and anti-tumor immune responses, can exhaust lymphocytes and prevent them from functioning ^[56,57].

In this model we consider activated CD8$^+$ T-cells and regulatory T-cells as the primary lymphocytes in the system. This is due to the firmly established role of CD8$^+$ T-cells in cancer therapies; research on NK and NKT-cells is still rapidly evolving ^[58,59]. Regulatory T-cells are the primary immune cells that produce TGF-$\beta$ to suppress inappropriate immune responses ^[60], and are thus essential to our inclusion of TGF-$\beta$ in the model.

As in ^[38], we assume that cancer stem cells have the potential to divide either symmetrically into two cancerous stem-cells, or asymmetrically into either a cancerous stem-cell and a non-stem cancer cell or two non-stem cancer cells, with respective probabilities $a_1, a_2, \text{ and }a_3$. We also assume that the typical population of cancer stem cells begins at the low end of ranges calculated through experimental data: approximately 1 percent of the total cancer cell population ^[14,15]. We denote the natural death rates of each population by $\delta$, while the kill-terms resulting from various population interactions are denoted by $\mu.$ The following system of ordinary differential equations describes the changing populations of each of the five state variables.

Equation 2.1 describes the changes occuring in the stem cell population, $S$. The stem cell population grows at a rate dictated by the quantity $a_1-a_3$, which governs the number of daughter stem-cells produced from both symmetric and asymmetric stem cell division. Interactions between cancer cells and CD8$^+$ T-cells result in cancer stem cell death, described by the $\mu_{ES}ES$ term. Finally, cancer stem cells die naturally at a rate of $\delta_S$-we note that due to stem cell longevity, the natural death rate for stem cells is almost negligible.

Equation 2.2 describes the population change of the effector CD8$^+$ T-cells, $E$. These cells have a basal growth rate of $b$ and are recruited to the tumor site as the tumor population begins to swell, as described in ^[39]. The inhibitory effects of the TGF-$\beta$ and the tumor on recruitment and proliferation are accounted for by the $(1+c_3TB)^{-1}$ term. The effector cells revert to regulatory T-cells at a rate of $r$ per unit time, and the shutdown of the activated effector cells by the regulatory T-cells is expressed by the $\mu_{RE}RE$ term.

The change in TGF-$\beta$ concentration over time is expressed in Equation 2.3. The production of TGF-$\beta$ is dependent on the tumor and Treg populations, as presented in ^[39]. We note that there are other existing models for describing TGF-$\beta$ population. A recent example is by Khatibi et al. in ^[53]. Upregulation of TGF-$\beta$ as a result of the increasing tumor cell population is given by the Michaelis-Menten term, $a\frac{T^2}{c_2+T^2}$. Conversely, the production of TGF-$\beta$ by the Tregs is linear, as described in the second term.

As defined in Equation 2.4, the increase in regulatory T-cells occurs via two methods. First, reversal of the activated effector CD8$^+$ T-cells contributes to an increasing Treg population. In addition, Tregs are activated by TGF-$\beta$; thus, an increasing concentration of TGF-$\beta$ results in an increase in the Treg population as well, described by the $\alpha_{BR}\frac{RB}{c_4+B}$ term.

The change in the tumor population over time is described in Equation 2.5. We assume that in the absence of immune cells, the tumor grows logistically. We note that another natural choice of model for tumor growth is the Gompertz Equation, but we found no significant change in the behavior of the populations when implementing this model and therefore opted to use the simpler logistic growth model. The cancerous stem cell division into non-stem cancer cells also contributes to the tumor population, governed by the $(a_2+2a_3)S$ term. Interactions between the effector CD8$^+$ T-cells and the tumor cells result in tumor apoptosis at a rate inversely proportional to the TGF-$\beta$ production.

Figure 2 shows the behavior of each of the five state variable populations in the no-treatment case. Parameter values used for this simulation can be found in Table 1.

2.2.   Model with chemotherapy

In this section, we incorporate the effects of chemotherapy in the dynamics between the immune system, cancer stem cells and tumor cells. We assume that chemotherapy has a negative impact on CD8$^+$ T-cells, tumor cells and regulatory T-cells, but that extermination rates may differ. We assume that chemotherapy does not affect the cancer stem cells; it is designed to target fast-cycling cells, while cancer stem cells are characterized by their slow division rates ^[14,15]. We let $M$ denote the concentration of chemotherapy drug in the blood. The following system of equations describes the model with this chemotherapy regimen incorporated.

In the updated model, the cancer stem cell and TGF-$\beta$ equations remain unchanged from their no-treatment counterparts. The equations for the effector CD8$^+$ T-cells, regulatory T-cells, and tumor cells each have an additional removal term to allow for the destruction of cells by the chemotherapy. Equation 2.11 describes the concentration of chemotherapy drug in the blood. We assume that drug concentration decays exponentially with time and include an impulse term, $\nu_b$, to describe the discrete dosage administrations.

3.   Model calibration and validation

To ensure model alignment with experimental data, the tumor state equations were fit to data from ^[40], using the Delayed Rejection Adaptive Metropolis (DRAM) algorithm for Bayesian inference ^[64]. This variation of the Metropolis-Hastings algorithm allows for periodic updating of the covariance matrix used in the proposal distribution via the added adaptation step. Additionally, use of this algorithm results in better mixing of posterior chains via the inclusion of a delayed rejection procedure (which allows for selection of an alternate sample to be drawn from a narrower proposal distribution in place of outright rejection).

For the non-treatment model, the parameters $\delta_T$, $\eta$ and $k$ were fit, with parameter samples drawn from the priors $\mathcal{U}(0, 1)$, $\mathcal{U}(0, 2)$ and $\mathcal{U}(10^8, 10^{10})$, respectively. Uniform prior bounds were selected to ensure biologically realistic posterior values. All remaining model parameters were fixed at their literature values, or estimated where literature values were not available. See Table 1 for fixed parameter values. The data used for calibration was from the first observed tumor (left untreated) on mouse CM37 in ^[40]. The final fit of Equation 2.5 to the data is shown in Figure 3.

Upon completion of the calibration for the no-treatment model, a further calibration was performed upon the chemotherapy model, using data from mouse CM41 in ^[40]. This second mouse was treated with 5-Fluorouracil with a dosage of 50 mg/kg every seven days, beginning at day 26 when the tumor first exceeded 3 mm in size. Of the three tumors for CM41, we used the data from the second tumor as this tumor most closely mimicked the growth rate of the tumor chosen for the no-treatment case. With the assumption that proliferation rate, carrying capacity, and tumor death rate should be very similar in this second tumor, only $K_T$-the kill rate for tumor cells due to chemotherapy-was estimated in this second calibration procedure. A uniform prior of $\mathcal{U}(0, 20)$ was used to initialize the DRAM procedure. The chemotherapy decay rate, denoted by $\gamma$, was fixed at $49.9066$ days$^{-1}$ according to the half-life function, $\gamma = \text{ln}(2)/t_{1/2}$, where the half-life of 5-Fluorouracil is $t_{1/2} = 20$ minutes ^[65]. The kill rates for CD8$^+$ and T-regulatory cells, $K_E$ and $K_R$, were chosen to be half of the determined $K_T$ value, as these rates should be smaller than $K_T$ but still on the same order of magnitude ^[63]. The final fit of Equation 2.10 to the data from ^[40] is shown in Figure 3.

4.   Treatment simulations

Figure 4 shows the dynamics of the six populations over the chemotherapy administration period outlined in ^[40], which is administering 50 mg/kg of 5-Fluorouracil (5-FU) per week to mice with skin carcinomas, beginning when the tumor reached 3-4 mm in size. It can be observed that while this treatment regimen reduces the rate at which the tumor develops, chemotherapy at this dosage level is still ultimately unsuccessful in eliminating the tumor. Another regimen Loizides et al experiment with is to administer mice with 100 mg/kg 5-FU weekly for 12 weeks. The result of this course of treatment can be seen in Figure 5.

5-FU is a very commonly used drug to treat various cancers including breast cancer, rectal cancer, cancers of head and neck and skin cancer. It is particularly effective for colorectal cancer ^[66]. Smith et al. find that continuous infusional 5-FU chemotherapy is more effective on colorectal cancer than conventional chemotherapy in ^[67]. In the trial that is addressed in the work of Smith et al., 5-FU is administered continuously by an ambulatory pump for 18 weeks at a dose of $200$ mg/m$^2.$ When this human dosage is converted to mice, we see the resulting cell populations in Figure 6. Smith et al. report that toxicity was well tolerated however, lethargy, diarrhea and vomiting occurred more frequently in the group treated with 5-FU than the group who received conventional chemotherapy.

We experimented with the two weekly doses, 50 mg/kg and 100 mg/kg of 5-FU as in the study of Loizides by administering the chemotherapy twice a week, half the dose each time. The 100/2 mg/kg twice weekly 5-FU keeps the tumor load lower until about day 100, even though the outcome remains the same. The behavior of the tumor in both of this regimens can be seen in Figures 7 and 8.

Despite experimenting with the dosage and the frequency of drug administration, we have not been able to eliminate the tumor using chemotherapy, only. More specifically, we have simulated: 100/7 mg/kg, daily for 12 weeks, 100/7 mg/kg daily for 24 weeks and 450 mg/m$^2$ daily for five days, as presented in ^[66]. Caroli et al discuss the serious adverse effects of a high dose of 5-FU, more specifically at a level of 2600 mg/m$^2$ which is equivalent to 70 mg/kg for humans, administered weekly ^[68]. Therefore, we did not exceed 100 mg/kg in our simulations.

5.   Sensitivity analysis

Sensitivity analysis was performed on both the non-treatment and chemotherapy models using Partial Rank Correlation Coefficients (PRCC) as a metric. Partial correlation metrics measure the effect of an input $x_j$ on an output $y$ after removing the effects of all other inputs on $y$. Specifically, the partial correlation coefficient between an input $x_j$ and an associated scalar output $y$ is computed by measuring the correlation between the residuals $x_j-\hat x_j$ and $y-\hat y$, where

We use the rank-transformed partial rank correlation coefficients as these tend to perform better in cases where the relationship between the parameter and output is nonlinear but monotonic ^[69]. Note that parameters to which the output is more sensitive have PRCCs close to 1 or -1, where a negative score indicates that the parameter and the quantity of interest are inversely related. The scalar quantity of interest (QoI) was defined to be the number of cancer stem cells ($S$) plus the number of tumor cells ($T$) at the final day of the data observation period. For the computation of the PRCC scores, 1000 parameter samples were generated by sampling each parameter uniformly from an interval ranging from $\pm 25\%$ of the nominal parameter value (listed in Table 1).

Figure 9 illustrates the results of the sensitivity analysis for both models. In both cases, the QoI is highly sensitive and inversely related to the tumor cell death rate ($\delta_T$), and also highly sensitive but positively related to the tumor cell proliferation rate ($\eta$). In comparison, the QoI has negligible sensitivity to all other parameters.

6.   Stability analysis

Using the no-treatment model (2.1)-(2.5), we analyze the stability of the system. With five highly-coupled states and 23 parameters, the Jacobian matrix cannot be determined symbolically, so we proceed using the parameter values from Table 1. Our system of nonlinear differential equations yields $18$ equilibrium points, only two of which are real, nonnegative, and biologically relevant. The first, $(S, E, B, R, T) = (0, 0, 0, 0, 0)$, represents the tumor-free equilibrium, while $(S, E, B, R, T) = (5.999998213, 0, 0, 0, 1.009648946\times 10^8)$ represents the equilibrium with the tumor near maximum capacity. The lack of additional biologically relevant equilibria suggests that there is no set of state variables at which a small, nonlethal tumor can be controlled by an effective immune system. We investigate the stability of our two equilibrium points below.

6.1.   The Tumor-free equilibrium $(0, 0, 0, 0, 0)$

Linearization of the system of equations (2.1)-(2.5) yields the Jacobian matrix

with eigenvalues

The existence of three negative and two positive eigenvalues at the origin indicates that this equilibrium is a saddle point. To gain insight into directions which are stable versus unstable, we investigate the eigenvectors corresponding to each of the eigenvalues given as columns of the following matrix:

Based on the eigenvectors, it can be seen that the two unstable directions align roughly with the $S$ and $T$ axes, while the three stable directions exist in $E-B-R$ space. That is, a solution trajectory approaching the origin will eventually veer away towards a non-zero tumor state.

6.2.   The maximum tumor capacity equilibrium $(5.999998213, 0, 0, 0, 1.009648946\times 10^8)$

This second equilibrium represents the system near maximum tumor capacity (note that tumor carrying capacity was defined to be $9.06865 \times 10^8$ cells in Section 3). The Jacobian matrix is given by

with eigenvalues

The eigenvalues indicate stability in three directions, with a spiral source in the remaining two-dimensional plane. The corresponding eigenvectors are the columns of the following matrix:

The dominance of the tumor state in the fourth eigenvector (corresponding to a negative eigenvalue), suggests that solution trajectories approaching the equilibrium along the $T$ axis will be drawn towards a maximum tumor capacity state.

Overall, the stability analysis of this system suggests that unless the cancer stem cells and tumor cells are eradicated completely ($S = T = 0$), simulation trajectories will inevitably converge to a maximum tumor capacity equilibrium.

7.   Results

In this work, we have modeled the role of cancer stem cells and TGF-$\beta$ with no-treatment and with chemotherapy using experimental data from ^[40] to calibrate and validate the models. With no-treatment, the tumor grows as indicated by experimental data in ^[40], as can be seen in Figure 2. In this case, the tumor reaches its carrying capacity while the CD8$^+$ T-cell population approaches zero. Although the CD8$^+$ T-cells will never actually vanish, this behavior indicates that the CD8$^+$ T-cell population is functionally zero, as they are shut-down by the TGF-$\beta$ whose concentration reaches a maximum. Stability analysis of the no-treatment model reveals two unstable equilibria. An analysis of the tumor-free equilibrium predicts a growing tumor in the presence of tumor or cancer stem cells. The close to maximum tumor capacity equilibrium is stable in the direction of the tumor and hence, yields a similar predicament. In the presence of tumor or cancer stem cells, the tumor cells reach their maximum population.

It is shown that under chemotherapy, which targets only tumor cells but negatively affects the immune cells, the cancer will relapse. As chemotherapy does not target the cancer stem cells, we conclude that the tumor cells recover due to the presence of the cancer stem cells. This matches current experimental data ^[40,70,71]. This behavior is illustrated in Figure 4. The CD8$^+$ T-cells become functionally zero, as in the no-treatment case, with maximum tumor capacity reached shortly after termination of the chemotherapy.

Due to the simplifying assumptions made, this model can only provide a theoretical framework for the actual physical system and reveal the impacts cancer stem cell populations have on the tumor and immune system from an abstract perspective.

Funding statement

This work was supported by NSF DMS $[1239280]$. EZ was supported by the St. Joseph's College Faculty Summer Research Grant.

Acknowledgements

The authors would like to thank Dr. Ami Radunskaya for valuable discussions.

Conflict of interest

All authors declare no conflicts of interest in this paper.