Personalized Radiotherapy Design for Glioblastoma: Integrating Mathematical Tumor Models, Multimodal Scans and Bayesian Inference

07/02/2018 ∙ by Jana Lipkova, et al. ∙ 0

Glioblastoma is a highly invasive brain tumor, whose cells infiltrate surrounding normal brain tissue beyond the lesion outlines visible in the current medical scans. These infiltrative cells are treated mainly by radiotherapy. Existing radiotherapy plans for brain tumors derive from population studies and scarcely account for patient-specific conditions. Here we provide a Bayesian machine learning framework for the rational design of improved, personalized radiotherapy plans using mathematical modeling and patient multimodal medical scans. Our method, for the first time, integrates complementary information from high resolution MRI scans and highly specific FET-PET metabolic maps to infer tumor cell density in glioblastoma patients. The Bayesian framework quantifies imaging and modeling uncertainties and predicts patient-specific tumor cell density with confidence intervals. The proposed methodology relies only on data acquired at a single time point and thus is applicable to standard clinical settings. An initial clinical population study shows that the radiotherapy plans generated from the inferred tumor cell infiltration maps spare more healthy tissue thereby reducing radiation toxicity while yielding comparable accuracy with standard radiotherapy protocols. Moreover, the inferred regions of high tumor cell densities coincide with the tumor radioresistant areas, providing guidance for personalized dose-escalation. The proposed integration of multimodal scans and mathematical modeling provides a robust, non-invasive tool to assist personalized radiotherapy design.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

Code Repositories

s3

Skull stripping and brain tissue segmentation


view repo
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

Glioblastoma (GBM) is the most aggressive and most common type of primary brain tumor, with a median survival of only 15 months despite intensive treatment [1]. The standard treatment consists of immediate tumor resection, followed by combined radio- and chemotherapy targeting the residual tumor. All treatment procedures are guided by magnetic resonance imaging (MRI). In contrast to most tumors, GBM infiltrates surrounding tissue, instead of forming a tumor with a well-defined boundary. The central tumor, which is visible on medical scans, is commonly resected. However, the distribution of the infiltrating residual tumor cells in the nearby healthy-appearing tissue, which are likely to contribute to tumor recurrence, is not known.
Current radiotherapy (RT) planning handles these uncertainties in a rather rudimentary fashion. Guided by population-level studies, standard-of-care RT plans uniformly irradiate the volume of the visible tumor extended by a uniform margin [1, 2, 3], which is referred as the clinical target volume (CTV). However, the extent of this margin varies by few centimeters even across the official RT guidelines [4]. Moreover, GBM infiltration is anisotropic and thus a uniform margin very likely does not provide an optimal dose distribution. In addition, GBM invasiveness is highly patient-specific, and thus not all patients benefit equally from the same margin, which impairs comparison and advancement of RT protocols.
Despite treatment almost all GBMs recur [5]. Biopsies [6] and post-mortem studies [7] show that tumor cells can invade beyond the CTV, which reduces RT efficiency, and is a possible cause of recurrence. At the same time radioresistance of tumor cells inside the CTV can also reduce RT efficiency. Radioresistance tends to occur in regions with complex microenviroment and hypoxia [8], both of which are commonly encountered in areas of high tumor cellularity. To address tumor radioresistance, several studies have suggested local dose-escalations [8, 9, 10]. In these approaches, a boosted dose is delivered into a single or multiple co-centered regions defined by adding uniform margins to the tumor outlines visible in MRI scans [11]. The Radiation Therapy Oncology Group (RTOG) phase-I-trial [10] showed an increase in median survival of 8 months with dose-escalation. However, no benefit in progression-free survival was observed, indicating a complex relationship between true progression of the underlying disease and the tumor extent visible in MRI scans.
Alternatively, positron emission tomography (PET) scans, which map tumor metabolic activities targeted by specific tracers, can be used to identify radioresistant regions. A promising tracer used in GBM imaging is 18F-fluoro-ethyl-tyrosine (FET) [12], whose uptake values have been shown to be proportional to tumor cell density, although the constant of proportionality is unknown and patient-specific [13, 14]. A prospective phase-II-study [5] demonstrated that dose-escalation based on FET-PET enhancement delineates the tumor structure better than uniform margins, thus leading to lower radiation toxicity. Still, FET-PET based dose-escalation did not increase progression-free survival. One possible explanation is that PET enhances mainly the tumor core, which is usually resected, while the PET uptake values in the remaining tumor periphery coincide with the baseline signal from the healthy tissue. This, together with a rather low resolution of PET scans, limits their ability to fully target radioresistant tumor residuals. This is also consistent with our results.

