## Abstract

Gliomas are well known for their potential for aggressive proliferation as well as their diffuse invasion of the normal-appearing parenchyma peripheral to the bulk lesion. This review presents a history of the use of mathematical modeling in the study of the proliferative-invasive growth of gliomas, illustrating the progress made in understanding the in vivo dynamics of invasion and proliferation of tumor cells. Mathematical modeling is based on a sequence of observation, speculation, development of hypotheses to be tested, and comparisons between theory and reality. These mathematical investigations, iteratively compared with experimental and clinical work, demonstrate the essential relationship between experimental and theoretical approaches. Together, these efforts have extended our knowledge and insight into in vivo brain tumor growth dynamics that should enhance current diagnoses and treatments.

## Introduction to Invasive Brain Tumors

Gliomas, the most common primary brain tumors, are thought to arise from the supporting glial cells of the brain or their precursors (1). Because they generally grow and invade extensively before the patient notes any symptoms, gliomas are almost impossible to cure. Glioblastomas are the most malignant and most common gliomas in adults, accounting for about 50% of all gliomas. Glioblastomas are distinguished by necrosis, which may be massive centrally or more irregular between vascular supplies, and peripheral cells that diffusely invade the surrounding tissue. The aggressive behavior of these tumors is reflected in their 100% fatality rate within approximately one year even after extensive surgery, radiotherapy, and chemotherapy.

Despite continual advances in imaging technology, glioma cells invade far beyond the abnormality shown on clinical imaging (e.g. CT, MRI, or PET) and even beyond gross and microscopic observations at autopsy. The extent is certainly beyond that guiding present-day radiotherapy of gliomas, which targets therapy to only an arbitrary 2 cm beyond the imaged bulk mass of the tumor. Clearly, it is the invasion into the normal-appearing surrounding tissue that is responsible for the tumor recurrence even in those tumors that are radio-sensitive.

## Discovering the In Vivo Dynamics of Cancer Cells, Specifically Glioma Cells

In the mid-1900s, major experimental and human studies converged on the idea that cancer cells divided at rates that, although varying widely among different cancers, were relatively constant over long periods of time in individual cases. Among the pioneers in the study of human cancers were Collins et al, who observed that metastases in the lungs visualized by plain x-rays of the chest grew at constant volume-doubling rates according to a simple exponential "law" (2). This formulation evoked a tremendous response in investigative human oncology, which has been extensively reviewed by Araujo and McElwain (3). However, this "law" fails a priori when one considers that an infiltrative primary neoplasm, such as a glioma, has an unknowable fraction that invisibly invades adjacent normal-appearing tissue, below the threshold of detection, leaving a correspondingly unknowable fraction that is visible and can be measured.

The failure of the exponential law was not immediately obvious as various relatively static models were developed (4-11), but recognition of the dynamic invasiveness of gliomas soon required inclusion of cellular motility in addition to proliferative growth. Unfortunately, without such a biomathematical model as an hypothesis, decades of potential serial observations without intervening treatment were lost in the early clinical therapy trials when treatments were generally ineffective and the only known way to evaluate the effectiveness of therapy was in double-blind comparisons involving large numbers of randomized patients, rather than a more detailed study of individual patients.

The pioneer in the theoretical analysis of the relation between cell kinetics and growth of the gross tumor was Steel (12). Steel's work made it quickly apparent that there was an order of magnitude difference between the times involved in the definitions of cellular and gross kinetics: hours to a few days for individual cells, many days and even months for gross tumors. His work formulated several equations. Each equation, although relatively simple arithmetically, involved a large number of unknowns: T_{p} (the potential doubling time), T_{d} (the actual gross volume-doubling time), T_{c} (the cell cycle time), T_{s} (the time during which DNA is being synthesized within T_{c}), γ (a factor defining the timing of T_{s} within T_{c} and the age distribution of the cells), LI (the labeling index), ϕ (the cell loss factor, the proportion of newly formed cells that are lost), and GF (the growth fraction, the proportion of the total cells that are actually cycling).

