Abstract
Glioblastoma multiforme (GBM) is a highly invasive primary brain tumour that has poor prognosis despite aggressive treatment. A hallmark of these tumours is diffuse invasion into the surrounding brain, necessitating a multimodal treatment approach, including surgery, radiation and chemotherapy. We have previously demonstrated the ability of our model to predict radiographic response immediately following radiation therapy in individual GBM patients using a simplified geometry of the brain and theoretical radiation dose. Using only two pretreatment magnetic resonance imaging scans, we calculate net rates of proliferation and invasion as well as radiation sensitivity for a patient's disease. Here, we present the application of our clinically targeted modelling approach to a single glioblastoma patient as a demonstration of our method. We apply our model in the full threedimensional architecture of the brain to quantify the effects of regional resistance to radiation owing to hypoxia in vivo determined by [^{18}F]fluoromisonidazole positron emission tomography (FMISOPET) and the patientspecific threedimensional radiation treatment plan. Incorporation of hypoxia into our model with FMISOPET increases the model–data agreement by an order of magnitude. This improvement was robust to our definition of hypoxia or the degree of radiation resistance quantified with the FMISOPET image and our computational model, respectively. This work demonstrates a useful application of patientspecific modelling in personalized medicine and how mathematical modelling has the potential to unify multimodality imaging and radiation treatment planning.
1. Introduction
Glioblastoma multiforme (GBM) is the most aggressive primary brain tumour and accounts for the majority of primary brain tumours [1]. Following diagnosis, GBM is often treated with surgical intervention followed by concurrent radiation and chemotherapy [2]. GBM is characterized by invasive tumour cells that can be found as far away as 4 cm from the tumour mass [3]. This tumour cell invasion is not revealed by magnetic resonance imaging (MRI), the principal means of monitoring GBM progression and response to therapy [4].
1.1. Routine magnetic resonance imaging for glioblastoma
Clinically, glioblastoma progression and response to therapy are monitored with MRI [4,5]. For the past two decades, routine clinical imaging protocols for glioblastoma have consisted of T1weighted and T2weighted MRI sequences, with the addition of the gadolinium contrast agent on T1 MRI (T1Gd) as well as fluid attenuation inversion recovery (FLAIR) based on the T2 MRI [6]. Regions of hyperintensity on the T1Gd MRI correlate with high tumour cell density and ‘bulk tumour’. In contrast, the T2 and FLAIR images reflect vasogenic oedema typically associated with inflammatory response to infiltrating tumour cells, at a much lower density than found on the T1Gd region.
Glioblastoma is a particularly hypoxic neoplasm and is, in part, histologically defined by the presence of endothelial proliferation, angiogenesis and necrosis. The cascade of events which initiate and propagate angiogenesis may involve both acute and chronic hypoxic events leading to a heterogeneously hypoxic neoplasm [7]. Although the human brain is oxygen rich in its native state, glioblastoma tumour cells consume oxygen through cooption of existing vasculature as well as stimulation of new vasculature [7]. Increased hypoxia is a hallmark of aggressive tumour growth [8–10] and is known to reduce the efficacy of radiation therapy (RT) [11] and is negatively correlated with prognosis, although this is debated [12]. In order to investigate the role of hypoxia in mediating response to RT, we studied an MRIbased patientspecific computational model of tumour growth and response to RT and combined it with hypoxia determined with [^{18}F]fluoromisonidazole (FMISO) positron emission tomography (PET) [13,14].
Our patientspecific computational model quantifies a prediction of response to RT using in vivo MRI data. Patientspecific growth rates are quantified by net rates of proliferation and invasion and are calculated using routinely available MRI obtained prior to treatment. Response to RT can be described by the linearquadratic dose–response model [15]. The motivation for this patientspecific computational modelling is the need for more quantitative and individualized medicine that unifies molecular and anatomical imaging modalities as well as incorporates tumour growth and response rates, which have been shown to vary from patient to patient [16–19]. We present an analysis for one GBM patient with MRI and FMISOPET prior to RT to document our method and demonstrate an approach to quantifying hypoxiamediated resistance to RT using patientspecific computational modelling.
2. Material and methods
2.1. Glioblastoma patient case study
We study a single glioblastoma patient with two MRIs and an FMISOPET study prior to RT. The patient received the standardofcare chemoRT following diagnostic biopsy and was followed serially with MRI throughout the disease course. Tumour size data were recorded prior to and following RT to be compared with predictions of a patientspecific computational model. The model was used to investigate and quantify the role of regional hypoxia in determining MRIdefined response to RT.
2.2. Natural history and diagnosis
The subject of this study is a 73yearold man who provided informed consent to participate in an observational study approved by the local institutional review board (IRB). An MRI obtained at the time of presentation demonstrated a left temporal lobe lesion surrounded by oedema. A needle biopsy via a bur hole procedure was performed of the left temporal lobe lesion. Multiple tissue sections revealed areas of glial neoplasm as well as areas of necrosis with associated pseudopalisading of neoplastic nuclei. On the basis of these pathological findings, the biopsy was most consistent with WHO grade IV GBM [7].
2.3. Imaging data
The patient's diagnostic and preoperative MRIs were performed 13 days apart. Three days after the biopsy procedure and 2 days prior to the first fraction of RT, the patient underwent an [^{18}F]FMISOPET study on an IRBapproved research protocol. MRI and PET protocols can be found in the electronic supplementary material. The patient underwent an MRI study 3 days following the completion of RT, and subsequent images were taken at twomonth intervals. The MRI and PET images were spatially coaligned to a common coordinate system using a rigid body transformation to the BrainWeb phantom [20] using the PFUS package within the PMOD software [21] and statistical parametric mapping, available through the Matlab software suite [22].
2.4. Tumour volume data
Tumour volumes were measured for the T1Gd and T2 sequences for MRI studies (table 1) using a semiautomated thresholdbased pixel intensity background subtraction software developed in Matlab. The accuracy and reproducibility of this method is comparable to manual tumour delineation [23]. Specifically, tumour volume V (cm^{3}) was calculated numerically using the formulawhere the image series is composed of N twodimensional slices, acquired in the axial plane, and where slice i has thickness dz_{i} and tumour area A_{i}. The tumour area A_{i} is determined by summing the number of pixels containing tumour on slice i. The tumour growth velocity, shown in table 1, was calculated based on changes in radius of a sphere of equivalent volume measured prior to biopsy on the T1Gd images.
2.5. Quantifying hypoxia in vivo with FMISOPET
The coregistered FMISOPET images were scaled to the average venous blood concentration of FMISO activity (see the electronic supplementary material for details) to produce a tumour/blood (T/B) ratio image seen in figure 1 [24]. A T/B ratio greater than or equal to 1.2 (T/B ≥ 1.2) was associated with regions of hypoxia and used to determine the total hypoxic volume (HV; table 2) [24]. PET imaging is inherently noisy, and there can be isolated voxels of FMISO uptake scattered throughout the brain; defining HV as T/B ≥ 1.2 and restricting the FMISO signal to the region of T2weighted MRI abnormality largely excludes this noise and isolates the FMISO signal to the tumour area [25]. The HV is distributed within the tumour, from the bulk tumour mass (T1Gd) to the invasive edge and tumour periphery, defined by the T2/FLAIR abnormality with a 2 cm margin (T2^{+}; table 2). The FMISO T/B values in the whole brain have a mean and standard deviation of 0.813 ± 0.223 and are not normally distributed (onesample Kolmogorov–Smirnov test p < 0.0001).
2.6. Radiation treatment
Five days after the diagnostic biopsy, the patient was treated with RT delivered with a threedimensional conformal treatment plan using 6 MV photons. The target volumes were defined as the T2defined abnormality with a 2.5 cm margin, which had a planned target dose of 54 Gy delivered in 30 daily fractions of 1.8 Gy per fraction. An additional dose was delivered to a smaller volume defined by the gadoliniumenhancing region plus a 2 cm margin of 7.2 Gy delivered in an additional four daily fractions of 1.8 Gy per fraction for a total of 61.2 Gy to this region. The target volumes are labelled as isodose curves of 61.2 and 54 Gy on the planning MRI in figure 1. Radiation was delivered with concurrent temozolomide (TMZ) 75 mg m^{−2} given daily during course of RT and continued adjuvantly for a total of 14 months [2].
2.7. A mathematical model of glioblastoma growth and response to radiation therapy
2.1and2.2The proliferation–invasion radiation therapy (PIRT) model is a partial differential equation that quantifies the spatial and temporal rates of change of glioblastoma cell density and incorporates the delivery and effect of RT. The model describes tumour cell density, denoted c = c(x, t) at time t and location x = (x,y,z) in units cells per mm^{3}, in terms of diffuse invasion and densitydependent logistic growth. Logistic growth relates the per capita growth rate to available space for the cells to grow, so that if there are few cells in a unit volume of tissue, the overall growth rate is higher than if there are many cells per unit volume. The maximum number of tumour cells that can fit in a cubic millimetre of tissue is known as the carrying capacity, denoted K, and is computed to be 1.91 × 10^{6} (cells mm^{−3}) assuming a 10 µm diameter tumour cell. Once the tumour cell density reaches the carrying capacity (c = K) at a particular spatial location, the tumour cells are spacerestricted and therefore do not proliferate. Because we do not incorporate the clearance of dead cells, when the tumour cell density reaches the tissue carrying capacity, we assume the tumour cells are dead or become quiescent in this region.
The PIRT model (equations (2.1) and (2.2)) quantifies glioblastoma growth in terms of two net rates: proliferation (ρ, per year) and invasion (D(x), mm^{2} per year). The invasion of malignant tumour cells into the brain parenchyma is influenced by the anatomy of the brain, codified in equation (2.2): tumour cells preferentially migrate along the myelinated axons of neurons composing the white matter and move more randomly and slowly through the dense grey matter which composes the cortical surface and some internal structures of the brain [26,27]. The last term in equation (2.1) represents the loss of tumour cells owing to RT and is based on the linearquadratic model and the clinical RT plan, discussed below and given in equations (2.3) and (2.4).
2.8. Differential motility
The invasive migration of tumour cells throughout the brain presents a challenge to understanding the true extent of the subclinical disease, as tumour dispersal speeds can vary up to 100fold between pioneering cells in white matter compared with the more random motion in the core of the tumour and in grey matter [7,27,28]. To model this behaviour, the tumour cell invasion rate (equation (2.2)) is a function of the spatial variable x, so that glioma cells migrate 100 times faster in the white matter than in the more dense grey matter, D_{w} = 100D_{g} [29,30]. The BrainWeb phantom was used for tissue classification, partitioning the brain into grey and white matter in addition to cerebrospinal fluid (CSF) on a 1 mm^{3} cubic grid. The phantom is used to define the invasion rate of tumour cells spatially in the brain D(x), figure 2 [20,31]. The simulated tumour is initiated as a single voxel of cells located at approximately the centre of mass of the T1Gddefined tumour volume. The model equations are solved using a numerical approach, with time and spatial grids determined to meet stability requirements (see the electronic supplementary material) [17].
2.9. Modelling the effect of radiation therapy on glioblastoma multiforme tumour cell density
2.3Equation (2.3) is the surviving fraction of cells (S) from the linearquadratic dose–response model. Equation (2.3) converts radiation dose in units of Gy defined in space and time (denoted dose(x,t)) as determined by the patient's clinical radiation treatment plan into a biological effect of each radiation treatment [11,15,32,33]. The linearquadratic model is the most widely used dose–response model in radiation oncology [11]. Firstorder repair kinetics are not included in this model because these effects are negligible at dose fractions of less than 2 Gy.2.4
The radiation loss term in the PIRT model (equation (2.1)) is nonlinear and follows the same densitydependent logistic formalism as the tumour growth model, so that the rate of tumour cell death is related to the rate of killing and the density of tumour cells. The probability of cell death (1 − S) from the linearquadratic model (equation (2.3)) determines the rate of cell killing during each treatment fraction time (equation (2.4)). For low cell densities, the effect of RT is linearly related to the fraction of tumour cells killed. This assumption is consistent with the common understanding that cells actively undergoing mitosis are more susceptible to DNA damage and are found at the rim and periphery of the tumour more than the dense core [11]. However, when the cell density is close to the carrying capacity of the tissue K, it is assumed that the effect saturates owing to increased interstitial pressure and decreased dose per cell. When the tumour cell density reaches carrying capacity (c = K), there is no radiation effect. This situation corresponds to an unresponsive necrotic core, where we assume that there is no radiationinduced cell killing, because the cells are already dead. A necrotic core is a histologic and radiographic hallmark of glioblastoma that develops when the tumour growth exceeds the tissue carrying capacity.
2.10. Patientspecific radiation sensitivity
In order to quantify radiation sensitivity on a patientspecific basis, we use the linearquadratic model with the ratio α/β fixed. This allows us to regard α as the single parameter to define radiation sensitivity. For each point in space and time, an effective dose and probability of cell survival can be calculated that corresponds uniquely to the individual patient's treatment plan and the linearquadratic dose–response model parameter α. Increasing α decreases the probability of cell survival, S, and therefore increases the probability of RTinduced cell death. Increasing values of α correspond to increasing treatment effect and deviation from untreated growth. With the ratio α/β fixed, the single parameter α can be uniquely determined using either the T1Gd or T2 postchemoradiation tumour size, using a bootstrap optimization technique, yielding a onetoone relationship between α and model prediction error, as described in Rockne et al. [17]. The first T2 MRI following chemoradiation was used to determine the radiation response parameter α owing to the localization of FMISOPET activity within the T2^{+} region (table 2).
2.11. Model of hypoxiamodulated radiation resistance
The ratio of the parameters α/β in the linearquadratic model provides a measure of tissue response to radiation exposure. In our treatment simulations, the ratio α/β is held constant at 10 Gy, which is a reasonable assumption for tumour tissue and consistent with previous work [11,25,34–36]. The relative contribution of the ratio α/β to radiation response is modulated by the presence of hypoxia. We use the scaling oxygen enhancement ratio (OER), which determines relative resistance to RT in regions of low oxygen [37–39]. The implementation of the OER to the linearquadratic model is to modify the radiobiological parameters α and β to be spatially defined, so that the OER is applied in regions of hypoxia, as follows:2.52.62.7The OER assumes values between one and three [40], corresponding to no hypoxic effects and maximum effect, respectively. Although the OER is often considered a function of oxygen tension measured in units of mm Hg, we associate the T/B FMISOPET image with a volume fraction of hypoxia and implement the OER in a binary, voxelwise manner as shown in equations (2.5)–(2.7).
2.12. Quantifying hypoxiamediated radiation resistance with the proliferation–invasion radiation therapy model
The PIRT model is used to investigate the role of FMISOPETdefined hypoxia in determining modelpredicted radiation response and to test the sensitivity of our definition of hypoxia in influencing radiographic response to RT. Each simulation was run using an anatomically accurate threedimensional brain phantom [20] and the patient's threedimensional conformal radiation dose prescription extracted from the treatment planning system (Philips Pinnacle). The following treatment scenarios are considered:

