#### Abstract

It is increasingly argued that cancer stem cells are not a cellular phenotype but rather a transient state that cells can acquire, either through intrinsic signaling cascades or in response to environmental cues. While cancer stem cell plasticity is generally associated with increased aggressiveness and treatment resistance, we set out to thoroughly investigate the impact of different rates of plasticity on early and late tumor growth dynamics and the response to therapy. We develop an agent-based model of cancer stem cell driven tumor growth, in which plasticity is defined as a spontaneous transition between stem and nonstem cancer cell states. Simulations of the model show that plasticity can substantially increase tumor growth rate and invasion. At high rates of plasticity, however, the cells get exhausted and the tumor will undergo spontaneous remission in the long term. In a series of *in silico* trials, we show that such remission can be facilitated through radiotherapy. The presented study suggests that stem cell plasticity has rather complex, nonintuitive implications on tumor growth and treatment response. Further theoretical, experimental, and integrated studies are needed to fully decipher cancer stem cell plasticity and how it can be harnessed for novel therapeutic approaches.

#### 1. Introduction

After stem cells have been discovered at the top of the hematopoietic system hierarchy [1], it became apparent that human acute myeloid leukemia is also organized hierarchically. Leukemia is initiated and fueled by a leukemic stem cell that gives rise to transit-amplifying progenitor cells and eventually differentiated cancer cells with limited lifespan [2]. A cellular hierarchy in a tumor and the accompanying stem cell hypothesis has been hailed as a significant breakthrough in the cancer research community, as its concept holds new promises for cancer therapy. If only CSCs are uniquely able to initiate, sustain, and propagate a tumor, then, selective eradication of CSCs, however difficult it might be to target them, would be sufficient to cure a cancer [3]. Despite its conceptual beauty, recent reports consolidated the skepticism that stemness might not be a prescribed cell phenotype but a transient state that cells can acquire and discard depending on the cellular environment and signaling context [4–7]. Then, the heterogeneous tumor population becomes a dynamic, moving target that is increasingly difficult to treat [8].

The complex biology of stem and nonstem cancer cells and their interactions with each other as well as with the intra- and extratumoral environment is yet to be fully deciphered experimentally. Inroads have been made to use mathematical and computational models to identify first-order principles and key biological mechanisms in cancer plasticity, from which new actionable hypotheses can be derived [7, 9–12]. Herein, we propose an* in silico *agent-based computational model to help decipher parts of the complexity that arises from the myriads of stem and nonstem cancer cell interactions and phenotypic plasticity. Agent-based models are increasingly utilized in theoretical oncology [13–20] to derive emerging population level dynamics from defined single cell properties and their perturbation. Such modeling approach has previously shown that the proliferation capacity (or telomere length [21, 22] or Hayflick limit [23, 24]) of nonstem cancer cell (CC) is a pivotal force in driving tumor evolution [25, 26]. Then, early progenitor cells that adapt a stem cell state confer different kinetic properties including proliferation potential to the new stem cell than dedifferentiating cells that are closer to their terminal phenotype.

The role of the cellular hierarchy in a tumor has been widely ignored in the modeling community, such that CCs that transition into a CSC have the same properties as the initial population-founding CSC. Thus, in addition to conferring the general CSC traits of (a)symmetric division and longevity, the new CSC gets also bestowed with the historic initial telomere length with the stupendous consequences of increased aggressiveness and treatment failure. Herein, we give explicit consideration of the degrading proliferative potential in the cellular hierarchy in a phenotypic plasticity model. We show that phenotypic plasticity promotes early tumor growth; in the long term, however, plasticity can impede tumor progression and ultimately lead to the collapse of the tumor population. We will use radiotherapy as an example of an external catalyst to eventually complete remission.

#### 2. Materials and Methods

##### 2.1. Mathematical Model

We adapt an* in silico* agent-based model [26, 27], in which each cancer cell occupies a 10 × 10 *μ*m grid point on a dynamically expanding 2D lattice [28]. A dynamic computational domain prevents boundary-imposed spatial constrains that may influence outcomes.

