Toward image-based personalization of glioblastoma therapy: A clinical and biological validation study of a novel, deep learning-driven tumor growth model

Abstract

Background

The diffuse growth pattern of glioblastoma is one of the main challenges for accurate treatment. Computational tumor growth modeling has emerged as a promising tool to guide personalized therapy. Here, we performed clinical and biological validation of a novel growth model, aiming to close the gap between the experimental state and clinical implementation.

Methods

One hundred and twenty-four patients from The Cancer Genome Archive (TCGA) and 397 patients from the UCSF Glioma Dataset were assessed for significant correlations between clinical data, genetic pathway activation maps (generated with PARADIGM; TCGA only), and infiltration (D_w) as well as proliferation (ρ) parameters stemming from a Fisher–Kolmogorov growth model. To further evaluate clinical potential, we performed the same growth modeling on preoperative magnetic resonance imaging data from 30 patients of our institution and compared model-derived tumor volume and recurrence coverage with standard radiotherapy plans.

Results

The parameter ratio D_w/ρ (P < .05 in TCGA) as well as the simulated tumor volume (P < .05 in TCGA/UCSF) were significantly inversely correlated with overall survival. Interestingly, we found a significant correlation between 11 proliferation pathways and the estimated proliferation parameter. Depending on the cutoff value for tumor cell density, we observed a significant improvement in recurrence coverage without significantly increased radiation volume utilizing model-derived target volumes instead of standard radiation plans.

Conclusions

Identifying a significant correlation between computed growth parameters and clinical and biological data, we highlight the potential of tumor growth modeling for individualized therapy of glioblastoma. This might improve the accuracy of radiation planning in the near future.

glioblastoma, tumor growth modeling, personalized therapy, deep learning

Key Points

• Tumor growth modeling promises individualized therapy for glioblastoma patients.

• Here, the model-derived parameters significantly correlate with clinical and genetic variables.

• With a novel neural solver, growth modeling becomes applicable in short computing time.

Importance of the Study

There is a crucial need to improve treatment in glioblastoma patients. Tumor growth modeling promises to individualize radiotherapy, providing a personalized dose distribution and thus potentially improving patient outcomes. We here present an in-depth validation of our deep learning-based model for inferring tumor-specific growth model parameters from preoperative MRI scans. In particular, our study revealed a significant correlation between the predicted patient-specific tumor parameters and the overall patient survival as well as activation of proliferation pathways stemming from mRNA data. Moreover, we showed that the radiotherapy plans derived from the patient-specific tumor predictions provide better coverage of tumor recurrence than standard plans. Our results highlight the potential of this approach to improve patient treatment in the future and constitute a starting point for further clinical implementation of tumor growth modeling.

Glioblastoma (GBM), the most aggressive primary brain tumor, is unique in several ways. Most characteristic, its spread is highly infiltrative and thus, standard magnetic resonance imaging (MRI) fails in distinguishing surrounding brain tissue and peritumoral edema from non-enhancing tumor infiltration. Moreover, GBM is the only cancer in which morbidity is almost solely due to local progression.¹ Therefore, a reliable prediction of individual tumor spread, even beyond the tumor outlines visible on MRI scans, would mark a tremendous advancement in the diagnosis and treatment of this fatal disease.

Standard clinical treatment of GBM consists of maximal safe tumor resection, followed by combined radio-chemotherapy targeting the tissue around the resection cavity to account for residual tumor cells.² Due to a lack of knowledge about the individual tumor spread, the clinical target volume (CTV) for irradiation is defined by a standard margin of usually 2 cm around the resection cavity (and the FLAIR-hyperintense region, depending on the specific target volume recommendations) sparing sensitive areas such as the brainstem and the chiasm.³ Besides local recurrences occurring directly at the resection cavity, the high recurrence rate of GBM can be possibly attributed to the microscopic invisible tumor infiltration beyond the margins of the CTV.⁴ Radiotherapy plans might therefore benefit from an accurate prediction of tumor infiltration pathways which is the reason for the increasing amount of research focused on this field.^5–7

