Modeling the efficacy of the extent of surgical resection in the setting of radiation therapy for glioblastoma.

Standard therapy for glioblastoma (GBM) includes maximal surgical resection and radiation therapy. While it is established that radiation therapy provides the greatest survival benefit of standard treatment modalities, the impact of the extent of surgical resection (EOR) on patient outcome remains highly controversial. While some studies describe no correlation between EOR and patient survival even up to total resection, others propose either qualitative (partial versus subtotal versus complete resection) or quantitative EOR thresholds, below which there is no correlation with survival. This work uses a mathematical model in the form of a reaction-diffusion partial differential equation to simulate tumor growth and treatment with radiation therapy and surgical resection based on tumor-specific rates of diffusion and proliferation. Simulation of 36 tumors across a wide spectrum of diffusion and proliferation rates suggests that while partial or subtotal resections generally do not provide a survival advantage, complete resection significantly improves patient outcomes. Furthermore, our model predicts a tumor-specific quantitative threshold below which EOR has no effect on patient survival and demonstrates that this threshold increases with tumor aggressiveness, particularly with the rate of proliferation. Thus, this model may serve as an aid for determining both when surgical resection is indicated as well as the surgical margins necessary to provide clinically significant improvements in patient survival. In addition, by assigning relative benefits to radiation and surgical resection based on tumor invasiveness and proliferation, this model confirms that (with the exception of the least aggressive tumors) the survival benefit of radiation therapy exceeds that of surgical resection.

Current standard therapy for glioblastoma (GBM) includes maximal surgical resection and radiation therapy with temozolomide.1 Yet, the impact of the extent of surgical resection (EOR) on patient outcome remains controversial.2, 3, 4, 5 Two large reviews examining the literature prior to 1990 and one later review in 1999 were unable to conclude that there was a relationship between EOR and patient survival.3, 4, 5 A more recent review of 28 studies by Sanai et al.6 spanning the period from 1990 to 2007 demonstrated mixed results but the weight of the evidence suggested a positive correlation between survival and EOR, particularly with gross total resection. Essentially all of the studies included radiation therapy, but analysis is complicated both by non‐volumetric assessment of tumor size and by the different reporting schemes used in the papers. Most of the papers used a three‐tier scheme of biopsy, subtotal resection and gross total resection. Of the 24 papers using non‐volumetric analysis, 14 papers found a positive correlation between EOR and survival. More standardization in measurement and reporting is available in the few papers that used volumetric analysis. However, no clear consensus exists in this group either. For instance, Pope et al.7 found no significant impact of EOR on patient survival, while Keles et al.8 and Lacroix et al.9 were positive studies.

Mathematical modeling of GBM treatment that incorporates the impact of both surgery and radiation therapy may provide further insight into the effect of standard therapy on patient tumor burden and survival. Accurate modeling may also provide an avenue to understanding some of the disparities that currently exist in the literature. This paper uses the well‐studied reaction–diffusion model of GBM to investigate the combined effects of standard radiation therapy and surgical resection on patient survival. Although there has been significant work with this model in the published literature, including both the effect of surgical resection,10, 11 and more recently the effect of radiation therapy,12, 13, 14 prior analysis assessing the impact of the combination of both surgery and radiation therapy to provide a full treatment model for GBM is limited.

Materials and Methods

Mathematical modeling

The proliferation and infiltration of glioblastoma (GBM) into adjacent tissues has previously been modeled using a partial differential equation reaction–diffusion model of tumor growth.11, 15, 16, 17 In one spatial dimension, this model can be written as

where the various terms are defined as follows: c(t,x) is the tumor cell density, in terms of cells/mm3, which is a function of position x and time t; D(x) is the diffusion term, in mm2/day, which models local tumor invasion of tumor cells; and is a logistic tumor growth term, where ρ is the tumor proliferation rate in units of (/day), governed by a tissue tumor carrying capacity K, in units of cells/mm3.

The carrying capacity for the tissue, K, can be considered a cell density of 105 cells per mm3.18 The solution domain is L = 200 mm, and we use the standard zero‐flux boundary conditions at x = 0 and at x = L.12 The partial differential equations are solved using the pdepe function in MATLAB.

The above reaction–diffusion model provides an evolution of the tumor cell density with time, where cell density gradually increases toward the carrying capacity K and where the tumor grows in space, with the tumor cell‐density curve moving to the right (Fig. 1).

Figure 1

The time evolution of the cell density profile as given by the standard reaction–diffusion model of Equation (1), with , K = 105 cells per mm3, D (0.4 mm2/day) and ρ (0.04/day).