To these unknowns, Alvord and Shaw added the amount of tumor (mostly invisible) left behind after surgical resection and developed a nomogram (Fig. 1) that approximated the interrelations between most of Steel's unknowns (13). By taking advantage of the simple arithmetic calculations that two 36.5 doublings of 10 μm cells would produce about 100 g of tumor and that 36.5 doublings every 10 days would require 365 days or 1 year to produce that 100 g of tumor, they were able to include survival time. Survival time could either be predicted or, more generally, be used retrospectively to estimate how much tumor that was all or mostly invisible to the neurosurgeon or neuroradiologist that had been left behind at the time of treatment. Based on what little evidence was present at the time, Wilson estimated that the gradient of glioblastoma cells might be approximately exponential, with the concentration decreasing to approximately 10% every 2 cm (14).

## Early Glioma Model Development: Growth and Diffusion

With this initial research on glioma growth dynamics, the foundation for a mathematical model began to take form. The formulation of the model was characterized by a series of unique model developments, each continually transformed and modified from the previous one. With each successive model came an intensified understanding of glioma growth at both the micro and macro level. The deficiency of each proposed model provided insight into the necessary improvements for the subsequent proposals.

In the early 1990s, Murray's group defined the basic spatio-temporal model, based on the classical definition of cancer, as uncontrolled proliferation of cells with the capacity to invade and metastasize (15). This model was simplified by taking advantage of the fact that gliomas practically never metastasize outside the brain, producing a conservation-diffusion equation, written in words in Equation 1.

Under the assumption of classical gradient-driven Fickian diffusion, this word equation could be quantified mathematically to produce Equation 2.

Equation 2 describes the dynamics of glioma cells where c(**x**,t) is the concentration of tumor cells at location **x** and time t. D is the diffusion coefficient representing the net motility of glioma cells and ρ represents the net proliferation rate of the glioma cells. The ▽ term is the spatial differentiation operator, which is in effect a gradient. Initial conditions for the model were c(x,0) = f(x), where f(x) defined the initial spatial distribution of malignant cells, presumably a point source at the center of tumorigenesis.

## Applications of Glioma Model to Clinical Cases

An extension of this mathematical formulation including two cell sub-populations was suitable for analyzing the first clinical case. Figure 2 illustrates the analysis of the growth of a recurrent malignant astrocytoma, treated and followed during the patient's last year with serial contrast-enhanced CTs (15). Optimizing the fit between these observations and the model provided the first estimates of the model parameters D and ρ; D = 10^{−2} cm^{2}/day (with 7%-9% "type 2" cells, defined as resistant to the first chemotherapy) (15) or 10^{−3} cm^{2}/day (no "type 2" cells) (16) and ρ = 10^{−2}/day.

Woodward et al speculated that these values of D and ρ derived for this one patient might be close to the average for all high-grade gliomas, with ±50% of each possibly encompassing the whole range of such tumors (16). This led to the development, at the level of CT resolution, of a two-dimensional model in homogeneous tissue bounded by the ventricles centrally and the skull peripherally. The model imposed no migration of cells beyond these brain boundaries. The tumor was proposed as spherically symmetrical and appearing grossly to be 3 cm in diameter at diagnosis and 6 cm in diameter at death. These sizes were averages derived from very wide ranges reported in the few publications available at the time (6,17,18), and similar to those in subsequent reports (19-21).

## Effects of Resection

Because gliomas are routinely operated upon with resection of "all" or variable portions of the tumor visible on imaging, creation of a mathematical model that could accurately mimic postoperative glioma dynamics was essential. Such a model was achieved by setting the cell concentration in the resection site to zero and allowing the surrounding tumor to continue to grow and diffuse at the same rates as before resection, continuing to the same grossly visible average fatal diameter of 6 cm. At this stage the model, crude as it was, was sufficiently elegant to be compared with actual glioma growth statistics. Kreth et al had reported survival curves for two groups of patients, one biopsied only and one subjected to "gross total resection," as defined by the neurosurgeon at operation (22). Both groups also received radiation therapy. The 7-week improvement by "gross total resection" over biopsy was predicted mathematically (16) and found clinically (22), but even 115 patients were not sufficient to prove this difference to be statistically significant.

## Low-Grade Gliomas