Mathematical tumor growth modeling aims to simulate tumor growth by using established mathematical concepts such as the reaction–diffusion equation that has first been applied by Murray et al. in the early 1990s.⁸ In the context of medical imaging-based modeling, current models for tumor growth describe gross biomechanical phenomena such as migration and proliferation of tumor cells, tumor interaction with brain tissue (mass effect),⁹ or formation of the necrotic tumor core.¹⁰

To identify parameters of the growth model that most exactly match the observed tumor anatomy on medical images, one has to solve the so-called inverse problem. This can be achieved using formulations such as constrained optimization^11,12 or Bayesian inference.^13,14 However, most of these approaches still require large amounts of computing time, reaching for instance from approximately 5 h for a scheme providing point parametric estimates,¹² up to several days for models that also provide uncertainty estimates for the inferred model parameters.¹³ Further, some models require specific input sequences (such as amino acid PET) to initialize the solver, limiting their broad clinical applicability.¹³

In this work, we implement a novel model that has recently been introduced by our working group.¹⁵ This model applies a neural network-based methodology that shows potential for clinical implementation by estimating the individual tumor growth parameters only from preoperative MR standard imaging in a matter of minutes. Here, we present validation and clinical application of this tumor growth model in 3 independent ways: First, in 2 large, independent glioblastoma cohorts, we investigate correlations between patient overall survival and the calculated tumor volume as well as the model-derived parameters. In the second part, we compare the derived model parameters, that is, proliferation and infiltration indices, with corresponding genetic data accessed from The Cancer Genome Archive Network (TCGA) looking for biological correlates of the estimated growth behavior.¹⁶

Finally, to assess the potential of the model to guide personalized radiotherapy plans, we compare model-based target volumes with standard radiotherapy plans from our institution with regard to the complete radiation volume and coverage of the following tumor recurrence. Thus, this study aims to facilitate translation of computational tumor modeling toward clinical applications.

Materials and Methods

Patient Characteristics and Study Design

We investigated 3 different datasets in this study with the same growth model. First, imaging, clinical, and genetic data from 124 adult glioblastoma patients from TCGA¹⁷ were assessed, as well as clinical and imaging data from 501 glioma patients from the UCSF Preoperative Diffuse Glioma MRI Dataset (UCSF).¹⁸N = 6 patients in the TCGA dataset had an IDH1 mutation (R 132 H) and therefore were excluded from our study. In the UCSF dataset, we also excluded all IDH mutant tumors (n = 104), leaving a total of 397 IDH wildtype glioblastomas.

In the second part, we evaluated 30 glioblastoma cases from our prospective glioma cohort (BraTUM) that have been approved by the local ethics committee. All of these underwent gross-total tumor resection and were confirmed by histopathological study as IDH wildtype glioblastoma as per WHO 2021 classification of CNS tumors.¹⁹ In our institutional cohort, all patients were scanned in a 3T whole-body MRI scanner (Achieva, Philips Medical Systems, Best, The Netherlands). The protocol consisted of T2 turbo spin echo (T2), T2-FLAIR, non-enhanced and contrast-enhanced T1 (CE-T1).

Image Processing

For our local patient cohort, we exported preoperative and postoperative MR images as well as the MR exam showing the first tumor recurrence from our PACS. Images were processed and segmented using the openly available BraTS Toolkit.^20–22 All postoperative scans from our patients were registered onto the preoperative scan using a 2-stage approach (rigid followed by a deformable SyN registration) using ANTs.²³ For TCGA and UCSF data, we downloaded preprocessed preoperative imaging data.¹⁷ All registrations and segmentations were checked by neuro-radiologist (M.M. with 5 years of experience in glioma imaging).

Preparation of Genetic Data

Genomic and clinical (meta)data for TCGA data were downloaded from the TCGA server. These encompass DNA methylation-based classification to exclude G-CIMP cases (carrying an IDH mutation) as well as PARADIGM pathway analysis data.^24,25 In brief, Pathway Representation and Analysis by Direct Inference on Graphical Models (PARADIGM) predicts the activity of a diverse set of molecular processes from a probabilistic belief propagation strategy that incorporates multimodal data such as copy number and gene expression estimates. The 433 pathways were filtered for keywords such as “proliferation” that seemed to be associated with cell growth and tumorigenesis. This resulted in 129 relevant pathways that were further assessed.

