Quantifying the Role of Angiogenesis in Malignant Progression of Gliomas: In Silico Modeling Integrates Imaging and Histology

One of the key outstanding questions in the clinical management of gliomas lies in the unpredictable pattern of progression from low-grade, with its associated modest survival time measured in decades, to high-grade, with its aggressive growth and short survival time measured in months. According to the World Health Organization (WHO) grading scheme for gliomas, low and high grade are differentiated by features of cellularity, mitoses, and vascular proliferation. The characteristic vascular morphology of gliomas, specifically high-grade glioblastomas, has led many to the hypothesis that the formation of new blood vessels, angiogenesis, is essential to their growth (1). The angiogenic property of high-grade gliomas was first revealed in a study by Brem (2) in which a glioblastoma fragment was implanted into a rabbit cornea (a well-known angiogenesis assay). The endothelial cell mitogen and permeability factor VEGF was found to be present in the hypoxic pseudopalisading tumor cells lying adjacent to necrotic regions and hyperplastic vessels, identifying their role in glioma angiogenesis (3). Angiogenic features therefore are considered a key to the differentiation between low- and high-grade gliomas.

We have previously developed a mathematical model for quantifying the heterogeneity of glioma proliferation and invasion kinetics across patients. This proliferation–invasion model has been successful in predicting untreated growth and invasion kinetics for each patient (4), providing predictions related to surgical resection (5), chemotherapy (6, 7), and radiation therapy (8). Furthermore, the model allowed for direct estimation of the extent of diffuse glioma invasion (beyond that visible on MRI; ref. 4). However, this model neglects the complex interplay between tumor and its environment (e.g., the blood vessels and stroma) and, thus, does not incorporate the angiogenic cascade.

In the present article, we develop the proliferation–invasion model to produce a new model that incorporates distinct cellular and microenvironmental changes related to the angiogenic cascade. Specifically, we incorporate 3 distinct glioma cellular compartments (normoxic, hypoxic, and necrotic) with a vascular compartment and diffusible angiogenic factors (e.g., VEGF). This new model will allow for an investigation of malignant progression by studying the dynamics of tumor grade as a function of these cellular and microenvironmental compartments. In summary, by comparing imaging, histologic, and molecular characteristics of gliomas with in silico results, we identify predictable patterns of malignant progression, extent of hypoxic burden, and the size of the necrotic region in 3 glioblastoma patients. Overall, we find that the accumulation of genetic mutations is not necessarily required for malignant progression but rather interactions of tumor cells with a fixed net invasion rate (D) and net proliferation rate (ρ) alone can stimulate malignant progression seen on histology.

One of the major limiting factors in the usefulness of mathematical models lies within the ability to parameterize them. Multiple parameters are often difficult to obtain and may only be “guesstimates” at best; however, the strength of models is that one can examine a wide range of parameter space that incorporates the extremes of real-world values, thus allowing a full understanding of the possible model outcomes. Our previous work has shown a marked quantitative heterogeneity in the net rates of proliferation and invasion, even within each grade of gliomas [e.g., glioblastoma multiforme (GBM); ref. 4]. We can use this previous work to provide reasonable bounds for exploration for D and ρ. The rest of the model parameters are estimated from literature sources or calculated to approximate physiologic conditions (Supplementary Table S1).

Because of the diffuse nature of this disease, the in silico grade of each simulated tumor will vary spatially (as is the case in vivo). To obtain an appropriate summarized grade, analogous to the WHO grading scheme for astrocytomas (22), for each case (D, ρ, and time combination), we sample virtual biopsies at random locations within the tumor region visible on imaging (i.e., the simulated T2-weighted tumor region). Following our previous work (4, 17), the simulated T2 region corresponds to a simulated total cell density greater than the T2 threshold of detection (i.e., cell density greater than a threshold, approximated to be 16% of the maximal tissue carrying capacity k). Each virtual biopsy is the size of a typical microscopic high-power field [(HPF) 40× HPF ∼ 0.067 mm²]. We consider the in silico pathologic specimen to consist of 10 random HPFs to capture the typical clinicopathologic review protocol. For example, on the basis of the WHO criteria, if within those, there is a region of hypercellularity (i.e., cell density greater than a threshold) and more than 2 mitoses encountered (but not a significant amount of neovasculature or necrosis), then the lesion would be considered a grade III anaplastic glioma.

