Single Cell Mathematical Model Successfully Replicates Key Features of GBM: Go-Or-Grow Is Not Necessary

Glioblastoma (GBM) is a malignant brain tumor that continues to be associated with neurological morbidity and poor survival times. Brain invasion is a fundamental property of malignant glioma cells. The Go-or-Grow (GoG) phenotype proposes that cancer cell motility and proliferation are mutually exclusive. Here, we construct and apply a single glioma cell mathematical model that includes motility and angiogenesis and lacks the GoG phenotype. Simulations replicate key features of GBM including its multilayer structure (i.e.edema, enhancement, and necrosis), its progression patterns associated with bevacizumab treatment, and replicate the survival times of GBM treated or untreated with bevacizumab. These results suggest that the GoG phenotype is not a necessary property for the formation of the multilayer structure, recurrence patterns, and the poor survival times of patients diagnosed with GBM.


Introduction
Glioblastoma is a malignant brain tumor that causes high morbidity and poor survival times. Brain invasion, motility, and rapid proliferation are characteristic features of malignant glioma cells. The hypothesis that cancer cell motility and proliferation are mutually exclusive first emerged in the work of Giese et al., who studied cell migratory and proliferative response to extracellular matrix proteins in vitro. Their findings suggested a dichotomy between the two behaviors [1]. Hatzikirou et al. later coined the phrase "Go-or-Grow" (GoG) to describe this dichotomy, suggesting that this phenomenon best explains the transition of malignant tumors to the invasive phenotypes in the presence of hypoxia. To test this hypothesis, they apply a lattice-gas cellular automata model, which we hereafter refer to as the Dresden model [2]. Finally, Dhruv et al. recently reported that the coordinated suppression and activation of certain transcription factors may explain the shift in glioma cells from the proliferating to migrating phenotype [3].
In previous work, we have developed a mathematical model of Glioblastoma multiforme (GBM), which not only replicates the known multi-layer structure of GBM, i.e.necrosis, enhancing ring and edema, but also reproduces specific tumor progression patterns and their a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 associated average patient survival times [4,5]. This model incorporates the GoG phenotype by including two specific cell types, Invasive (I) cells and Proliferative (P) cells, which can switch back and forth between one phenotype or another depending on the nutrient environment in the brain. As suggested by their names, I cells move throughout the brain but do not divide, and P cells divide but do not move; under the assumption of the GoG phenotype, a single tumor cell will never move and divide at the same time. The equations in [4,5] also include an angiogenic term.
Our previous work uncovered motility as a possible predictor of GBM progression patterns. Unique to our model is a hypoxia-driven motility term, which is different from diffusion, or concentration-driven dispersion, in that it allows tumor cells to react to low-oxygen environments with increased motility. Through titrations of these motility terms, we have identified three distinct patterns of motility in GBM tumors: (1) highly dispersive, (2) moderately dispersive, and (3) hypoxia driven. These motility phenotypes are associated with different survival times as well as the way a GBM reacts to anti-angiogenic (AA) treatment. The math model includes parameters that model fundamental properties of the tumor including the motility rate, which we found to be a distinguishing factor in the progression types (PT) [6]. We have identified these three types as Expanding FLAIR, Expanding FLAIR + Necrosis, and Expanding Necrosis.
Baker et al. recently developed an agent-based model of GBM, simulating perivascular brain tumor growth and invasion, that does not include the GoG phenotype [7]. In their model, which simulates the displacement of each cell, cancer cells in contact with blood vessels move slower due to adhesion and divide at a fixed mitotic rate. Experimental data from Baker et al. supports the idea that some glioma cells retain the ability to both divide and migrate on blood vessels [7]. In an effort to better understand the GoG phenotype, here, we develop a new mathematical model at the scale of MRI which includes all the features of our two-cell model (angiogenesis, hypoxia-driven motility, concentration-driven motility, and proliferation) with the exception of the GoG phenotype. Rather than include two specific cell types, one that divides and one that moves, our new model includes one cell type, Glioma cells (G), which can both move and divide at the same time.
In this investigation, our goals were to use the new single glioma cell mathematical model to (1) replicate the multilayer structure of GBM for untreated tumors in all three motility phenotypes of GBM (highly dispersive, moderately dispersive and hypoxia-driven), (2) replicate the three progression patterns associated with bevacizumab-treated tumors (Expanding FLAIR, Expanding FLAIR + Necrosis, Expanding Necrosis), and (3) replicate the survival times associated with these three progression patterns. Our results suggest that the GoG phenotype may not, in fact, be necessary.
In the sections that follow, we first present the equations governing the single glioma cell model, explaining the key differences between this model and our original two-cell model. We then show the results of simulations using our new model, including a simulated clinical trial. We conclude by proposing that the GoG phenotype is not an absolute necessity for the formation of the multilayer structure of GBM and its recurrence patterns.