Computational Tumor Growth Modeling

To estimate tumor cell growth and infiltration, we model the tumor using the established reaction–diffusion model:

$${\displaystyle \begin{array}{l}{\displaystyle \frac{\mathrm{\partial }c}{\mathrm{\partial }t}=\mathrm{\nabla }\cdot (D\mathrm{\nabla }c)+\rho c(1-c),}\end{array}}$$

where ∂c/∂t is the evolution of tumor cell density c over time, D is the diffusion tensor, and ρ denotes the tumor proliferation rate in (1/s). In particular, the tensor D is defined through D_w and D_g, where D_w is the diffusion coefficient in white matter, and D_g is the diffusion coefficient in gray matter, both with the unit (m²/s). Notably, those diffusion tensors are not associated with or derived from diffusion-weighted imaging, but rather mathematically describe the tendency (and direction) of tumor cells to infiltrate the surrounding brain. Both tensors in diffusion-tensor imaging and in our models are mathematical tensors that describe varying spatially propagation of cells. In the case of diffusion tensor imaging, it is a migration of water molecules, while in the case of our tumor model, it is a propagation of tumor cells. It is typically believed that tumor cells migrate faster in white matter compared to gray matter. Also, no cell diffusion takes place in the cerebrospinal fluid. Thus, in our model, we define the diffusion coefficient in white matter to be 10 times larger than in gray matter. Furthermore, tumor cell migration into the cerebrospinal fluid is constrained in the model.

Traditionally, these partial differential equations are solved numerically, as in Lipkova et al.²⁶ Given that these solvers have high demands in computational power and time, deep learning-based models are an attractive alternative.¹⁴ In recent work, we have trained a deep neural network to solve the inverse problem, that is, estimate the model parameters best describing an individual patient’s tumor.¹⁵ In brief, a patient’s tumor segmentation is morphed into the SRI24 atlas space²⁷ using ANTs. Next, a deep neural network that has been trained on 80,000 tumor simulations (that have been generated by randomly sampling parameter combinations within plausible ranges²⁸) in this atlas space (and has been validated on 8000, and tested on 12,000 simulations, respectively) is used to find the parameters of the above growth model which best fit the individual tumor. In other words, the exact degree of the cell diffusion (D) and proliferation (ρ) are the parameters that are inferred by connecting the tumor model with the tumor information available on medical scans. Using these parameters as input, we run a single forward simulation to obtain tumor cell distribution and warp this result back to patient space using the inverse of the registration transform.

The training was done on Nvidia P8000 GPU with 48 GB of VRAM but could also be done on CPU without any relevant prolongation of the total pipeline running time (as the network prediction both on CPU and on GPU is a lot faster than the registration step in the pipeline). Detailed technical information about the applied growth model, its architecture, and training process can be found in our recent technical publication.¹⁵

Comparison with Standard Radiotherapy Planning

For the local patient cohort, we exported CTV maps from the radiation oncology planning system. Planning was performed according to the ESTRO-ACROP recommendations.³ These maps (planned on the postoperative MR image) were mapped to the preoperative space as described above. To binarize the continuous tumor cell distributions from our computational growth model, we re-scaled densities by the 99th percentile (to calculate relative tumor cell densities) and tested various relative tumor cell density cutoffs (c_t = {0.33, 0.5, 0.66, 0.75}) to compare both the coverage of the later recurrence as well as the total radiation volume.

Statistical Analysis

To assess the prognostic significance of growth model parameters (in particular the ratio D_w/ρ, which is associated with infiltration length,²⁹ and the simulated tumor volume), we calculated multivariable Cox regression models, including established prognostic factors, that is, patient age at diagnosis and volume of the contrast-enhancing tumor. For reference, single-variable Cox regression was performed as well. Further, to avoid misleading results. we checked for collinearity of the predictors in the multivariable Cox regression models. For correlation of PARADIGM pathway activations and model parameters, we employed Spearman correlation testing. To account for multiple tests, we performed 1000 label permutations and used the resulting distribution to perform significance testing. Lastly, to check if there are significant differences between volume and recurrence coverage of model-derived CTV and standard CTV from the clinical routine we performed Wilcoxon signed-rank test. All analyses were done in Python 3.8.