Figure 2 illustrates simulations of the PIHNA model varying only the net proliferation rate (ρ), revealing stereotypical differences between glioma grades. Simulations are represented as plots of the density or concentration of the model variables (normoxic, hypoxic, and necrotic tissue, vasculature, and VEGF) with respect to the distance from the center of the in silico tumor. Note the increasing cellularity, vasculature, hypoxia, and necrosis with increasing grade consistent with the currently used histopathologic grading scheme. The simulated pattern showing hypoxia as a wave at the leading edge of the hypercellular tumor was recently observed by Swanson and colleagues when comparing MRI observations of human glioblastoma with imaging of hypoxia with [18F]fluoromisonidazole positron emission tomography (18).

To explore the prognostic possibilities of knowing D and ρ, we created 225 virtual patients representing 15 × 15 combinations of D and ρ. The range of D and ρ was chosen to cover those observed in our experience with glioblastoma patients, for whom we have estimated these parameters from serial MRIs (4, 19, 20). We anticipate a wider range (specifically lower D and ρ values) to encompass a wider spectrum of low-grade gliomas. On this log–log scaled plot of ρ (y-axis) and D (x-axis), lines of constant velocity (for velocities of a form similar to the Fisher velocity of |$2\sqrt {D\irho }$|⁠) would be parallel diagonal lines. Assuming a fatal tumor burden (20) equivalent to a T2 radius of 4 cm, we calculated the in silico survival time, beginning at (probably) the smallest detectable size (T2 radius = 1 cm), as shown in Fig. 3. The pattern clearly is a diagonal one, related to the velocity of radial expansion, with short survival to be expected with high D/high ρ and long survival with low D/low ρ. Note, however, that as D gets large, the tumor is dominated by its diffuse extent so that it takes longer for the detection threshold to be met; however, once it is met, the in silico tumor immediately fills the virtual brain, albeit at a low cell density—similar to gliomatosis cerebri. Such a constant size for the fatal tumor burden may represent the average but does not seem realistic in view of the great variability of glioma size at autopsy (21). If we assume a fatal tumor burden equivalent to a total number of glioma cells (our current best estimate is about 1.1 × 10¹¹ cells; ref. 20), the pattern of survival is quite different (results not shown), more vertical, with short survival expected with high ρ and long survival with low ρ, but again, biased by the invisible invasive cells included distantly because of high D.

Motivated by the clinical paradigm, we considered the effect of both tumor size and time on the evolution of grade. For each pair of ρ and D (and tumor size or time point), as shown in Fig. 3 for mapping predicted survival times for a spectrum of in silico patients, we apply the in silico grading scheme (refer to Materials and Methods) to each simulated tumor to create Fig. 4A. Because of the diffuse nature of this disease, the in silico grade of each simulated tumor will vary spatially. We therefore, as before, took virtual biopsies at random locations within the tumor region to obtain an appropriate summarized grade (expressed as a probability of observing a given grade on virtual biopsy) for each case (D, ρ, and time/size combination) generating Fig. 4A:

1. As a function of tumor size, the progression wave of the in silico grading map appears to descend obliquely (diagonally) from the top left to the bottom right. Thus, for small tumors to attain the histologic features of GBM, there must be both low D and high ρ; however, for larger tumors, there is a wider range of D and ρ that will generate the features of GBM.

2. As a function of time, the progression wave of the in silico grading map quickly spreads obliquely (diagonally) from high ρ/low D (top left corner) but then more vertically, such that practically all are grade IV. Again, this is evidence of the biasing effect of high D and the invisible low densities of invading glioma cells.

The color gradient between each grade is a reflection of the probability that a given grade is observed on virtual biopsy. Movies comparing the progression of these 2 simulated patients (cases * and +) are available in Supplementary Material.