The tumor population is divided into cancer stem cells (CSCs) and nonstem cancer cells (CCs) with an individual proliferation capacity, *ρ*, representative of the telomere length [21, 29]. Telomeres are shortened during mitosis [22, 30] reducing the proliferation capacity in each daughter cell (), which is a visualization of the Hayflick limit [23, 24]. CSCs are believed to upregulate telomerase which rebuilds telomeric DNA and thus prevents telomere erosion and confers longevity to the cell [31–34]. We assume a cell with exhausted proliferation capacity () will undergo cell death in the next mitotic attempt as previously assumed [25, 35] without explicit consideration of replicative senescence [18, 36]. In addition to replicative cell death, we consider spontaneous cell death in CC with probability *α*, which is prevented in CSC. Without considering cell plasticity, CC and CSC populations are only connected through asymmetric division of a CSC, when one of the progeny adopts a CC phenotype (Figure 1(a)). We denote the probability of symmetric CSC division by . Plasticity may occur at successful proliferation with probabilities (CSC differentiation) and (CC dedifferentiation). If , plasticity is averted and cells have a persistent phenotype. With and , cell phenotypes are plastic and stemness becomes a transient state. We note that stemness is defined by the ability to divide (a)symmetrically and prevention of telomere erosion [33, 37]. Therefore, a CC that adopts a CSC state through dedifferentiation will be equipped with current telomere length, that is, , which will be bequeathed to subsequent daughter cells (Figure 1(b)). This is in stark contrast to previous modeling attempts that in addition to conferring general CSC traits also reset time and equip new CSC with historic uniform initial telomere length.

**(a)**

**(b)**

**(c)**

At discrete simulation time steps representative of hour, cells are randomly selected and updated. In case of a CC, spontaneous cell death is considered with probability . Proliferation and migration of surviving cells are mutually exclusive (, i.e., once per day; , i.e., 150 *μ*m per day) and subject to available space in the immediate cell neighborhood. All model parameters are summarized in Table 1, and the simulation procedure is visualized in a flowchart in Figure 1(c).

##### 2.2. Tumor Morphology Analysis

In case of frequent migration events, defining tumor periphery on a two-dimensional lattice is not straightforward, as cells may separate from the main tumor mass. In order to define the tumor periphery, we first transform the lattice into binary information (cell present or not). Then, we substitute the value of each pixel by the average number of positive pixels in the immediate 8 neighbors’ Moore neighborhood and apply an image intensity threshold of 3/8 (this includes vacant sites with more than 2 cancer cells in the neighborhood) and select regions containing more than 10^{4} cells. We define the tumor periphery as the periphery of the selected regions. All transformations are performed using MATLAB and Image Processing Toolbox Release 2014a, The MathWorks, Inc., Natick, Massachusetts, United States.

Tumor circularity is calculated aswhere is the tumor area and is the perimeter length. For a perfect circle, the circularity is equal to 1. Both and are calculated using the* regionprops* function from MATLAB Image Processing Toolbox.

In order to calculate the CSC fraction in the proximity of the tumor boundary, we dilate the tumor periphery by 20 pixels in four orthogonal directions and select all cells within the dilated area. This procedure selects cells that are within 200 *μ*m from the tumor boundary.

##### 2.3. Metastatic Potential

We estimate the potential for metastatic spread using a virtual 5-well plate, where wells are connected sequentially by a narrow canal of 100 *μ*m width and 400 *μ*m length. Each well is circular with a 500 *μ*m radius. Each simulation is initiated with a single CSC in the geometric center of the first well. Metastatic spread is simulated until the first CSC enters the last blind-ended canal.

##### 2.4. Radiotherapy

We simulate fractionated radiotherapy of 30 fractions of dose Gy applied every 24 h hours. The dose-dependent surviving fraction is calculated using the linear-quadratic formalismwhere describes radioresistance of quiescent cell, that is, for a cell with no available space in the neighborhood, and otherwise and describes increased radioresistance of CSC [38, 39]. For CC, . In line with previous estimates, we set , , , and [39].