Standard RT plans can be improved by incorporating information from computational tumor models. These models, calibrated against patient medical scans, provide estimates of tumor infiltration that extend the information available in medical images and can guide personalized RT design. Despite extensive development of tumor growth models

[15, 16, 17, 18, 19, 20] and calibration strategies [21, 22, 23, 24, 25, 26, 27, 28], their translation into clinical practice remains very limited. We postulate that there are (at least) three translational weakness: 1) Most model calibrations rely on data not commonly available in clinical practice. For example, in [22, 23, 24, 25, 26, 27, 28] medical scans with visible tumor progression from at least two time points are used for the model calibration. However, for GBM patients only scans acquired at single preoperative time point are available. 2) Models are based on simplified assumptions motivated by in-vitro studies. For instance, it is frequently assumed that the tumor cell density is constant along the tumor borders visible on MRI scans (e.g.,[21, 22, 23, 24, 25, 26, 27]). However, the tumor cell density varies significantly along the visible lesion borders due to anatomical restrictions and anisotropic tumor growth. 3) Even if advanced calibration techniques as in [28] are used, it is not clear how robust the model predictions are and what benefits they offer over the standard treatment protocols.
Here, we address these translational issues and provide clinically relevant patient-specific tumor predictions to improve personalized RT design. We present a Bayesian machine learning framework to calibrate tumor growth models from multimodal medical scans. We show that an integration of information from complementary structural MRI and functional FET-PET metabolic maps enables the robust inference of the tumor cell densities from scans acquired at single time point. To the best of our knowledge, this is the first study making joint use of FET-PET and MRI scans for the patient-specific calibration of a tumor growth model. Our Bayesian approach infers modeling and imaging parameters under uncertainties arising from measurement and modeling errors. We propagate these uncertainties through the computational tumor model to obtain robust estimates of the tumor cell density together with credible intervals that can be used for personalized RT design. The patient-specific tumor estimates offer an advantage in determining margins of CTV as well as regions for dose-escalation. A clinical study is used to assess benefits of the personalized RT design over standard treatment protocols.
In the remainder of the paper, Section II introduces the Bayesian framework for model calibration, including the tumor growth and imaging models. The results are presented in Section III where the framework is applied to synthetic and clinical data, followed by a personalized RT study. Conclusions are presented in Section IV. Additional technical details are given in the Supplementary Materials (SM) available in http://ieeexplore.ieee.org.

Ii Bayesian model calibration

The Bayesian framework we develop combines a deterministic model for tumor growth with a stochastic imaging model

relating model predictions with tumor observations available from patient medical scans. Bayes theorem is used to estimate the probability distribution of the unknown parameters of both models, accounting for modeling and measurement uncertainties. Identified parametric uncertainties are propagated to obtain robust patient-specific tumor predictions. An overview of the framework is given in

Fig. 1.

Fig. 1: Overview of the inference framework. I) Medical images show preoperative patient scans: (I.A) 3D reconstruction of T1Gd images and (I.B) slices across all modalities. Tumor observations (e.g., segmentations) extracted from each modality are illustrated in (I.C). II) Computational medicine includes a tumor growth model (II.B), which simulates tumor evolution in the patient anatomy (II.A), and imaging model (II.E) which relates the image observations with the modeled tumor cell density u. Subplot (II.D) shows a schematic representation of the actual tumor cell density along the dashed line shown in (I.C) and its relation to each tumor observation available from the medical scans. The unknown, patient-specific parameters for each model are listed in Tables (II.C,F). III) Bayesian inference is used to identify the probability distributions of the unknown parameters, accounting for the modeling and measurement uncertainties. Parametric uncertainties are then propagated to obtain robust predictions about patient-specific tumor cell densities. The most probable tumor cell density, given by the maximum a posterior (MAP) estimate, is shown in (III.A,B)

, while the mean and standard deviation (std) are shown in

(III.B). IV) Personalized Radiotherapy uses the patient-specific predictions to improve dose-distribution and escalation design. A comparison of a standard and personalized dose distributions is shown in (IV.A-C). The regions of estimated high tumor cell densities are marked by the orange isocontour in (IV.D,F). Observed tumor recurrences in T1Gd and FLAIR, marked by the pink and yellow curves in (IV.E,F), are used to compare the treatment plans. The personalized plan spares more healthy tissue, while achieving tumor coverages comparable to the standard protocol and guides the design of personalized dose escalation plans.

Ii-a Tumor growth model

Many tumor growth models are based on the Fisher-Kolmogrov (FK) equation [23], which captures the main tumor behaviour: proliferation and infiltration. We use FK equation to describe the tumor model . The equation is solved in a patient-specific brain anatomy reconstructed from MRI scans, where each voxel corresponds to one simulation grid point. Let be the brain anatomy consisting of white and grey matter and be normalized tumor cell density at time and voxel at location , where is index across all voxels. The dynamics of the tumor cell density is modeled as:

(1)
(2)

The term

denotes proliferation rate. The tumor infiltration into the surrounding tissues is modeled by the tensor

where is a identity matrix and

(3)

The terms and denote percentage of white and grey matter at voxel , while and stand for tumor infiltration in the corresponding matter. We assume [28]. The skull and ventricles are not infiltrated by the tumor cells and act as a domain boundary with an imposed no-flux boundary condition Eq. 2, where is the outward unit normal to . The tumor is initialized at voxel and its growth is modeled from time until detection time . The parameters are considered unknown and patient-specific. The model is implemented in a 3D extension of the multi-resolution adapted grid solver [29] with a typical simulation time of 1-3 minutes using 2 cores. An overview of the model and its parameters is shown in Fig. 1 (II. A-C).

Ii-B Multimodal imaging model

We consider routinely acquired T1 gadolinium enhanced (T1Gd) and fluid attenuation inversion recovery (FLAIR) MRI scans in combination with FET-PET maps. A stochastic imaging model is designed to relate model predictions of the tumor cell density u and the tumor observations available from the medical scans. Here,

denotes a vector of tumor observations obtained from a certain image modality,

is an entry in and enumerates all voxels in the given image. A voxel corresponds to the same location in each scan and the simulation domain . See Fig. 1 for an overview of all the imaging modalities (I.B), corresponding tumor observations (I.C) and their relation to tumor cell density u (II.D).
The MRI scans provide morphological information about the visible tumor in the form of binary segmentations. The segmentation assigns a label to each voxel with visible tumor and otherwise. The probability of observing a segmentation with a simulated tumor cell density u

is modeled by a Bernoulli distribution 

[28]:

(4)

Here is the probability of observing the tumor in the MRI scan and it is assumed to be a double logistic sigmoid:

(5)

where denotes an unknown cell density threshold below which tumor cells are not visible in the MRI scan, while the term represents uncertainty in . The parameters are assumed unknown and patient-specific.
The FET-PET signal is proportional to tumor cell density with an unknown constant of proportionality [13, 14]. Let be the normalized FET-PET signal after subtracting the patient-specific baseline signal from healthy tissue i.e., and the corresponding constant of proportionality. We assume that can be related with the modeled tumor cell density as

(6)

where

is prediction error accounting for modeling and measurement uncertainties. Because of the noisy nature of the PET scan, the error term is assumed to be a normal distribution

. The probability of observing the PET signal with the simulated tumor cell density u is then modeled as

(7)

The PET scan, acquired at resolution, is registered to the MRI scans with resolution. To justify the product in Eq. 7 only voxels separated by distance are used. The parameters are unknown and patient-specific. An overview of the imaging model and its parameters is shown Fig. 1 (II. D-F).

Ii-C Parameters estimation and uncertainty propagation

The parameters of the model where , are assumed unknown and a probability distribution function (PDF) is used to quantify their plausible values. A prior PDF is used to incorporate any prior information about . Bayesian model calibration updates this prior information based on the available data . The updated posterior PDF is computed by the Bayes theorem:

(8)

where is the likelihood of observing data from the model for a given value of . Since each of the medical scans captures a different physiological process, the tumor observations are assumed independent and the likelihood function can be expressed as:

The prior PDF is assumed uniform with details specified in the SM. Since an analytical expression for Eq. 8 is not available, sampling algorithms are used to obtain samples from the posterior

. We use Transitional Markov Chain Monte Carlo (TMCMC) algorithm

[30] which iteratively constructs series of intermediate PDFs:

(9)

where and is a generation index. The term controls the convergence of the sampling procedure and is computed automatically by the TMCMC algorithm. TMCMC method constructs a large number of independent chains that explore parameter space more efficiently than traditional sampling methods [30] and allow parallel execution. We use a highly parallel implementation of the TMCMC algorithm provided by the 4U framework [31].
The inferred parametric uncertainties are propagated through the model to obtain robust predictions about u given by:

(10)

or by simplified measures such as the mean

and variance

derived from the first two moments

, :

where is the space of all unknown parameters. The most probable tumor cell density estimate, is given by the maximum a posteriori (MAP) defined as .

Iii Results

The Bayesian framework described in the previous section is first applied to synthetic data to test sensitivity of the inference and to show the role of multimodal image information on the model calibration. Afterwards, clinical data are used to infer patient-specific tumor cell densities and to design personalized RT plans. Tumor recurrence patterns are used to compare the proposed and standard RT plans. The software and data used in this paper are publicly available111https://github.com/JanaLipkova/GliomaSolver.