As can be seen from Figure 1, at any time point, the model predicts a gradient of cell density, decreasing with distance from the center of the tumor. These mathematical predictions mirror direct histological observations of cell density gradients in GBM.19 The model estimates this gradient by setting cell density thresholds for tumor visibility on both T1 post‐contrast images and FLAIR/T2 images (Fig. 2).14

Figure 2

Schematic showing cell density versus distance from tumor center, illustrating the use of cell density thresholds to delineate calculated tumor radius on T1 post‐contrast and T2 weighted images.

The T1 detection threshold, used to simulate the enhancing T1 tumor radius, is set at 0.8 K, while the T2 threshold is 0.16 K, maintaining consistency with values from prior work.14

Rockne et al.12, 13 add a radiation therapy term to this equation to develop the proliferation–invasion–radiation therapy (PIRT) version of the model:

The PIRT model allows the modeling of a tumor's response to radiation therapy as the tumor both proliferates and spreads in space, differentiating this model from other models in which there is no motion of tumor cells. This model is also able to take into account the behavior of tumor cells below the MR detection threshold.

The R term quantifies the loss of tumor cells due to radiation therapy, which is delivered in discrete doses, and is, hence, amenable to modeling different dosing schedules. R is defined as a function of S, the fraction of cells surviving a radiation dose, using the well‐known linear‐quadratic dose‐response model: , where α (in units of per Gy) and β (in units of per Gy2) reflect type A (single ionizing event) and type B (pairwise interaction of ionizing events) tissue damage. Because tissues can be somewhat characterized by an α/β ratio, and to simplify the model to a single radiation parameter, the α/β is held fixed, as per Rockne et al.12, 13 In this work, similar to prior works, this ratio is held at 10 Gy.12, 13, 14 Thus, S can be written as: .

If S is the surviving cell fraction for a given dose, then (1 − S) is the probability of cell death. In this formulation, α can be regarded as the parameter defining radiation sensitivity. A higher α corresponds to a decreased probability of survival. In addition, d reflects the given radiation dose, which is a function of both space and time, d(x,t), allowing specific radiation therapy protocols to be modeled.

In the PIRT model, the effect of radiation therapy is also cast as a function of cell density, using the same logistic formulation as the tumor growth model as follows: .13

Hence, at locations of high cell density, when the tumor cell proliferation is decreased according to the logistic configuration, the effect of radiation therapy is likewise decreased. This is consistent with the understanding that radiation is most effective in regions of high mitotic activity and high cellular turnover. Conversely, at low cell densities, the effect of radiation therapy is essentially linearly related to the fraction of cells killed. The above formulation is equivalent to the Taylor series approximation of having the logistic term as part of the exponent of S.

This mathematical model can be further modified by the addition of a surgical term, G(x,t):

This allows surgeries of differing radii to be simulated, and even allows the simulation of repeat surgical resection.

Methods

The above model was used to simulate growth and treatment of GBM. Radiation therapy was modeled as a total of 61.2 Gy/34 doses, using the University of Washington protocol. The first 28 doses were given in a field using the T2 tumor boundary + 2.5 cm, while the last 6 doses were given as a more spatially limited booster field using the enhanced T1 radius + 2 cm.12

Surgical radii of 0, 5, 10 and 20 mm were applied across a range of model tumors in four groups covering the combinations of high and low invasiveness and proliferation (designated as high and low D and high and low ρ), using published ranges for D and ρ for GBM.12, 17 This is done to mirror the heterogeneity of biologic behavior of GBM, including the newly discovered genetic heterogeneity, with factors such as MGMT promoter status and IDH1 mutations, which significantly affect biologic behavior.

The high D–high ρ group (HH) consisted of a center value set of D = 0.2 cm2/day and ρ = 0.1/day. The low D–low ρ group (LL) consisted of a center value set of D = 0.04 cm2/day and ρ = 0.02/day. Using these values, high D– low ρ (HL) and low D–high ρ (LH) tumors were also modeled. For each group, nine theoretical patients were generated by taking the combinations of D ± 50% and ρ ± 50%, following the method of Woodward et al.10 Theoretical survival curves generated for each group using a preset survival endpoint of enhanced T1 tumor radius of 3.5 cm were used to allow inter‐group and intra‐group comparison. This lethal radius is the same as that used in prior modeling work by Rockne et al.13 and Swanson et al. Survival curves were generated for no therapy, radiation therapy alone (i.e. surgical radius = 0), and surgical radii of 5, 10 and 20 mm, followed by radiation therapy. A surgical radius of 5 mm is termed partial resection, a radius of 10 mm is designated as subtotal resection, while a radius of 20 mm, which actually goes beyond the 1.5‐cm enhanced tumor margin, is termed gross total resection.