Results

Clinical Characteristics of the Study Samples

In the TCGA-dataset (n = 124), mean patient age at initial diagnosis was 57.6 years (SD = 14.3). There were 76 male and 48 female patients. In the UCSF cohort (n = 397), 235 were male and 162 were female with a mean age of 61.5 years (SD = 11.9). In our institutional cohort (BraTUM, n = 30), mean age at diagnosis was 62.2 (SD = 9.1). There was also a male predominance in the dataset (18 vs. 12).

Correlation Between Growth Model Parameters and Clinical Metrics

For the TCGA and UCSF cohorts, we performed multivariable Cox regression analysis with the estimated parameter ratio D_w/ρ as well as the estimated tumor volume (using a relative tumor cell density cutoff of 0.5). Since we found a fairly strong collinearity between those predictors (Pearson’s r = 0.5), we did not integrate both factors in the same Cox regression model. For reference, we included the volume of the contrast-enhancing tumor as well as patient age at diagnosis in both models which has been shown to correlate with survival in prior studies.^30,31 Of note, between real contrast-enhancing tumor volume and simulated tumor volume as well as D_w/ρ, the correlation was weaker (Pearson’s r = 0.31 and −0.2, respectively). All results of the correlation analysis can be found in the Supplementary material (Additional file 1).

In the TCGA, D_w/ρ and the estimated tumor volume significantly correlated inversely with survival, whereas the contrast-enhancing tumor volume did not (Table 1A). More precisely, the parameter ratio D_w/ρ had a hazard ratio (HR) of 1.9 (P < .05). Patient age at the time of diagnosis was significantly correlated with survival in both models. Figure 1 depicts the survival curves of the 5th, 25th, 50th, 75th, and 95th percentile of D_w/ρ in TCGA that show shorter survival with an increasing ratio of D_w/ρ, that is, with a higher infiltration of tumors. The baseline is the median value.

Table 1.

Multivariable Comparison Analysis of Clinical Parameters and Growth Model-Derived Parameters.

Variables	Multivariable Analysis I: (Cox Proportional Hazard Model)		Multivariable Analysis II: (Cox Proportional Hazard Model)
	Hazard Ratio	P-value	Hazard Ratio	P-value
(A) TCGA (n = 124 patients)
Volume of contrast-enhancing tumor	1	0.07	1	0.85
Diagnosis age	1.03	<0.05	1.02	0.01
Dw/ρ	1.9	<0.05	–	–
Simulated tumor volume	–	–	1.16	0.04
(B) UCSF (n = 397 patients)
Volume of contrast-enhancing tumo	1.01	0.05	1	0.14
Diagnosis age	1.03	<0.05	1.03	0.01
Dw/ρ	1.12	0.2	–	–
Simulated tumor volume	–	–	1.08	0.04

Figure 1.

Survival probability curves show shorter survival times of glioblastoma patients with increasing percentile of the computed parameter ratio D_w/ρ, which is a surrogate for tumor infiltration, derived from the tumor growth model.

In the UCSF cohort, we again observed a significant association of simulated tumor volume with overall survival when also including age and contrast-enhancing tumor volume (HR 1.08, P = .049). Also here, contrast-enhancing tumor volume was not significantly associated with survival after inclusion of simulated tumor volume. For the parameter ratio D_w/ρ we found no significant association with survival (Table 1B). For reference, we also performed a single-variable Cox regression analysis and added the results to the Supplementary Material (Additional file 2).

Biological Validation Utilizing Pathway Activation Maps

Eleven pathways showed a significant correlation with D_w/ρ after permutation-based adjustment for multiple comparisons. Intriguingly, we observed a clear association with the proliferation parameter ρ, as most correlation parameters were negative with D_w/ρ. For example, patients whose genetic alterations led to an over-activity of a specific cellular pathway for cell proliferation (eg, pathway “119_Cell Proliferation”) showed a significantly lower parameter ratio of D_w/ρ in our computational growth model. A lower parameter ratio might be due to higher values of the proliferation parameter ρ, and/or due to lower values of the infiltration parameter D_w. Therefore, our model’s estimation matches the genetic observations of the tumor’s proliferation activity. Table 2 gives an overview of those signaling pathways and the results of the correlation analysis. The complete table of 129 relevant pathways can be found in the Supplementary Material (Additional file 3a) as well as the complete, unfiltered table of 433 pathways (Additional file 3b).