#### 3. Results and Discussion

We initiate each simulation with a single CSC with proliferation capacity and probability of symmetric division and simulate tumor growth for 720 days. We consider bidirectional plasticity with equal probabilities (static phenotypes), 0.01%, 0.1%, 1%, and 10%. For each set of parameters, we performed 100 simulations, and statistical analyses were performed using Student’s -test.

##### 3.1. Impact of Plasticity on Tumor Growth Characteristics

Phenotypic plasticity increases initial tumor growth rate yielding larger tumors after 720 days than the tumors with static phenotypes (Figure 2(a)). However, for larger plasticity (), tumor growth saturates around day 320 keeping the tumor in a dormant state followed by a decrease in total cell number. In contrast with initial growth that is favored by higher plasticity rates, increasing phenotypic plasticity inhibits tumor growth later. This population level behavior mimics the evolution of CSC number (Figure 2(b)) and ratio (Figure 2(c)). CSC fraction in phenotypic plasticity tumors appears to saturate around the value of 50% for all considered transition rates, but this ratio is only achieved in the considered time frame with .

**(a)**

**(b)**

**(c)**

**(d)**

In addition to increasing the average tumor size during the growth phase, plasticity reduces the amount of variation in tumor size (Figure 2(d)). For , standard deviation is only about 20% of the mean at days, compared to about 3 times larger S.D./mean for tumors with static phenotypes. The source of variation in tumor size across independent simulation is opportunistic competition for available space between CSC and CC in the tumor interior where most CSCs are located without plasticity [19, 25, 35, 40]. Phenotypic plasticity partially averts intratumoral CSC inhibition and prolonged phases of population level dormancy [41] as new CSCs are continuously created at the boundary with fewer spatial constrains.

Interestingly, in case of , one simulated tumor died spontaneously early during development. In that simulation, the initial CSC differentiated before the first stochastic symmetric division event, and all of its CC progenies died before a dedifferentiation event. In the Appendix, we show analytically that the probability of that event is approximately 0.5% for all considered plasticity rates and, thus, is large enough to manifest in numerical simulations.

To analyze analytically if tumors with phenotypic plasticity can undergo spontaneous remission, we first consider the possible division fates of a CSC with proliferation capacity (). With probability , the CSC will differentiate before a symmetric division and it will be lost. In case of a dedifferentiation event of a CC, the new CSC will have a proliferation capacity (Figure 3(a)). If the CSC divides symmetrically, a new will be created. In case of asymmetric division, we need to follow the fate of the CC progeny. At each iteration, the CC can die spontaneously with probability *α*. Let us denote by the probability that the newly created CC will die before a proliferation attempt (). Then, no new will be created with probability . If the CC divides, a new is only created in case of a dedifferentiation event with probability . If after dedifferentiation the undergoes symmetric division, then two new are to be created. This theoretical consideration has four possible outcomes: the number of (1) decreases with probability , (2) remains the same with probability , (3) increases by one with probability , and (4) increases by two with probability . Hence, the probability that all will die after the initial seeding of one fulfills the relation , where and , from which we obtain thatObviously, for , that is, no possible differentiation event, we have , and for the problem reduces to . Importantly, if , then all CSCs will die off regardless of the initial number of CSCs. Equation contains only one unknown parameter, , which depends on the extent of spatial inhibition and the proliferation probability . However, the probability, , is an increasing function of , and . Thus, if a tumor dies spontaneously with probability 1 for , then it also dies spontaneously with probability 1 for any larger . Similarly, if for , then the tumor will not die spontaneously in each single simulation iteration. It is worth mentioning that a nonzero value of parameter is crucial for relating the agent-based model to probability , as the average time between the CSC proliferation events can grow without bound if , due to intratumoral spatial inhibition.

**(a)**

**(b)**

**(c)**

**(d)**