Relying on previously published estimates, including the variation of α with tumor proliferation, an α of 0.03/Gy was used for the low p group and an α of 0.09/Gy was used for the high p group. The α/β ratio was held fixed at 10 Gy.12, 13, 14, 20

Results

The results of the radiation therapy model can be displayed in several ways. An instructive approach is to look at the cell density in space, at given time points, such as just prior to therapy and immediately following the conclusion of therapy (Fig. 3).

Figure 3

(a) Cell density profile in space at two time points: just prior to (red curve) and just after radiation therapy (blue curve) for a model tumor (D = 0.2, ρ = 0.02, α = 0.04). It is noted that radiation therapy causes a decrease in the cell density curve. (b) Effect on the cell density profile of a partial (5‐mm) surgical resection, followed by radiation therapy (blue cruve).

A second way to look at the model data is tumor radius versus time plots, using the thresholding technique described above. Figure 4 shows the radius versus time plot of a tumor undergoing therapy when the enhanced T1 radius (T1C) reaches 1.5 cm. As can be seen from the plots, the tumor regrows following therapy, until the pre‐set endpoint of patient death (T1C radius of 3.5 cm).

Figure 4

(a) T1C radius versus time plots for a model tumor (D = 0.2, ρ = 0.03, α = 0.04). The blue curve represents tumor growth with no therapy; death occurs at 357 days. The red curve represents the changes caused by radiation therapy, with a reduction in radius during therapy, followed by tumor rebound; death now occurs at 422 days. It is noted that the two curves overlap prior to therapy (initial blue segment). (b) Compares radiation therapy alone (blue then red curve) to gross total surgical resection followed by radiation therapy (blue curve); death now occurs at 468 days. It is noted that with surgery followed by radiation therapy, the model tumor disappears on T1 post‐contrast images, with a long lag time (approximately 200 days), followed by tumor recurrence.

Applying the model in this way, theoretical survival curves can be generated for the different tumor groups (Fig. 5).

Figure 5

Theoretical survival curves for nine patient cohorts in the HH, HL, LH and LL tumor groups. The first curve is without therapy. The second curve represents the improved survival with radiation therapy. Overlapped on this curve are the curves for partial and subtotal resection, showing no added survival benefit. The last curve represents gross total resection followed by radiation therapy, showing an additional benefit for total resection.

Average untreated patient survival time across all tumor groups was found to be 7.8 months, ranging from 3.0 months in the HH group to 14.8 in the LL group. Modeling of radiation therapy alone showed a significant increase in patient survival across all tumor groups, increasing average survival time from 7.8 to 10.1 months. This constitutes an average of a 29% increase in patient survival time. Radiation therapy alone was found to be most effective in the HH group, increasing average survival time by 54%, and least effective in the LL group for which average survival time increased by 19%. Incorporation of surgical resection into the simulated treatment regimen demonstrated vastly different effects depending on EOR. Gross total resection (20 mm) increased average survival across all tumor groups from 10.1 months with radiation alone to 12 months, showing an independent additive effect to radiation therapy. The combination of radiation therapy and gross total resection was found to be most effective in the HL group in which survival time was increased by 83% as compared to untreated tumor progression. Contrary to gross total resection, partial (5 mm) and subtotal resection (10 mm) showed no significant improvement in survival in any of the tumor groups.

Discussion

For most cancers the extent of surgical resection would seem likely to correlate with patient survival. However, the relationship between these two variables for GBM remains controversial. The spread of GBM cells beyond visible tumor margins and the inability to perform broad resections while preserving critical brain functions may underlie this disconnect. Indeed, GBM recurs even after total hemispherectomy.2 Earlier literature reviews could not reach a conclusion as to the relationship between EOR and survival.3, 4, 5 Interpretation of published data is complicated by lack of standardization in both the reporting and measurement of EOR. For example, many studies did not use volumetric analysis for optimal measurement of EOR, and oftentimes EOR was not confirmed by postoperative imaging.2, 21, 22

Given these difficulties, mathematical modeling of GBM growth and treatment accounting for both surgery and radiation therapy might add insight beyond that provided by empirical evidence, particularly as to whether the extent of surgical resection adds a survival advantage to radiotherapy alone.4