Mathematical Model
and 3), the brain tissue is taken to be homogenous so that the rate of diffusion is constant throughout the brain. The new system of equations includes a single PDE for glioma cells (see Table 1 and Fig 1b). These cells can multiply and migrate using two modes of motility: concentration-driven (passive transport) and hypoxia-driven (active transport) motility. A description of the major differences in these two modes of motility can be found in [4,5]. Most notably, passive diffusion is driven by glioma cell concentration, and active transport is driven by hypoxia, or low nutrient conditions, which varies inversely with total cell concentration. Passive diffusion is blind to hypoxia, whereas active transport causes cancer cells to move in bulk away from necrosis and into healthy brain tissue.
By reducing the two glioma cell model to a single glioma cell model, we were able to eliminate two parameters: the rate of transition from P cells to I cells during hypoxia and the rate of transition from I cells to P cells during normoxia. All other parameter values used in the new model were either identical or within a similar range as those values used in the previous GoG model ( Table 2). Fig 1 shows an interactive cell diagram comparing the two models.
The system of equations for the two-cell GoG model is also reprinted below (Table 3). Most notably, our previous model includes two thresholds, one for hypoxia (C hyp ) and one for death (C ltm ). These thresholds are functions of total cell concentration. The difference in the two thresholds (C ltm − C hyp = F) created a transition period for the switch of P cells to I cells, and Assumptions regarding angiogenesis, hypoxia, mitosis, and necrosis: Rate of necrosis : Necrotic threshold : where σ = 1.0 to simulate angiogenesis, and σ = 0 to simulate anti-angiogenesis.
A description of the parameters and their values/units may be found in Table 2.

Numerical Methods
The original two-cell mathematical model is discussed in detail in [4,5]. Specific parameter choices, including those differentiating the three motility phenotypes, are also detailed in Table 2 and [4,5]. These papers also review the numerical methods used to solve the system of equations, which is identical to the methods used to simulate results for the new single glioma cell model. The computations were performed at the Alabama Supercomputer Authority (www.asc.edu).

The Necrotic and Angiogenic Tumor Threshold
The initial brain concentration in our model is 1.0. The initial necrotic threshold, C ltm , is 1.1, which means the tumor mass has reached ten percent of the initial brain concentration, the maximum amount of additional tissue the brain can sustain in our model without further vascularization. Without angiogenesis, the tumor and brain tissue will begin to die once the total cell concentration reaches the initial necrotic threshold. The tumor cells also stop dividing at this point. With angiogenesis, C ltm begins at the initial hypoxia threshold and increases as a Assumptions regarding angiogenesis, hypoxia, mitosis, and necrosis:

Measure of Local Hypoxia for P=I Conversion
Mitotic Rate : Rate of necrosis : Hypoxic threshold : where σ = 1.5 to simulate angiogenesis, and σ = 0 to simulate anti-angiogenesis.
Necrotic threshold : function of glioma cells (G) at a fixed angiogenic rate (per hour), which may vary from tumor to tumor. When the total cell concentration reaches C ltm , division slows and cells begin to die. Without angiogenesis, a primary characteristic of GBM, tumor cells cannot exceed a concentration of 10% of the brain (i.e. total cell concentration cannot exceed the initial necrotic threshold of 1.1). Hence, one of our two methods of detection for GBM is the presence of total cell concentrations above the threshold for the detection of vascular proliferation, C vas . This threshold is set at 1.12, meaning local vascularization has allowed tumors cells to reach a concentration of 12% or more of the brain, which exceeds the initial death threshold by 20%. The plots in Fig 2 show (a) the initial (in the absence of angiogenesis) rate of mitosis and onset of cell death as a function of total cell concentration and C ltm and (b) the effect of tumor angiogenesis on this rate and death threshold, which allows total cell concentrations to exceed C vas .