Table 2

Relevant Proliferation Pathways That Revealed Significant Correlation With Growth Model-Derived Parameters.

Pathway	Dw/ρ
	Spearman Correlation Coefficient	P-value (Spearman)
119_Cell proliferation	−0.251	0.009
119_Apoptosis	0.252	0.007
71_Regulation of calcium-dependent cell–cell adhesion	−0.234	0.012
45_Regulation of cell–cell adhesion	−0.233	0.007
71_Regulation of cell–cell adhesion	−0.234	0.008
22_Re-entry into mitotic cell cycle	−0.234	0.011
38_Homophilic cell adhesion	−0.222	0.014
109_Apoptosis	−0.198	0.034
43_Positive regulation of cell migration	−0.159	0.031
118_Cell proliferation	−0.178	0.058
51_Mitosis	−0.181	0.047
108_Cell migration	0.172	0.05

Comparison With Standard Radiotherapy Plans

In the third part of the study, we evaluated the growth model for its clinical potential on data from 30 patients from our institutional glioma cohort. Therefore, we compared the simulated tumor infiltration with the corresponding CTV as planned for the clinical routine. We investigated: (1) how much the radiation volumes differ and (2) if coverage of future contrast-enhancing recurrent tumors might be improved when using the predicted tumor spread instead of standard radiotherapy plans. The results for the different cutoff values (c_t) of tumor cell density are shown in Table 3 (with the raw data accessible in Supplementary Material Additional file 4). Naturally, choosing lower values for c_t, leads to larger simulated tumor volumes and hence, to an improvement of recurrence coverage. On the other hand, this of course also yields larger radiation volumes. Figure 2 illustrates improvement of recurrence coverage for each patient for 4 different cutoff values of the tumor cell density (c_t = {0.33, 0.5, 0.66, 0.75}). To approach clinical utility, we were looking for a reasonable tradeoff between improved recurrence coverage and a radiation volume as low as possible. A c_t-value of 0.5 leads to a significant mean improvement of recurrence coverage of 4.567% (P = .044) without a significant increase in tumor volume compared to the CTV (P = .926), thus showing the most promising results for clinical application. The absolute and relative numbers of the plots can be found in the Supplementary Material (Additional file 5). For example, with a model-derived radiation plan defined by a c_-value of 0.5 we would have improved recurrence coverage in 37% of the patients, while there would have been no difference in 43%, and even less recurrence coverage in 20% of the cases, compared to standard radiotherapy planning.

Table 3.

Results of the Clinical Validation Study. Comparison of Clinical Target Volume and Recurrence Coverage Between Standard Radiotherapy (RT) Plan and Tumor Growth Model (TGM)—Derived Target Delineations With Different Isolines Accounting for CutOff Values for Tumor Cell Density.

Cutoff Value for Tumor Cell Density	Clinical Target Volume			Recurrence Coverage
	Mean Volume Ratio (TGM vs. RT Plan)	T-Value (Wilcoxon Signed Rank Test)	P-Value (Wilcoxon Signed Rank Test)	Mean Improvement of Recurrence Coverage (TGM vs. RT Plan)	T-Value (Wilcoxon Signed Rank Test)	P-Value (Wilcoxon Signed Rank Test)
0.33	1.493	90.0	0.003	5.160%	5.0	0.005
0.5	1.220	228.0	0.093	4.567%	34.0	0.044
0.66	1.011	151.0	0.094	2.461%	77.0	0.469
0.75	0.907	106.0	0.009	−0.626%	92.0	0.627

Figure 2.

Waterfall plots illustrate the change in recurrence coverage when utilizing isolines of cutoff values for tumor cell density for target delineation instead of standard radiotherapy planning. Subplots A‐D represent cutoff values for tumor cell density of 0.33, 0.5, 0.66, and 0.75, respectively. In general, the lower the cutoff value is chosen, the larger the estimated target volume becomes, and naturally, the more recurrence area is covered.