GT 0.1300 0.0250 302.00 80.64 171.52 128.00 0.0230 0.8800 0.7000 0.2500 -
MAP 0.1226 0.0296 271.13 80.61 171.55 127.92 0.0280 0.9884 0.8440 0.4874 0.0502
mean 0.1543 0.0361 253.85 80.64 171.47 127.80 0.0294 0.9775 0.8325 0.4841 0.0505
std 0.0530 0.0125 107.59 0.1280 0.1024 0.2048 0.0022 0.0139 0.0132 0.0076 0.0004
TABLE I: Results of the Bayesian calibration for the synthetic case generated with the ground truth (GT) values. Reported is maximum a posteriori (MAP), mean and standard deviation (std). The units are ; ; ; and , , .
Fig. 2: Synthetic test case. Orange box: A 2D slice of the synthetic ground truth (GT) tumor cell density (A) and corresponding image observations (B): the normalized FET-PET signal with additive noise (red-blue color scale) and the outlines of the T1Gd (yellow) and FLAIR (pink) binary tumor segmentations. Blue box: Results of the Bayesian calibration with multimodal data. The results are in close agreement with GT data. Green box: Calibration results using only the MRI data, which do not provide enough information to recover the tumor cell density profile correctly.
Fig. 3: The results of the Bayesian calibration for the synthetic case. Above the diagonal: Projection of the TMCMC samples of the posterior distribution in 2D space of the indicated parameters. The colors indicate likelihood values of the samples. The number in each plot shows the Pearson correlation coefficient between the parameter pairs. The colored triangles mark the four selected parameters used in Fig. 4. Diagonal: Marginal distributions obtained with Gaussian kernel estimates. Boxplot whiskers demarcate the 95% percentiles. Below the diagonal: Projected densities in 2D parameter space constructed by 2D Gaussian kernel estimates. The black dots mark the values used to generate the synthetic data.
Fig. 4: Insensitivity of the tumor cell density to the speed of the growth. Shown are slices of the tumor cell densities computed with different combinations of parameters as listed at the bottom of each plot. These correspond to the colored triangles in Fig. 3. Despite significant variation in the parameter values, all combinations lead to very similarly-appearing tumors. In the absence of temporal information, the time dependent parameters are not identifiable, since the model calibration cannot distinguish between compensating effects among the parameters that affect the dynamics. As shown here, tumors with similar and values appear very similar to one another (here , ). Hence, the Bayesian calibration identifies the probability distribution of all the plausible values.

Iii-a Sensitivity study

The model is used to generate a 3D synthetic tumor in a brain anatomy obtained from [32] using the parameters reported in Table I. A 2D slice of the simulated ’ground-truth’ (GT) tumor cell density is shown in Fig. 2 (A). The synthetic T1Gd and FLAIR tumor segmentations are constructed by thresholding the GT tumor cell density at and . The FET-PET signal is designed by taking the GT tumor cell density within the T1Gd and FLAIR segmentations, adding Gaussian noise with zero mean and standard deviation (std) , and normalizing the result. The value of is chosen as average std of the FET signal from the healthy brain tissue. The generated synthetic image observations are shown in Fig. 2 (B).
A sensitivity study for the number of samples is performed, indicating that 6000 samples is adequate for the model. The manifold of the inferred probability distribution is presented in Fig. 3 and the calibrated parameters are given in Table I. As seen from the probability distribution manifold, tumor observations from a single time point do not contain enough information to infer time dependent parameters exactly, since different combinations of these parameters can generate the same tumor cell density as shown in Fig. 4. The lack of identifiability of poses a challenge for calibration approaches searching only for a single value of . Instead, Bayesian calibration provides fairer estimate; the inferred probability distribution shows a strong correlation between the parameters , while the high std values imply low confidence in these parameters. On the other hand, parameters that affect the tumor spatial pattern, e.g. , are identified with high accuracy, which is reflected by their low std. The image related parameters are slightly overestimated due to the assumed correlation length and the effect of the complex brain anatomy. The role of the anatomy is discussed further in the SM.

The inferred parametric uncertainties are propagated to obtain robust posterior predictions about the tumor cell density shown in