Figure 4 shows the time course of the radial expansion of the simulated abnormality seen on T2 MRI for 2 different choices of the model parameters—a high ρ and low D case (case *; Fig. 4B) and a lower ρ and higher D case (case +; Fig. 4C). These 2 simulated cases differ in their time course for appearance on imaging, as well as their in silico grade at the time of first detection (assumed here to be 1 cm in radius on T2 MRI). In fact, these 2 cases represent in silico “primary” and “secondary” glioblastomas (pGBM and sGBM, respectively) in the sense that from the time of detectability on T2 MRI (assumed to be 1 cm in radius), case * maintained a grade IV in silico grade, whereas case + was considered a lower grade until it was over 2 cm in radius on T2 MRI (22). To further explore this, we illustrate the portion of parameter space in which pGBMs and sGBMs live in Fig. 5, denoting the location of the case * and + examples.

As reviewed by Harpold and colleagues (4), our studies with the proliferation–invasion model have culminated in a characterization of pretreatment glioma kinetics with a net proliferation rate (ρ) and a net invasion rate (D). An important consequence of the prior proliferation–invasion model is an asymptotically constant rate of radial expansion visible on imaging (23). This constant (untreated) radial growth has been confirmed in both low-grade (24) and high-grade (19) gliomas. Furthermore, this rate of radial expansion, which is mathematically related to the model parameters as |$2\sqrt {D\irho }$|⁠, is prognostically significant in gliomas (8, 17, 20, 25, 26). We explored this in the simulations of the expanded PIHNA model (c.f., Fig. 4B and C) and found a similar pattern of approaching a constant speed of radial expansion (dependent predominantly upon the D and ρ).

By randomly sampling from the regions of D–ρ space containing in silico pGBMs and sGBMs, we noticed several similarities to published observations. Specifically, we confirmed the in vivo observation that pGBMs tend to have increased expression of VEGF when compared with grade III and sGBMs (22, 27), as shown in Fig. 6A. Regions of the D–ρ space (as shown in Fig. 5) were categorized as pGBM, sGBM, or grade III depending on their in silico grade status at rT2 = 3 cm. Then, from each category, 10 simulations were selected at random, and the ratio of total VEGF to total tumor cells (obtained by integrating concentrations over the whole domain) was calculated. The error bars represent the full range of data in each category, whereas the gray boxes indicate the 25th to 75th percentile. Similarly, with increasing grade, an increasing number of these virtual biopsy cells express VEGF or hypoxia-inducible factor (HIF)-1α consistent with the in vivo observations of Korkolopoulou and colleagues (28) and Flynn and colleagues (ref. 29; Fig. 6B).

Using patient-specific estimates of net proliferation rate (ρ) and net invasion rate (D; refs. 8, 25, 30, 31), it is possible to assess the validity of patient-specific PIHNA simulations bridging the elements interpretable through the model such as histology, immunohistochemistry (e.g., as in Fig. 6), and other imaging features. Specifically, we considered 3 patients who had a similar size on T2 MRI of approximately 2 cm but who had different hypoxic burdens (percentage of the cells on histologic specimens expressing HIF-1α), different necrotic core sizes seen on gadolinium-enhanced T1-weighted MRI (T1Gd MRI), and different D–ρ combinations. Figure 7 shows the following: (i) a map of histologic grade, (ii) a map of the maximum of the hypoxic cell population (h), and (iii) the size of the necrotic core (n, assuming a cutoff value of 5% of the maximum cell density defining the edge of the necrotic core) all at simulated time point for which the simulated T2 radius is 2 cm. Figure 7 shows a direct mapping of parameters of net rates of proliferation (ρ) and invasion (D) for 3 glioblastoma patients (stars), allowing a direct correlation between model-predicted and actual observations of histologic grade (Fig. 7A) and hypoxia measured histologically by HIF-1α expression (Fig. 7B). Moreover, those same model simulations also accurately predict the volume of necrosis seen in these patients as assessed by regions of hypointensity on contrast-enhanced T1-weighted MRI (Fig. 7C). All patients were accurately predicted to be of histologic grade IV (red in Fig. 7A). Patient 1 with the lowest ρ was accurately predicted as having a smaller hypoxic fraction (30%) and smaller necrotic radius (8 mm) than others. Patients 2 and 3 were predicted to have similar levels of hypoxia (red/orange on Fig. 7B), which was observed on histology (40% HIF-1α⁺ in both cases). Patient 2 had a 16-mm necrotic radius on T1Gd MRI, whereas patient 3 had a slightly larger 18-mm necrotic core, which is consistent with the model predictions seen in Fig. 7C. In summary, these patient-specific simulations are compared with patient-specific information from imaging and histology, and we find that the simulations provide consistent patterns between these 5 distinct data points: overall T1Gd radius, T2 radius, tumor grade, central necrosis visible T1Gd MRI, and HIF-1α immunohistochemistry. Furthermore, Gu and colleagues (32) recently reported, using the PIHNA model, to generate patient-specific simulated [18F]fluoromisonidazole positron emission tomography images of hypoxia.