Figure 3 presents an exemplary case of a patient with recurrent GBM following the fiber tracts of the corpus callosum. Of note, on the preoperative images, no visible alterations were visible in most of the area of later recurrence. In this case, standard CTV planning on postoperative MRI captured only a small proportion of the later recurrent site. In contrast, the model-derived target volumes based on three different isolines for estimated tumor cell density showed significant improvement in the recurrence coverage, as this tumor was simulated to have a highly infiltrative nature (high D_w), proven by the later recurrence.

Figure 3.

Comparison of standard clinical target volume (CTV) and computed target delineations derived from isolines of different estimated tumor cell densities by the tumor growth model (TGM). Underlying images are contrast-enhanced T1 (CE-T1).

Exemplary case of a patient with recurrent glioblastoma. (A) Depicts initial size of the tumor in CE-T1 imaging. (B) Shows the clinical target volume (CTV) planned on postoperative imaging and registered to the recurrence image that is visible underneath the overlay (C). Subparts D‐F illustrate target volumes derived from isolines of estimated tumor cell density by the growth model (TGM) with cutoff values of c_t = 0.33 (D), c_t = 0.5 (E), and c_t = 0.66 (F). Note that percentage of recurrence coverage improves with lower cutoff value, and is, in this case, worst with traditional radiotherapy planning. However, there was also an increasing tumor volume with a relative volume ratio of 1.03, 1.25, and 1.52. Therefore, a reasonable trade-off needs to be found between improved recurrence coverage and a radiation volume as small as possible.

Discussion

Important Aspects of Growth Modeling

To aid clinicians in individualizing the therapy of GBM patients, tumor growth models need to be (1) reliable and sufficiently validated, (2) applicable with widely available computational resources, and (3) time-efficient and feasible in daily radiological practice. Here, we introduce a growth model that can be trained in a supervised fashion with a data set of numerical simulations. Hence, unlike in other successful studies,³² there is no need for a large data set of longitudinal brain tumor images which can be difficult to obtain for smaller institutions and may fail to capture the diffuse, multi-faceted growth patterns seen in this disease. Furthermore, the model is able to predict patient-specific spatial distribution of the tumor on data from a single time point (preoperatively) and commonly acquired MR sequences, namely CE-T1 and T2/FLAIR, which also makes it easily accessible. A major advantage of our model is also its short computing time of 4–7 min. in total, which in this case includes registration, morphing the input images to atlas space, inference, tumor simulation, and morphing the growth simulation back to patient space. Importantly, the computing time is stable across tumor models of different complexity. Thus, this model outclasses prior approaches in terms of clinical applicability.^12,26

Different Growth Model Parameters and Survival

The 2 most important parameter outputs of our growth model are the proliferation coefficient ρ and the diffusion coefficient D_w, which can be interpreted as a marker of infiltration. Further, from these individually estimated parameters we are able to infer tumor cell distributions, even beyond the visibly contrast-enhancing tumor. To investigate the biological transferability of those, we performed survival analysis in a multivariable approach, including established parameters, namely patient age and volume of the contrast-enhancing tumor at initial diagnosis. While patient age remained an important prognostic marker in our study, the contrast-enhancing tumor volume failed to be of prognostic significance when including the simulated tumor volume in 2 large, independent cohorts (TCGA and UCSF, combined > 500 glioblastoma patients). Of note, by itself, contrast-enhancing tumor volume was significantly associated with survival in both cohorts. On the other hand, the simulated tumor volume with a cutoff value of 0.5 tumor cell density correlated significantly with survival in the multivariable setting. This might be due to inclusion of presumably relevantly tumor-infiltrated areas beyond the contrast-enhancing tumor margins. Effectively, this assumption is supported by studies that confirmed the non-enhancing, FLAIR-hyperintense infiltration volume on preoperative MRI as an important independent prognostic factor.^33,34