Fig. 2 (C-E). Despite the large parametric uncertainties, the MAP and mean tumor cell density estimates are almost indistinguishable from the GT tumor. The low std values imply that, using our Bayesian formulation, the information contained in multimodal data is sufficient to infer tumor cell density from single time point scans.
For comparison, if the model calibration is performed only with the MRI data, i.e. , the estimated tumor cell densities shown in Fig. 2 (F-G) deviate from the GT tumor mainly in the central part of the lesion, which is also consistent with the regions of high std shown in Fig. 2 (H). Nonetheless, the outlines of the predicted tumor are similar to those of the GT tumor. This is because the tumor morphology is mainly constrained by the MRI data, since the FET-PET signal coincides with the baseline signal of the healthy tissue in the regions of lower tumor infiltration. On the other hand, the FET-PET signal constrains the tumor cell density profile in the regions of high tumor infiltration. This highlights the importance of integrating structural and functional image information for the model calibration when dealing with single time point data.

Fig. 5: Results of the Bayesian calibration for patients P1-P8. Orange box: Preoperative scans showing (A) T1Gd; (B) FLAIR; (C) Tumor segmentations: T1Gd (yellow) and FLAIR (pink); (D) FET-PET. Blue box: (E) MAP and (F) std of the inferred tumor cell densities shown in the preoperative FLAIR scans. The margin is shown as the green curves. (G) The 3D reconstructions show outlines of the MAP tumor (blue) together with tumor extent visible on the FET-PET scans (orange) in the preoperative anatomy (white). Green box: Scans of the first detected tumor recurrence. (H) FLAIR tumor recurrence (pink), (green), and (blue) margins. (I) T1Gd tumor recurrence (yellow) and the dose-escalation outlines proposed by FET enhancement (magenta) and MAP estimates (orange). (J) Multilevel dose-escalation designed by MAP estimates.

Iii-B Patient study

A retrospective clinical study is conducted on 8 patients diagnosed with GBM. Scans of the patients P1-P8 are shown Fig. 5 and the details about acquisition protocols and image processing are reported in the SM. All patients received the standard treatment, surgery followed by combined radio- and chemotherapy [2]. There was no visible tumor after the treatment and patients were regularly monitored for recurrence. The preoperative scans shown in Fig. 5 (A-D) are used for the Bayesian inference. The calibrated parameters are reported in Table 1 in the SM and the posterior, patient-specific predictions for the tumor cell densities are shown in Fig. 5 (E-G). These patient-specific predictions provide estimates about the possible tumor cell migration pathways in the surrounding of the visible tumor, constrained by the patient anatomy and the available tumor observations. The predicted tumor infiltration pathways can be validated by the patterns of the first detected tumor recurrence shown in Fig. 5 (H), where the outlines of the predicted infiltrations (blue) and recurrence tumors (pink) are depicted. For patients P5,P7,P8 the model accurately predicts tumor infiltration also inside the healthy-appearing collateral hemisphere, whereas for cases P1-P4 the tumor predictions are correctly restricted only to one hemisphere. Moreover, despite a similar appearance of the preoperative tumors in patients P1 and P2, the model correctly predicts more infiltrative behaviour for the patient P1, which is consistent with the recurrence pattern. The high confidence in the predictions is reflected by low std shown in Fig. 5 (F).

Iii-C Personalized radiotherapy design

The patient-specific tumor cell density predictions can be used to design margins of the CTV and to identify high cellularity regions that could mark areas of increased radioresistance. The personalized RT plan can be based either on the most probable scenario given by MAP estimate or the worst case scenario given as a sum of the mean and std of the tumor cell density. Since in the present study, the mean and MAP estimates are very similar, and the std values are small, the MAP estimates are used. An overview of the proposed personalized RT design is shown in Fig. 1 (IV), while the details are described in the following subsections. The tumors recurrence patterns are used to assess the benefits of the proposed RT plan over the standard treatment protocol. For evaluation purposes, all recurrence scans are registered to the preoperative anatomy. To prevent registration errors arising from mapping the anatomy with the resection cavity to the preoperative brain, rigid registration is used. This provides a sufficient mapping for most cases, however it cannot capture the post-treatment tissue displacement around the ventricles in patients P7-P8, making the mapping less accurate in these regions. The design of methods that provide robust registration between pre- and post-operative brain anatomies is still an open problem.

Iii-C1 Dose distribution

An ideal CTV covers all the residual tumor, including infiltrating tumor cells that are invisible on the pretreatment imaging scans, while sparing healthy tissue. We use the tumor recurrence pattern to evaluate the efficiency () of the CTV, defined here as the relative volume of the recurrence tumor () contained within the CTV:

(11)