Up to this point, model development had focused on describing the growth of high-grade gliomas because of their rapid advancement and ability to be viewed on CT and other imaging modalities. In order to increase the applicability of the model to all gliomas, Woodward et al proposed that low-grade gliomas might grow one tenth as fast as the high-grade gliomas (again, with ±50% each of D and ρ to encompass the whole range) (16). This was tested against all of the available published results (Fig. 3A) (23), and established that in silico the model was able to predict all possible actual survivals (Fig. 3B).

## The Detectable Edge as a Traveling Wave

One relationship derived from equation 2 is the Fisher approximation, which equates the velocity v of the detectable tumor margin to twice the square root of the product ρD as follows:

or D = *v*^{2}/4ρ. This approximation was originally derived from the observation that a population governed by growth and diffusion alone expands at the constant velocity

for large time, thereby expanding linearly for any given ρ and D.

Mandonnet et al were able to show that low-grade gliomas did indeed grow both slowly and linearly (24). Figure 4 illustrates that in their first 27 patients it appeared that the average velocity of the diameter was about 4 mm/year. We have found it more convenient to speak of radial expansion, which would be half of this, or approximately 2 mm/year. Concurrently, the constancy of velocity, albeit much faster, was being confirmed in a single rare patient with a glioblastoma that was followed with repeated MRIs for a year without intervening treatment (25). This case produced a radial velocity of 12 mm/year, six times that of low-grade gliomas but close to the 10 times that previously predicted by Woodward et al (16).

This comparison with actual patient data led to the concept of linear radial velocity of the traveling wave of the detectable "edge" of all gliomas. A very different geometric formulation by Mayneord (26) led to the same concept of linear radial growth of Jensen's sarcoma, a particular rat tumor, but this tumor was characterized as having a markedly necrotic center with the viable tumor cells growing essentially in two dimensions as a thin surface of viable tumor cells surrounding the necrotic core. With solid but infiltrating gliomas, however, the linear radial growth results in a cubic growth of its visible volume. While all of the cells can potentially be proliferating exponentially, the radius of the visible bulk is increasing linearly and its volume is growing cubically. The major difference between cubic and exponential growth is that the "volume-doubling time" is not constant. Instead, the apparent volume-doubling time decreases progressively, as though any treatment is being effective, or as though the tumor is obeying Gompertzian or logistic growth and outgrowing its blood supply. Furthermore, consecutive volumes of an infiltrating glioma cannot be converted to a meaningful volume-doubling time according to the classical exponential hypothesis. Conversely, calculated "apparent volume-doubling times" cannot be converted to linear radial velocities if the initial volume is not specified for each case. If this conclusion is true for neoplasms generally, a lot of work of the last 50 years on growth rates of potentially infiltrating cancers has become irretrievably damaged.

## Quantifying in vitro Behavior

In addition to the analysis of in vivo patient data, the modeling approach has been successfully used to quantify in vitro studies of glioma invasiveness. Palfi et al studied the invasiveness in rodent brain slice assays of 42 low-and high-grade human glioma specimens obtained at biopsy (27). The most significant result of this analysis was the linking of phenotypic behavior of the glioma specimens with their genotypic expression: those glioma cells with 1p and 19q loss were less invasive than their 1p/19q-positive counterparts. This correlates with the general observation that in vivo low-grade gliomas with 1p and 19q loss generally have a better prognosis, related to therapy (28) or not (29). Of interest also, if still unexplained, is the observation that most such tumors occur in non-temporal sites (30,31), especially bifrontal (30).

## Brain Heterogeneity: Isotropic Migration

The original analyses of the mathematical model assumed homogenous brain tissue so that the diffusion coefficient D, defining random motility of glioma cells, was constant and uniform throughout the brain (15,16,32). Recognizing that the model had to be improved to accommodate the advances in MRI technology that were coming along in parallel, Swanson et al reformulated the model to accept different diffusion rates in grey and white matter (33). This modified model introduced the complex geometry of the brain and presented diffusion (motility) as a function of the spatial variable **x** to accommodate the observation that glioma cells demonstrate greater motility in white matter than in grey matter (33). The original word equation, Equation 1, continued to apply, but the mathematics changed to involve a spatially varying diffusion parameter, D(**x**), as shown in Equation 3.

