 Research article
 Open Access
 Published:
Modeling treatmentdependent glioma growth including a dormant tumor cell subpopulation
BMC Cancer volume 18, Article number: 376 (2018)
Abstract
Background
Tumors comprise a variety of specialized cell phenotypes adapted to different ecological niches that massively influence the tumor growth and its response to treatment.
Methods
In the background of glioblastoma multiforme, a highly malignant brain tumor, we consider a rapid proliferating phenotype that appears susceptible to treatment, and a dormant phenotype which lacks this pronounced proliferative ability and is not affected by standard therapeutic strategies. To gain insight in the dynamically changing proportions of different tumor cell phenotypes under different treatment conditions, we develop a mathematical model and underline our assumptions with experimental data.
Results
We show that both cell phenotypes contribute to the distinct composition of the tumor, especially in cycling low and high dose treatment, and therefore may influence the tumor growth in a phenotype specific way.
Conclusion
Our model of the dynamic proportions of dormant and rapidly growing glioblastoma cells in different therapy settings suggests that phenotypically different cells should be considered to plan dose and duration of treatment schedules.
Background
Gliomas are the most common type of primary brain tumors including their highly malignant form, the glioblastoma multiforme (GBM), which accounts for about 15% of all brain tumors [1]. Despite current standard treatment of GBM by surgical resection and adjuvant radio and chemotherapy, the median survival time for GBM patients is still poor, approximating 12–15 months [2], mostly due to unsatisfactory response of the tumor to treatment strategies. Additionally, combined aggressive radio−/chemotherapy is causing severe side effects frequently necessitating interruptions of the therapy due to e.g. blood toxicity [3]. GBMs and also many other tumors are heterogeneous tumors, being composed of cells with different, partly specialized phenotypes [4]. Besides e.g. rapidly proliferating tumors cells, invading immune cells, endothelial cells and (tumor) stem cells, also a subpopulation of so called dormant tumor cells exists in the heterogeneous tumor mass. These cells enter a quiescent state driven by cellintrinsic or extrinsic factors, including permanent competition for nutrients, oxygen, and space (“cellular dormancy”) [5,6,7,8]. In several tumors and metastases, dormant cells have been shown to be not proliferative or only very slowly cycling [9,10,11,12]. Linking dormancy and effects of chemotherapy, studies on glioma cells showed that cells underwent a prolonged cell cycle arrest upon treatment with temozolomide (TMZ), the most common chemotherapeutic in GBM therapy [13].
Evolutionary forces, such as competition and selection, shape the growth of the tumor and therefore the progression of the cancer. These forces create different ecological niches within the tumor encouraging the adaption of specialized tumor cell phenotypes. Accordingly, the proportional balance between different tumor cellular phenotypes can drastically change with treatment conditions. Indeed, compared to rapidly proliferating tumor cells, especially dormant cells exhibit a much higher robustness against chemotherapeutic drugs [5]. This dormant state seems to be reversible [13], so that the conversion to dormancy and the exit from dormancy may be a mechanism that facilitates tumor survival and progression even upon adverse or changing conditions. Hence, a better understanding of the proportional dynamics of different cell phenotypes within gliomas under chemotherapeutic treatment may improve further therapeutic approaches.
Mathematical models are beneficial resources to gain insight into key mechanisms of cancer development, growth, and evolution and to help identifying potential therapeutic targets [14]. Among these approaches, evolutionary game theory [15, 16] models the interactions between different individuals as a game between agents playing different strategies and relates the payoff from this game to the reproductive fitness of the corresponding agent [17,18,19,20,21].
Here, we use evolutionary game theory to model the proportions of two different phenotypes of GBM cells in a variety of different treatment conditions, see Basanta and Deutsch [18] for a related approach in GBM. Defining the fitness of the different cell types as growth rate in comparison to cells of the respective other phenotype, we focus especially on the balance between the rapidly proliferating and the cellular dormant phenotype and describe the corresponding payoffs in a payoff matrix which also includes the effect of treatment. Then, we use a special form of the replicatormutator equation [22, 23], which takes into account that conversion from dormant to rapidly proliferating phenotype and vice versa is possible. To strengthen our theoretical assumptions, we analyzed cell numbers and the cellular expression of a dormancy marker under different chemotherapy dosages and the phenotypic conversion modalities in cultured GBM cells in vitro. Taken together, the aim of our study was to develop a simple theoretical model which describes the dynamically changing proportions of two different GBM cell phenotypes, rapidly proliferating and dormant cells, under different treatment conditions. Showing this, we suggest that different properties of cell phenotypes should be taken into account for the development of more efficient, less toxic treatment schedules in order to improve patient’s prognosis and quality of life.
Methods
Theoretical model
We analyze the proportions of two different GBM cell phenotypes, dormant (D, please refer to Table 1 for symbols used in the equations) and rapidly proliferating (P) cells, in a mathematical model including the influence of different treatment conditions. In the following, we characterize the cells in terms of their fitness, which we define as the growth rate in comparison to cells of the other phenotype. Dormant cells always have a very low or even zero growth rate ε, which we assume to be independent of the exact composition of the population and the treatment condition. Rapidly proliferating P cells, on the other hand, have a very large fitness advantage compared to dormant cells, which means they proliferate much faster, but they also compete with each other for space and resources. Facing another P cell, a focal P cell has an intermediate fitness, which we assume to be still much larger than the growth rate of D cells ε. Their fitness therefore depends on the relative fraction of D vs. P cells. Due to the very slow growth of D cells, P cells will represent the vast majority of glioblastoma cells in the absence of treatment.
Under treatment conditions, however, the population composition changes. Even though D cells still have the same very low (or zero) growth rate ε, P cells experience a fitness cost λ due to treatment. The reduction of the fitness due to treatment only applies to P cells, because cytotoxic drugs mostly affect rapidly dividing cells. The fitness cost parameter λ can be adjusted to account for the strength of the applied treatment. In principle we can continuously vary this parameter. However, for simplicity we focus on two different treatment strategies: In high dosage (HD) chemotherapy the treatment strength parameter λ is large compared to the growth rate of the P type. Since high dosage chemotherapy has strong side effects for the whole organism (for GBM: [3]), in reality this treatment strategy cannot be maintained for extended time periods. Therefore, strong treatment needs to be applied in turns with weaker or no treatment. For low dosage (LD) chemotherapy, λ means only a small reduction of the growth rate of the P cells. As the sideeffect stress to the organism should also be lower, this treatment regime could be applied for longer time spans.
Dormant (D) and rapidly proliferating (P) phenotypes in glioblastoma and their aforementioned interactions can be described by the following payoff matrix [18]:
This matrix gives the fitness for each type if confronted with any of the two other types. Here, we find for example that the fitness of a focal P cell interacting with a D cell is 1 − ε − λ, which includes both the small or zero growth rate of D cells ε and the fitness cost for P cells under treatment λ.
As the phenomenon of dormancy is presumably a reversible process that also occurs without any treatment, we assume that conversion between both phenotypes is possible with a small rate σ. Thus, P cells may enter a dormant phenotype, and D cells may exit from their quiescent state, converting into a P phenotype at any time point.
In the following, we include these fitness effects and phenotypic conversion into a set of ordinary differential equations. In general, the growth of a whole cell population can be explained in terms of a differential equation that describes the change in the number of individuals over time
Here n is the number of individuals, t is the time and r(n, t) is the growth rate, which can itself depend on the number of cells and the time.
At first, we focus on the number of D cells, n_{ D }, in the population over time, which have a very small but constant growth rate ε
For P cells on the other hand, the growth rate of n_{ P, } given by the average fitness from the payoff matrix (weighted to the cell fractions), changes with the composition of the population
Since the system under consideration is constrained, both in terms of nutrients and space, in reality the cell population only grows exponentially as indicated by the growth equations in the very beginning of the process where the constraints regarding space or nutrients are negligible. However, we are mainly interested in the fraction of D cells \( {x}_D=\frac{n_D}{n_D+{n}_{\mathrm{P}}} \) in the population and vice versa the fraction of P cells \( {x}_{\mathrm{P}}=1{x}_D=\frac{n_{\mathrm{P}}}{n_D+{n}_{\mathrm{P}}} \). To obtain the change in fractions for both types, we subtract the average growth rate \( \overline{f} \) of the population from both individual growth rates,
From this we obtain two differential equations for the fractions of D and P cells,
Next, we include the spontaneous conversion between phenotypes with a constant rate σ, which is independent of the cellular growth. This leads to an additional term to the differential equation of both phenotypes
These equations have the important difference to the usual replicatormutator equation [15] that phenotype conversion is a spontaneous process with a constant rate and is independent of the growth in the population. This allows conversion from D to P even if D cells do not grow at all.
Using these equations, we model different therapy schedules combining different treatment strengths in different cycling time plans. Since the equations are nonlinear, we use numerical integration with Odeint of the Python library Scipy^{Footnote 1} to examine the temporal dynamics of the system under different treatment regimes. Additionally we analytically determine the fixed points of the system and their stability.
Experimental model
Cell culture and cell number determination
The GBM cell line LN229 was purchased from ATCC/LGC Standards (Middlesex, UK, ATCCCRL 2611) and cultured in Dulbecco’s modified eagle medium (DMEM) plus 10% fetal calf serum (FCS, PAN Biotech, Aidenbach, Germany). Mycoplasma contaminations were routinely excluded by bisbenzimide staining. The GBM cell line identity was proven routinely by STR (Short Tandem Repeat) profiling at the Department of Forensic Medicine (Kiel, Germany) using the Powerplex HS Genotyping Kit (Promega, Madison, WC). Briefly, DNA was amplified with a STR multiplex PCR, electrophoretic separation was performed with the 3500 Genetic Analyser (Thermo Fisher Scientific, Waltham, MA, USA), and evaluated using the Software GeneMapper IDX (Thermo Fisher Scientific). For determination of cell numbers after low and high dose chemotherapy treatment, 25,000 cells/well were seeded in 6 well plates (Greiner Bioone, Frickenhausen, Germany). Cells were grown for 24 h, then washed with phosphate buffered saline (PBS), supplemented with fresh DMEM + 10% FCS and temozolomide concentrations (SigmaAldrich, St. Louis, MO, USA; dissolved in dimethyl sulfoxide DMSO) as indicated in Fig. 2a (5, 50 or 100 μg/ml for 10 days). Temozolomide (TMZ) is a DNA alkylating drug causing apoptotic cell death and the most commonly used chemotherapeutic in GBM therapy. Control cells were supplemented with 0.5% DMSO, which corresponds to the solvent concentrations of each TMZ stimulated sample. Cells were stimulated for 10 days with TMZ, while media were changed every 2–3 days. After 10 days, cells were detached by trypsination and total cell numbers per well counted using trypan blue exclusion and a Neubauer chamber (Brand, Wertheim, Germany). DMSO stimulated control cells were already detached after 6 days of stimulation, split 1:10 and seeded again to exclude limitations of growth due to space and nutrient limitations. This splitting factor (1:10) was considered when relative cell numbers of TMZ treated samples in comparison to DMSO controls were determined for n = 5–6 independent experiments.
Immunocytochemistry
For immunocytochemistry, 50,000 cells were seeded onto polyDlysine coated glass cover slips, grown for 24 h and supplemented with indicated TMZ or DMSO concentrations as described above. From day 6, growth media were additionally supplemented with 10 μM 5bromo2′deoxyuridine (BrdU, SigmaAldrich, St. Louis, MO) to allow for incorporation in the DNA in the S phase of the cell cycle. After 10 days, cover slips were fixed with an icecold mixture of methanol and acetone (1:1) for 10 min, rinsed with 0.1% Tween / PBS (3 × 5 min), incubated with 1 M HCl for 30 min, neutralized with 0.1 M sodium borate buffer (pH 8.5), and rinsed again with 0.1% Tween/PBS. Afterwards, cells were blocked for unspecific bindings with 0.5% bovine serum albumin (BSA) / 0.5% glycine in PBS (1 h) and incubated over night with the primary antibody against H2BK (1:300, Biorbyt, Cambridge, UK), a marker of glioma dormancy [24, 25] and the primary antibody against BrdU (1:200, Abcam, Cambridge, UK). Then cover slips were incubated with the secondary antibodies (donkey antirabbit IgG, labelled with Alexa Fluor 488, and donkey antisheep labelled with Alexa Fluor 555, both Invitrogen, Carlsbad, CA, USA) for 1 h at 37°, and 4′, 6diamidino2phenylindole (DAPI; Sigma Aldrich, St. Louis, MO, USA; 1 mg/ml, 1:30,000, 30 min at room temperature) to stain nuclei. Cover slips were embedded using ImmuMount (Thermo Fisher Scientific, Rockford, IL, USA), and analysed with equal exposure times using an Axiovert microscope and digital camera (Zeiss, Jena, Germany). H2BKimmunopositive, BrdUpositive and double positive cells were counted and normalized to total cell numbers in 6 (DMSO controls) to 10 (TMZ samples) fields of view for n = 4 independent experiments.
DiO retention and cell countings on phenotype conversion
To monitor the conversion to and from dormancy we used the green fluorescent vital dye DiO (Invitrogen), as rapidly proliferating cells lose the dye due to repeated divisions, while resting, dormant (or very slowly cycling) cells retain the dye and can be detected by fluorescence microscopy. Investigating the conversion to dormancy, 150,000 LN229 cells were seeded into 6wellplates, stained with Vybrant® DiO CellLabeling Solution (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s instructions and stimulated with 100 μg/ml TMZ (or equal volume of the solvent DMSO) for 10–12 days. Cells were photographed combining transmittedlight microscopy and fluorescence microscopy with equal exposure times for TMZ and control treated cells, and green fluorescent cell portions were determined in comparison to total cell counts. To determine the influence of different cell densities on the incidence of conversion, 50,000 and 150,000 cells were seeded, respectively, into 6wellplates and treated with 100 μg/ml TMZ (or equal volumes of the solvent DMSO) for 10 days. As the DMSO control treated cells rapidly proliferate, cells were detached at day 6 (50,000) or day 3 and 6 (150,000), cell numbers counted using a Neubauer chamber to determine the growth rate over this time period, and seeded again at initial density, to allow for cell growth without limitation of space and nutrients. After 10 days, TMZ and control treated cells were detached and counted. To extrapolate the total cell numbers of control cells, growth rates determined at day 3, 6 and 10 were used, and TMZ surviving cells were calculated as percentage of extrapolated total cells.
Statistical analysis
Statistical analysis and graphical presentation of experimental data were performed with Graph Pad Prism using a twotailed ttest (*** p < 0.001).
Results
Modelling the dynamics of cell frequencies
The temporal dynamics of the proportion of D against P cells in GBM strongly depends on the treatment conditions. Therefore, we first analyze the fixed points of the dynamical system and how they change for different treatment strengths λ, without considering possible conversions of phenotypes. The fixed points mark a stable equilibrium between the portions of P and D cells under certain, predefined conditions and are found by setting Eq. 1 to zero. Of particular interest are stable fixed points, as the system returns into this state again after a small perturbation [26].
For our system, there is only one stable fixed point for each treatment condition (Fig. 1a). If we consider the case of no phenotype conversion, σ = 0, we can give the exact position of this point for each treatment condition λ. The fraction of D cells at the fixed point is then given by
For small treatment strengths λ the fraction of D cells in the population at the stable fixed point is zero, but after reaching a threshold, the fraction of dormant cells increases linearly with λ until the whole population consists of dormant cells at very high treatment strengths.
With conversion between the two phenotypes (σ > 0), the analytical calculation for the stable fixed points is more difficult and the result is not instructive. By contrast to the previous results without phenotype conversion, there is always a small proportion of dormant cells in the population, even at very low treatment strengths. The proportion of dormant cells at the fixed points increases immediately with increasing treatment strength until it approaches a maximum at high treatment strengths. For both D cell growth rates ε (ε = 0, orange lines; and ε = 0.1, blue lines) the population composition is very similar or even the same without phenotype conversion at very small or large treatment strengths λ. In contrast, the largest effect of ε on the population composition is at intermediate values of treatment strength.
The average fitness \( \overline{f} \) for the whole tumor cell population including P and D cells decreases linearly from the maximum at treatment strength λ = 0 until it reaches the minimum of \( \overline{f}=\varepsilon \) at the point where the fraction of dormant cells in the population starts to increase (Fig. 1b). Interestingly, with spontaneous conversion σ > 0, the average fitness at the fixed point can become smaller than ε and even negative for high treatment strengths, potentially leading to a shrinking tumor. This is caused by conversion of D cells into P cells which are then susceptible to treatment.
Comparison to experimental data
To test our mathematical model of phenotype composition upon treatment, we used LN229 cells as an experimental in vitro model. We treated these cells for 10 days with temozolomide (TMZ), the most commonly used cytotoxic drug in glioma therapy. In a first step, we focused on different treatment strength and analysed the portions of surviving cells in comparison to control cultures and the percentage of cells expressing H2BK (histone cluster 1), a marker of glioma dormancy [24, 25], alongside with incorporation of BrdU in the late treatment phase (day 6–10). In general, after 10 days of treatment, samples stimulated with 5, 50 and 100 μg/ml TMZ had significantly less total cell numbers than control treated cells (Fig. 2a). By immunocytochemistry of H2BK, we could detect and quantify the fraction of dormant cells within the cultures, and by adding BrdU to the cells from day 6 of treatment and immunocytochemical staining of BrdU, we could in parallel mark cells that incorporate BrdU in the DNA (examples of microscopic pictures in Fig. 2b). While DMSOtreated control cells showed a low fraction of H2BKpositive cells (mean: 9.7% ± 3.5), TMZ treatment yielded increased numbers of dormant cells reaching a plateau at high concentrations (5 μg/ml: mean 26.8% ± 9.0, 50 μg/ml: mean: 82.8% ± 5.3, 100 μg/ml: 87.7% ± 8.0, compare Fig. 2b, grey graph portions). In parallel, we investigated the incorporation of BrdU in the DNA and determined a high portion (66.0% ± 7.8) of BrdU positive cells in the control cultures and lower portions upon TMZ treatment (5 μg/ml: 53.6% ± 14.5; 50 μg/ml: 33.4% ± 5.3; 100 μg/ml: 33.7% ± 10.1, compare Fig. 2b, hatched graph portions). Interestingly, BrdU incorporation also took place in TMZ treated cultures, so that staining for the dormancy marker H2BK and for BrdU could be observed in the very same cells (compare examples of microscopic photographs in Fig. 2b) indicating that cell cycle arrest may occur after the S phase of the cell cycle. Together with our experiments described in the following section and Fig. 2c, showing that dormant cells hardly divide within our experimental time frame, these observations suggest that dormant glioma cells are not or only very slowly cycling. Furthermore, taking into account that we use a clonal cell line, the occurrence of dormant cells needs to be a phenotypic adaption to the environmental conditions as all cells are genetically homogenous (as proven by routinely STR profiling, compare Materials and Methods section).
To investigate if the conversion to a dormant phenotype depended on the cell density, initially, we determined in a DiO retention assay that nearly all cells (98.3% ± 1.2) retain the green fluorescent dye when treated with 100 μg/ml (“high dose”) TMZ for 10–12 days, while in control cultures (treated with equal volumes of the solvent DMSO) only 2.9% ± 1.7 retained the dye (Fig. 2c, left part). The vital dye is included in every cell at the moment of staining, and is transferred to every daughter cell upon cell division. However, this means the staining is diminishing after several divisions of rapidly proliferating cells, but retained in nonproliferative or very slowly cycling dormant cells. Thus, assuming that nearly all cells that survive treatment with 100 μg/ml TMZ are dormant in our particular setting, we determined the relative incidence of phenotype conversion and the influence of the cell density on this conversion factor by determination of TMZ surviving cells in relation to (extrapolated) total cell numbers of control (DMSO treated) cultures. In our experimental setting, a portion of 0.68% ± 0.13 cells of initially seeded 50,000 LN229 glioblastoma cells survived this high dose treatment, while in cultures of initially seeded 150,000 cells, the portion of surviving cells was nearly similar (0.66% ± 0.13; Fig. 2c, right part) underlining the assumptions of our theoretical model.
Thus, treatment with TMZ significantly reduced total cell numbers of LN229 cells, while the share of dormant cells within the culture, as detected by the dormancy marker H2BK, was drastically elevated. The incidence of conversion to dormancy did not depend on cell densities in our particular experimental setting.
Treatment schedules
Next, we use our model to analyze the dynamics of the population composition for periodically changing treatment conditions. One example trajectory for a growth rate for D cells ε = 0 and a conversion rate between phenotypes σ =0.1 is depicted in Fig. 3. The fraction of D and P cells in the population alternates between the fixed points corresponding to the momentary treatment condition. The trajectory starts with a phase of no treatment, which is characterized by a high average growth rate and a cell population composition of mostly P cells and only very few D cells. After the first phase of unconstrained growth large parts of the tumor are removed (e.g. by surgery), leaving only a small number of cancer cells. Under the following high dosage treatment conditions, the dormant phenotype has the highest fitness and takes over the population. The relative fraction of D cells will increase until the steady state under high treatment conditions is reached. The impact of treatment on P cells leads to a strong initial decline in average growth rate, until the population has a significant proportion of dormant cells and the growth rate starts to recover slightly.
Under the following low dosage treatment conditions, P cells (making up a small fraction of the whole population at the end of the high dosage treatment) are less affected by the treatment and now grow faster. The average growth rate will have a maximum when then relative fraction of P cells in the population is still low, since they have a competitive advantage over D cells, and then declines afterwards towards an equilibrium well above the high dosage growth rate. Accordingly, the total number of cells increases strongly again in this regime.
Switching the order of high dosage and low dosage treatment only has a small effect on total number of cells: If treatment starts with low dosage, the system will go into a state with a slightly higher fraction of dormant cells, which makes it less susceptible to the following high dosage treatment. Starting with low dosage therefore does not help to reduce the tumor size.
In Fig. 4 we compare three different treatment schedules: just one switch from initial high (H) dose to low (L) dose treatment (HHHHHHLLLLLL, Fig. 4a, each instance of the letter H or L corresponds to the same time interval), slow cyclic switching (HHHLLL, Fig. 4b), and fast cyclic switching (HLHL, Fig. 4c) for two different growth rates of D cells (left panels ε = 0 and right panels ε = 0.1). In case of only switching once, the fixed points for each treatment are quickly reached. At high dosage treatment the number of cells increases very slowly or even decreases. In the following low treatment phase, however, P cells take over growing particularly fast and jeopardizing any positive effect from the previous strong treatment. This is true for both treatment strengths of the high dosage phase.
For the fast switching treatment schedule (HLHL, Fig. 4c), the fixed point of the population proportion is not reached before the treatment changes again. Therefore the population dynamics stays between the two stable fixed points for the two treatment regimes, but does not reach them. By contrast, in the slow treatment switching regime (HHHLLL, Fig. 4b) the fixed points for both high and low dosage treatment are reached such that the composition of the cell population essentially resembles the case of just one switching event (Fig. 4a). However the time spent at these fixed points is still significantly reduced compared to only a single switch.
The bottom panels of Fig. 4a, b and c show the total number of cells based on the average fitness of the population under the assumption of exponential growth. When the growth rate is positive the cell population grows, otherwise it shrinks. Interestingly, the average growth rate of the population is well below zero only for a short period during the high dosage treatment and only if the share of P cells is still very high and the fraction of D cells in the population is small. However, in this regime the fitness recovers fast and approaches equilibrium with an average fitness close to zero, such that the total number of cells does not change anymore. The strongly negative growth rate directly after switching to the high dosage treatment is therefore the reason why the number of cells for quickly changing treatment regimes is significantly smaller than for slowly changing treatment cases.
Population growth
To systematically examine the effect of switching treatment cycles on the growth rate of the population, we analyze the temporal dynamics of the population size for varying treatment cycle durations and two different growth rates of D cells ε (Fig. 5). Unsurprisingly, a lower growth rate of D cells has a diminishing effect on the overall growth of the cells. For increasing treatment cycle length cancer cells have an increasing overall growth, while maintaining the same total high and low dosage treatment durations. The overall growth rate approaches a maximum with increasing treatment cycle length when the dynamics reaches equilibrium in each cycle.
Taken together, our mathematical model allows us to theoretically predict the fitness and proportions of rapidly proliferating and dormant tumor cells under different treatment conditions. Strengthening our theoretical assumptions we could exemplarily show the effect of high and low chemotherapy doses on the cell numbers and the proportion of a dormant cell phenotype in cultured GBM cells in vitro. Simulating different therapy schedules, we observed that fast switching of low and high treatment doses yields a lower total tumor cell number at equal total drug dose in comparison to low switching schedules.
Discussion
In this study, we established a mathematical model and analyzed the proportions of two different cell phenotypes occurring in GBM, rapidly proliferating and dormant cells. Corroborated by experimental data obtained from in vitro experiments with cultured GBM cells, we observed that treatment strength influences the balance between both phenotypes which in turn influences the growth of the whole tumor population. Sequential switching of treatment strength may thus drastically influence the proportion of dormant and rapidly proliferating cells, especially if switching to the next condition takes place before the population dynamics reaches a steady state.
Dormancy in GBM has been shown by the existence of distinct fraction of temporarily nonproliferative cells in murine models [27], as well as by the identification of clones which were able to generate indolent dormant tumors both in subcutaneous and orthotropic intracranial sites [28]. Additionally, dormancy seems to be characterized by specific features in GBM, such as a nonangiogenic phenotype [24, 28, 29], and is influenced by the (micro)environment e.g. hypoxia [30] and coagulation [31,32,33]. However, cellular dormancy in tumors is not only regarded as a state to overcome times of adverse conditions but has also been assigned to DNA repair mechanisms [34]. Interestingly, dormant GBM cells are hallmarked by the upregulation of specific genes like angiomotin, ephrin typeA receptor 5 (EphA5), insulinlike growth factorbinding protein 5 (IGFBP5), and histone cluster 1 (H2BK) [24, 25]. We used the latter as a marker to detect dormant cells in our in vitro experiments and show that the proportion of dormant cells increases with increasing chemotherapy concentrations.
As the fitness in the competition for space and resources depends on the proportions of phenotypically different cell subpopulations, we used evolutionary game theory as a framework for our mathematical model. Previous studies discussed game theoretic interactions with more phenotypes for many different types of cancer, including glioma [35], prostate cancer [36], and multiple myeloma [19, 37] or general tumors [38]. Also, evolutionary game theory is often used in spatially structured populations to answer questions about the effect of environmental constraints on tumor composition and invasiveness of cancer cells [20, 39, 40]. However, including spatial structure in order to increase the realism of the model leads to a large number of additional assumptions and potential pitfalls [41, 42]. Other modeling approaches for dormancy in cancer focus on the interaction between the immune system and the tumor [43,44,45], the effect of angiogenesis [46], or spatial competition between cells [47]. However, these approaches do not explicitly model the conversion between phenotypes and its consequences under therapy with varying strength.
Thus, we decided to simplify our model on several levels: (i) We do not take spatial structure into account. (ii) We abstract from the interaction with the immune system. (iii) We concentrate on two tumor cell phenotypes – rapidly proliferating and dormant cells – although other cell phenotypes, such as fast migrating cells, cells mimicking vasculature, (cancer) stem cells and invading immune cells (e.g. [48,49,50]), also contribute to the whole tumor mass. (iv) We focus on the fitness of the respective phenotypes rather than the potentially underlying reasons for phenotypical changes (e.g. genetic or epigenetic changes).
An important aspect of our model is the conversion between the different cell phenotypes. Recent studies suggest that dormant cells may originate from “normal” tumor cells by currently intensively investigated mechanisms (e.g. [51,52,53]). As a fundamental criterion for tumor dormancy, dormant cells need to be able to reawake and start growing again, so that they then phenotypically resemble the rapidly growing cell phenotype. Thus, we introduced a conversion factor σ into our model capturing these phenotypical transformation processes. Whether such conversions occur spontaneously or can be induced specifically or randomly by extrinsic or intrinsic mechanisms is poorly understood. We thus assumed a spontaneous event which can be modeled by a constant rate.
Using our theoretical approach we showed that dependent on the applied treatment strength an equilibrium balance between rapidly proliferating cells and dormant cells is eventually reached. At this fixed point and with low dosage or no treatment, mostly rapidly proliferating cells dominate the population, similar to the findings of Basanta and colleagues [35]. At stronger treatment, the fraction of dormant cells becomes successively larger, yielding a lower growth rate of the whole tumor. However, high dosage treatment cannot be applied for longer time periods, as it causes severe side effects (for GBM: [3]). Hence, we focused on alternative treatment schedules.
Several previous models discuss the effect different treatment schedules on various aspects of the cancer growth like angiogenesis [44, 54] or evolution of resistance [55]. Here, either the dosage and timing or the type of the chemotherapeutic drug is varied, which can have a massive effect on the growth of the tumor. Accordingly, we combine sequential cycles of low and high dose treatment with different durations. Thereby, we observed that the total growth of the tumor is considerably lower for fast switching compared to a slow switching scheme.
Conclusion
In this study, we have developed a theoretical model to predict the tumor growth kinetics under different treatment strengths including a dormant cell phenotype and underlined our theoretical approach with experimental data. Using our model which allows for phenotypic conversion, we could simulate how different tumor cell phenotypes proportionally contribute to the growing tumor mass in cycling treatment schedules. Additionally, we could observe that switching between high and low dosage treatment (with equal total treatment amounts) remarkably affects tumor growth in a frequency dependent way.
Thus, the dynamic proportions between cell phenotypes should be taken into account in the optimization of treatment schedules in order to control tumor growth.
Notes
Abbreviations
 D:

Dormant cells
 DMSO:

Dimethyl sulfoxide
 FCS:

Fetal calf serum
 GBM:

Glioblastoma multiforme
 P:

Rapidly proliferating cells
 SD:

Standard deviation
 TMZ:

Temozolomide
References
 1.
Ohgaki H, Kleihues P. Epidemiology and etiology of gliomas. Acta Neuropathol. 2005;109:93–108.
 2.
Stupp R, Mason WP, van den Bent MJ, Weller M, Fisher B, Taphoorn MJ, Belanger K, Brandes AA, Marosi C, Bogdahn U, Curschmann J, Janzer RC, Ludwin SK, Allgeier A, Lacombe D, Cairncross JG, Eisenhauer E, Mirimanoff RO. Radiotherapy plus concomitant and adjuvant temozolomide for glioblastomas. N Engl J Med. 2005;35:987–96.
 3.
Niewald M, Berdel C, Fleckenstein J, Licht N, Ketter R, Ruebe C. Toxicity after radiochemotherapy for glioblastoma using temozolomide  a retrospective evaluation. Rad Oncol. 2011;6:141.
 4.
Gatenby RA, Gillies RJ, Brown JS. Of cancer and cave fish. Nat Rev Cancer. 2011;11:237–8.
 5.
Almog N. Molecular mechanisms underlying tumor dormancy. Cancer Lett. 2010;294:139–46.
 6.
Wikman H, Vessella R, Pantel K. Cancer micrometastasis and tumour dormancy. APMIS. 2008;116:754–70.
 7.
Bragado P, Sosa MS, Keely P, Condeelis J, AguirreGhiso JA. Microenvironments dictating tumor cell dormancy. Recent Results Cancer Res. 2012;195:25–39.
 8.
Yeh AC, Ramaswamy S. Mechanisms of Cancer cell dormancyanother Hallmark of Cancer? Cancer Res. 2015;75:5014–22.
 9.
Zhang XH, Giuliano M, Trivedi MV, Schiff R, Osborne CK. Metastasis dormancy in estrogen receptorpositive breast cancer. Clin Cancer Res. 2013;19:6389–97.
 10.
Sosnoski DM, Norgard RJ, Grove CD, Foster SJ, Mastro AM. Dormancy and growth of metastatic breast cancer cells in a bonelike microenvironment. Clin Exp Metastasis. 2015;32:335–44.
 11.
Linde N, Fluegen G, AguirreGhiso JA. The relationship between dormant Cancer cells and their microenvironment. Adv Cancer Res. 2016;132:45–71.
 12.
Sertil AR. Hypoxia regulation in cellular dormancy. In: Hayat MA, editor. Tumor dormancy, quiescence and senescence, Vol. 3. Dordrecht: Springer Netherlands; 2014. p. 13–24.
 13.
Hirose Y, Berger MS, Pieper RO. p53 effects both the duration of G2/M arrest and the fate of temozolomidetreated human glioblastoma cells. Cancer Res. 2001;61:1957–63.
 14.
Altrock PM, Lin LL, Michor F. The mathematics of cancer: integrating quantitative models. Nat Rev Cancer. 2015;15:730–45.
 15.
Nowak MA. Evolutionary dynamics. Cambridge: Harvard University Press; 2006.
 16.
Hofbauer J, Sigmund K. Evolutionary games and population dynamics. Cambridge: Cambridge University Press; 1998.
 17.
Basanta D, Anderson ARA. Exploiting ecological principles to better understand cancer progression and treatment. Interface Focus. 2013;3:20130020.
 18.
Basanta D, Deutsch A. A game theoretical perspective on the somatic evolution of Cancer. In: Bellomo N, Angelis E, editors. Selected topics in Cancer modeling. Basel: Springer/ Birkhäuser; 2008. p. 97–112.
 19.
Dingli D, Chalub FACC, Santos FC, Van Segbroeck S, Pacheco JM. Cancer phenotype as the outcome of an evolutionary game between normal and malignant cells. Brit J Cancer. 2009;101:1130–6.
 20.
Kaznatcheev A, Scott JG, Basanta D. Edge effects in gametheoretic dynamics of spatially structured tumours. J R Soc Interface. 2015;12:20150154.
 21.
Orlando PA, Gatenby RA, Brown JS. Cancer treatment as a game: integrating evolutionary game theory into the optimal control of chemotherapy. Physical Biol. 2012;9:065007.
 22.
Bomze I, Bürger R. Stability by mutation in evolutionary games. Games Econ Behav. 1995;11:146–72.
 23.
Page KM, Nowak MA. Unifying evolutionary dynamics. J Theor Biol. 2002;219:93–8.
 24.
Almog N, Ma L, Raychowdhury R, Schwager C, Erber R, Short S, Hlatky L, Vajkoczy P, Huber PE, Folkman J, Abdollahi A. Transcriptional switch of dormant tumors to fastgrowing angiogenic phenotype. Cancer Res. 2009;69:836–44.
 25.
Adamski V, Hempelmann A, Flüh C, Lucius R, Synowitz M, Hattermann K, HeldFeindt J. Dormant human glioblastoma cells acquire stem cell characteristics and are differentially affected by Temozolomide and AT101 treatment strategies. Onctarget. 2017;8(64):108064–78.
 26.
Strogatz S. Nonlinear dynamics and Chaos: with applications to physics, biology, chemistry, and engineering (studies in nonlinearity). New York: Westview Press; 2000.
 27.
Endaya BB, Lam PY, Meedeniya AC, Neuzil J. Transcriptional profiling of dividing tumor cells detects intratumor heterogeneity linked to cell proliferation in a brain tumor model. Mol Oncol. 2016;10:126–37.
 28.
SatchiFainaro R, Ferber S, Segal E, Ma L, Dixit N, Ijaz A, Hlatky L, Abdollahi A, Almog N. Prospective identification of glioblastoma cells generating dormant tumors. PLoS One. 2012;7:e44395.
 29.
Naumov GN, Bender E, Zurakowski D, Kang SY, Sampson D, Flynn E, Watnick RS, Straume O, Akslen LA, Folkman J, Almog N. A model of human tumor dormancy: an angiogenic switch from the nonangiogenic phenotype. J Natl Cancer Inst. 2006;98:316–25.
 30.
Hofstetter CP, Burkhardt JK, Shin BJ, Gürsel DB, Mubita L, Gorrepati R, Brennan C, Holland EC, Boockvar JA. Protein phosphatase 2A mediates dormancy of glioblastoma multiformederived tumor stemlike cells during hypoxia. PLoS One. 2012;7:e30059.
 31.
Magnus N, Gerges N, Jabado N, Rak J. Coagulationrelated gene expression profile in glioblastoma is defined by molecular disease subtype. J Thromb Haemost. 2013;11:1197–200.
 32.
Magnus N, Garnier D, Meehan B, McGraw S, Lee TH, Caron M, Bourque G, Milsom C, Jabado N, Trasler J, Pawlinski R, Mackman N, Rak J. Tissue factor expression provokes escape from tumor dormancy and leads to genomic alterations. Proc Natl Acad Sci U S A. 2014;111:3544–9.
 33.
Magnus N, D'Asti E, Meehan B, Garnier D, Rak J. Oncogenes and the coagulation systemforces that modulate dormant and aggressive states in cancer. Thromb Res. 2014;133(Suppl 2):S1–9.
 34.
Evans EB, Lin SY. New insights into tumor dormancy: targeting DNA repair pathways. World J Clin Oncol. 2015;6:80–8.
 35.
Basanta D, Simon M, Hatzikirou H, Deutsch A. Evolutionary game theory elucidates the role of glycolysis in glioma progression and invasion. Cell Prolif. 2008;41:980–7.
 36.
Basanta D, Scott JG, Fishman MN, Ayala G, Hayward SW, Anderson ARA. Investigating prostate cancer tumour–stroma interactions: clinical and biological insights from an evolutionary game. Brit J Cancer. 2012;106:174–81.
 37.
Pacheco JM, Santos FC, Dingli D. The ecology of cancer from an evolutionary game theory perspective. Interface Focus. 2014;4:20140019.
 38.
Kaznatcheev A, Vander Velde R, Scott JG, Basanta D. Cancer treatment scheduling and dynamic heterogeneity in social dilemmas of tumour acidity and vasculature. Brit J Cancer. 2017;9:1–24.
 39.
Gerlee P, Anderson ARA. An evolutionary hybrid cellular automaton model of solid tumour growth. J Theoretical Biol. 2007;246:583–603.
 40.
Anderson ARA, Hassanein M, Branch KM, Lu J, Lobdell NA, Maier J, Basanta D, Weidow B, Narasanna A, Arteaga CL, Reynolds AB, Quaranta V, Estrada L, Weaver AM. Microenvironmental independence associated with tumor progression. Cancer Res. 2009;69:8797–806.
 41.
Zukewich J, Kurella V, Doebeli M, Hauert C. Consolidating birthdeath and deathbirth processes in structured populations. PLoS One. 2013;8:e54639.
 42.
Hindersin L, Traulsen A. Most undirected random graphs are amplifiers of selection for birthdeath dynamics, but suppressors of selection for deathbirth dynamics. PLoS Comput Biol. 2015;11:e1004437.
 43.
Wilkie KP, Hahnfeldt P. Tumorimmune dynamics regulated in the microenvironment inform the transient nature of immuneinduced tumor dormancy. Cancer Res. 2013;73:3534–44.
 44.
Hahnfeldt P, Folkman J, Hlatky L. Minimizing longterm tumor burden: the logic for metronomic chemotherapeutic dosing and its antiangiogenic basis. J Theoretical Biol. 2003;220:545–54.
 45.
Page K, Uhr J. Mathematical models of cancer dormancy. Leuk Lymphoma. 2005;46:313–27.
 46.
Kareva I. Escape from tumor dormancy and time to angiogenic switch as mitigated by tumorinduced stimulation of stroma. J Theor Biol. 2016;395:11–22.
 47.
Enderling H, Anderson ARA, Chaplain MAJ, Beheshti A, Hlatky L, Hahnfeldt P. Paradoxical dependencies of tumor dormancy and progression on basic cell kinetics. Cancer Res. 2009;69:8814–21.
 48.
Adamski V, Schmitt AD, Flüh C, Synowitz M, Hattermann K, HeldFeindt J. Isolation and characterization of fast migrating human glioma cells in the progression of malignant gliomas. Oncol Res. 2017;25:341–53.
 49.
Hattermann K, Flüh C, Engel D, Mehdorn HM, Synowitz M, Mentlein R, HeldFeindt J. Stem cell markers in glioma progression and recurrence. Int J Oncol. 2016;49:1899–910.
 50.
HeldFeindt J, Hattermann K, Müerköster SS, Wedderkopp H, KnerlichLukoschus F, Ungefroren H, Mehdorn HM, Mentlein R. CX3CR1 promotes recruitment of human gliomainfiltrating microglia/macrophages (GIMs). Exp Cell Res. 2010;316:1553–66.
 51.
HoppeSeyler K, Bossler F, Lohrey C, Bulkescher J, Rösl F, Jansen L, Mayer A, Vaupel P, Dürst M, HoppeSeyler F. Induction of dormancy in hypoxic human papillomaviruspositive cancer cells. Proc Natl Acad Sci U S A. 2017;114:E990–8.
 52.
Sosa MS, Parikh F, Maia AG, Estrada Y, Bosch A, Bragado P, Ekpin E, George A, Zheng Y, Lam HM, Morrissey C, Chung CY, Farias EF, Bernstein E, AguirreGhiso JA. NR2F1 controls tumour cell dormancy via SOX9 and RARβdriven quiescence programmes. Nat Commun. 2015;6:6170.
 53.
Ranganathan AC, Adam AP, AguirreGhiso JA. Opposing roles of mitogenic and stress signaling pathways in the induction of cancer dormancy. Cell Cycle. 2006;5:1799–807.
 54.
Schättler H, Ledzewicz U, Amini B. Dynamical properties of a minimally parameterized mathematical model for metronomic chemotherapy. J Math Biol. 2016;72:1255–80.
 55.
Dhawan A, Nichol D, Kinose F, Abazeed ME, Marusyk A, Haura EB, Scott JG. Collateral sensitivity networks reveal evolutionary instability and novel treatment strategies in ALK mutated nonsmall cell lung cancer. Sci Rep. 2017;7:1232.
Acknowledgements
We thank Judith Becker, Martina Burmester, Fereshteh Ebrahim and Brigitte Rehmke for expert technical assistance.
Funding
This work was supported by a sponsorship of the Medical Faculty of the University Kiel (“Forschungsförderung 2016” to JHF and KH) and the Deutsche Forschungsgemeinschaft RTG2154 [project 7 (KH) and project 8 (JHF)]. The funding body did not have any influence in the design of the study, collection, analysis and interpretation of data and in writing of the manuscript.
Availability of data and materials
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request. The python code to reproduce Figs. 1, 3, 4, and 5 is open source and available under https://github.com/marvinboe/TreatCycle.
Author information
Affiliations
Contributions
MB, JHF, AT and KH designed experimental and theoretical models and prepared the manuscript; MB and AT calculated and visualized the theoretical model; JHF, MS, RL and KH contributed to the experimental model; all authors read and approved the final manuscript.
Corresponding author
Correspondence to Kirsten Hattermann.
Ethics declarations
Ethics approval and consent to participate
Not applicable
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Böttcher, M.A., HeldFeindt, J., Synowitz, M. et al. Modeling treatmentdependent glioma growth including a dormant tumor cell subpopulation. BMC Cancer 18, 376 (2018). https://doi.org/10.1186/s1288501842811
Received:
Accepted:
Published:
Keywords
 Evolutionary game theory
 Glioma
 Dormancy