Figure 3(b) shows the probability of a CSC population vanishing for and if . All simulated tumors will eventually die off spontaneously for values larger than ≈9.09%. This explains the previously observed decrease in the average tumor size for (cf., Figure 2(a)). For in the initially presented tumor growth simulations, for , and thus about 90% of tumors will grow successfully.

Let us consider 100 independent simulations initialized with CSC_{3} and . As predicted by the above theory, all 100 tumors will undergo remission as in that case (Figure 3(c)). Increasing the probability of symmetric division from to , only 17 of 100 tumors die out (Figure 3(d)). However, the above theory cannot conclusively predict if all of these tumors will eventually die off.

##### 3.2. Impact of Plasticity on Tumor Morphology

Figure 4(a) shows the largest tumors after being simulated for 720 days of 100 independent simulations for different plasticity probabilities. Although the average tumor size for is larger than for (cf., Figure 1(a)), the biggest simulated tumor is smaller (114,950 versus 283,504 cells). This is a manifestation of the larger coefficient of variation in the static phenotype cohort (cf., Figure 2(d)).

**(a)**

**(b)**

**(c)**

**(d)**

The morphology of tumors with static phenotypes is described as self-metastatic [35] with low circularity [42]. Tumors with intermediate plasticity probabilities, and , feature more circular morphology at day 720 as the self-metastatic morphology caused by the spatial inhibition of CSCs is averted by spontaneous CCs dedifferentiation at the tumor periphery (Figure 4(b)). For larger plasticity probabilities (), clusters of cells at the tumor periphery die off stochastically, which yields a less regular boundary and decreasing tumor circularity.

Low values of plasticity probabilities ( and ) are associated with a fivefold increase in the fraction of quiescent cells compared to tumors with static phenotypes, arguably at least in part due to the larger size and thus reduced surface-to-volume ratio (Figure 4(c)). For , the exhaustion and spontaneous death of CC are not effectively filled and thus a larger proportion of cells is actively proliferating, which is comparable to the one in the nonplastic tumor.

A monotonic dependence on plasticity probability is the fraction of cancer stem cells in the proximity of the tumor boundary (Figure 4(d)). For , more than 50% of the CSCs are close to the tumor boundary, which represents almost the entire cancer stem cell fraction (cf. Figure 2(c)). The prevalence of CSC in the tumor periphery suggests that plasticity may lead to increased potential for metastatic spread.

##### 3.3. Impact of Plasticity on Invasiveness

For each considered value of plasticity probability, we simulated 100 virtual 5-well plate experiments. For , 0.1%, and 1%, in all simulations a CSC successfully reached the right boundary of the experimental setup. However, for , as many as 97 out of 100 simulations ended with spontaneous death of all cancer cells before invading all wells. Time to reach the right boundary was significantly lower in all plastic tumors, with almost 1.5-fold and 2-fold reduction for and , respectively, when compared to tumors with static phenotypes (Figure 5(a)). However, the invasion speed follows a nonmonotonic behavior in plastic populations, as the time needed to reach the right boundary is reduced for compared to both and , indicating exhaustion of the population as discussed above.

**(a)**

**(b)**

Visualizations of successful simulations show a plasticity-dependent gradient of cancer stem cell fraction values in subsequent wells (Figure 5(b)). With lower values of plasticity, the difference between the composition of population in the first and the last wells increases, with about 16-fold change for nonplastic tumor. For the largest considered value of , however, there is no significant difference in composition between the wells.

##### 3.4. Impact of Plasticity on Radiation Outcome

We simulated the impact of radiotherapy on tumors that consisted of 250,000 cells generated for different plasticity probabilities. As expected, the tumor without plasticity responds best to radiotherapy (Figure 6(a)), as it has the smallest fraction of CSC with decreased radiosensitivity (cf., Figure 2(c)) and small proportion of quiescent cells with decreased radiosensitivity (cf., Figure 4(c)). The tumors with relatively low plasticity, and , show similar response to radiation, with the number of cells reduced by about 20 times compared to the pretreatment state. Tumor with the highest plasticity event probability, , shows an intermediate response, which can be explained by lower fraction of quiescent cells compared to other plastic tumors (Figure 4(c)) and larger fraction of stem cells compared to the plasticity-free tumor (Figure 2(c)). However, in none of the simulations was the tumor eradicated during the course of treatment.

