Prognostic Significance of Growth Kinetics in Newly Diagnosed Glioblastomas Revealed by Combining Serial Imaging with a Novel Biomathematical Model

Glioblastomas are the most aggressive primary brain tumors, characterized by their rapid and diffuse infiltration of the adjacent normal-appearing brain tissue. Several clinicopathologic factors have been found to allow stratification of glioma patients according to probable survival times: (a) for the tumor: type and grade, size, and site; (b) for the patient: age and neurologic functioning, and history of previous treatment (extent of surgery, X-ray irradiation, chemotherapy; ref. 16). None of these factors concerns patient-specific measures of the in vivo kinetics of the tumor itself—partially a result of the lack of tools to quantify overall growth kinetics of a tumor capable of both extensive invasion peripheral to the imaging abnormality and aggressive rates of cellular proliferation. We have recently developed a novel biomathematical model for glioma growth and invasion that allows the translation of routinely available, serial pretreatment magnetic resonance images (MRI) into patient-specific net rates of proliferation and dispersal of glioma cells (5, 17). It is the purpose of this presentation to show statistically that these two new measures of biological aggressiveness offer prognostic information independent of the classic clinicopathologic features and thus allow a more accurate prediction of the duration of survival of glioblastoma patients.

The cell density gradient implicit to these invasive tumors suggests that imaging reveals only the “tip of the iceberg” of each glioma, with a significant portion of the glioma cells already invaded peripheral to the imaging abnormalities as shown in Fig. 1A. The biomathematical model, combining the two key net rates, proliferation (ρ) and invasion (D), has accurately quantified the diffuse invasion as well as the overall growth kinetics of gliomas (5, 17). Specifically, based on this model, we have developed a technique for translating changes seen on serial MRIs pretreatment into measures of net proliferation and net invasion of glioma cells (17). Figure 1B shows the model-predicted diffuse gradient of glioma cells predicted by simulation of the biomathematical model with patient-specific model parameters (D and ρ) shown in Table 1 for patient 11. The color gradient in Fig. 1B represents the model-simulation of glioma cellular density and reveals the diffuse nature of the lesion peripheral to the nodule seen on T1Gd MRI (dark gray contour). The ratio ρ/D can be interpreted as the relative steepness of the gradient of invading cells peripheral to the frank abnormality seen on T1Gd imaging (gross pathology or histology). Alternatively, the inverse of this ratio (D/ρ) has been termed the “invisibility index” and can be used to estimate the number of the diffusely invaded tumor cells peripheral to (and therefore invisible to) the imaging abnormality (10). Autopsy studies have revealed the remarkable accuracy of these model-predicted gradients of invasion in other glioblastoma cases for which autopsy was available (13).

In addition to quantifying the diffuse extent of gliomas, a novel prediction of the model's combination of an exponential net proliferation rate (for low cell densities) and a diffuse invasion pattern is an asymptotic overall linear pattern of radial expansion of the detectable abnormality on imaging, dependent on both dispersal, D, and proliferation, ρ, with an overall velocity defined by Fisher's approximation: $v=2D⁢ρ$⁠. This consistent growth pattern has been validated in a variety of gliomas (13, 14) and has already been shown to be a prognostic indicator in low-grade gliomas (18, 19). In addition to relating the overall velocity of radial expansion on MRI to a net combination of the invasiveness and net proliferation of the glioma cells, the individual parameter values, D and ρ, can be calculated from serial pretreatment imaging for each individual contrast-enhancing glioma (5, 13, 20, 21). These parameters are characteristic for each patient's tumor and allow for a dynamic characterization of disease response to any given treatment by providing an untreated virtual control (UVC) against which to compare actual tumor behavior. We will refer to this ratio as the therapeutic response index (TRI), which is described as a fraction of the time to fatal tumor burden (FTB).