We modeled four GBM patient groups with various combinations of high and low invasiveness and proliferation in order to account for heterogeneity of the biologic behavior of GBM. In the aggregate, radiation therapy alone showed a significant increase in patient survival from 7.8 to 10.1 months (29% increase). These results, with a 2.3‐month increase, are comparable to published data on the effect of radiation therapy alone in a recent trial of elderly patients, where survival was prolonged by 2.9 months (from 3.9 to 6.8 months).23 The overall survival times in that trial are shorter than the averages for our model, potentially reflecting the more aggressive biologic behavior of GBM thought to be associated with elderly patients. If only high proliferation tumor cohorts in our model are considered, the average survival time with radiation therapy alone is 6.2 months, comparable to this published data.23

In our model, gross total resection (20 mm) further increased overall survival from 10.1 to 12.0 months, showing an independent additive effect to radiation therapy. However, partial and subtotal resection showed no additional benefit. These conclusions seem to match more recently published evidence showing a significant survival advantage for gross total resection when added to radiation therapy, in contrast to little survival advantage achieved with partial or subtotal resection. For instance, meta‐analysis indicates that patients undergoing subtotal resection have an average survival of 11.3 months, compared to 14.2 months for those with gross total resection; that is, a survival advantage of slightly less than 3 months.6 This association is particularly evident in papers using volumetric analysis,9, 24, 25, 26, 27 and correlates well with our theoretical results of 10.1 months for partial or subtotal resection and 12 months for gross total resection.

One advantage of mathematical modeling is that it may provide additional understanding of empirical data from clinical trials. For instance, Sanai et al. found a threshold value for EOR of 78%, below which there was no correlation with survival, but above which there was a progressive increase in survival with EOR.25 Chaichana et al.24 found a difference between the <70% and >70% groups (10.5 vs 14.4 months average survival, respectively), without graded improvement, while Grabowski et al.27, similar to Lacroix et al., found a cutoff at 98% EOR, with average survival of 14 months for <98% EOR and 16 months for >98% EOR. Our model may help understanding some of these somewhat discrepant findings. For instance, our model supports a threshold effect for EOR. We found that no tumor group showed survival benefit between 0, 5 and 10‐mm radii of resection, but all groups showed a significant improvement at 20 mm. An important additional insight achieved with our model is that the threshold value above which EOR correlates with survival could be highly dependent on tumor aggressiveness, and, in particular, on the proliferation rate ρ. For example, in our model a highly proliferative tumor from the HH group, with D = 0.2 mm2/day and ρ = 0.1/day shows no improvement in survival with a resection radius of 15 mm. This would correlate with a 100% resection of the enhancing tumor, because treatment in the model begins when the T1C is 1.5 cm. There is a very minimal survival advantage beginning at a resection radius of 16 mm, which increases progressively with EOR, becoming significant at 20 mm.

Thus, our model could help explain why some papers find no correlation between EOR and survival, even up to total resection. Our modeling of highly proliferative tumors shows a survival advantage only when surgical resection margins extend beyond the enhancing tumor. In this case, baseline survival with radiation is 16.7 weeks and is not improved by resections of 5, 10 or 15 mm. At 20 mm (i.e. resection with margins beyond T1C), survival improves to 19.1 weeks. If even wider resection margins are possible (e.g. a 30‐mm resection radius), survival increases dramatically to 26.4 weeks. This is because the wide resection reduces the residual cells from which the tumor recurs. However, in a tumor with a mid‐range proliferation index (e.g. D = 0.2 mm2/day and ρ = 0.05/day) a minimal discernable survival advantage (defined as 2 days) begins with a surgical resection radius of 14 mm (e.g. at 81% of enhanced tumor volume). Therefore, the different threshold values between the various papers examining this relationship may correspond to chance variations in the growth characteristics of tumors in their respective patient populations.

Thus, conclusions based on our model are that radiation therapy alone shows a significant increase in calculated patient survival across all tumor groups. Partial and subtotal resection in general show no additional benefit, but gross total surgical resection has an additional added survival benefit, which is less than that of radiation therapy in all groups except the least aggressive tumors. Our model predicts a definite threshold below which EOR has no effect on patient survival and further predicts that this threshold varies with tumor aggressiveness. For tumors with a relatively low proliferation coefficient, there may be progressive benefit with EOR beginning at volumes of even less than 50%. For very highly aggressive tumors, however, only a resection with margins beyond the enhanced T1 tumor radius can prolong survival, and this benefit increases significantly with the size of the surgical margin. For moderately proliferative tumors, survival advantages begin to be attained with resection volumes of approximately 80% of the enhanced tumor volume.