Definitions of FLAIR, Vascular Proliferation, and Necrosis
Note that a single pixel is approximately 3 mm 2 , which is at the resolution of MRI. The following describes definitions for Angiogenic Tumor, Necrosis, and Fluid Attenuated Inversion Recovery (FLAIR) in our simulations: 1. FLAIR (C > 1.003): Our definition of FLAIR is meant to capture any area of the brain that has been invaded by a small concentration of cancer cells. We assume that cell concentrations above 1.003, that is when tumor cells exceed 0.3 percent of the initial brain concentration, generate high signal in FLAIR MRI sequences. > 1.12): Vascular Proliferation is defined as any area (or pixel) of the brain where the total concentration of cells exceeds the threshold for detection of Vascular Proliferation (C vas ) of 1.12. We assume that at this point, the tumor cell density causes a disruption in and increased permeability of the blood-brain barrier, resulting in contrast enhancement on MRI.

Radiological Necrosis (B < 0.30):
Radiological necrosis (that which can be detected at the resolution of an MRI) is defined as any area (or pixel) of the brain where 70% or more of the brain cells in that area are dead.

Pathological Necrosis (B < 0.9995):
Pathological necrosis (that which can be detected under a microscope following tumor resection) is defined as any area (or pixel) of the brain where 0.05% or more of the brain cells in that area are dead.

Simulated Patient Diagnosis and Death: Clinical Trials
In this investigation, we monitor tumor growth, FLAIR, gadolinium enhancement and necrosis; the overall survival time is taken to be the difference in the time of death and the time of diagnosis. We define the time of diagnosis to be the time in the simulation when either (1) the area of tumor with vascular proliferation or (2) the area of radiological necrosis reaches a certain critical threshold (see Table 4). Likewise, the time of death is defined as the time when the area of (1) tumor with vascular proliferation, (2) FLAIR, or (3) radiological necrosis in the brain reaches a critical size. As in the time of diagnosis, the time of death is triggered when the first of these three thresholds is reached. As in our previous work with the GoG model, we use the single-cell model to simulate clinical trials with "patients" exhibiting the three identified tumor progression patterns: (1) Expanding FLAIR (PP1), (2) Expanding FLAIR + Necrosis (PP2), and (3) Expanding Necrosis (PP3). We also run simulations of untreated tumors in an attempt to replicate the known average survival time of untreated GBM. We selected five titrations of each diagnostic and death criteria to produce a total of 25 "patients" in each of the three tumor progression categories. Once one of the five diagnostic criteria was reached, we simulated anti-angiogenic tumor treatment by setting the angiogenic parameter σ to zero. We also let each of these simulations run untreated for a total of 75 untreated tumor simulations. Table 4 displays the titrations used to simulated death and diagnosis in these in silico clinical trials.

Results
In the absence of the GoG phenotype, the single-cell mathematical model was designed to replicate: 1) the multilayer structure of GBM in untreated tumors (vascular proliferation along with a necrotic core and a proliferating cancer cell ring), 2) the three tumor progression types under anti-angiogenic treatment, and 3) the known survival times associated with each progression type. In our previous two-cell model, which incorporates the GoG phenotype, we achieved all three goals by varying motility alone and fixing all other parameters. Specifically, highly dispersive tumors progressed by Expanding FLAIR, moderately dispersive tumors progressed by Expanding FLAIR and Necrosis, and hypoxia-drive tumors progressed by Expanding Necrosis. In this model, we achieve the same results by varying motility alone but without the GoG phenotype. Table 5 summarizes the conclusions of our simulations, while also providing a comparative analysis to the two-cell model and the other mathematical models referenced in Table 6.