Our hypothesis is that our biomathematical model–based net measures of proliferation and invasion rates provide novel insight into the overall glioma phenotype for each patient in vivo and are thereby prognostically significant. In this study, we considered 32 glioblastoma patients with serial pretreatment imaging to allow quantification of a net proliferation rate and a net invasion rate for each glioma. Using proportional hazard survival analyses, we first assessed the prognostic significance of these biomathematical model parameters relative to standard clinicopathologic parameters. Next, using the model-predicted survival of each patient's UVC, we assessed the net increase in actual survival (with treatment) for each patient. Complementary to the proportional hazard survival analysis, we found that those patients with the highest net proliferation rate were most likely to benefit significantly from therapy, consistent with our recent finding about radiotherapy (22). That some glioblastomas do not respond significantly to cytotoxic therapies merely confirms the wide range of growth rates, extending into the low ranges, where X-ray therapy (XRT) and chemotherapy are generally thought be less effective; the only problem with that thought being that net proliferation rates of gliomas could not previously be accurately measured.

The two main criteria for inclusion in the patient population were (a) a diagnosis of WHO grade 4 glioblastoma and (b) at least two gadolinium-enhanced T1-weighted (T1Gd) and T2-weighted (T2) MRIs before any surgical or treatment intervention to calculate the velocity of radial expansion and the biomathematical model parameters, D and ρ. To control for the role of heterogeneity in the primary treatment as it affects survival, we accepted only patients receiving standard conformal radiation therapy to a dose of at least 48 Gy (one patient; for others, 59.4–69 Gy). Of the 32 patients identified, 17 received concurrent radiation and temozolomide chemotherapy as their primary therapy (summarized in Table 1). Twenty-one of the patients were diagnosed and treated at the University of Washington Medical Center, six at the SFC Brain Imaging Research Centre at the University of Edinburgh (20), and the remaining five at the University of California at Los Angeles Medical Center. Patient data were accrued following Institutional Review Board–approved protocols.

As outlined in Table 1, 25 of the 32 patients had Karnofsky performance status (KPS) scores at diagnosis ranging from 60 to 100 (median, 80), and 32 had ages ranging from 22 to 76 (median, 55); 22 received either a biopsy or subtotal resection (STR) and the other 10 received a gross total resection (GTR). The Radiation Therapy Oncology Group (RTOG) recursive partitioning analysis (RPA) ranged from class III (17.1 mo expected median survival) to class IV (11.2 mo) and class V (7.5 mo; ref. 23). Due to the fairly small number of patients in each RPA class, the survival analysis was done treating RPA as a dichotomous variable, dividing at RPA III and IV (n = 19) versus RPA V (n = 10).

Model parameters for tumor growth kinetics (D and ρ) were calculated from two pretreatment MRIs as reviewed by Harpold and colleagues (10) and applied by Swanson and colleagues (5, 13, 21) to defined tumor volumes seen on T1Gd and T2. Velocities of radial expansion were computed as the average linear rate of increase in mean radius between two pretreatment observations on T1Gd. The ratio ρ/D was determined from comparison of the tumor volumes seen on T1Gd and on T2 at a single time point. Details of the sensitivity of model parameters to uncertainty in the magnetic resonance tumor volume measurements are provided in Supplementary text and Supplementary Table S1.

Following the protocol described in ref. 19, the untreated survival time for the UVC was estimated by developing the concept of a fatal tumor burden (FTB), either a size (35 mm radius) of the “solid tumor” as seen on T1Gd MRI or a total number of cells visible or invisible (1.1 × 10¹¹; ref. 20)—see Table 1, Supplementary Table S1, and Supplementary Fig. S1. The Munro-Kellie hypothesis (24) explains why brain tumors might actually have a FTB with a fixed volume of the skull. Both FTB end points yield practically the same answer (see Supplementary Fig. S1) and lead naturally to the TRI. The TRI is the ratio of the actual survival of the patient to the model-predicted survival of that patient's UVC.

Survival analysis was done using Cox proportional hazard regression (25) implemented in S-PLUS (26). For patients still alive, standard censoring procedures were taken using survival in days from the date of operation or biopsy. Significance was determined by log-rank P values.