Our model is also able to assign relative benefits to radiation and surgical resection based upon tumor invasiveness and proliferation. We found that the relative benefit of radiotherapy is greatest in the most aggressive tumor group (high D, high ρ), and that radiation therapy is least effective in the least aggressive tumor group. The greatest relative benefit of radiation plus total surgical resection was in the high D–low ρ group (expected to have the greatest tumor infiltration beyond visible MR boundaries, but with a relatively low proliferation rate).

It is important to put our model in context relative to similar efforts, as several other radiation therapy models have been proposed in the literature. For example, the model of Wein et al.28 is a continuum model that uses an expanding spherical shell model of tumor growth with a centrally necrotic core to estimate optimized radiation dosing schedules. This model, although quite useful in suggesting possible optimal radiation fractionation schedules, is only accurate for a prevascular tumor stage, and, as stated by the authors, is not expected to fully model actual tumor growth. Furthermore, because it only maintains a thin shell of viable tumor, it cannot be easily used to model surgery. Our model is ideally adapted to simulating surgery, and allows regrowth of tumor in the surgical bed, mirroring actual tumor recurrence patterns. Another extensively‐used approach to modeling has been that of stochastic models, using Monte Carlo techniques at the level of individual cells or cell clusters. As pointed out by Dionysiou et al. in his review of such models, most are useful only for small in vitro spheroids or very small prevascular tumors, and cannot be used for accurate modeling of large in vivo tumors.29 Moreover, these models cannot be easily “individualized” using parameters that can be measured for each patient. One exception is the four‐dimensional stochastic model proposed by Dionysiou and Stamatakos wherein the 3‐dimensional tumor lattice can be adapted to an individual patient's tumor as outlined on an MRI, and can also incorporate the temporal motion of cells through the cell cycle.29 However, the model proposed here has the advantage of simplicity, because it uses a continuum approach. Similar to the model of Dionysiou and Stamatakos, it can incorporate genetic information to adjust radiation sensitivity parameters. Moreover, it also has the capacity of using serial MRI scans to estimate D and p for an individual tumor, probably more critical factors than radiation sensitivity in patient survival. It would be interesting, in future papers, to compare the simulations of the continuum and discrete stochastic models.

The limitations of our model and similar types of analyses are many, and generally have to do with the simplifications inherent in mathematical modeling. For example, in our model, tumors are modeled as each having a uniform D and ρ, and growing in a spherically symmetric fashion. Tumor invasiveness is modeled as a diffusion process. Clearly, these are gross simplifications. Just as clearly, tumors are not uniformly treated exactly when the T1C radius reaches 1.5 cm, and patients do not uniformly die when that radius reaches 3.5 cm. However, the 3.5‐cm lethal radius is consistent with empirical data on lethal glioma diameters from which Alvord derived a range of tumor diameters from 3.1 to 8.9 cm.30 The hope, therefore, is as previously stated, that a good model “catches the essential truth behind the thing it is trying to depict (p. 32).”31

Future directions seek to extend the model by incorporating other significant variables, such as the effect of residual tumor volume on patient survival, independent of EOR, an association which has support in the literature.24, 27 The impact of temozolomide could also impact our model. Although temozolamide has become a standard part of glioblastoma treatment, how to incorporate the effects of temozolamide into the model remains unclear. Initial results by Barazzuol et al.20 suggest that temozolamide enhances radiation sensitivity, and that this can be modeled as a decrease in the α/β ratio. If so, that can be easily incorporated into the model, but, again, will require close correlation to empirical data. It is hoped that this will be the focus of a future project. However, it is also noted that evidence suggests that temozolamide is less effective for the subset of GBM with unmethylated MGMT promoter, and, thus, the impact of excluding temozolamide from the current treatment model may not adversely affect its predictive capacity over a large class of tumors.32

Overall, our very preliminary results suggest that at least for tumors with moderate proliferation indices, there is no impact on the threshold where EOR begins to improve survival, but that the contribution of radiation therapy in relationship to EOR towards extending longevity becomes more predominant. Finally, as genetic profiling of GBM becomes more precisely correlated with phenomena such as invasiveness and proliferation, additional modifications of our model may further improve its accuracy, achieving greater ability to tailor the model to individual patients with further potential optimization of treatment decisions.

Disclosure Statement

The authors have no conflict of interest to declare.