D(**x**) is still defined as the diffusion coefficient defining the net motility of the glioma cells but with D(**x**) = D_{G} or D_{W}, different constants for **x** in grey and white matter, respectively.

In order to accurately analyze the dynamics of the model in heterogeneous brain tissue, a detailed description of grey and white matter distribution in the brain was required. Fortunately, the neuro-anatomical atlas from the BrainWeb database was available, providing a spatial distribution of grey and white matter for the entire brain at a resolution of 1 mm^{3} voxels (34,35). Swanson et al applied resection to this model and found that the spatial heterogeneity of grey and white matter contributed to the patterns of recurrence following treatment, ranging from the typical observed ring recurrence at the boundary of the resection bed to the less common multifocal recurrence pattern (36).

## Brain Heterogeneity: Anisotropic Migration

Glioma cells are commonly thought not only to migrate more quickly in white matter but also to migrate preferentially along blood vessels and fiber tracts in the white matter of the brain. Magnetic resonance diffusion tensor imaging (MR DTI) has recently emerged as a powerful tool for analyzing the 3D geometry and producing an in vivo reconstruction of significant white matter tracts. Jbabdi et al incorporated MR DTI with the model's Equation 3 augmented to treat the spatial diffusion coefficient D(**x**) as a tensor **D**(**x**) (37). They simulated the extended model for the case of low-grade gliomas located within significantly anisotropic portions of the brain, such as the uncinate fasciculus between the frontal and temporal lobes (37). Such a model could easily be applied to other paralimbic low-grade gliomas invading the arcuate (superior longitudinal) fasciculus and spreading into the external capsule and lateral temporal lobe (38).

## Imaging with Limited Thresholds of Detection

In addition to the CT observations of Lewander et al (39) and Greene et al (40), Kelly et al (21,41) and Dalrymple et al (20) reported that the histologic "edge" of the "solid tumor" coincided approximately with the circumference of the tumor visualized by enhanced CT or MRI (T1-Gd) and that "isolated tumor cells" could be seen microscopically inside and less commonly outside of the circumference of the tumor visualized in the T2 image. Indeed, tumors can be cultured from normal-appearing tissue 4 cm from the edge of a glioblastoma (42). Further consideration of the possibility of a relationship between specific MRI sequences (particularly T1-Gd and T2) suggested a measure of the gradient of cellular invasion beyond the visible "edge" of enhancing malignant gliomas.

Figure 5 illustrates the T1 and T2 detection thresholds and their failure to image the entirety of the invasive glioma cells. By setting the T1-Gd concentration at 80% of maximum and T2 concentration at 2%, Swanson was able to approximate a gradient of glioma cells that corresponded to the ratio D/ρ (43). Varying ρ and D with D/ρ fixed allows the geometry of the tumor growth and invasion to stay the same, while the time scale on which the growth and invasion occur changes. For example, two different tumors with D/ρ fixed could appear exactly the same at a single time point but reach their ultimately fatal size in very different lengths of time.

The availability of 70 patients with glioblastomas and pre-operative volumes of T1-Gd and T2 images (Prof. Robert Rostomily, University of Washington, personal communication) provided an opportunity to test the limitations of these imaging modalities. Analysis of these patients improved the model by replacing the average 3 cm diameter at diagnosis with individual sizes and "gradients" at diagnosis (Fig. 6). Unfortunately, none of the patients had a second MRI before treatment, allowing no way to obtain individual velocities in this series of patients. There has also been no advance in the definition of the "fatal tumor burden," 6 cm on average, but with a very wide range (17, 18). The cause of death is obviously multi-factorial (44) and not every patient dies of the tumor itself. Nevertheless, about two thirds of all patients show herniation at death, indicating that increasing mass is a major cause of death in most patients and suggesting that there may be a reasonably constant "fatal tumor burden" for most patients, but this still remains to be proven.

## Creating Virtual Control Patients

Despite these challenges, using an average velocity allowed D and ρ to be calculated for each of the 70 patients and, with the same criterion of 6 cm diameter at death, allowed calculation of the predicted survival of each patient following any extent of resection (biopsy, subtotal resection, or gross total resection proven by post-operative enhanced CT). Although the prognosis of individual patients could not be accurately predicted, the median survivals and the shapes of the curves in silico were quite satisfactorily close to the actual, and it was strongly suspected that these would improve if real velocities could have been determined for each patient.