Interestingly, in the multivariable Cox regression, the parameter ratio D_w/ρ had a hazard ratio of 1.9 (P < .05) in TCGA, indicating that an increase in D_w/ρ (either through a larger D_w, or lower ρ) translates into shorter survival. In UCSF, we observed a similar trend (i.e. HR > 1), although this was not significant (which might be due to different input images or different therapeutic approaches that affect survival). While the negative correlation with the proliferation coefficient might seem counterintuitive at first sight, it is in accordance with some studies that find no significant correlation between the histopathological proliferation index KI-67 and survival in glioblastoma cases.^35–37

Importantly, since there is an expectable correlation between the simulated tumor volume and the parameter ratio D_w/ρ we analyzed those predictors in separate Cox regression models. Furthermore, there is a weaker correlation between the actual contrast-enhancing tumor volume and the simulated tumor volume as well as D_w/ρ (0.31 and ^_0.2) which might affect the Cox proportional hazard models to a small degree.

Connecting Growth Model Parameters with Tumor Genetics

In recent years, a pathway-level understanding of genomic perturbation has become a promising approach to assessing the underlying genetic mechanisms in tumorigenesis. Pilot studies, such as by TCGA, made clear that in glioblastoma, different genomic alterations, or aberrant expression in different genes often result in alteration of the same pathways.¹⁶ Here, we used PARADIGM which can integrate any number of genomic and functional genomic datasets to infer the molecular pathways altered in a patient sample and present them in the form of factor graphs.²⁵ Utilizing mRNA data from the TCGA dataset we created pathway activation maps for 118 glioblastoma patients and compared them with our parameter ratio D_w/ρ. Interestingly, we observed a significant, negative correlation between D_w/ρ and typical pathways that are associated with cell proliferation, regulation of cell‐cell adhesion, cell migration, or mitosis. This indicates that the proliferation coefficient ρ seems to have a strong correlation with typically altered pathways in tumorigenesis of GBM. To the best of our knowledge, this is the first study, that found a significant correlation between real genetic data and MRI-derived parameters of a mathematical tumor growth model, serving as a promising biological validation of this model.

Toward Personalized Radiotherapy Design

Previous studies have indicated that the common CTV in radiotherapy design of glioblastoma might not sufficiently cover all relevant tumor infiltrated areas.^5,38 One reason for that might be that current radiotherapy planning utilizes a uniform, isotropic margin around the resection cavity to define the CTV. However, the growth pattern of GBM is known to be anisotropic due to anatomical boundaries and the fact that its cells infiltrate gray matter much less than white matter.^6,39 Unfortunately, to this point, standard MRI sequences are not capable of capturing areas of lower tumor infiltration, thus limiting an exact target delineation for radiotherapy planning. Therefore, a growing branch of radiological research is focused on finding means to detect microscopic tumor infiltration by implementing advanced imaging modalities, such as amino acid PET, chemical exchange saturation transfer MRI, or whole-brain MR spectroscopy.^40–43 In addition to this, several studies suggest using mathematical models such as Fisher–Kolmogorov model, to estimate the tumor infiltration pathways beyond the regions of contrast-enhancing tumor, and use this information to guide radiotherapy dose distribution.^6,26,44,45

All these models rely on isolines of the simulated tumor cell density to define the margins of the target volume. However, these isolines can be chosen differently. For example, Unkelbach et al. decided on an isoline that encompasses the same total volume as the CTV used in clinical practice.⁶ For our study, we decided to test a few different isolines of estimated tumor cell density (c_t = {0.33, 0.5, 0.66, 0.75}) to find the best recurrence coverage with a total volume comparable to the manually delineated CTV, or even smaller, in order to keep radiation toxicity small. The results indicated that a cutoff value of 0.5 tumor cell density provides a reasonable tradeoff between a significant improvement of recurrence coverage while having a comparable volume. However, even with a not statistically significant enlarged target volume of 1.2 times the standard CTV, those model-based radiotherapy plans would need to be quality-proofed for practicability in terms of meeting dose constraints to the brain and organs at risk. Notably, as Figure 2 (and additional file 3 in the supplements) illustrates, there are many cases in which there was no improvement of recurrence coverage at all. This is due to the fact that most recurrences occur close to the resection cavity and are therefore already covered by the standard CTV. The most promising advantage of model-derived radiation plans might therefore be to prevent more distant recurrences by warping the plan toward areas with higher estimated tumor cell density.