“Histologic grading is a means of predicting the biologic behavior of a neoplasm” (22). This definition consists of a challenge, namely, to explain what is assumed to be true, that histology has some innate ability to explain why it can, if it indeed does, predict the future. For the past 15 years, we have been exploiting a biomathematical model of gliomas (ref. 23; the proliferation invasion model) based on only 2 key cellular characteristics, net rates of proliferation (ρ) and of diffusion/invasion (D) of glioma cells, revealing predictive insights related to surgical resection (5), chemotherapy (6, 7), and radiation therapy.

By directly extending the proliferation invasion model, we produced the PIHNA model, which encompasses the angiogenic cascade and allowed us to define tumor grade as a combination of both cellular and microenvironmental model–specific variables. Critically, this allowed us to summarize intuitively a complete parameter study using a graphical representation of changing grade. This highlighted a disparity between the current in vivo grading schemes and observed tumor behavior and suggested that malignant “progression” may be less related to acquisition of new genetic abnormalities than a natural consequence of a relatively small set of initial key molecular genetic changes that “turn on” the glioma cells to proliferate and invade. That is, the net rates of proliferation (ρ) and diffusion (D) of any particular glioma may ultimately be what drive the in silico behavior of glioma (including the whole spectrum from slow to fast and from localized to diffuse) and histologic appearances (including increasing degrees of cellularity, mitoses, hypoxia-driven neoangiogenesis, and necrosis), provided the microenvironment responds appropriately.

The goal of a good grading scheme is to capture prognosis implied by cellular behavior; however, the in silico grade maps (Fig. 4) suggest that the WHO grading scheme, independent of size and time, is preferentially biased toward assessing proliferative capacity, although Fig. 4 suggests that invasiveness has a significant impact on survival as well. A key difference between the grade maps in terms of size and time is in the ability to capture the impact of invasion rates on overall tumor growth kinetics. One would presume that rapid growth (as a result of a combination of invasion and proliferation) would be associated with shorter survival times. To illustrate this, we map the time from a typical diagnostic size (say 1 cm in radius on T2) to a typical fatal size (say 4 cm in radius on T2) to approximate a survival map. A comparison of this survival map with the grade map (for the same initial size of a T2 radius of 1 cm) shows that the shortest survival times are for those in silico tumors with both high D and ρ, whereas the WHO grading scheme applied to the simulated tumors would suggest that only high ρ is important. Furthermore, these model results suggest that exploration of invasion-specific markers [e.g., matrix metalloproteinase 2 (ref. 33) and vimentin (ref. 34)] may enhance the prognostic accuracy of the current grading scheme. Admittedly, assessing invasion from fixed histologic sections of dead tissue is a true challenge. However, this challenge may be mediated by dynamic modeling tools like that presented here.