Single-Cell model Replicates Progression Patterns of Bevacizumabtreated GBM
Our simulations also show that variations in the motility phenotypes can affect tumor response to anti-angiogenic therapy (see Fig 3(d)-3(f)). As in the two-cell model, highly-dispersive tumors (those governed by a high concentration-driven motility parameter) progress by Expanding FLAIR (Fig 3(b)). Treatment effectively reduces the proliferating tumor mass, halting further vascular proliferation; however, the tumor continues to spread at low densities throughout the brain, eventually killing the patient when 65% of the brain has been invaded by FLAIR (green arrow). For a moderate concentration-driven motility parameter or moderately-dispersive tumor, simulated treatment results in progression by Expanding FLAIR + Expanding Necrosis (see [5]). Fig 3(e) displays this progression pattern. When this moderately-dispersive tumor is treated (first columns), there is a decrease in the proliferating tumor mass, or gadolinium enhancement, but necrosis and FLAIR continue to expand, as evidenced by the growing hole in the brain (pink arrow) and the presence of brain invasion or FLAIR (green arrow) beyond the site of necrosis.
In the presence of very low concentration-driven motility, a high parameter choice for hypoxia-driven motility generates a GBM tumor that progresses by Expanding Necrosis. This treated hypoxia-driven tumor results in an aggressively expanding area of necrosis, as indicated by the pink arrow in Fig 3(f). Note that both the treated and untreated hypoxia-driven tumor simulations (Fig 3c and 3f), glioma cells remain in close proximity to the site of necrosis, which is consistent with the magnetic resonance imaging indicators distinguishing the Expanding Necrosis progression pattern from the Expanding FLAIR + Necrosis progression pattern [5,6]. The median survival times for each group were close to those found in [6] and [4]. For our computational trial, we found Log-Rank p < 0.0001 for comparisons between all trial  1% of the brain with detectable vascular proliferation (purple arrows) or radiological necrosis (pink arrows) served as triggers for detection of each GBM. In all three cases, the tumors progressed to the appearance of a necrotic core (pink arrows) surrounded by a proliferating ring (orange arrows). Death was triggered by either 20% tumor mass or 5% radiological necrosis. Treatment of these tumors reproduced progression by Expanding FLAIR (d), Expanding FLAIR + Necrosis (e), and Expanding Necrosis (f). For the treated tumors, the first time shot (treatment) is taken immediately prior to anti-angiogenesis treatment, the second time shot shows the 2-month or 3-month follow-up, and the final time shot displays tumor appearance at the simulated time of death. In progression by Expanding FLAIR (d), treatment effectively eliminates the spread of tumor with vascular proliferation (purple arrows) while low-density tumor cells, or FLAIR (green arrows) continue to invade the brain. Both (e) and (f) also show a reduction in the spread of vascular tumor. However, for moderately dispersive tumors (e), treatment results in progression by Expanding FLAIR (green arrow) + Expanding Necrosis (pink arrow). Treated hypoxia-driven tumors (f) progress by Expanding Necrosis alone (pink arrow), where the area of lowdensity invasive tumor hovers just beyond the periphery of the necrotic core.

Single-Cell model Replicates Survival Times of GBM
doi:10.1371/journal.pone.0169434.g003 Table 6 lists a few well-known mathematical models and summarizes their specific differences; only some incorporate the GoG phenotype. The Swanson-PI (proliferation-invasion) model was among the first models of GBM to appear in the literature and includes a single cell phenotype that can both proliferate and diffuse throughout the brain [8]; however, this model does not include angiogenesis or hyopxia-driven motility. Of those listed, note that our two-cell model is one of only two models that incorporate the GoG phenotype. The Swanson-PIHNA (proliferation-invasion-hypoxia-necrosis-angiogenesis) model also includes two distinct cell phenotypes (normoxic and hypoxic) [9]. However, unlike the Dresden or our Two-Cell GoG models, both cell phenotypes in the Swanson-PIHNA model have the ability to diffuse throughout the brain, and the normoxic phenotype can both proliferate and diffuse. Our twocell and single-cell models are the only models that include a hypoxia-driven motility term, which is different from diffusion. Both our models also include an angiogenic term. In summary, the single-cell model reproduces all the same biological milestones of GBM as the two-cell model (Multilayer Structure, Progression Patterns, and Survival Times), as summarized in Table 5, suggesting that GoG phenotype may not be a necessity in GBM. These results corroborate the findings of Beker et al. [7]. We conclude by proposing the hypothesis that the presence or absence of the GoG phenotype may be cell-specific and/or a function of the local environment.