(1) Clinical radiation dose delivered with spatially uniform treatment response. This corresponds to a value for the OER equal to one (OER = 1) at all locations in the brain, and assumes that all tumour cells are equally susceptible to RT damage.

(2) Clinical radiation dose with localized radioresistance owing to hypoxia via a binary relationship. This quantifies focal radioresistance defined by the patient's FMISOPET scan, where the tissue to blood value is greater than or equal to 1.2 [24], defines the HV where the OER is modified. Four simulations were performed with the OER equal to 1.5, 2.0, 2.5 and 3.0 to characterize a range of hypoxiamediated responses. The spatially defined OER is given by2.8
2.13. A patientspecific proliferation–invasion radiotherapy model
To establish a histopathological correlation of cellular density to imaging, we associate enhancing features on T1Gd and T2weighted MRIs with high and low tumour cell densities, respectively [41,42]. One corollary of this estimate is the modelling of a diffuse gradient of tumour cells invisible to imaging. The relative proportion of occult disease is characterized by an invisibility index (the ratio of the model parameters D/ρ [43]) which has been inversely related to the volume of hypoxia within the tumour measured on FMISOPET relative to the overall tumour mass on MRI, so that more nodular tumours (small ratio D/ρ) are more hypoxic than diffuse tumours (large ratio D/ρ) [9]. This can be interpreted as the difference between vascular cooption and tumourdriven angiogenesis in diffuse versus nodular neoplasms, respectively.
The PIRT model (equation (2.1)) predicts a nearly linear radial growth of the abnormality seen on imaging, which approaches a constant velocity defined by Fisher's approximation [44,45]. A constant velocity of imageable growth on MRI has been demonstrated for 27 untreated lowgrade gliomas [46] computed from serial MRI observations. The invisibility index and velocity of growth are computed from gross tumour volumes and combined to calculate patientspecific net rates of diffusion (D) and proliferation (ρ) [17,42].
2.14. Patientspecific model parameters
Patientspecific model parameters D and ρ were calculated based on baseline tumour growth kinetics using prebiopsy MRI. Specifically, within the spectrum of dynamics observed in glioblastoma, the patient's net rates of invasion and proliferation that characterize the tumour growth lie within 1 standard deviation of the population mean observed in 63 patients [16]. Similarly for the patientspecific radiosensitivity parameter α, the calculated value of α = 0.055 (per Gy) reflects neither an exceptionally resistant nor sensitive response (table 3). The relative invasiveness or ‘invisibility index’ (D/ρ = 0.93 mm^{2}) reflects a tumour growth pattern which is balanced between cooptive and angiogenic vascularity and therefore predicts a modest hypoxic burden.
2.15. Spatial metrics of model accuracy
In order to assess the accuracy of model predictions, we compare the simulated and actual tumour regions on T1Gd and T2weighted MRI using spatially defined similarity metrics. Each metric returns a value indicating the quality of agreement between model and observed tumour growth. This analysis was performed on the two MRIs prior to RT and on the first MRI following RT on both the T1Gd and T2 tumour volumes. Modelpredicted RT response is defined with the first MRI performed postRT [17]. Because we do not explicitly model chemotherapy, we cannot confidently apply this model or analysis to MRI observations beyond the first two postradiation scans, when adjuvant chemotherapy is often administered. Similarity metrics using algebraic combinations of truepositive (TP), falsepositive (FP), falsenegative (FN) and truenegative (TN) are computed on a voxelwise basis and include the positive predictive value = TP/(TP + FP), sensitivity = TP/(TP + FN), specificity = TN/(FP + TN), Jaccard index = TP/(FN + TP + FP) and volume similarity = 1 − (FPFN/(FP + 2TP + FN)). These metrics take values from 0 to 1, with one indicating exact agreement. Modelpredicted and observed tumour radii were also compared.
In addition to voxelwise metrics of concordance, we quantified morphological similarity by computing the distance between the predicted and observed tumour surfaces in three dimensions. This measure returns a distribution of distances, so that a value of zero in this distribution indicates intersection of the simulated and actual tumour surfaces. We report the median and standard deviation of this distribution—the closer to zero and smaller the variance, the better the model prediction. Reporting the median and standard deviation allows us to evaluate both variance and bias in our model predictions and avoid a cancelling effect of including both positive and negative distances that would be reflected in the mean of the distribution.
3. Results
Using our patientspecific model for glioma growth and response to RT, we find the incorporation of hypoxiamediated radiation resistance defined with FMISOPET leads to an order of magnitude decrease in relative volumetric error, from 14.6% to 1.1% (table 4). Incorporation of hypoxiamediated resistance provided better qualitative and quantitative predictive value to the model in regions of high cellular density where the hypoxia was localized, despite the relatively small volume of hypoxia within the tumour region, representing only 13% of the bulk tumour mass. Relative error on the T2weighted MRIbased tumour size was improved more modestly, from 0.5% to 0.2%. Absolute differences between volumebased tumour radius ranged from 0.04 to 2.63 mm. Interobserver error in gross tumour radius has been estimated as ±1 mm [47], indicating that simulation predictions are comparable to estimated uncertainty in measurable tumour size, summarized in table 4 and illustrated in figure 3. Additionally, when comparing the PIRT model simulation with and without incorporating the OER, 238% more tumour cells survived owing to local hypoxia resistance effects and 24% of the modelpredicted tumour cells prior to treatment were killed by the radiation, indicating a 10fold overestimate of RT effect without hypoxic resistance.
Spatial concordance between the modelpredicted disease burden prior to and following RT was measured through voxelwise similarity metrics. These metrics are challenged by the spatial alignment of the clinical MRI with the atlas and by a large number of true negatives for pixels far outside the gross tumour region. With this limitation in mind, we observed an increase in voxelwise concordance by incorporation of localized hypoxiamodulated radiation resistance for most metrics. Volume similarity was increased from 0.844 to 0.91, and the Jaccard metric was increased from 0.72 to 0.83. Sensitivity increased from 99.75 to 99.87, although specificity decreased from 30.25 to 15.61.
The distance between surfaces analysis was performed for all three MRI studies, two prior to biopsy and one immediately following the full course of RT on T1Gd and T2 sequences are illustrated in figure 4. The median distance (±standard deviation) between measured and predicted surfaces, not including overlapping voxels with a distance measure of zero, decreased on T1Gd from 2.8 ± 2.0 to 2.2 ± 2.2 mm and were unchanged on T2, from 2.8 (±3.4). These results indicate an improved model–data agreement in the tumour core (T1Gd) where the hypoxia is localized and a stable model prediction in the tumour periphery (T2), with an overall bias of overestimating the outward growth of the tumour.
3.1. Sensitivity analysis
To investigate the sensitivity of our predictions to definitions of hypoxic regions and parameter selection, we considered a range of OER values in the HV ranging from 1.5 to 3.0. This produced a parabolic relationship between relative volumetric error and OER using T1Gd size such that the minimum model–data error occurred at approximately OER = 2.5. OER values between 2.0 and 3.0 yielded relative volumetric errors at or below approximately 1%, and all resulted in an order of magnitude improvement to the model prediction. This analysis demonstrates the robustness of the model predictions to variations in OER, and by extension, the volume fraction of hypoxic cells the T/B FMISOPET image represents.
To further investigate the role of spatial localization of radiation resistance, simulations were also performed in isotropic threedimensional spherical symmetry without anatomical complexity, using the T1Gd isodensity as a lower bound threshold to determine the hypoxic region. In these simulations, we observed a similar outcome, that is, that model–data agreement is substantially improved with the inclusion of hypoxiamediated radiation resistance, where HV is within the bulk tumour region. However, the model predictions were not as accurate with this simplification. This suggests that the contrastenhancing region is not a suitable surrogate for hypoxic tumour, and that a more specific indicator such as FMISOPET in three dimensions not only improves the model prediction, but is also needed in a spatially heterogeneous tumour such as GBM. Sensitivity of modelpredicted tumour burden following RT based on variability in the tumour growth (D, ρ) or response rates (α, α/β) is discussed in detail in the supplement to Rockne et al. [17] and is not sufficient to account for the improvement in model–data agreement.
3.2. Limitations and considerations
This study focuses on one patient to document our method and demonstrate an approach to quantifying hypoxiamediated resistance to RT using patientspecific computational modelling. More patients must be studied in order to support our findings that incorporation of FMISOPETdefined hypoxia improves the predictions of the computational model.
As in our previous work [17], we assume the delivery and effect of RT to be an instantaneous, deterministic event using the linearquadratic equation and its corresponding probability of cell survival/death. Concurrent chemotherapy is assumed to be included in the net effect of RT and is not modelled explicitly. The survival benefit of adding TMZ chemotherapy to radiation treatment was demonstrated in the landmark study by Stupp et al. [2] and established a ubiquitous standard of care for newly diagnosed glioblastoma. However, it remains an open question of how best to translate the additive therapeutic benefit of TMZ into a mathematical model. The assumption that the effect of TMZ is included in the response rate parameter α does restrict the interpretation of the model predictions for this study.
To further study the additive effect of TMZ to radiation response, one could identify patient cohorts that received RT alone and compare the radiation sensitivity parameters in that cohort with those who received both RT and TMZ. A paired analysis could be performed to control for variations in the tumour growth rates D and ρ, as is done in the recent study by Adair et al. [48]. For intravenously delivered therapy such as TMZ, the blood–brain barrier will result in heterogeneous drug delivery spatially within the tumour. Advanced imaging such as dynamic contrast enhanced or perfusion MRI could be used to infer a drug concentration gradient within the tumour, and a similar mathematical formalism for the loss of tumour cells owing to RT could be used to model the effects of TMZ. A study of this kind would aid in the development of a TMZ model that could be used to study the effects of TMZ that is administered alone following chemoradiation.
The patient's stereotactic needle biopsy was not modelled, as the volume of tissue removed was at the core of the tumour and not a significant portion of the overall tumour mass. The inflammatory response of the tumour and normal appearing surrounding brain owing to surgery may impact MRI abnormalities, particularly on T2/FLAIR. We assume any such effects arising from only a needle biopsy did not impact response to RT, presentation of disease on followup imaging or tumour progression. Most glioblastoma patients receive extensive surgical removal of their lesion, which can be modelled by setting that portion of the computational domain and tumour model to zero.
This study assumes no changes to hypoxia through the course of radiotherapy, which does not include the effects of reoxygenation that could change the distribution of hypoxia within the tumour [49]. The linearquadratic model does not reflect myriad repair processes and microenvironmental changes induced by radiotherapy [50], although extensions to the L–Q model exist to approximate the effects of reoxygenation and hypoxia [51]. Hypoxia dynamics can be studied with kinetic FMISOPET imaging [52]; however, this approach remains relegated to the imaging time point. Moreover, we consider hypoxia a binary variable, which does not account for intravoxel heterogeneity. A more detailed mathematical model of the angiogenic process as it relates to hypoxia could be implemented to explore this heterogeneity, as proposed in [53,54].
Finally, this approach is driven by imaging modalities and does not account for molecular or genetic heterogeneity that is known to exist within glioblastoma, based upon data generated by the TCGA Research Network: http://cancergenome.nih.gov/ [55]. It has been shown that the molecular subtypes defined by the TCGA vary spatially within a single patient's tumour [56], suggesting that a single molecular subtype may not be sufficient to characterize a patient's disease. This emphasizes the need for targeted and imagelocalized biopsies within various regions of the tumour, as shown by Gill et al. [57]. To account for genetic and molecular heterogeneity influence in imageable response to treatment, one would need a multiscale, patientspecific model framework. Such an undertaking is not attempted here, although the authors recognize that molecular alterations likely play an important role in determining tumour evolution and response to treatment.
4. Discussion
We have investigated a computational model of human glioblastoma growth that incorporates hypoxiamediated resistance to RT based on FMISOPET in the complex architecture of the brain on a patientspecific basis. Incorporation of focal radiation resistance improves model–data agreement, measured with a variety of metrics, despite the small hypoxic burden (13%) of the tumour. The model predictions are improved by the incorporation of an OER to approximate the degree of radioresistance created by hypoxic conditions within the tumour. The large improvements in model–data agreement attributable to a modest volume of hypoxiamediated resistance to radiation effect underscores the significance of spatial heterogeneity in delivery and response to RT in glioblastoma and the complexity of a threedimensional model.
Without the incorporation of OER and localized radiation resistance, the model is unable to fully capture postRT tumour size in regions of high cellular density, motivating the addition of the OER parameter derived from in vivo clinical data. Ideally, a more patientspecific, biologically driven model of the process of angiogenesis, hypoxia and necrosis using a modelling approach which integrates multimodality imaging is desired [53,54]. In conjunction, prospective interventional imaging studies which capture changes in intratumoural hypoxia throughout the course of therapy are needed in order to improve models of hypoxiadriven resistance to RT to account for changes in the hypoxic state of the tumour. Moreover, daily fraction radiotherapy likely introduces phenotypic and genotypic selection pressures which eliminate the sensitive tumour cells and leave the most aggressive, resistant clones to repopulate following therapy completion [58].
The goal of modelling biological response to therapy has been long sought after, but there have been no demonstrable successes with the potential for impacting individual treatment planning [59]. Other efforts to incorporate biological effect into radiation treatment plans have either relied on static features of the tumour [60] or are not capable of being truly patientspecific because of the large number of parameters to be estimated [61–65]. We have developed a technique to incorporate proliferation, invasion and response to RT of the tumour over the time course of treatment within the threedimensional anatomy of the brain. Granting all the assumptions and limitations of the simple proliferation–invasion tumour growth model, this framework has provided a methodology to investigate the role of FMISOPETdefined hypoxia in modulating radiation response in vivo quantitatively. By providing quantitative metrics of a patient's response to radiotherapy, our model has the potential to unify multimodality imaging and treatment planning and establish a useful application of patientspecific modelling in personalized medicine.
Funding statement
We gratefully acknowledge the support of the James S. McDonnell Foundation, the Academic Pathology Fund and NIH grant nos. R01 NS060752, R01 CA164371 and P01 CA042045.
Acknowledgements
We thank and acknowledge Jeanne Link and Mark Muzi for the collection and processing of the PET study data. We also express profound gratitude to Dr Ellsworth C. Alvord Jr for his mentorship and guidance.
 Received October 24, 2014.
 Accepted December 5, 2014.
© 2014 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.