Obviously, histopathologic grading is based on a static snapshot of the tumor tissue at a certain size (and time in the tumor's natural history). Applying our in silico WHO grading scheme, a virtual patient presenting with a 1-cm tumor with a relatively high D and low ρ (case + in Fig. 4), would appear histologically as a grade II. However, according to our model, following this patient in time would show a rapid progression to GBM, whereas a patient with a relatively high ρ and low D (case *) would have an overall slower growth rate but would have histologic characteristics of a GBM from diagnosis. This is consistent with clinical observation showing wide variability in survival and growth patterns even among histologically similar tumors and suggests that morphologically based grading is naturally limited (22). Such cases of malignant progression, as illustrated in Fig. 4C (case +), are common among low-grade (grade II or III) gliomas but are currently considered relatively unpredictable events associated with dramatic decrease in median survival time (22).

Figure 6 provides preliminary confirmation that the PIHNA model maps directly to, and connects with, clinical metrics of malignant progression, medical imaging, and histology in individual patients. These data suggest that patient-specific simulations of the PIHNA model may be able to recapitulate and, in fact, predict characteristics of each patient's tumor onto an orthogonal data set that is not obtainable directly from the original patient imaging used as input to the model. Although a similar analysis on a broader patient set is needed to confirm these preliminary findings, these data support an integrative role of the PIHNA model in bridging imaging, histology, and immunohistochemistry in individual patients.

Traditional molecular genetic characterization of gliomas shows a marked heterogeneity within and across grades (22). We propose that our model parameters may be simply consequences (downstream) of these molecular signatures. Our model simulation results suggest that one can think of each glioma as having an implicit invasive and proliferative potential, which manifests quantitatively in terms of our model parameters D and ρ. As the cells grow and invade according to their preprogrammed D–ρ combination, the local microenvironment may become harsh (e.g., hypoxic), which may lead to the production of angiogenic factors that induce vascular growth (via HIF-1α and similar pathways). Other than those cells within the hypoxic tumor becoming less proliferative and the vascular endothelial cells beginning to proliferate, the kinetics of the cell population (c) driving the growth of the tumor overall have not changed (i.e., D and ρ remain the same), whereas depending on the size and time of observation, the morphologic features (seen on histologic specimen) can be strikingly different. These results are consistent with recent findings (35), suggesting that microenvironment and/or tumor–stromal interactions are the key to malignant progression and changes in cell kinetics may not be necessary to see these pathologic changes.

It is important to note that this conclusion does not preclude changes in cellular behavior as a result of gene expression changes (caused by any number of mechanisms including environmental influences as well as mutational events), nor does it exclude the possibility that early in glioma development (and perhaps posttreatment), the tumor consisted of a more heterogeneous population with a range of D and ρ. Rather, these results simply illustrate that changes in the dominant cellular kinetics (e.g., D and ρ) are not necessary to obtain the range of morphologic patterns seen in vivo within individual patients. This means that the wildly different morphologic changes on imaging and histopathology observed in some patients through their disease course may not actually be a result of an underlying change in the dominant tumor cellular phenotype (D and ρ) but simply a manifestation of the evolution of that cellular phenotype interacting with its microenvironment.

We realize, of course, that we are proposing a simplified model for the dynamic interactions between the tumor cells and the microenvironment and therefore miss several key characteristics that can be readily explored in future models. Overall, our proposed minimal PIHNA model for glioma–microenvironment interactions, which explicitly incorporates the interactions of normoxic glioma cells, hypoxic glioma cells, vascular endothelial cells, diffusible angiogenic factors, and the formation of necrosis, is sufficient to describe the spectrum of in vivo dynamics of gliomas visualized with medical imaging and may be an invaluable research analysis and clinical assessment tool. An important area of future exploration for this model is to understand the changing tumor microenvironment in the context of treatment. This is particularly relevant to increasing our understanding of the effects of antiangiogenics on disease progression.

No potential conflicts of interest were disclosed.

K.R. Swanson thanks the generous and timely support of the McDonnell Foundation, the Dana Foundation, the Academic Pathology Fund, NIH/NINDS R01 NS60752, and the NIH/NCI Moffitt-UW Physical Sciences Oncology Center U54 CA143970.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.