Figure 7A gives a comparison between the model prediction and the actual patient data, plotting survival time against percentage of survivors. It is important to note that the virtual untreated controls differ for each group of patients, being modeled as individual patients with exactly the same size of tumor (both T1-Gd and T2) at the time of diagnosis but not treated. It is as though each treated patient had an identical twin that could have been treated differently. The result is that the duration of survival of the matched virtual controls for the gross total resection (GTR) and biopsy/subtotal resection (BX/STR) groups differed by 12 weeks (Fig. 7B), which is about half the difference of the 25 weeks between the actual GTR and (unmatched) actual BX/STR groups. GTR could not be defined more precisely clinically except as having resected all detectable tumor by postoperative enhanced CT. GTR, of course, could be precisely defined mathematically and we predicted that the real GTR would probably be between 100% and 125% resection, the model removing either exactly the entire hyper-intensity on T1-Gd scan or a 25% margin surrounding the T1-Gd image. Indeed, the prediction proved correct.

Notably, all of these developments have used routine clinical MRIs (5-mm slices with T1-Gd and T2 sequences). Further advances can easily be incorporated, increasing the resolution, but current technology is already adequate for many practical purposes.

## Brain Heterogeneity: Modeling Chemotherapy

Chemotherapy would seem to be the most likely way to reach all tumor cells, both the local bulk and the diffusely invading cells. However, the drugs must be able to penetrate the normal blood-brain barrier. Even then, the heterogeneity of the vascular density within the grey and white matter can affect the results. Swanson et al showed that, although the total number of tumor cells within the brain may be decreasing with chemotherapy, the extent of invasion of the tumor remains practically unaffected due to the continuing motility of the tumor cells within the white matter (36). Even though the tumor appears to be regressing on MRI, extensively invaded tumor cells remain occult, below the detection abilities of the MRI, primarily throughout the white matter. Once treatment is stopped, diffuse recurrences seem inevitable. This suggests a potential difficulty with the design of clinical trials relying solely on MRI data as a measure of success of treatment.

## Visualizing The Relationships of D and ρ

Building on the mathematical model of the behavior of gliomas based solely on their net rates of proliferation (ρ) and diffusion (D), one can construct a log-log graph of D vs ρ (Fig. 8B). One of the advantages of a log-log graph is that both the products and ratios of D and ρ can be expressed as straight lines. In Figure 8B the diagonal between the point D = 100 mm^{2}/year and ρ = 1/year and the point D = 10 mm^{2}/year and ρ = 10/year represents all of the points where the product Dρ = 100 mm^{2}/year^{2} and the square root of the product Dρ = 10 mm/year. Thus, using the Fisher approximation of the velocity being twice this square root, v = 20 mm/year. One can then add parallel lines representing any desired velocities (Fig. 8B).

Similarly, the other diagonal between the point D = 10 mm^{2}/year and ρ = 1/year and the point D = 100 mm^{2}/year and ρ = 10/years represents all of the points where D/ρ = 10 mm^{2} (Fig. 8A). Parallel lines can be added representing any desired ratios of D/ρ (Fig. 8A). Swanson has made computerized simulations of glioma growth revealing that at the average time of diagnosis, when the visible radius is 15 mm, the ratios D/ρ of 1, 10, and 100 mm^{2} represent the invisible portion to be about 4%, 50%, and 87% of the total tumor cells (43). Thus, we can speak of the ratio D/ρ as an approximation of the "Invisibility Index."

These graphs could, of course, have been constructed on a simple linear scale (Fig. 8D), but the products would have been curved (hyperbolic) and the ratios diverging straight lines, so that the relationships would have been obscured rather than revealed.