Table 1 summarizes the pretreatment clinical characteristics and treatment details of the patients as well as the MRI measurements and related biomathematical model parameter estimates of proliferation and migration kinetics for each patient characterizing their UVC. The median net proliferation rate, ρ, of 15/y corresponds to a net cellular doubling time of 17 days and is consistent with reported estimates (27). The net migration rates, D (3.2–468.9 mm²/d; median, 29.7 mm²/d), are similar to the random motility rate estimates of glioma cells observed under a variety of experimental conditions (9).

The results of univariate Cox proportional hazard analysis are summarized in Table 2. In univariate analysis, the only significant predictors found were the UVC estimates for tumor net proliferation and invasion kinetics, ρ and ρ/D (P < 0.04 and P < 0.02)_. Among the standard clinicopathologic factors, RPA was the closest to reaching significance in univariate analysis (P < 0.06).

Although there are a variety of multivariate models to consider, given the limited size of our population, we retained only the terms that were closest to significance in the univariate analysis. To control for these known prognostic factors, multivariate survival analysis (Table 3) was done controlling specifically for RPA classification and the model-defined UVC measures. Consistent with the finding of Stupp and colleagues (28) showing improved survival (17 patients who received temozolomide concurrent with primary radiotherapy), a chemotherapy flag was added to the multivariate analysis (28). Even when controlling for these known clinical factors in multivariate analysis, the biomathematical model parameters for tumor kinetics, ρ and ρ/D, remained significant (P < 0.02), whereas neither RPA nor the concurrent chemotherapy flag attained further significance. The log hazard ratio coefficients resulting from the multivariate Cox model indicate that a 19/y increase in ρ (equivalent to a decrease of 13 days in the glioma cell doubling time) or a 1.3/mm² increase in ρ/D would be associated with a 50% decrease in survival. Comparison between the ρ and ρ/D multivariate models (in terms of the logarithm of the partial likelihood) shows no clear advantage for either measure of glioma growth kinetics. Further analysis is required to reveal the relative importance of these two means of quantifying glioma aggressiveness.

Next, we asked the question: Does the baseline untreated growth and invasion kinetics of each patient's UVC differentiate patients that will survive longer or shorter than that would be expected by their RPA classification (Fig. 2)? We used the RPA-associated median survival times to divide each group into two subgroups, short and long survivors relative to their RPA-associated median survival. A comparison of the ratio of the patient's actual survival time to the median survival associated with each patient's RPA classification compared with each patient's UVC velocity (Fig. 2A) or proliferation rate ρ (Fig. 2B) reveals a striking pattern, suggesting that only tumors with low velocity of radial expansion on T1Gd (v < 20 mm/y in Fig. 2A) or low net proliferation rate (ρ < 10/y in Fig. 2B) survive longer than the median RPA prognosis. Further, only patients with high velocity (v > 20 mm/y in Fig. 2A) or high net proliferation rate (ρ > 10/y in Fig. 2B) survived less than the median RPA prognosis.

Each patient's UVC also allows for estimation of a model-predicted untreated survival time for the UVC. This allows creation of a TRI defined as the ratio of the patient's actual survival to the survival of their UVC defined as the time required to attain an appropriate fatal tumor burden (Table 1). Using our TRI, we begin our analysis of individual patients. A quick glance at Fig. 3A suggests a random distribution of the actual (treated) survival versus predicted (untreated) survival of the UVCs, with 8 patients practically on the line of equality of predicted and actual (TRI = 1.0), 6 below the line, and 18 above the line. By contrast, Fig. 2C and D shows that those tumors that are most likely to have a TRI significantly greater than 1 are those with a high velocity (v > 30 mm/y in Fig. 2C) and an elevated net proliferation rate (ρ > 10/y in Fig. 2D).

