Skull stripping and brain tissue segmentation
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
Glioblastoma is a highly invasive brain tumor, whose cells infiltrate
Existing radiotherapy (RT) plans for brain tumors derive from population...
Modeling of brain tumor dynamics has the potential to advance therapeuti...
Tumor growth prediction, a highly challenging task, has long been viewed...
The emerging field of precision oncology relies on the accurate pinpoint...
Currently self-report pain ratings are the gold standard in clinical pai...
With the advance of imaging technology, digital pathology imaging of tum...
Skull stripping and brain tissue segmentation
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 . 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 . 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 . Biopsies  and post-mortem studies  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 , 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 . The Radiation Therapy Oncology Group (RTOG) phase-I-trial  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) , 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  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  are used, it is not clear how robust the model predictions are and what benefits they offer over the standard treatment protocols.
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 inFig. 1.
Many tumor growth models are based on the Fisher-Kolmogrov (FK) equation , 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:
denotes proliferation rate. The tumor infiltration into the surrounding tissues is modeled by the tensorwhere is a identity matrix and
The terms and denote percentage of white and grey matter at voxel , while and stand for tumor infiltration in the corresponding matter. We assume . 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  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).
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,
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
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).
is modeled by a Bernoulli distribution:
Here is the probability of observing the tumor in the MRI scan and it is assumed to be a double logistic sigmoid:
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
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
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).
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:
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 which iteratively constructs series of intermediate PDFs:
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  and allow parallel execution. We use a highly parallel implementation of the TMCMC algorithm provided by the 4U framework .
The inferred parametric uncertainties are propagated through the model to obtain robust predictions about u given by:
or by simplified measures such as the mean
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 .
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.
The model is used to generate a 3D synthetic tumor in a brain anatomy obtained from  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 inFig. 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.
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 . 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).
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.
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:
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.
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  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 . 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.
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 , 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 , 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 , 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.
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.