One is now in a position to see where any particular glioma fits in this scheme of things, where all of the lines are superimposed, using a 3-cycle log-log format to include all of the values so far found for low-and high-grade gliomas (Fig. 9). Low-grade gliomas are in the bottom left, their velocities averaging 2 mm/year (24); and high-grade gliomas are in the top-right, enclosed within a rectangle defined by D/ρ of 2 to 20 mm^{2} and velocities of 10 to 200 mm/year. Gliomatosis without detectable mass would be in the lower right corner, above which would come gliomatosis with small masses in grey matter (e.g. "hypertrophy" of the pons or thalamus). In the upper left would appear those gliomas that are curable by resection (e.g. pilocytic astrocytomas). In between these extremes are the other gliomas, with "slow growing," "rapidly growing," "slightly infiltrative," "markedly infiltrative," etc. being favorite expressions of most grading schemes. However, a glance at Figure 9 should convince almost anyone that the overlapping of the two major behavioral characteristics at right angles to each other makes a single-digit grading scheme unlikely to be successful.

Although the ratio of the volumes of T1-Gd and T2 images is not the correct calculation, larger volume ratios would generally correlate with larger values of D/ρ and, therefore, with more diffuse tumors that would be less susceptible to GTR (45). That such tumors are more likely to be EGFR-amplified (45) with a poorer prognosis (46) are interesting observations, but it would be even more interesting to recalculate the volume ratios to D/ρ, the prediction being that EGFR-amplified gliomas should have higher than average D/ρ.

Now that the specific behavioral characteristics can be defined for each case, we believe that the time is ripe for the correlation of morphologic features of individual cases, not just the average or range within groups of cases. Beyond this, one can look forward to defining when and why contrast enhancement occurs in MRIs and the probable concomitant neo-angiogenesis induced by hypoxia/hypoglycemia resulting from local hypercellularity as a "low-grade" glioma "progresses" to a higher grade. There may be an increase in velocity-an unfavorable ratio of D/ρ-which in the "natural history" of that glioma leads to one or more foci of greater metabolic demand than the native blood vessels can provide. The increase in vascularity may or may not be adequate, with inadequacy leading to necrosis, a hallmark of glioblastoma. But here too, an unfavorable D/ρ allows the range of velocities to be quite marked, so that patients may survive only a few months or even a few years without the neuropathologist being able to recognize the difference in D and/or ρ characteristic of each case. There may be some mathematical relationship between proliferation (e.g. MIB-1) and apoptosis (e.g. caspase-3), but it is not likely to be a simple subtraction or division in view of the different time courses of the processes involved.

## Summary

We have presented an historical account of the iterative comparisons of theory and reality, which have allowed the progressive improvement of a relatively simple bio-mathematical model of the visible and invisible expansion of an infiltrating neoplasm, specifically a glioma. The spatio-temporal model began at the relatively crude level of resolution available on CT scans in homogeneous tissue bounded by the ventricles and skull. It has advanced to accommodate heterogeneous tissue with differences in grey and white matter with an anatomic accuracy of 1 mm^{3} and can use current routinely available MRIs to evaluate the effects of specific treatments in individual patients. The present model has the novel, if not unique, capability to distinguish different mobilities in grey and white matter, with an estimated factor of 5 but an expected range of 10 to 100. Mobility in white matter is not likely to be isotropic, but more likely to be differentiated into specific tracts along which malignant cells may migrate even more easily. This has recently been modeled by anisotropic diffusion (37). This grey-white novelty is currently useful in two areas: 1) the theoretical construction of virtual examples that mimic real MRIs, and 2) the comparison with autopsy specimens. The virtual examples provide insight into the extent of the invisible infiltrating portion of gliomas, especially in explaining why they are essentially incurable regardless of the extent of resection and why chemotherapy must enter through normal blood vessels well in advance of the infiltrating tumor cells. As resolution increases in routine clinical MRIs, these virtual examples will undoubtedly begin to be approached in vivo.

The progress of mathematical modeling in the field of human oncology has greatly expanded from the major contributions of Collins et al in 1956 (2) and Steel in 1977 (12). This article has discussed only gliomas, but there is certainly a significant overlap with current mathematical modeling efforts concerning other cancers. With the inevitable and continuous advancements in imaging, it is clear that the progress in modeling will continue to transform our understanding of in vivo tumor dynamics. We believe that increased understanding of tumor growth dynamics will lead to improvements in the diagnosis and treatment of these diseases.