A closer analysis of the individual results allows estimates of the specificities and sensitivities of potential predictors of good response to therapy (TRI > 1.5, those surviving 50% longer then their virtual control if untreated). These estimates are plotted in Fig. 3B. That is, we hypothesize that because both radiotherapy and chemotherapy predominantly target proliferating cells, those tumors with high net rates of proliferation (high ρ) will have the largest TRI (largest deviation from their UVC), although they may not survive the longest overall because tumors with high ρ are aggressive overall (Tables 2 and 3). A similar analysis was done with patients with KPS with a cutoff of 70, RPA with a cutoff of V, GTR versus biopsy/STR, and whether the patient received the Stupp protocol of XRT concurrent with temozolomide chemotherapy as their primary therapy. This analysis (Fig. 3B) reveals that using only a single parameter, ρ with a cutoff of 14/y, is the most sensitive and specific in predicting those patients who will have a TRI >1.5. We explored the benefit of combining this threshold for ρ with the standard clinical parameters [KPS, RPA, extent of resection (EOR), and Stupp protocol] to find any synergistic effect on combining these markers. Figure 3B shows that the most favorable combination is a threshold of 14/y for ρ and RPA < V.

The real world can be divided into two unequal populations: a relatively large number of patients and a relatively small number of therapists. Most likely, from time to time, many individuals in both groups probably ask if there are any controls that have taken this treatment and gotten better and other controls who have gotten worse. Both “controls” being exactly like the patient about to be treated, not just an average individual who might not even resemble the present patient as to age, gender, etc., and even more importantly, with exactly the same disease, at the same location and at the same stage, etc., wouldn't it be nice to have an UVC (not a real live twin who could not be ethically treated or not)? Such has been the target of our research for the last several years: to develop such a control for each of our patients with a glioma, especially a glioblastoma, and to develop a virtual control that has all the necessary and sufficient characteristics of the glioma that could be treated or not to compare with the real patient at hand. We believe we have such a model and have given it a trial in the present article: Thirty-two real glioblastomas each matched exactly with those characteristics that should repeat the real course from its origin through diagnosis and treatment to death. In summary, all we need are two sets of MRIs (pretreatment) from which to measure the velocity of radial expansion and to calculate the net rates of dispersal and proliferation of the UVC.

Because the concept of a virtual control (our UVC) is new, perhaps even unique to our investigative group, we began the presentation of our results with classic statistics with which the reader might be more familiar to show that our criteria have some comparability in the standard scheme of survival analysis as well as having some value in searching for other interesting details. These details on the evaluation of the individual patients can be approached by our novel biomathematical model, but not by ordinary statistical approaches.

Stratification of glioblastoma patients according to probable prognosis has relied predominantly on clinicopathologic factors including age, neurologic status, KPS, EOR (29–32), and RPA (23). The relatively few studies that have sought to estimate the proliferative aggressiveness of glioblastomas in vivo have focused on tissue-based histologic analysis, such as MIB-1 labeling index (33, 34), or novel imaging agents, such as [¹⁸F]fluorothymidine positron emission tomography (PET; refs. 35–37). Histologic estimates for proliferation rate with the Ki-67 (MIB-1) labeling index, bromodeoxyuridine (BrdUrd) labeling index, and flow cytometry have predominantly suggested differences in these measures across grades of gliomas but not necessarily among WHO grade 4 glioblastomas (33, 34, 36). Imaging cellular proliferation using fluorothymidine PET (36) has been found to result in prognostically significant estimates of proliferation across grades of gliomas. In addition to their limited application to glioblastomas, each of these measures is difficult to quantify in individual cases because they are highly spatially variable and representative only of the positive component (mitosis) of the proliferation rate (ρ) without counting the negative component resulting from apoptosis or other cell death mechanisms. Further, none of these measures directly considers the role of invasiveness in the overall growth pattern.

Based on injections of i.v. iododeoxyuridine and BrdUrd into patients before tumor resection, Hoshino and colleagues (27) reported a potential doubling time of gliomas ranging from 2 days to 1 month. Our patient-specific estimates of the net proliferation rate ρ corresponding to doubling times ranging from 1 to 87 days with a median of 17 days are thus consistent with these ranges. Furthermore, we have shown a similar agreement between our patient-specific estimates of D and that estimated from in vitro data (9).