Limitations of the Study

Radiotherapy personalization based on a reaction-diffusion model relies on multiple assumptions, for example, an initial hypothesis about the tumor cell density distribution over the whole brain at diagnosis time. Since histopathological validation is impossible for that purpose, assumptions about the relation between tumor cell density and MRI features have to be utilized. It is common practice to use thresholds for the tumor cell density following the outlines of the enhancing tumor core, and of the FLAIR-hyperintense tumor portions comprising of non-enhancing tumor and peritumoral edema and gliosis. Upon agreement with experienced neuropathologists from our institution, we decided on cutoff values of c_enhancing = 0.6, and c_edema = 0.3, for our model. However, in a recent postmortem histopathological validation study, Martens et al. invalidated the assumption of a tumor cell density iso-value at the edema outlines.⁴⁶ Although our validation results already support the biological accuracy of our model, some assumptions might need to be adapted in the future to facilitate more accurate results.

In terms of methodology, the registration steps from the postoperatively planned CTV to the recurrence image might be a source of error in accurately evaluating recurrence coverage. In general, comparing standard and model-based CTVs by utilizing recurrence coverage as a proxy is just one vivid way to assess their potential superiority.

Further, since our growth model relies on the initial tumor morphology in the preoperative setting to estimate tumor cell density in the whole brain, it cannot be applied to postoperative images. Therefore, this method of defining the CTV is conceptually different from the standard radiotherapy planning on postoperative images.

Importantly, due to the retrospective study design, we have to acknowledge the fact that all patients have been irradiated according to the standard radiotherapy planning which of course might affect the recurrence pattern. This potential bias could only be solved by a prospective study design.

Furthermore, inclusion of advanced imaging methods, such as DTI or amino acid PET to account for regional differences in the tumor cell density maps might be of great interest. Lastly, a prospective follow-up study needs to prove the benefit of our growth model for patient survival in the real clinical setting.

Conclusions

In this work, we tested a novel, deep learning-based model for radiotherapy personalization for its clinical applicability in the advanced diagnostics and treatment of glioblastoma. Of note, all patients used in this study were independent test data that the model had not seen during training. Our model is easy to deploy in clinical routine due to its few and readily available input sequences of CE-T1 and T2/FLAIR from preoperative images only, and a convincingly fast computing time on widely available hardware. Furthermore, we found significant evidence that the derived model parameters reflect underlying tumor biology as shown by our extensive validation with clinical and biological data. Ultimately, those models could be utilized for improved target delineation in radiotherapy of glioblastoma in the future as our first clinical application study indicates. Our results highlight the potential of this approach to improve patient treatment in the future and constitute a starting point for further clinical implementation of tumor growth modeling.

Funding

German Research Foundation (D.F.G.), Collaborative Research Center (S.F..B) 824 (subproject B12 to F.K., B.M., and B.W). D.F.G., TUM International Graduate School of Science and Engineering, GSC 81, and German Excellence Initiative, Institute for Advanced Studies (to F.K., D.W., B.M., B.W.). EU Marie Sklodowska-Curie program, Translational Brain Imaging Training Network (Grant ID: 765148 to I.E.). DCoMEX (Grant ID 956201 to I.E.). Helmut Horten Foundation (to BM).

Conflict of interest statement

M.M. and F.S. serve as part-time consultants for Novocure GmbH.

Authorship statement

Conceptualization: B.W., B.M. Software: I.E., J.L., F.K., D.W., B.M. Data curation: I.E., B.W., M.M., J.P., J.B. Formal analysis: M.M., B.W. Resources: J.P., J.B., Cl.D., Ch.D., D.B., F.S., J.G. Visualization: M.M., B.W. Writing—original draft: M.M. Writing—review and editing: B.W., I.E. Supervision: B.W., B.M., S.E.C., C.Z.

Author notes

Bjoern Menze and Benedikt Wiestler senior authors contributed equally.

Marie-Christin Metz and Ivan Ezhov contributed equally to this work.