**(a)**

**(b)**

**(c)**

Simulations of tumor growth following the treatment show that tumors without and with low plasticity show regrowth to pretreatment cell counts within 200 days after the onset of treatment (Figure 6(b)). However, almost all of the simulated tumors with undergo spontaneous remission (Figure 6(c)) within two years after transient regrowth. Radiotherapy on highly plastic tumors may serve as an accelerator of CSC exhaustion due to frequent differentiation events with subsequent loss of proliferation potential.

#### 4. Conclusions

The observation of cancer stem cell plasticity has led to the general belief of more aggressive tumor growth and reduced treatment response [12]. We introduced an* in silico* agent-based model of cancer stem cell driven tumor growth to study the impact of different rates of phenotypic plasticity on tumor growth, morphology, invasion, and treatment response. Simulations of our model show that plasticity accelerates early tumor growth, as prolonged phases of tumor dormancy due to intratumoral competition [41] are averted by cells on the periphery acquiring stem cell traits and propagating the tumor population. Whilst tumors with fixed phenotypes show a gradual increase in CSC over time, phenotypic plasticity yields relatively constant population fractions, with CSC predominantly located at the tumor boundary, thereby facilitating invasion and metastatic spread. The explicit consideration of the cellular hierarchy within a tumor and the accompanying reduction of proliferation potential, however, offer a previously unappreciated aspect of CSC plasticity. Transitions between stem and nonstem cancer cell states at high frequencies yield a reduction of proliferation potential in each cell, thereby reducing the lifespan of each daughter cell and inevitably cell death. This may lead to population level dormancy and ultimately collapse of the tumor with complete remission.

Mathematical and computational models, by virtue of their very purpose, are subject to gross oversimplifications of reality [43]. However, they may provide interesting and nonintuitive insights into the tumor progression dynamics, including novel discussion points in understanding cancer stem cell plasticity. Further theoretical, experimental, and integrated studies are needed to fully decipher cancer stem cell plasticity and how it can be harnessed for novel therapeutic approaches. In a simple experimental approach, cancer stem cell driven tumor cell lines can be propagated* in vitro*, and at different passages cells may be sorted for cancer stem cell marker expression. The negative population can be exposed to radiation or hypoxic and acidic environmental conditions, which have been shown to reprogram cells towards a cancer stem cell phenotype [44–47]. Repeated sorting of the previously marker-negative population should yield marker-positive cells, arguably plastic new cancer stem cells, which can be further propagated* in vitro*. With increasing passages, our model suggests a significant reduction in telomere length of marker-positive cells, which may translate into different* in vivo *tumor growth dynamics.

#### Appendix

#### Probability of Early Tumor Death Event

We set out to approximate analytically the probability that the initial CSC with proliferation capacity *ρ* will differentiate before first symmetric division and all of its nonstem progenies will die before any dedifferentiation event, which we denote as . We consider a very early stage of tumor development and, thus, we can neglect space inhibition (we have a large migration probability of 15/24). The probability that a nonstem cell with given proliferation potential, *ρ*, and all its progenies will not dedifferentiate before dying off, , follows the recurrence relationwith . The probability that the CSC will undergo differentiation before first symmetric division is equal to and the number of nonsymmetric divisions follows geometric distribution with probability of success . Thus, we obtain the following expression for :In Figure 7, we plotted under the assumption that and for parameters considered in the main text. We see that for all considered plasticity event probabilities the probability of early tumor death is close to 0.5%. It is worth noticing that if spatial inhibition will increase value more for lower probabilities of plasticity event, on average more nonstem progenies are created before differentiation event.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgment

Jan Poleszczuk would like to thank the Foundation for Polish Science for partial support.