Figure 5 (H) shows the FLAIR scans with the first detected tumor recurrence outlined by the pink lines. The margin of the administered , designed by the standard RTOG protocol with a margin around the visible tumor, is marked by the green lines in (E,F,H). The personalized CTV, referred as , is constructed by thresholding the MAP tumor cell density at for all patients. (This value was chosen so that the efficiency of the is comparable to that of the ). The outlines of the proposed are shown as the blue isocontours in Fig. 5 (H-I). A visual comparison of the CTVs shown in Fig. 5 (H) and Fig. 1 (IV), and a quantitative comparison presented in Fig. 6, show that the proposed personalized plans spare more healthy tissue, hence reducing radiation toxicity, while maintaining the efficiency of the standard RTOG protocol. Both plans show reduced efficiency for patients P7-P8, mainly around the ventricles, which may be caused by misalignment between the preoperative and recurrence anatomies.
These preliminary results imply that the regions predicted as a tumor-free by the model, remain tumor-free and thus the model predictions have potential to guide personalized CTV design. The standard or hospital-specific protocols can be updated by the model predictions to spare brain tissue not infiltrated by the tumor. This can lead to significant savings in the healthy tissue, especially in the cases of large lesions or lesions close to hemispheres separation and other anatomical constraints.

Fig. 6: A comparison of the RT plan based on the RTOG protocol (green) and MAP estimates (blue). (A) The overall irradiated volume and (B) the corresponding efficiency . The uses a smaller irradiation volume while having a comparable efficiency as the .

Iii-C2 Dose-escalation

No dose-escalation plan for GBM patients has been yet approved by phase-III-clinical trials. Here, we present a theoretical comparison of two escalation plans targeting high tumor cellularity regions identified by: 1) FET-PET enhancement as proposed in [5] and 2) MAP estimates. We evaluate the efficiency of an escalation plan by its capability of targeting T1Gd-enhanced tumor recurrence regions. In these regions, the recurrent tumor has high cellularity, despite having received the full radiation dose, suggesting tumor radioresistance. Figure 5 (I) shows the T1Gd scans with the first detected tumor recurrence. The margins of the T1Gd-enhanced tumor recurrence are marked by the yellow lines, while the outlines of the dose-escalation plans designed by the FET-PET enhancements are shown in magenta. The FET enhancements do not fully cover the T1Gd recurrent tumor in patients P4-P7, providing a possible explanation for why improvements in progression-free survival have not been observed in [5]. In comparison, the MAP estimates, calibrated by the FET-PET signal, extend the information about the tumor cell density in the periphery of the visible lesion. Figure 5 shows two possible dose-escalation plans based on the inferred MAP tumor cell density: (I) a single-level dose-escalation based on the thresholded MAP solution with threshold marked by the orange lines and (J) a cascaded four-level escalation plan constructed by thresholding the MAP tumor cell density at . The optimal design of a personalized dose-escalation plan would require more extensive studies. However, these preliminary results show that the inferred high cellularity regions coincide with the areas of tumor recurrence better than those suggested by the FET-PET enhancement alone. The Bayesian inference framework developed here thus provide a promising tool for a rational dose-escalation design.

Iv Conclusion

We have demonstrated that patient-specific, data-driven modeling can extend the capabilities of personalized RT design for infiltrative brain lesions. We combined patient structural and metabolic scans from a single time point with a computational tumor growth model through a Bayesian inference framework and predicted the tumor distribution beyond the outlines visible in medical scans. The patient-specific tumor estimates can be used to design personalized RT plans, targeting shortcomings of standard RT protocols. The software and data used in this work are publicly released222https://github.com/JanaLipkova/GliomaSolver to facilitate translation to clinical practice and to encourage future improvements. In the future, the Bayesian framework developed here could also be extended to predict individual patient responses to RT by incorporating data obtained during the course of treatment as done in [33], in which non-spatial tumor models were used. In this way, the treatment can be further improved by adaptively refining the RT plans based on the predicted patient responses. Moreover, the basic FK tumor growth model could be replaced by a Fokker-Planck diffusion model [16], which would not increase the number of unknown parameters or affect the computational complexity significantly, but might provide a better description of biological diffusion. Future work could also incorporate more advanced models, such as [34], that account for cancer stem cells, their progeny and nonlinear coupling between the tumor and the neovascular network. However, it remains to be seen whether scans acquired at a single time point would provide enough information to calibrate the advanced models sufficiently. If not, then simpler, well-calibrated models may prove to be more informative. Finally, in future studies, the computational framework developed here will be tested on a larger patient cohort and prospective clinical trials will be performed. In summary, the results presented here provide a proof-of-concept that multimodal Bayesian model calibration holds a great promise to assist the development of personalized RT protocols.

Acknowledgment

