Abstract
Glioblastoma, the most aggressive form of primary brain tumor, is predominantly assessed with gadolinium-enhanced T1-weighted (T1Gd) and T2-weighted magnetic resonance imaging (MRI). Pixel intensity enhancement on the T1Gd image is understood to correspond to the gadolinium contrast agent leaking from the tumor-induced neovasculature, while hyperintensity on the T2/FLAIR images corresponds with edema and infiltrated tumor cells. None of these modalities directly show tumor cells; rather, they capture abnormalities in the microenvironment caused by the presence of tumor cells. Thus, assessing disease response after treatments impacting the microenvironment remains challenging through the obscuring lens of MR imaging. Anti-angiogenic therapies have been used in the treatment of gliomas with spurious results ranging from no apparent response to significant imaging improvement with the potential for extremely diffuse patterns of tumor recurrence on imaging and autopsy. Anti-angiogenic treatment normalizes the vasculature, effectively decreasing vessel permeability and thus reducing tumor-induced edema, drastically altering T2-weighted MRI. We extend a previously developed mathematical model of glioma growth to explicitly incorporate edema formation allowing us to directly characterize and potentially predict the effects of anti-angiogenics on imageable tumor growth. A comparison of simulated glioma growth and imaging enhancement with and without bevacizumab supports the current understanding that anti-angiogenic treatment can serve as a surrogate for steroids and the clinically driven hypothesis that anti-angiogenic treatment may not have any significant effect on the growth dynamics of the overall tumor cell populations. However, the simulations do illustrate a potentially large impact on the level of edematous extracellular fluid, and thus on what would be imageable on T2/FLAIR MR. Additionally, by evaluating virtual tumors with varying growth kinetics, we see tumors with lower proliferation rates will have the most reduction in swelling from such treatments.
Introduction
Glioblastoma Multiforme (GBM) is a highly aggressive and invasive primary brain tumor. The standard treatment protocol is to surgically remove as much of the tumor as is reasonably safe, followed by a combination of chemotherapy with radiation. Despite aggressive treatment, the prognosis remains poor with a median survival time of 14 months (Stupp et al., ). The inability to accurately determine the extent of diffuse tumor cell infiltration of the normal brain affects the ability to assess response to treatment through clinical imaging, confounding clinical progress. Currently, clinicians rely primarily on three magnetic resonance imaging (MRI) modalities to monitor the development of the tumor, the T2 weighting, FLAIR, and T1 weighting with gadolinium contrast enhancement (T1Gd) sequences illustrated in Figure 1. However, it is known that none of these sequences are able to show the entire extent of the malignant cells (Silbergeld and Chicoine, ), since the abnormal regions highlighted in the MR images are as dependent on the microenvironment around the disease, particularly the vasculature, as on the tumor cells themselves.
Figure 1
Increased vasculature is a primary hallmark of GBM, and while angiogenesis is a hallmark of cancer in general, there are some important factors that separate GBMs from other tumors. First, the glioma cells inhabit an organ that is highly vascularized in its native state. Second, glioma cells are diffusely invasive and are known to co-opt the existing vasculature and migrate and grow along the vessels (Holash et al., ; Leenders et al., ). Nevertheless, GBMs can form hypoxic regions, often leading to regions of necrosis, and thus, downstream of this hypoxic signaling, emit an abnormally large amount of angiogenic factors such as vascular endothelial growth factor (VEGF) for the recruitment of additional vasculature, analogous to observations in solid tumors (Kerbel, ). This process results in the vasculature developing abnormally large vessel radial sizes and, unique to the brain, results in a breakdown of the blood brain barrier in the tumor region.
Features characterizing MR imaging observation
In the case of GBM, the enhancing abnormalities on all of the primary MR imaging modalities, T1Gd, T2, and FLAIR primarily result from a compromised blood brain barrier. The T1Gd image signal is enhanced where the contrast agent has been able to leak into the extracellular space through breakdowns in the blood brain barrier due to tumor-induced neo-angiogenesis. The T2-weighted and FLAIR images show edema or swelling; FLAIR is different from T2 in that the signal from the cerebral spinal fluid (CSF) is inverted. In the case of GBM, the non-CSF T2/FLAIR hyperintense signal is primarily vasogenic edema, defined as fluid originating from blood vessels that accumulates around cells (Marmarou, ). The fluid leaves the vessels due to pressure and osmotic gradients induced by the breakdown in the blood brain barrier. Thus, whilst changes in any primary imaging modality (T1Gd or T2/FLAIR) are often interpreted as corresponding to changes in tumor cell density, they may be artifacts of MR imaging.
The role of anti-angiogenics in GBM
The concept of anti-angiogenic treatment for cancer has been popular ever since the landmark paper by Folkman () stating that malignant tumors were angiogenesis-dependent and has been used with some success for other solid tumors in combination with chemotherapy (Hurwitz et al., ; Sandler et al., ). Since a defining hallmark of GBM is increased vasculature through endothelial cell proliferation (Louis et al., ), this disease seems like an obvious candidate for vascular targeting treatment. However, the differences between the vasculature in GBMs and other solid tumors produce different treatment effects. In GBMs, one of the effects of anti-angiogenic treatments is to, at least transiently, repair the blood brain barrier and allow the vessels to return to their normal radial size, increasing their efficiency (Jain, ; Batchelor et al., ) – referred to as vascular normalization. Ostensibly, this improved efficiency of the vasculature is not the desired impact, though it may help in delivery of other therapeutic agents.
More concerning, however, is that this normalization may directly impact the efficacy of the MR imaging. It is possible for glioma patients with enhancing lesions on T1Gd and T2/FLAIR imaging to have decreased enhancement within a day of anti-angiogenic treatment (Batchelor et al., ; Norden et al., ), as illustrated by patients 1 and 2 in Figure 2, but upon stopping treatment, the imageable lesion is even larger and more disperse than before (Iwamoto et al., ). However, responses are varied; a patient may see no deflection in growth but faster growth after treatment, such as the third patient in Figure 2, or see stabilized disease returning to the previous growth patterns after treatment, as illustrated by the fourth patient in Figure 2. These patients were consented to this study with approval by the local institutional review board at either the University of Washington or the University of California, Los Angeles, and their relevant demographic and therapeutic information is given in Table 1. These conundrums have led to two hypotheses: first, anti-angiogenic treatment has minimal cytotoxic effect but does influence the imaging so that the tumor cannot be effectively visualized, and second, that the treatment may be selecting for a more aggressively invasive phenotype (Verhoeff et al., 2009; Keunen et al., ).
Figure 2
Table 1
| Age | Sex | Grade | Race | XRT dose (cGy) | Concurrent TMZ | Bev given at recurrence | Concurrent therapies with Bev | |
|---|---|---|---|---|---|---|---|---|
| Patient 1 | 48 | M | III | Caucasian | Given, but unknown dosage | Y | Y | Irinotecan, dexamethasone |
| Patient 2 | 55 | M | IV | Caucasian | 6000 | Y | Y | Carboplatin |
| Patient 3 | 61 | M | IV | Unknown | 6000 | Y | Y | Irinotecan |
| Patient 4 | 66 | F | IV | Caucasian | 6120 | Y | Y | Irinotecan |
Demographic and treatment information corresponding to patients in Figure 2.
Previous studies have shown the use of anti-angiogenic drugs, specifically bevacizumab (Avastin), tends to increase progression-free survival and reduce symptoms of recurrent GBMs, but they have failed to consistently show a significant increase in overall survival and there is concern that the measured radiographic responses do not reflect changes in tumor cell counts (Verhoeff et al., 2009; Deming, ). Additionally, animal studies have revealed that treatment with anti-angiogenic drugs may be creating an environment favorable for local invasion and metastasis (Ebos et al., ; Pàez-Ribes et al., ). Even though there is a lack of significant evidence for increased overall survival after bevacizumab, and it is possible that treatment selects for a more invasive phenotype, the increase in quality of life for some of the patients, due to the relief from edema-related symptoms, means anti-angiogenic therapy is an attractive and relevant treatment option. However, the inability to determine a priori which patients will receive more benefit than harm from anti-angiogenic therapy ultimately keeps clinicians wary (Deming, ).
In this paper, we aim to illustrate how a previous mathematical model of glioma growth can be extended to explicitly incorporate edema formation allowing us to directly characterize and potentially predict the effects of anti-angiogenics on imageable tumor growth. The ultimate goal of this model is to help the treatment planning process by identifying exactly those patients that would receive the most benefit from anti-angiogenic treatment.
Materials and Methods
The proliferation-invasion-hypoxia-necrosis-angiogenesis-edema model
Over the last decade we have made a significant effort toward the development of patient specific mathematical models of GBM that are able to capture the growth kinetics of individual patients (Swanson, ; Swanson et al., , ,; Harpold et al., ; Rockne et al., ). The simplest form of the model, referred to as the Proliferation-Invasion (PI) model is based on patient specific net rates of proliferation and invasion and has been successful in predicting untreated growth rates for individual patients (Harpold et al., ) and providing predictions of outcomes following surgical resections (Swanson et al., 2003), chemotherapy (Swanson et al., ,, 2003), and radiation (Rockne et al., ), while also providing insight into glioma ontogeny (Bohman et al., ).
A mathematical model of the angiogenic cascade in glioblastoma
The Proliferation-Invasion-Hypoxia-Necrosis-Angiogenesis (PIHNA) model first discussed in Swanson et al. (2011) incorporates the angiogenic cascade and characterizes malignant gliomas with relative proportions of well-oxygenated “normoxic” tumor cells, (c), poorly oxygenated hypoxic tumor cells, (h), necrotic cells, (n), and vascular, or endothelial cells, (v), along with a generic population of angiogenic factors, (a) (Swanson et al., 2011). In words, it assumes the level of nutrients present in the local microenvironment, as inferred from the number of vasculature cells, determines whether the present tumor cells will exhibit normoxic or hypoxic phenotypes. That is, if there is a sufficient level of nutrients present, the cells will remain normoxic, but if the nutrient level falls below a given threshold, the cells will become hypoxic. If the nutrients provided by the vasculature fall below an even lower threshold value, the hypoxic cells will undergo necrosis, at a rate of αh (1/year) and remain in the necrotic cell population. Normoxic tumor cells are allowed to move (invade) and divide while, due to restricted amounts of nutrients, the hypoxic cells are only allowed to move. The hypoxic cells produce a large amount of angiogenic factors which ultimately cause an increase in the number of vasculature cells. The system is described with a mathematical model composed of the following five coupled reaction-diffusion equations:
In these equations, D(x) is the net rate of invasion (mm2/year). Glioma cells migrate faster along myelinated axons in the white matter than in the dense and less structured cortical gray matter. For this reason, we consider the net rate of invasion as piecewise constant, with non-zero values in the gray and white matter, Dg and Dw, respectively, with Dw > Dg, and zero in the regions of cerebral spinal fluid. Additionally, ρ (1/year) is the net proliferation rate of the normoxic cells, γ (1/year) and β (1/year) are the maximum conversion rates between the hypoxic and normoxic cell populations, αn (1/year) is the rate at which cells undergo necrosis when in contact with necrotic cells (contact necrosis), αh (1/year) is the rate of conversion of hypoxic cells to necrotic cells when nutrient levels fall too low, Dv (mm2/year) is the rate of dispersal of vasculature cells, estimated from Sherratt and Murray (), Levine et al. () μ (1/year) is the vasculature proliferation rate, estimated from Xiu et al. (2006), T = (c + h + n + v)/K (dimensionless), where K is the carrying capacity (cells/mm3), and V = v/(v + c + h) (dimensionless) and is a surrogate for the local vasculature efficiency. Angiogenic factors are produced by both normoxic and hypoxic cells with rates δc (1/year) and δh (1/year) respectively, with δh > δc and are consumed by the vasculature for both regular vasculature maintenance [with rate ω (1/year)] and for vasculature proliferation [with rate q (1/year)]. Finally, the angiogenic factors are assumed to decay over time with rate λ (1/year) and disperse with rate Dv (mm2/year). Values for parameters related to the angiogenic factors were derived in part from work done in Levine et al. (), Serini et al. (), Mac Gabhann and Popel (). The reader is referred to (Swanson et al., 2011) for further details.
It is known that GBM tumors are extremely genetically heterogeneous both within a single tumor and between different tumors (Dunn et al., ). A large effort has been put forth to identify subtypes of GBMs by their dominating genotype (Verhaak et al., 2010). While the model parameters do not directly try to capture effects of single mutations, it is our belief that the dominating genotypes characterizing subtypes of GBMs ultimately result in different net rates of proliferation and invasion which would be used in our model. For example, pro-neural tumors are more likely to have the IDH-1 mutation and be secondary GBMs. In our model, this would manifest as a low-D, low-ρ tumor which begins as low grade and progresses into higher grade/malignancy.
In this model, there are different cell populations, normoxic, hypoxic, and necrotic competing for space and each with differing phenotypes: normoxic cells proliferating and invading, hypoxic cells only invading, and necrotic cells which are dead and just taking up space. Since each cell population is evolving in space and time, there is an effective spatial heterogeneity of predicted proliferative activity across space and time which could be analogized to heterogeneous Ki67 labeling across glioma specimens. Thus, while this model is attempting to capture the overarching phenotype of different tumors and assumes global constants for individual tumors, spatial heterogeneity in behavior is possible due to regional levels of vasculature which may result in hypoxia and/or necrosis. We also remark the aim of this model is not to predict cell-level behaviors (ex. Ki67), rather, use information obtained from routine imaging to quantify and explain imaging scale behavior and evolution.
The PIHNA model only captures cellular species and angiogenic factors. While these all have an impact on what is ultimately seen on MR imaging, in and of themselves they are not sufficient to describe enhancing regions of T2 and T1Gd MR images. Here we extend the PIHNA model to capture the imaging responses post anti-angiogenic treatment. To achieve this, we add to the model one additional element: edema, (l), to create the merged proliferation-invasion-hypoxia-necrosis-angiogenesis-edema (PIHNA-E) model. A schematic of the six species interactions is shown in Figure 3.
Figure 3
A mathematical model for tumor-induced edema formation
Generically, edema refers to a swelling phenomenon. While there are different types of cerebral edema, in the case of GBM, it is almost exclusively vasogenic edema which results from fluid and protein leakage from the breakdown of the blood brain barrier (Marmarou, ). Over the last few decades, there have been quite a few attempts to model vasogenic edema (Rapoport, ; Kumagai, ; Nagashima et al., ) and its associated interstitial pressure and interstitial fluid velocity (Baxter and Jain, , ). These models were primarily based on Starling’s equation which describes fluid exchange between compartments due to pressure and osmotic gradients. These models are very detailed and are generally solved on shorter time scales, i.e., a few days versus months. In our efforts, while we do have interest in the specific mechanisms, we will take a coarser grained approach allowing us to approximate the phenomena over longer time scales relevant to tumor growth kinetics.
To begin, we make the simplifying assumption that the edema is only composed of fluid which has leaked into the extra-cellular space and has not yet been reabsorbed into the system. This fluid is assumed to leak into the extracellular space where the blood brain barrier has been compromised. From the PIHNA model, we can approximate these regions along with the degree of permeability from the local levels of present angiogenic factors. Once the fluid is in the extra-cellular space it moves via diffusion and is reabsorbed into the system at a constant rate. This process is written in the form of a partial differential equation as:
Here l is the concentration of edematous fluid, Dl (mm2/year) is the diffusion rate of the edematous liquid which would be analogous to an Apparent Diffusion Coefficient, ADC, value derived from diffusion-weighted MRI (Moritani, ), lv (fluid/mm3) is the normal level of fluid in the vasculature, δl (1/year) is the reabsorption rate, and Ktrans (1/year) is the transmission rate [analogous to the value Ktrans measured on dynamic contrast enhanced MRI (DCE-MRI)], capturing the permeability surface area product per unit volume of tissue (Tofts, 1991) and is assumed to depend on the level of angiogenic factors, a, present. Homogeneous Neumann boundary conditions are assumed at the boundary of the brain to ensure no fluid leaves the brain.
The dependence of the Ktrans coefficient on the angiogenic factors is assumed to take a Michaelis–Menten type I form: to reflect that VEGF (also known as vascular permeability factor, VPF) strongly influences vascular permeability (Bates, ). Here Kmax (1/year) is the maximum possible value of Ktrans, the value of which is calibrated to Grade IV gliomas (Patankar et al., ), and Khalf (angiogenic factors/mm3) is the concentration of a at which Ktrans reaches half of its maximum value. Thus, as the concentration of the angiogenic factors increases, the degree of vessel permeability will also increase until saturated.
Modeling anti-angiogenic treatment
The PIHNA-E model describes the evolution of the tumor and its microenvironment in an untreated context. By understanding the premise of how specific therapies are meant to alter the system, one can also model the effects of various treatments. Here we are interested in anti-angiogenic treatment and, while there are many different types of drugs for this action, will focus on the drug bevacizumab.
Bevacizumab is a drug specifically targeted at the molecule vascular endothelial growth factor A (VEGF A). This particular angiogenic factor stimulates the growth of new vessels by binding with the vascular endothelial growth factor receptor (VEGFR2) on endothelial cells. Bevacizumab inhibits angiogenesis by binding to the free molecules of VEGF A and preventing them from binding to VEGFR2. An unintended consequence of this drug, however, is that beyond preventing the growth of new vessels, it also “normalizes” pre-existing vasculature (Jain, ; Verhoeff et al., 2009). That is, once the levels of stimulating angiogenic factors are reduced, the vessels are able to repair their leakiness and return to a normal size – making them more efficient nutrient deliverers. In our model, both of these phenomena can be captured by requiring higher levels of angiogenic factor to be present to have the same level of “action” in the contexts of both vessel proliferation and vessel permeability. Additionally, since the treatment is making the vessels more efficient, the level of vasculature needed for a cell to be normoxic will decrease, which we can capture by modifying the cell conversion rates from hypoxic to normoxic and from normoxic to hypoxic. Treatment is approximated by decreasing the parameter for conversion from normoxic to hypoxic (β) by a factor of 10, increasing the parameter for conversion from hypoxic and normoxic (γ) by a factor of 10, and increasing the required levels of angiogenic factors for inducing vascular growth and vessel permeability by 2 as supported by the studies in Desjardins et al. (), Zhang et al. (2009). The treatment modification of β and γ is representative of a dramatic increase in the efficiency of the blood vessels, though exact changes are not available from experimental data.
Simulations of glioblastoma growth and response to anti-angiogenic therapy
For simplicity, we consider in all simulations here a two-dimensional tumor growing on one axial slice of the brain, with the brain geometry defined from the BrainWeb atlas (Cocosco et al., ). The brain is primarily composed of three different types of matter, CSF, gray matter, and white matter. Glioblastomas originate in gray or white matter and due to physical barriers will not enter into the regions of CSF. New mass will often deform the barriers, a phenomena called mass effect, and while there are some models that attempt to capture this (Clatz et al., ; Mohamed and Davatzikos, ; Hogea et al., ), here the brain is considered a stationary domain.
In all simulations, the domain is taken to be a slice of human brain embedded in a grid [0, 147] mm × [0, 185] mm and the equations are spatially discretized on a grid with resolution of 1 mm × 1 mm using first order accurate finite volumes. Time integration is done with an operator splitter technique utilizing backward Euler for the diffusion terms and the TR-BDF2 algorithm (Leveque, ) for the reaction terms with a time step size of 1 day. The simulations were initiated with a small amount of normoxic cells distributed as where (x0, y0) = (103, 83). Vasculature is set at 3% of the cell carrying capacity in all gray and white matter based on estimates from Blinkov and Glezer (), and all other quantities in the PIHNA-E model are initiated to zero. Unless otherwise stated, parameter values used in simulations for the edema equation are in Table 2, the additional parameter values are taken as specified in Swanson et al. (2011).
Table 2
| Parameter | Definition | Value | Reference |
|---|---|---|---|
| Khalf | Michaelis–Menten half-max of response of ECs to angiogenic factors | 5.75e−7 (mmol/cc tissue) | Derived from Mac Gabhann and Popel () |
| Dl | Net rate of edematous fluid diffusion | 0.77e−3 (mm2/s) | Chosen as average ADC value in normal brain tissue as given in Moritani () |
| Kmax | Maximum Ktrans value in response to angiogenic factors | 36 (1/day) | Taken to match the maximum Ktrans value observed in Grade IV gliomas in Patankar et al. () |
| δI | Edematous fluid reabsorption rate | 0.3 × Kmax (1/day) | Assumed proportional to vessel permeability |
Parameter values for the edema equation used in all simulations unless stated otherwise in the text.
All other parameter values are taken as described in Swanson et al. (2011).
Results
Decoupling imaging changes from tumor response
To highlight the real impact of anti-angiogenic treatment as captured by our model, we chose parameter values that represent a patient with an aggressive GBM (net invasion rates Dw = 53 mm2/year and Dg = 0.53 mm2/year and net proliferation rate ρ = 75 1/year) and simulate tumor growth without treatment (Figure 4) and then compare to tumor growth with treatment (Figure 5). For comparison to what was done in previous work (Swanson et al., 2008a, 2011; Rockne et al., ), we refer to the region with total cell density summing to 80% of the cell carrying capacity (K) to correspond to what would enhance on the T1Gd, and started treatment when the T1Gd spherically symmetric equivalent radius was equal to 1 cm and was terminated 100 days later. This is consistent with a typical size of an abnormality seen clinically for consideration of anti-angiogenic treatment. Although treatment length can vary, 100 days is representative of the length of a typical course of treatment with anti-angiogenics in human GBM.
Figure 4
Figure 5
Snapshots of the untreated case are shown in Figure 4 with the analogous snapshots corresponding to the same time points of the treated tumor being shown in Figure 5. The top row in both Figures 4 and 5 shows the density of the bulk tumor (the summed density of all the cell populations: normoxic, hypoxic, necrotic, and vasculature). The second row shows what the simulated FLAIR corresponding to the microenvironmental levels of edematous extracellular fluid. These figures also contain radial growth plots showing the equivalent spherical radii for the regions of interest corresponding to the tumor and the edema. The line representing the bulk tumor is calculated from the volume of tissue containing abnormal cells, normoxic, hypoxic, and necrotic, at levels greater than or equal to 16% of the carrying capacity, i.e., a density five times lower than what can be visualized on T1Gd (Swanson et al., 2008b). The edema radius was defined by considering the volume containing edematous fluid above 50% of the fluid level in the capillaries. There is no literature to guide the choice of the cutoffs for fluid volume constituting T2/FL visible edema. Thus, cutoffs were chosen to roughly match clinically observed behavior.
In the untreated case (Figure 4) the radial plot shows the size of the region impacted by edema evolves very similarly to the size of the region occupied by the bulk tumor throughout the entire course of growth. In the treated case (Figure 5), the edema grows at the same rate as the bulk tumor until the treatment begins at which point the edema begins to decline. Edema begins to increase again once the angiogenic factors have been able to accumulate at levels which overcome the impact of the anti-angiogenic drug. Once treatment is terminated the edema levels rise to again occupy a region of the same size as the bulk tumor.
These simulations support the hypothesis that anti-angiogenic treatment may not have a significant effect on the growth dynamics of the overall cell populations, while having a large impact on the level of edematous extracellular fluid and thus on what would be imageable on T2/FLAIR MRI. This is also in agreement with the current understanding that anti-angiogenic treatment serves as a surrogate to steroids for reducing swelling.
Exploring response across tumor kinetics
The virtual control experiment illustrated in Figures 4 and 5 is only providing insight into tumor/edema response for the case of one set of tumor growth kinetics. However, the range of radiographic response patterns seen clinically is broad, a few examples of which are shown in Figure 2. Previous work has shown patient specific values of net proliferation and invasion range over many orders of magnitude (Tracqui et al., 1995; Harpold et al., ). To investigate if the different types and extents of radiographic responses could be explained by different underlying tumor growth kinetics, we simulated tumor growth and the associated edema under treatment for many different combinations of net proliferation rates, ρ, and net invasion rates, Dw and Dg.
For all cases, treatment was started when the tumor reached a simulated 1 cm T1Gd radius for a total of 100 continuous days. Treatment was implemented in the same manner as the first case. Illustrative results are shown in Figure 6 from simulations with ρ = [5, 75, 125] (1/year) and Dw = [5.3, 53] (mm2/year) and Dg = Dw/10 consistent with the observed range hypothesized in human gliomas (Harpold et al., ) and since observed (Wang et al., 2009; Rockne et al., ).
Figure 6
From these six scenarios, summarized in Figure 6, we were able to observe a few trends. First, none of the simulations showed a significant change in the bulk tumor growth rate after treatment had begun. However, levels of edema were impacted and by different degrees in each of the simulations. In general, tumors with higher proliferative capacities (higher ρ’s), due to their higher metabolic needs, have larger regions of hypoxia and thus produce greater levels of angiogenic factors. Treatment of these tumors initially reduces the level of edema, however, the tumor quickly produces enough angiogenic factors to continue progressing under imaging. In contrast, the slower growing tumors (low ρ) produce low levels of angiogenic factors and treatment may result in a complete disappearance of the abnormality on imaging. Additionally, higher dispersion rates (high Dw and Dg) reduce the local metabolic needs and the production of angiogenic factors. Thus, the imaging of tumors with high Dw and high Dg improve for a longer time, however, ultimately the tumor does still produce enough angiogenic factors to be seen progressing on imaging. For all cases, once treatment is discontinued, the edema levels quickly rise to again reflect the underlying tumor burden.
What is particularly encouraging is that many of the different behaviors observed for the four patients illustrated in Figure 2 can be connected to different simulation predictions. For example, the first patient’s behavior is analogous to the moderate proliferation rate with a low invasion rate in that the imageable lesion initially decreased, stabilized, but after treatment dramatically increased. Additionally, the second patient can be compared to either the moderate or highly proliferative rate with a high invasion rate where the hyperintensity is seen to decrease at the beginning of treatment, but even while treatment is continuing, start growing again.
Discussion
Anti-angiogenics remain a controversial form of treatment for GBM due to the difficulty in assessing tumor response using MR imaging. The resulting reduction in swelling and related symptoms for a subset of patients keeps it an attractive option despite the lack of evidence of an increase in overall survival and the possibility of the treatment selecting for a more aggressive phenotype (Ebos et al., ; Pàez-Ribes et al., ; Verhoeff et al., 2009; Keunen et al., ). It is unclear how anti-angiogenic treatment impacts the cell phenotypes present and there is not yet a deep enough understanding or a unifying theory to provide explanation for all the different response patterns. Thus, a priori identification of patients who will receive a significant symptom-reduction benefit remains difficult. As there are other side effects and consequences from anti-angiogenic therapy, being able to make this early distinction would help remove the controversial nature of this therapy.
The model developed in this work, built on the PIHNA model for glioma proliferation and invasion (Swanson et al., 2011; Gu et al., ), is meant to illustrate a first step toward the creation of a tool for identifying patients who will receive the greatest benefit from anti-angiogenic treatment. It captures the formation of edema caused by leaky vasculature, and is thus able to decouple what would be seen on the T2/FLAIR MRI from the true underlying disease burden. In effect, this would help the clinicians to “turn the light back on” by being able to infer the disease burden that lies beyond what is captured by imaging alone. While this model is not meant to capture individual cell behavior, it does provide a map between overall tumor growth kinetics and treatment response on the imaging/continuum scale.
Many simplifying assumptions have been made in the creation of this model such as ignoring the possible direct impact of anti-angiogenic therapy on cell proliferation, the likely presence of thrombosis (Tehrani et al., 2008), and that the hyperintensity on the T2/FLAIR image is entirely a result of fluid leaking from the vasculature. The model is clearly incomplete, and future modifications of the model will need to consider these other phenomena as well as other possible factors on T2/FLAIR hyperintensity, such as higher cell density and additional cytoplasm. However, even in its current state, it has been able to exhibit many of the types of response patterns observed clinically. It is particularly encouraging because the modeling effects of treatment were held constant and only varying the net dispersal and net proliferation rates, D and ρ, respectively, was sufficient to produce a wide range of imaging responses analogous to what is seen clinically. That is, by modeling the treatment in the exact same way for different values of proliferation and diffusion in the tumor growth model, the visible levels of edema are seen to respond in different ways. In general, the simulations predict edema (swelling) to decrease, supporting the role of these drugs as surrogates for steroids for reduction of symptoms, analogous to the most current understanding (Deming, ). Though, the model clearly illustrates that not all patients would receive the same benefit.
Another interesting implication from these simulations is that while different imaging responses to treatment were achieved for the same treatment conditions, in all cases the bulk tumor is seen to progress with little deflection in overall tumor growth rates. This result could be considered evidence against the cytotoxic effects of anti-angiogenic treatment when administered exclusively, also in agreement with current clinical understanding (Verhoeff et al., 2009). While this work is not directly speaking to survival outcomes, we remark that these results highlight the potential for the mathematical model paradigm to serve in evaluating clinical trial outcomes by analyzing relative benefit from anti-angiogenics, especially in the case of low-N trials. As a particular example, such models may have the potential to be applied to patient cohorts for exploring how differential effects of anti-angiogenics on imaging may or may not relate to overall outcomes.
A drawback of this model is the large number of parameters required. In this document, we assumed the primary influential factor were the net rates of invasion and proliferation and thus held all other parameters constant. These other parameters are likely different patient to patient, however, as demonstrated here, changes in the small number of parameters are sufficient to produce a wide variety of imaging responses. While future sensitivity analysis is required, we believe the work is here is evidence that a complicated explanation for the different imaging responses to anti-angiogenic therapy may not be needed. Major next steps of this work will involve developing techniques for obtaining patient specific growth parameters from pre-treatment images which we believe will be successful from previous accomplishments with a simpler model capturing just the proliferation and invasion tumor characteristics PI (Swanson, , ; Swanson et al., , 2003, 2004, 2008a; Szeto et al., 2009; Wang et al., 2009; Rockne et al., ).
Also, modification of the model to capture the pressure induced from the vasogenic edema and possible herniation would allow for deeper understanding of the steroid-like reduction in swelling. For each new feature eventually added validation tests will be required, however, we believe the results presented here in and of themselves represent a significant step in overcoming clinical imaging restrictions with mathematical models.
Statements
Acknowledgments
The authors would like to greatly acknowledge the funding of this work which was sponsored in part by NIH Grants: R01 CA 16437 and U54 CA 143970, along with the NSF Fellowship 1122322.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
References
1
BatchelorT. T.SorensenA. G.di TomasoE.ZhangW. T.DudaD. G.CohenK. S.et al (2007). AZD2171, a pan-VEGF receptor tyrosine kinase inhibitor, normalizes tumor vasculature and alleviates edema in glioblastoma patients. Cancer Cell11, 83–95.10.1016/j.ccr.2006.11.021
2
BatesD. O. (2010). Vascular endothelial growth factors and vascular permeability. Cardiovasc. Res.87, 262–271.10.1093/cvr/cvq105
3
BaxterL. T.JainR. K. (1989). Transport of fluid and macromolecules in tumors. I. Role of interstitial pressure and convection. Microvasc. Res.37, 77–104.10.1016/0026-2862(89)90074-5
4
BaxterL. T.JainR. K. (1990). Transport of fluid and macromolecules in tumors. II. Role of heterogeneous perfusion and lymphatics. Microvasc. Res.40, 246–263.10.1016/0026-2862(90)90023-K
5
BlinkovS. M.GlezerI. I. (1968). The Human Brain in Figures and Tables: A Quantitative Handbook. New York: Basic Books.
6
BohmanL. E.SwansonK. R.MooreJ. L.RockneR.MandigoC.HankinsonT.et al (2010). Preoperative MRI characteristics of glioblastoma multiforme: implications for understanding glioma ontogeny. Neurosurgery67, 1319–1327.10.1227/NEU.0b013e3181f556ab
7
ClatzO.SermesantM.BondiauP. Y.DelingetteH.WarfieldS. K.MalandainG.et al (2005). Realistic simulation of the 3-D growth of brain tumors in MR images coupling diffusion with biomechanical deformation. IEEE Trans. Med. Imaging24, 1334–1346.10.1109/TMI.2005.856734
8
CocoscoC.KollokianV.KwanR. S.EvansA. (1997). “Brainweb: online interface to a 3D MRI simulated brain database,” in Neuroimage, Proceedings of the Third International Conference on the Functional Mapping of the Human Brain, Copenhagen.
9
DemingS. (2012). Bevacizumab for glioblastoma: symptom relief proven, overall survival benefit being studied. OncoLog57. Available at: http://www2.mdanderson.org/depts/oncolog/articles/pf/12/6-jun/6-12-2-pf.html
10
DesjardinsA.BarboriakD. P.HerndonJ. E.IIReardonD. A.QuinnJ. A.RichJ. N.et al (2007). Dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) evaluation in glioblastoma (GBM) patients treated with bevacizumab (BEV) and irinotecan (CPT-11). J. Clin. Oncol.25(Suppl.), abstr. 2029.
11
DunnG. P.RinneM. L.WykoskyJ.GenoveseG.QuayleS. N.DunnI. F.et al (2012). Emerging insights into the molecular and cellular basis of glioblastoma. Genes Dev.26, 756–784.10.1101/gad.187922.112
12
EbosJ. M.LeeC. R.Cruz-MunozW.BjarnasonG. A.ChristensenJ. G.KerbelR. S. (2009). Accelerated metastasis after short-term treatment with a potent inhibitor of tumor angiogenesis. Cancer Cell15, 232–239.10.1016/j.ccr.2009.01.021
13
FolkmanJ. (1971). Tumor angiogenesis: therapeutic implications. N. Engl. J. Med.285, 1182–1186.10.1056/NEJM197111182852108
14
GuS.ChakrabortyG.ChampleyK.AlessioA. M.ClaridgeJ.RockneR.et al (2012). Applying a patient-specific bio-mathematical model of glioma growth to develop virtual [18F]-FMISO-PET images. Math. Med. Biol.29, 31–48.10.1093/imammb/dqr002
15
HarpoldH. L.AlvordE. C.Jr.SwansonK. R. (2007). The evolution of mathematical modeling of glioma proliferation and invasion. J. Neuropathol. Exp. Neurol.66, 1–9.10.1097/nen.0b013e31802d9000
16
HogeaC.DavatzikosC.BirosG. (2008). An image-driven parameter estimation problem for a reaction-diffusion glioma growth model with mass effects. J. Math. Biol.56, 793–825.10.1007/s00285-007-0139-x
17
HolashJ.MaisonpierreP. C.ComptonD.BolandP.AlexanderC. R.ZagzagD.et al (1999). Vessel cooption, regression, and growth in tumors mediated by angiopoietins and VEGF. Science284, 1994–1998.10.1126/science.284.5422.1994
18
HurwitzH.FehrenbacherL.NovotnyW.CartwrightT.HainsworthJ.HeimW.et al (2004). Bevacizumab plus irinotecan, fluorouracil, and leucovorin for metastatic colorectal cancer. N. Engl. J. Med.350, 2335–2342.10.1056/NEJMoa032691
19
IwamotoF. M.AbreyL. E.BealK.GutinP. H.RosenblumM. K.ReuterV. E.et al (2009). Patterns of relapse and prognosis after bevacizumab failure in recurrent glioblastoma. Neurology73, 1200–1206.10.1212/WNL.0b013e3181bc0184
20
JainR. K. (2005). Normalization of tumor vasculature: an emerging concept in antiangiogenic therapy. Science307, 58–62.10.1126/science.1104819
21
KerbelR. S. (2000). Tumor angiogenesis: past, present and the near future. Carcinogenesis21, 505–515.10.1093/carcin/21.3.505
22
KeunenO.JohanssonM.OudinA.SanzeyM.RahimS. A.FackF.et al (2011). Anti-VEGF treatment reduces blood supply and increases tumor cell invasion in glioblastoma. Proc. Natl. Acad. Sci. U.S.A.108, 3749–3754.10.1073/pnas.1014480108
23
KumagaiY. (1986). Mathematical model of brain edema and optimal control of intracranial pressure. No To Shinkei38, 885–891.
24
LeendersW. P.KüstersB.de WaalR. M. (2002). Vessel co-option: how tumors obtain blood supply in the absence of sprouting angiogenesis. Endothelium9, 83–87.10.1080/10623320212006
25
LeVequeR. (2005). Numerical Methods for Conservation Laws, 2nd Edn. Basel: Birkhäuser.
26
LevineH. A.PamukS.SleemanB. D.Nilsen-HamiltonM. (2001). Mathematical modeling of capillary formation and development in tumor angiogenesis: penetration into the stroma. Bull. Math. Biol.63, 801–863.10.1006/bulm.2001.0240
27
LouisD. N.OhgakiH.WiestlerO. D.CaveneeW. K.BurgerP. C.JouvetA.et al (2007). The 2007 WHO classification of tumours of the central nervous system. Acta Neuropathol.114, 97–109.10.1007/s00401-007-0278-6
28
Mac GabhannF.PopelA. S. (2004). Model of competitive binding of vascular endothelial growth factor and placental growth factor to VEGF receptors on endothelial cells. Am. J. Physiol. Heart Circ. Physiol.286, H153–H164.10.1152/ajpheart.00254.2003
29
MarmarouA. (2007). A review of progress in understanding the pathophysiology and treatment of brain edema. Neurosurg. Focus22, E1.10.3171/foc.2007.22.3.1
30
MohamedA.DavatzikosC. (2005). Finite element modeling of brain tumor mass-effect from 3D medical images. Med. Image Comput. Comput. Assist. Interv.2005(Pt 13749), 400–408.
31
MoritaniT. (2009). Diffusion-Weighted MR Imaging of the Brain. New York: Springer.
32
NagashimaT.HorwitzB.RapoportS. I. (1990). A mathematical model for vasogenic brain edema. Adv. Neurol.52, 317–326.
33
NordenA. D.YoungG. S.SetayeshK.MuzikanskyA.KlufasR.RossG. L.et al (2008). Bevacizumab for recurrent malignant gliomas: efficacy, toxicity, and patterns of recurrence. Neurology70, 779–787.10.1212/01.wnl.0000304121.57857.38
34
Pàez-RibesM.AllenE.HudockJ.TakedaT.OkuyamaH.ViñalsF.et al (2009). Antiangiogenic therapy elicits malignant progression of tumors to increased local invasion and distant metastasis. Cancer Cell15, 220–231.10.1016/j.ccr.2009.01.027
35
PatankarT. F.HaroonH. A.MillsS. J.BalériauxD.BuckleyD. L.ParkerG. J.et al (2005). Is volume transfer coefficient (K(trans)) related to histologic grade in human gliomas?AJNR Am. J. Neuroradiol.26, 2455–2465.
36
RapoportS. I. (1978). A mathematical model for vasogenic brain edema. J. Theor. Biol.74, 439–467.10.1016/0022-5193(78)90224-2
37
RockneR.RockhillJ. K.MrugalaM.SpenceA. M.KaletI.HendricksonK.et al (2010). Predicting the efficacy of radiotherapy in individual glioblastoma patients in vivo: a mathematical modeling approach. Phys. Med. Biol.55, 3271–3285.10.1088/0031-9155/55/12/001
38
SandlerA.GrayR.PerryM. C.BrahmerJ.SchillerJ. H.DowlatiA.et al (2006). Paclitaxel-carboplatin alone or with bevacizumab for non-small-cell lung cancer. N. Engl. J. Med.355, 2542–2550.10.1056/NEJMoa061884
39
SeriniG.AmbrosiD.GiraudoE.GambaA.PreziosiL.BussolinoF. (2003). Modeling the early stages of vascular network assembly. EMBO J.22, 1771–1779.10.1093/emboj/cdg176
40
SherrattJ. A.MurrayJ. D. (1990). Models of epidermal wound-healing. Proc. Biol. Sci.241, 29–36.10.1098/rspb.1990.0061
41
SilbergeldD. L.ChicoineM. R. (1997). Isolation and characterization of human malignant glioma cells from histologically normal brain. J. Neurosurg.86, 525–531.10.3171/jns.1997.86.3.0525
42
StuppR.MasonW. P.van den BentM. J.WellerM.FisherB.TaphoornM. J.et al (2005). Radiotherapy plus concomitant and adjuvant temozolomide for glioblastoma. N. Engl. J. Med.352, 987–991.10.1056/NEJMoa043330
43
SwansonK. R. (1999). Mathematical Modeling of the Growth and Control of Tumors. Applied Mathematics. Ph.D. thesis, University of Washington, Seattle.
44
SwansonK. R. (2002). Serial imaging observations and postmortem examination of an untreated glioblastoma: a traveling wave of glioma growth and invasion. Neuro-oncology4, 340.
45
SwansonK. R.AlvordE. C.Jr.MurrayJ. D. (2000). A quantitative model for differential motility of gliomas in grey and white matter. Cell Prolif.33, 317–329.10.1046/j.1365-2184.2000.00177.x
46
SwansonK. R.AlvordE. C.Jr.MurrayJ. D. (2002a). Quantifying efficacy of chemotherapy of brain tumors with homogeneous and heterogeneous drug delivery. Acta Biotheor.50, 223–237.10.1023/A:1022644031905
47
SwansonK. R.AlvordE. C.Jr.MurrayJ. D. (2002b). Virtual brain tumours (gliomas) enhance the reality of medical imaging and highlight inadequacies of current therapy. Br. J. Cancer86, 14–18.10.1038/sj.bjc.6600021
48
SwansonK. R.AlvordE. C.MurrayJ. D. (2004). Dynamics of a model for brain tumors reveals a small window for therapeutic intervention. Discrete Continuous Dyn. Syst. Ser. B4, 289–295.
49
SwansonK. R.BridgeC.MurrayJ. D.AlvordE. C.Jr. (2003). Virtual and real brain tumors: using mathematical modeling’ to quantify glioma growth and invasion. J. Neurol. Sci.216, 1–10.10.1016/j.jns.2003.06.001
50
SwansonK. R.HarpoldH. L.PeacockD. L.RockneR.PenningtonC.KilbrideL.et al (2008a). Velocity of radial expansion of contrast-enhancing gliomas and the effectiveness of radiotherapy in individual patients: a proof of principle. Clin. Oncol.20, 301–308.10.1016/j.clon.2008.01.006
51
SwansonK. R.RostomilyR. C.AlvordE. C.Jr. (2008b). A mathematical modelling tool for predicting survival of individual patients following resection of glioblastoma: a proof of principle. Br. J. Cancer98, 113–119.10.1038/sj.bjc.6604125
52
SwansonK. R.RockneR. C.ClaridgeJ.ChaplainM. A.AlvordE. C.Jr.AndersonA. R. (2011). Quantifying the role of angiogenesis in malignant progression of gliomas: in silico modeling integrates imaging and histology. Cancer Res.71, 7366–7375.10.1158/0008-5472.CAN-11-1399
53
SzetoM. D.ChakrabortyG.HadleyJ.RockneR.MuziM.AlvordE. C.Jr.et al (2009). Quantitative metrics of net proliferation and invasion link biological aggressiveness assessed by MRI with hypoxia assessed by FMISO-PET in newly diagnosed glioblastomas. Cancer Res.69, 4502–4509.10.1158/0008-5472.CAN-08-3884
54
TehraniM.FriedmanT. M.OlsonJ. J.BratD. J. (2008). Intravascular thrombosis in central nervous system malignancies: a potential role in astrocytoma progression to glioblastoma. Brain Pathol.18, 164–171.10.1111/j.1750-3639.2007.00108.x
55
ToftsP. S. (1991). Measurement of the blood-brain barrier permeability and leakage space using dynamic MR imaging. 1. Fundamental concepts. Magn. Reson. Med.17, 357–367.10.1002/mrm.1910170208
56
TracquiP.CruywagenG. C.WoodwardD. E.BartooG. T.MurrayJ. D.AlvordE. C.Jr. (1995). A mathematical model of glioma growth: the effect of chemotherapy on spatio-temporal growth. Cell Prolif.28, 17–31.10.1111/j.1365-2184.1995.tb00036.x
57
VerhaakR. G.HoadleyK. A.PurdomE.WangV.QiY.WilkersonM. D.et al (2010). Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1. Cancer Cell17, 98–110.10.1016/j.ccr.2009.12.020
58
VerhoeffJ. J.van TellingenO.ClaesA.StalpersL. J.van LindeM. E.RichelD. J.et al (2009). Concerns about anti-angiogenic treatment in patients with glioblastoma multiforme. BMC Cancer9:444.10.1186/1471-2407-9-444
59
WangC. H.RockhillJ. K.MrugalaM.PeacockD. L.LaiA.JuseniusK.et al (2009). Prognostic significance of growth kinetics in newly diagnosed glioblastomas revealed by combining serial imaging with a novel biomathematical model. Cancer Res.69, 9133–9140.10.1158/0008-5472.SABCS-09-2028
60
XiuM.TurnerS. M.BuschR.GeeT. A.HellersteinM. K. (2006). Measurement of endothelial cell proliferation rate in vivo using 2H20 labeling: a kinetics biomarker of angiogenesis. FASEB J.20, A718–A718.
61
ZhangW.KreislT. N.SolomonJ.ReynoldsR. C.GlenD. R.CoxR. W.et al (2009). Acute effects of bevacizumab on glioblastoma vascularity assessed with DCE-MRI and relation to patient survival. Proc. Intl. Soc. Mag. Reson. Med.17, 282.
Summary
Keywords
glioma, edema, mathematical model, anti-angiogenic therapy
Citation
Hawkins-Daarud A, Rockne RC, Anderson ARA and Swanson KR (2013) Modeling Tumor-Associated Edema in Gliomas during Anti-Angiogenic Therapy and Its Impact on Imageable Tumor. Front. Oncol. 3:66. doi: 10.3389/fonc.2013.00066
Received
20 December 2012
Accepted
14 March 2013
Published
04 April 2013
Volume
3 - 2013
Edited by
Heiko Enderling, Tufts University School of Medicine, USA
Reviewed by
Janusz Rak, McGill University, Canada; Kerri-Ann Norton, Johns Hopkins University, USA
Copyright
© 2013 Hawkins-Daarud, Rockne, Anderson and Swanson.
This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.
*Correspondence: Kristin R. Swanson, Department of Neurological Surgery, Feinberg School of Medicine, Northwestern University, 680 N Lake Shore Drive, Suite 1500, Chicago, IL 60611, USA. e-mail: kristin.swanson@northwestern.edu
This article was submitted to Frontiers in Molecular and Cellular Oncology, a specialty of Frontiers in Oncology.
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.