We used proportional hazard analysis of newly diagnosed glioblastoma patients considering the prognostic effect of tumor measures such as MRI-detectable tumor radius, tumor radial expansion velocity, and the model-defined net rates of dispersal and proliferation (10, 18–21). We found that, in our patient population, model parameters for these rates of net dispersal and proliferation of glioma cells were strong predictors of survival, with only the RTOG RPA classification (among the standard clinicopathologic parameters considered) nearing significance (a marginal P value of <0.06) in univariate analysis. Multivariate analysis reveals that the significance of the biomathematical model biological aggressiveness parameters remains even when controlling for age, RPA classification, as well as the existence of concurrent chemotherapy during primary radiotherapy. Further study on a larger patient population is required to enforce this significance and any role these measures may play in the routine clinical management of glioblastoma patients.

For this patient population, the net invasiveness, D, only reaches prognostic significance when considered in combination with ρ. The ratio ρ/D gives a measure not only of the steepness of the gradient of the invading cells peripheral to the T1Gd edge but also of the net proliferation rate normalized to the net diffuse invasion, which has been shown to be associated with increased hypoxic burden in glioblastomas (38). Thus, tumors with high ρ/D would be expected to show not only rather distinct borders on histology but also more necrosis associated with hypoxia due to the crowding and other microenvironmental constraints inherent to a highly proliferative tumor with a relatively low invasion. The present analysis suggests that prognostically significant measures of biological aggressiveness can be quantified from routinely available pretreatment MR imaging of glioblastomas.

It is difficult to tease out subpopulations of patients with varying degrees of responsiveness in standard clinical study approaches. We have presented a method to predict the patient's course without any treatment and to compare how the predicted result compares with the actual outcome. Our hypothesis is that ρ and D are sufficient to create such an UVC, but as we gain more experience in making such comparisons, we are willing to admit that there may be other factors that will improve the comparability of the real patient with the virtual control. At the current state, we have a tool that we believe represents an advance over the current techniques.

Following the vein of our early work (1, 2, 15, 39, 40), a number of model extensions and novel approaches to modeling gliomas have been explored by us (41) and others (42–45). Most others have focused on modeling the experimental setting in which, for instance, tumor spheroid data are analyzed (43, 46), but some have considered realistic patient geometries (42, 44, 47). Most of these are impractical for direct application to individual patients because of the need to estimate numerous patient-specific parameter values from scant data (42, 44), and none have applied their techniques to any significantly sized patient data sets to yield patient-specific predictions that can be tested, validated, and used to improve the care of glioma patients.

Other than our recent report on a small subset of eight of these patients (20), to our knowledge, this is the first time in which routine clinical MRIs pretreatment have been translated to prognostically significant net rates of glioblastoma proliferation and invasion. Such a novel tool for quantifying glioma growth kinetics has numerous potential applications focusing on the creation of an UVC for each glioma against which treatment effects in each patient can be measured (6, 48, 49) and in which novel treatment approaches can be optimized.

We have developed a novel tool to potentially identify each individual patient who did or did not respond to treatment (through the TRI). Because all patients received XRT, what radiotherapist would like to learn is whether his/her therapy (even with all the others combined) did or did not work in particular patients who can now be specifically identified; whether the addition of chemotherapy did or did not matter very much; and whether high ρ was associated with longer survival and low ρ with only average survival. If the treatment(s) cannot provide a longer survival then the UVC, then that treatment can be discarded or retained to better suit the dynamics of the tumor growth. These general finding are consistent with our more detailed analysis of an extended model including a term to quantify the effect of radiation therapy and patient-specific radiation dose plans (48, 49).

This study highlights the novel role for a biomathematical model that takes into account both net tumor cell proliferation and net tumor cell invasion to allow for translation of information provided by routine serial (pretreatment) MR imaging to overall tumor growth kinetics or phenotype in vivo. Thus, dynamic insight from routinely available pretreatment imaging may be useful both in predicting and in assessing individual patients. Further, these results suggest that these measures of kinetics are important for the accurate prediction of disease course independent of the routine clinical parameters currently relied upon.

No potential conflicts of interest were disclosed.

Grant support: K.R. Swanson gratefully acknowledges the timely and generous support of the McDonnell Foundation and the Academic Pathology Fund.

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.