JL acknowledges partial funding from the National Science Foundation-Division of Mathematical Sciences (NSF-DMS) through grant DMS-1714973 and the Center for Multiscale Cell Fate Research at UC Irvine, which is supported by NSF-DMS (DMS-1763272) and the Simons Foundation (594598, QN). JL additionally acknowledges partial funding from the National Institutes of Health (NIH) through grant 1U54CA217378-01A1 for a National Center in Cancer Systems Biology at the University of California, Irvine, and NIH grant P30CA062203 for the Chao Comprehensive Cancer Center at the University of California, Irvine.
JL acknowledges the partial funding from the Bavaria California Technology Center (BaCaTec) through grant 6090142.

References

  • [1] R. Stupp, M. Brada, M. Van Den Bent, J.-C. Tonn, and G. Pentheroudakis, “High-grade glioma: ESMO Clinical Practice Guidelines for diagnosis, treatment and follow-up,” Annals of oncology, 2014.
  • [2] R. Stupp, W. P. Mason, M. J. Van Den Bent, M. Weller, B. Fisher, M. J. Taphoorn, K. Belanger, A. A. Brandes, C. Marosi, U. Bogdahn et al., “Radiotherapy plus concomitant and adjuvant temozolomide for glioblastoma,” New England Journal of Medicine, 2005.
  • [3] N. G. Burnet, S. J. Thomas, K. E. Burton, and S. J. Jefferies, “Defining the tumour and target volumes for radiotherapy,” Cancer Imaging, 2004.
  • [4] A. K. Paulsson, K. P. McMullen, A. M. Peiffer, W. H. Hinson, W. T. Kearns, A. J. Johnson, G. J. Lesser et al., “Limited margins using modern radiotherapy techniques does not increase marginal failure rate of glioblastoma,” American journal of clinical oncology, 2014.
  • [5] M. Piroth, M. Pinkawa, R. Holy, J. Klotz, S. Schaar, G. Stoffels et al., “Integrated boost IMRT with FET-PET-adapted local dose escalation in glioblastomas,” Strahlentherapie und Onkologie, 2012.
  • [6] L. Souhami, W. Seiferheld, D. Brachman, E. B. Podgorsak, M. Werner-Wasik, R. Lustig et al., “Randomized comparison of stereotactic radiosurgery followed by conventional radiotherapy with carmustine to conventional radiotherapy with carmustine for patients with glioblastoma multiforme: report of radiation therapy oncology group 93-05 protocol,” International Journal of Radiation Oncology* Biology* Physics, 2004.
  • [7] E. C. Halperin, G. Bentel, E. R. Heinz, and P. C. Burger, “Radiation therapy treatment planning in supratentorial glioblastoma multiforme: an analysis based on post mortem topographic anatomy with CT correlations,” International Journal of Radiation Oncology* Biology* Physics, 1989.
  • [8] S. Rockwell, I. T. Dobrucki, E. Y. Kim, S. T. Marrison, and V. T. Vu, “Hypoxia and radiation therapy: past history, ongoing research, and future promise,” Current molecular medicine, vol. 9, no. 4, 2009.
  • [9] S. E. Combs, I. Burkholder, L. Edler, S. Rieken, D. Habermehl, O. Jäkel et al., “Randomised phase I/II study to evaluate carbonion radiotherapy versus fractionated stereotactic radiotherapy in patients with recurrent or progressive gliomas: The CINDERELLA trial,” BMC cancer, 2010.
  • [10] C. Tsien, J. Moughan, J. M. Michalski, M. Gilbert et al., “Phase I three-dimensional conformal radiation dose escalation study in newly diagnosed glioblastoma: Radiation Therapy Oncology Group Trial 98-03,” International Journal of Radiation Oncology* Biology* Physics, 2009.
  • [11] M. Hingorani, W. P. Colley, S. Dixit, and A. M. Beavis, “Hypofractionated radiotherapy for glioblastoma: strategy for poor-risk patients or hope for the future?” The British Journal of Radiology, Sep 2012.
  • [12] S. Rieken, D. Habermehl, F. L. Giesel, C. Hoffmann et al., “Analysis of FET-PET imaging for target volume definition in patients with gliomas treated with conformal radiotherapy,” Radiotherapy and Oncology, 2013.
  • [13] F. Stockhammer, M. Plotkin, H. Amthauer, F. K. H. van Landeghem et al., “Correlation of F-18-fluoro-ethyl-tyrosin uptake with vascular and cell density in non-contrast-enhancing gliomas,” J Neurooncol, 2008.
  • [14] M. Hutterer, M. Nowosielski, D. Putzer, N. L. Jansen, M. Seiz, M. Schocke, M. McCoy, G. Göbel, C. la Fougère, I. J. Virgolini et al., “[18F]-fluoro-ethyl-L-tyrosine PET: a valuable diagnostic tool in neuro-oncology, but not all that glitters is glioma,” Neuro-oncology, vol. 15, no. 3, pp. 341–351, 2013.
  • [15] J. Alfonso, K. Talkenberger, M. Siefert, B. Klink et al., “The biology and mathematical modelling of glioma invasion: A review,” J. R. Soc. Interface, 2017.
  • [16] A. Swan, T. Hillen, J. C. Bowman, and A. D. Murtha, “A patient-specific anisotropic diffusion model for brain tumour spread,” Bulletin of mathematical biology, 2018.
  • [17] A. Hodgkinson, M. A. Chaplain, P. Domschke, and D. Trucu, “Computational approaches and analysis for a spatio-structural-temporal invasive carcinoma model,” Bulletin of mathematical biology, 2018.
  • [18] V. Cristini and J. Lowengrub, Multiscale modeling of cancer.   Cambridge University Press, 2010.
  • [19] A. Mohamed and C. Davatzikos, “Finite element modeling of brain tumor mass-effect from 3D medical images,” in International Conference on Medical Image Computing and Computer-Assisted Intervention.   Springer, 2005, pp. 400–408.
  • [20] Hogea C, et al., “A robust framework for soft tissue simulations with application to modeling brain tumor mass effect in 3D MR images,” PHYSICS IN MEDICINE AND BIOLOGY, 2007.
  • [21] J. Unkelbach, B. Menze, E. Konukoglu, F. Dittmann et al., “Radiotherapy planning for glioblastoma based on a tumor growth model: improving target volume delineation,” Physics in medicine and biology, 2014.
  • [22] E. Konukoglu, O. Clatz, H. Delingette, and N. Ayache, “Personalization of reaction-diffusion tumor growth models in MR images: application to brain gliomas characterization and radiotherapy planning,” 2010.
  • [23] K. R. S. Hana L.P. Harpold, Ellsworth C. Alvord, “The evolution of mathematical modeling of glioma proliferation and invasion,” J Neuropathol Exp Neurol, 2007.
  • [24] P. R. Jackson, J. Juliano, A. Hawkins-Daarud, R. C. Rockne, K. R. Swanson et al., “Patient-specific mathematical neuro-oncology: using a simple proliferation and invasion tumor model to inform clinical practice,” Bulletin of mathematical biology, 2015.
  • [25] R. C. Rockne, A. D. Trister, J. Jacobs, A. J. Hawkins-Daarud et al., “A patient-specific computational model of hypoxia-modulated radiation resistance in glioblastoma using 18 F-FMISO-PET,” Journal of the Royal Society Interface, 2015.
  • [26] M. Lê, H. Delingette, J. Kalpathy-Cramer, E. R. Gerstner et al., “Bayesian personalization of brain tumor growth model,” in International Conference on Medical Image Computing and Computer-Assisted Intervention.   Springer, 2015.
  • [27] A. Hawkins-Daarud, S. K. Johnston, and K. R. Swanson, “Quantifying Uncertainty and Robustness in a Biomathematical Model Based Patient-Specific Response Metric for Glioblastoma,” bioRxiv, 2018.
  • [28] B. Menze, K. Van Leemput, A. Honkela, E. Konukoglu, M.-A. Weber et al., “A generative approach for image-based modeling of tumor growth,” in Information Processing in Medical Imaging.   Springer, 2011.
  • [29] D. Rossinelli, B. Hejazialhosseini, W. van Rees, M. Gazzola et al., “MRAG-I2D: Multi-resolution adapted grids for remeshed vortex methods on multicore architectures,” Journal of Computational Physics, 2015.
  • [30] J. Ching and Y.-C. Chen, “Transitional markov chain monte carlo method for bayesian model updating, model class selection, and model averaging,” Journal of engineering mechanics, 2007.
  • [31] P. E. Hadjidoukas, P. Angelikopoulos, C. Papadimitriou, and P. Koumoutsakos, “4u: A high performance computing framework for bayesian uncertainty quantification of complex models,” Journal of Computational Physics, 2015.
  • [32] C. A. Cocosco, V. Kollokian, R. K.-S. Kwan, G. B. Pike et al., “Brainweb: Online interface to a 3D MRI simulated brain database,” in NeuroImage.   Citeseer, 1997.
  • [33] I. Tariq, T. Chen, N. F. Kirkby, and R. Jena, “Modelling and bayesian adaptive prediction of individual patients’ tumour volume change during radiotherapy,” Physics in Medicine & Biology, vol. 61, no. 5, p. 2145, 2016.
  • [34] H. Yan, M. Romero-López, L. I. Benitez, K. Di et al., “3D mathematical modeling of glioblastoma suggests that transdifferentiated vascular endothelial cells mediate resistance to current standard-of-care therapy,” Cancer research, 2017.