# Personalized Radiotherapy Planning Based on a Computational Tumor Growth Model

Matthieu Lê, Hervé Delingette, Jayashree Kalpathy-Cramer, Elizabeth R

Gerstner, Tracy Batchelor, Jan Unkelbach, Nicholas Ayache

► To cite this version:

Matthieu Lê, Hervé Delingette, Jayashree Kalpathy-Cramer, Elizabeth R Gerstner, Tracy Batchelor, et al.. Personalized Radiotherapy Planning Based on a Computational Tumor Growth Model. IEEE Transactions on Medical Imaging, 2016, pp.11. 10.1109/TMI.2016.2626443. hal-01403847

HAL Id: hal-01403847

https://inria.hal.science/hal-01403847

Submitted on 27 Nov 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Personalized Radiotherapy Planning Based on a Computational Tumor Growth Model

Matthieu Lê<fnref n="1" />, Hervé Delingette<fnref n="1" />, Jayashree Kalpathy-Cramer<fnref n="2" />, Elizabeth R. Gerstner<fnref n="3" />, Tracy Batchelor<fnref n="3" />, Jan Unkelbach<fnref n="4" />, Nicholas Ayache<fnref n="1" />

**Abstract**—In this article, we propose a proof of concept for the automatic planning of personalized radiotherapy for brain tumors. A computational model of glioblastoma growth is combined with an exponential cell survival model to describe the effect of radiotherapy. The model is personalized to the magnetic resonance images (MRIs) of a given patient. It takes into account the uncertainty in the model parameters, together with the uncertainty in the MRI segmentations. The computed probability distribution over tumor cell densities, together with the cell survival model, is used to define the prescription dose distribution, which is the basis for subsequent Intensity Modulated Radiation Therapy (IMRT) planning. Depending on the clinical data available, we compare three different scenarios to personalize the model. First, we consider a single MRI acquisition before therapy, as it would usually be the case in clinical routine. Second, we use two MRI acquisitions at two distinct time points in order to personalize the model and plan radiotherapy. Third, we include the uncertainty in the segmentation process. We present the application of our approach on two patients diagnosed with high grade glioma. We introduce two methods to derive the radiotherapy prescription dose distribution, which are based on minimizing integral tumor cell survival using the maximum a posteriori or the expected tumor cell density. We show how our method allows the user to compute a patient specific radiotherapy planning conformal to the tumor infiltration. We further present extensions of the method in order to spare adjacent organs at risk by re-distributing the dose. The presented approach and its proof of concept may help in the future to better target the tumor and spare organs at risk.

**Index Terms**—Radiotherapy planning, computational tumor growth model, personalization, uncertainty, segmentation, glioblastoma

## I. INTRODUCTION

HIGH grade glioma is one of the most common and aggressive types of primary brain tumors. The treatment of high grade glioma usually involves resection when possible, followed by concurrent chemotherapy and radiotherapy.

Previous works on computational growth models for gliomas have focused on reaction-diffusion equations to model cell proliferation and infiltration into surrounding brain tissue [1]. The model has been extended to model response to chemotherapy, surgical resection, and radiotherapy. For instance, a sink term can be added to the reaction-diffusion

<figure><img src="image_1.png" /><figcaption>Figure 1. The clinical segmentation of the T1Gd abnormality (Top, orange line) is used to define the clinical target volume (CTV, white dashed line) as a 2 cm expansion of the segmentation. In clinical settings, 60 Gy is prescribed to the CTV. We propose to personalize the prescription dose (Bottom) to account for tumor infiltration and segmentation uncertainty.</figcaption></figure>

equation in order to model the impact of chemo or radiotherapy [2], [3]. The resection of a brain tumor can also be modeled by deleting the tumor cells in the resected region [4], [5]. More advanced therapy schedules using for instance anti-angiogenic drugs can also be studied with more complex models [6], [7], [8].

In this article, we provide proof of concept of a method for the automatic planning of personalized radiotherapy for glioblastoma (Figure 1). The beneficial impact of radiotherapy for glioblastoma patients has been clearly demonstrated [9], [10]. However, its planning is made difficult by the infiltrative nature of the disease, and the uncertainty in delineating the abnormality in Magnetic Resonance Images (MRI). To account for the tumor infiltration, a margin of 1 to 3 cm is added to the abnormality visible on MRI to define the clinical target volume (CTV) [11] (Figure 1). The exact extent of this margin is left at the discretion of the clinician.

<footnote n="1">¹Asclepios Project, Inria Sophia Antipolis, France.</footnote>

<footnote n="2">²Martinos Center for Biomedical Imaging, Harvard-MIT Division of Health Sciences and Technology, Charlestown, MA, USA.</footnote>

<footnote n="3">³Department of Neurology, Massachusetts General Hospital, Boston, MA, USA.</footnote>

<footnote n="4">⁴Department of Radiation Oncology, Massachusetts General Hospital and Harvard Medical School, Boston, MA, USA.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 2. Summary of the method: the segmentation of the tumor on the different MRIs is used to personalize the tumor growth model. This is combined with a dose response model to define the prescription dose. Finally, the delivered dose is optimized using 9 equally spaced coplanar photon beams. The color code indicates which data is used for the different scenarios: one or two MRI acquisition at two different time points, the clinical segmentations or plausible samples to take into account the segmentation uncertainty.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 3. First time point on the left, second time point on the right. (Top) The proliferative rim is outlined in orange on the T1Gd MRI. (Middle Top) The edema is outlined in red on the T2-FLAIR MRI. The edema encloses the proliferative rim. (Middle Bottom) Tumor cell density computed with the reaction-diffusion model. The black (resp. white) line is the threshold values τ₁ (resp. τ₂ ) corresponding to the T1Gd (resp. T2-FLAIR) abnormality. (Bottom) Comparison between the clinician segmentation and the contours from the model.</figcaption></figure>

In order to account for the infiltrative nature of the tumor,
several studies recently proposed to personalize radiotherapy
planning based on a computational growth model. Corwin *et* al. [12], [13] personalized spherically symmetric doses based
on a 1D reaction-diffusion tumor growth model using the
T1Gd and T2-FLAIR abnormalities radius as observations
[14], [15]. In this framework, they showed that personalizing
the delivered dose could improve therapy in terms of days
gained by the patients. However, this spherically symmetric
assumption prevents taking into account boundaries of the
tumor progression such as the ventricles. Unkelbach *et* al.
[16], [17] studied the optimization of the radiotherapy planning
based on a tumor growth model in order to automatically
define realistic 3D prescription dose distributions, taking into
account the natural boundaries and privileged pathways of the
tumor progression. The proposed planning was personalized
to the patients geometry, but without personalizing the tumor
growth model parameters.

In this article, we extend previous works by personalizing
a 3D tumor growth model in order to define radiotherapy
prescription doses. This allows one to automatically compute
realistic 3D prescription doses conformal to the tumor infiltra-
tion (see Figure 1). Moreover, we study the impact of taking
into account the uncertainty in the different inputs of the model
(segmentations and model parameters). We use a tumor growth
model based on a reaction diffusion equation, which models
the infiltrative spread of tumor cells in the surrounding white
and gray matter. A Bayesian approach is taken to estimate
the posterior distribution over the model parameters based
on the MRIs of the patient. A recently proposed method to
sample plausible image segmentations is used to incorporate
uncertainty in the segmentation of the tumor in the MR
images [18]. The tumor growth model is then combined
with an exponential cell survival model to describe the effect
of radiotherapy. The probability distribution over tumor cell
densities, together with the cell survival model, is used to
define the prescription dose distribution, which is the basis
for subsequent Intensity Modulated Radiation Therapy (IMRT)
planning. The scope of this paper is the personalization of
radiotherapy planning. As such, we focus on patients which
were not treated with surgical resection. The proposed model
could however be extended in order to included the impact of
such therapy following the developments done in [4], [5].
In this article, we consider three different scenarios. In the first one, we only consider a single MRI acquisition of the T1Gd and T2-FLAIR MRI before therapy planning. This scenario is the closest to the clinical setting where radiotherapy planning is usually based on a single MRI acquisition. In the second, we consider two MRI acquisition at two time points for a total of four MRIs: the T1Gd and T2-FLAIR at the first and second time point (see Figure 3). In the third scenario, we include the uncertainty in the segmentation of the abnormality visible on the different MRIs to the personalization strategy. The second and third scenarios are proofs of concept of a method to include additional information to the personalized therapy pipeline. We acknowledge that patients are usually subject to therapy between the two time points, and as such, the growth model personalization is biased by the impact of therapy. Note however that if the therapy does not result in a decrease of the tumor volume, its impact is implicitly taken into account in the personalization of the growth parameter.

Based on those different scenarios, we propose three principled approaches to compute the prescription dose. First, we minimize the surviving fraction of tumor cells after irradiation for the most probable tumor cell density. Second, we minimize the expected survival fraction tumor cells after irradiation. Third, we present an approach to correct the prescription dose to take into account the presence of adjacent organs at risk.

The generation of different plausible segmentations based on the clinical ones is presented in Section II. The forward model of tumor growth is presented in Section III. The personalization method for the three different scenarios is presented in Section IV. The three principled approach for the personalization of the dose response model to define the prescription dose and the IMRT is detailed in Section V. A summary of the method is illustrated in Figure 2. To our knowledge, this is the first work that uses a personalized model of brain tumor growth taking into account the uncertainty in tumor growth parameters and the clinician's segmentations in order to optimize radiotherapy planning.

## II. SEGMENTATION SAMPLES

The T1Gd abnormality, which is the active part of the tumor, and the larger T2-FLAIR abnormality, which is usually called the edema, were segmented by a clinician. In order to take into account the uncertainty in the segmentation, we propose to randomly modify the original clinician segmentations. The method is based on [18], where samples of such segmentations are generated from a high dimensional Gaussian process, as the zero crossing of a level function. The samples are efficiently produced on the regular grid using the separability and stationary properties of the squared exponential covariance function (see [18] for details). The samples take into account the image intensity information using the signed geodesic distance as the mean of the Gaussian process.

Segmentation samples for the T1Gd and T2-FLAIR abnormalities at the first and second time points are generated. Let $S_i^0$ denote the clinical segmentations for the T1Gd and T2-FLAIR abnormalities at the first and second time points, where the index $i = 1, ..., 4$ refers to the 4 available images (see

Figure 3). Let $S_i = \{S_i^k\}_{k=1,...,K}$ denote sets of K plausible segmentations per modality and time point, where each $S_i^k$ is a plausible sample from $S_i^0$, the i-th clinician segmentation.

Figure 4 shows examples of such samples for $K = 5$. The samples automatically respect the boundaries of the tumor progression such as the ventricles, because of the presence of large intensity gradients. The five presented samples per abnormality correspond to an average DICE of 87%, which is comparable to the inter-expert DICE measured in the BraTS Challenge for brain tumors delineation [19]. Comparing the output of the forward tumor growth model with these plausible noisy segmentations allows to include the uncertainty of the original clinician segmentations.

Note that other approaches could allow the handling of segmentation uncertainty. For instance, one could compare the output of the tumor growth model with probabilistic segmentation approaches which have been proposed for glioblastoma [20].

## III. TUMOR GROWTH MODEL

The tumor growth model is based on the reaction-diffusion equation,

$$ \frac{\partial u}{\partial t} = \underbrace{\nabla(D \cdot \nabla u)}_{\text{Diffusion}} + \underbrace{\rho u(1-u)}_{\text{Logistic Proliferation}} \quad (1) $$

$$ D \nabla u \cdot \vec{n}_{\partial \Omega} = 0 \quad (2) $$

Equation (1) describes the spatio-temporal evolution of the tumor cell density $u$, which infiltrates neighboring tissues with a diffusion tensor $D$, and proliferates with a net proliferation rate $\rho$. Equation (2) enforces Neumann boundary conditions on the brain domain $\Omega$. Following [21], we define the diffusion tensor as $D = d_w \mathbb{I}$ in the white matter, and $D = d_w/10 \mathbb{I}$ in the gray matter, where $\mathbb{I}$ is the 3x3 identity matrix. Below, we identify the scalar parameter $d_w$ with $D$.

The solution of the reaction-diffusion equation (1) is a tumor cell density $u$ computed over the whole brain domain. However, parts of the brain that glioblastomas usually do not invade were excluded from the tumor simulation such as the CSF or the cerebellum. In order to relate the tumor cell density $u$ to the MRIs, the frontier of the visible abnormalities is assumed to correspond to a threshold value of the tumor cell density $u$. We note $\tau_1$ the value of the tumor cell density $u$ corresponding to the frontier of the T1Gd abnormality, and $\tau_2$ the value corresponding to the frontier of the T2-FLAIR abnormality (see Figure 3).

The initialization of the tumor cell density $u(t = t_1, x)$ at the time of the first acquisition is of particular importance, as it impacts the rest of the simulation. In this work, the tumor tail extrapolation algorithm described in [22] is used. The method is based on the assumption that the solution of equation (1) at the first time point has converged to its asymptotic, traveling wave type solution. Thereby, the tumor cell density is propagated outward (and inward), starting from the T1Gd segmentation, and drops approximately exponentially with distance. The steepness of the falloff, i.e. the distance at which the cell density drops by a factor $1/e$ is given by the invisibility
index $\lambda = \sqrt{D/\rho}$. By construction of the initialization, the T1Gd abnormality falls exactly on the threshold $\tau_1$ of the tumor cell density at the first time point.

The reaction-diffusion equation is solved using the Lattice Boltzmann Method [21], [23], [24] which allows for easy parallelization and fast computations. On a 1mm×1mm×1mm resampled MRI, simulating 30 days of growth takes approximately 50 seconds on a 2.3Ghz 50 core machine.

Note that this model is an approximation of the complex growth of the disease. For instance, it could be extended in order to include mass effect [25], or a more detailed description of the disease [6]. In other works, this model has been extended to model different types of therapy such as resection [26], [5], chemotherapy [2], or anti-angiogenic therapy [7]. The common approach taken in these works is to add a death term to the reaction-diffusion equation, which allows to model the shrinkage of the tumor due to the therapy. It was also shown in [27] that the personalized parameters of a reaction-diffusion model were good predictors of certain mutations status of the patient.

## IV. PERSONALIZATION

The personalization of the tumor growth model is combined with a dose response model in order to define the radiotherapy planning. We compare three different scenarios. First we only use a single time point (the second acquisition) to personalize the model such that the radiotherapy plan will be defined using a single acquisition, similarly to what is being done in clinic. Second we use two time points in order to personalize the model. The radiotherapy plan will then be defined on the latest acquisition. Third, we use two time points and include the uncertainty in the segmentation.

### A. Scenario 1: One time point only

In this section, we are interested in the posterior probability of the model parameter $\theta = (D, \rho)$, knowing the clinical segmentations $S_3^0$ on the T1Gd and $S_4^0$ on the T2-FLAIR at the second time point. To cast the problem in a probabilistic framework, we follow the Bayes rule: $P(\theta|S_3^0, S_4^0) \propto P(S_3^0, S_4^0|\theta) P(\theta)$. The likelihood is modeled as

$$P(S_3^0, S_4^0|\theta) \propto \exp \left( - \frac{H(D, \rho, S_3^0, S_4^0)^2}{\sigma^2} \right) \quad (3)$$

where $H(D, \rho, S_3^0, S_4^0)$ is the 95th percentile of the symmetric Hausdorff distance between the border of the segmentation $S_4^0$, and the isoline at $\tau_2$ of the simulated tumor cell density $u$ using $(D, \rho)$, and initialized with the segmentation $S_3^0$. We further model the prior as log-uniform and independent between the parameters,

$$P(\theta) = P(D)P(\rho) \quad (4)$$

We sample from the posterior distribution using a Metropolis-Hasting algorithm. Note that this section only uses the initialization algorithm (see Section III) which only depends on the invisibility index $\lambda = \sqrt{D/\rho}$. Note that this

section can be related to the method described in [16], where a single time point is used to propose a dose planning. However, Unkelbach et al. [16] use a nominal value of the invisibility index whereas it is personalized in this scenario. Moreover, the Bayesian methodology allows to take into account the uncertainty in the personalization.

### B. Scenario 2: Two time points

In this section, we are interested in the posterior probability of the model parameter $\theta = (D, \rho)$, knowing the clinical segmentations $S_i^0$ for $i = 1, 2, 3, 4$ on the T1Gd and T2-FLAIR at the first and second time point respectively. In this case, the likelihood is model as

$$P(\{S_i^0\}_{i=1,2,3,4}|\theta) \propto \exp \left( - \frac{1}{\sigma^2} \left( \frac{\sum_{i=2}^4 H_i(D, \rho, S_1^0, S_i^0)}{3} \right)^2 \right) \quad (5)$$

where $H_i(D, \rho, S_1^0, S_i^0)$ is the 95th percentile of the symmetric Hausdorff distance between the border of the segmentation $S_i^0$ for $i = 2, 3, 4$, and the isoline of the simulated tumor cell density $u$ using $(D, \rho)$, and initialized with the segmentation $S_1^0$. We model the prior as described in Section IV-A.

We sample from the posterior distribution using the Gaussian Process Hamiltonian Monte Carlo (GPHMC) algorithm first described by [28], and used for tumor growth personalization in [21].

### C. Scenario 3: Two time points and segmentation uncertainty

In this section, we want to include the uncertainty in the segmentation to the personalization process. We denote the set of plausible segmentations by $S = \{S_i\}_{i=1,2,3,4}$ (see Section II). We introduce the random variables $Z_i = (Z_{i1}, ..., Z_{iK})$ for $i = 1, 2, 3, 4$, which are one-hot binary vectors where $P(Z_{ij} = 1|S) \propto P(S_i^j)$, and $Z_{il} = 0$ for $l \neq j$ when $Z_{ij} = 1$. The random variable $Z_i$ is a measure of the plausibility of the samples: $P(Z_i) = \prod_{i=1}^K P(Z_{ij} = 1)^{Z_{ij}}$. We are interested in the posterior probability of the model parameter $\theta = (D, \rho, Z_1, Z_2, Z_3, Z_4)$, knowing the observations $S$. We model the likelihood as

$$P(S|\theta) \propto \exp \left( - \frac{1}{\sigma^2} \left( \frac{\sum_{i=2}^4 H_i(D, \rho, Z_1, Z_i)}{3} \right)^2 \right) \quad (6)$$

where $H_i(D, \rho, Z_1, Z_i)$ is the 95th percentile of the symmetric Hausdorff distance between the border of the segmentation indexed by $Z_i$, and the isolines of the simulated tumor cell density $u$ using $(D, \rho)$, and initialized with the contour selected with $Z_1$. We model the prior independent between the parameters, log-uniform for $D$ and $\rho$, and uniform for $Z_i$ (i.e. $P(Z_{ij} = 1) = 1/K$),

$$P(\theta) = P(D)P(\rho) \prod_{i=1}^4 P(Z_i) \quad (7)$$
We sample from the posterior distribution using the GPHMC like in Section IV-B. The only difference is that at each iteration, we randomly sample segmentations from the prior $P(Z_i)$.

## V. RADIOTHERAPY PLANNING

In this section, we detail how we use the personalization of the tumor growth model in order to define the best radio-therapy plan at the time of the second acquisition. We start by coupling the growth model with a cell survival model (Section V-A). We then detail how to compute the **prescription doses** in Section V-B, and how to compute the **delivered dose** in Section V-C.

### A. Cell survival

Cell survival after irradiation is often modeled using the linear-quadratic model. In this article, we follow the derivations of [16], and consider the linear approximation of the linear-quadratic model. In this framework, the density of surviving tumor cells $s$ after irradiation with a cumulative dose $d$ in Gray (Gy = Joules / kg) is given by

$$s = u \exp(-\bar{\alpha}d) \quad (8)$$

where $u$ is the tumor cell density before irradiation, and $\bar{\alpha}$ is the radiosensitivity parameter, corrected for the fact that we consider a linear approximation of the linear-quadratic model.

### B. Prescription Dose Optimization

A prescription dose can be defined as the dose minimizing the surviving fraction of tumor cells. This is formally defined as the dose solving the following optimization problem [16],

$$\underset{d}{\text{minimize}} \quad f(d, u) = \sum_{j \in I} u_j \exp(-\bar{\alpha} d_j) \quad (9)$$

$$\text{subject to} \quad \sum_{j \in I} d_j \le d^{\text{int}} \quad (10)$$

$$d_j \ge 0 \quad (11)$$

where $I$ is the set of voxels in the image. Equation (9) aims at minimizing the number of surviving tumor cells. Equation (10) constrains the integral dose to be lower or equal to a user defined value $d^{\text{int}}$, in order to avoid the trivial solution of delivering an infinite dose. The parameter $d^{\text{int}}$ can be defined based on clinical considerations related to the total dose a brain can tolerate. Equation (11) constrains the dose to be positive.

The optimal prescription dose can be found by setting the derivative of the corresponding Lagrangian to zero, resulting in

$$d_i = \max \left[ 0, \frac{1}{\bar{\alpha}} \ln \left( \frac{u_i \bar{\alpha}}{\mu} \right) \right] \quad (12)$$

where $\mu$ is the Lagrange multiplier for the constraint (10). This solution leads to a surviving tumor cell density $s = \mu/\bar{\alpha}$ where the dose is strictly positive, and $s < \mu/\bar{\alpha}$ elsewhere.

<figure><img src="image_4.png" /><figcaption>Figure 4. From top to bottom: segmentation samples for the T1Gd at the first and second time points, and for the T2-FLAIR at the first and second time points. The sample are generated independently for the different time points and modalities. The different colors correspond to the different samples. The original clinical segmentation $S_i^0$ for $i = 1, 2, 3, 4$ is in red on the different modalities.</figcaption></figure>

A local maximum dose constraint of 60 Gy following clinical recommendation can also be included.

The personalization of the tumor growth model provides samples $\{\theta_l\}$ from the posterior distribution $P(\theta|S)$. We propose three principled methods to compute prescription doses based on the computed samples.

1) **MAP Dose:** The MAP (Maximum A Posterior) dose is defined as the dose minimizing the surviving fraction of the most probable tumor cell density denoted as $u(\theta^{\text{MAP}})$. This dose does not take into account the uncertainty in the personalization.

2) **Probabilistic Dose:** The probabilistic dose is defined as the dose minimizing the expectation of the surviving fraction of tumor cells. This expectation can be estimated using samples from the posterior distribution as follows,

$$\mathbb{E}_{\theta} [f(d, u(\boldsymbol{\theta}))] = \mathbb{E}_{\theta} \left[ \sum_{i \in I} u_i(\boldsymbol{\theta}) \exp(-\bar{\alpha} d_i) \right] \quad (13)$$

$$\approx \sum_{i \in I} \frac{1}{N} \sum_{\theta} u_i(\theta) \exp(-\bar{\alpha} d_i) \quad (14)$$

$$\approx \sum_{i \in I} \hat{u}_i \exp(-\bar{\alpha} d_i) \quad (15)$$

$$\approx f(d, \hat{u}) \quad (16)$$

where $\hat{u} = \frac{1}{N} \sum_{\theta} u(\theta)$ is the empirical mean of the tumor cell density. Computing the probabilistic dose is then
<figure><img src="image_5.png" /><figcaption>Figure 5. MAP (Top), mean (Middle) and standard deviation (Bottom) of the tumor cell density at the second time point computed with 100 random samples of the posterior, when taking into account two time points and the uncertainty in the segmentations. From left to right: axial, coronal, and sagittal views. The brainstem is outlined in white, and the confined target volume T is outlined with a dashed white line. The arrows indicate regions of varying uncertainty above the brainstem (see Figure 14).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 6. Normalized histogram of the distribution of the invisibility index $\lambda = \sqrt{D/\rho}$. The distribution using a single time point $t_2$ is more peaked (in black) than using two time points (in red), or two time points and the segmentation uncertainty (in blue).</figcaption></figure>

equivalent to minimizing the fraction of surviving tumor cells using the empirical mean tumor cell density $\hat{u}$.

3) Corrected Dose: The corrected dose is defined as the prescription dose corrected for the presence of neighboring organs at risk (OARs). We minimize the surviving fraction of tumor cell density minus the surviving fraction of the OARs cell density (i.e. we penalize the death of OAR cells) as follows,

$$ \underset{d}{\text{minimize}} \quad f(d, \hat{u}) - \delta f(d, \beta c) \quad (17) $$

where $\beta$ is the standard deviation map of the tumor cell density (Figure 5), $c$ is the cell density of the OARs, and $\delta$ is a factor which weighs the impact of the correction. The term $\beta c$ translates the fact that we only consider the impact of the OARs in the regions of high uncertainty in the tumor cell density. Note that $f(d, \hat{u}) - \delta f(d, \beta c) = f(d, \hat{u} - \delta \beta c)$. Hence taking into account the OARs is equivalent to minimizing the original problem using the corrected tumor cell density $\hat{u} - \delta \beta c$.

### C. IMRT Planning

We optimize an Intensity Modulated Radiation Therapy (IMRT) plan using 9 equally spaced coplanar 6MV photon beams and a piece-wise quadratic objective function, as detailed in [16], [17]. Dose-calculation is performed using the software CERR [29].

## VI. RESULTS

We first present the results for one high grade glioma patient. This patient was not subject to surgical resection, but was under a complex treatment of concurrent chemo- and

radiotherapy. We picked two time points separated by 28 days which revealed a visible growth large enough to conduct our experiments.

The threshold for the T1Gd and T2-FLAIR abnormalities is set to $\tau_1 = 80\%$ and $\tau_2 = 16\%$ respectively [26]. The log-uniform prior is bounded such that $D \in [10^{-4}, 10]$ mm$^2$/days, and $\rho \in [10^{-5}, 10]$ days$^{-1}$, and we use a noise level of $\sigma = 5$mm for the likelihood.

For the scenario 1, 4000 thousand samples are drawn from the posterior with a normal distribution with standard deviation 0.3 for the proposal function, leading to an acceptance rate of 30%. For scenarios 2 and 3, 2000 samples are drawn from the posterior distribution, leading to an acceptance rate of 60%.

Figure 6 shows the histograms of the invisibility index $\lambda = \sqrt{D/\rho}$. Including the second time point, and the uncertainty in the segmentation, increases the uncertainty in the invisibility index. Indeed, the histogram is more peaked when using a single time point than when including the second time point. It also results in a larger invisibility index. Figure 7 shows the samples from the posterior density of the parameters $D$ and $\rho$ knowing the considered segmentations when using two time points (Left) and when including the uncertainty in the segmentation (Middle). We can see that the presence of two close modes in the region of high probability disappears when the uncertainty in the segmentation is considered. The samples reveal an asymmetric posterior distribution where the mode and mean are different, suggesting that the probabilistic and MAP dose distributions will be different as well. Moreover, the histograms of the random variable $Z_i$ for $i = 1, 2, 3, 4$ (Figure 7 Right) show that all the segmentation samples are equally probable.

Figure 3 (Right) shows the most probable tumor cell density taking into account the uncertainty in the segmentation, along with the extracted contours for the T1Gd and T2-FLAIR
abnormalities. Note that the tumor infiltration is approximated in the regions of tumor necrosis. Indeed, clinician sometimes excluded the necrotic core from the segmentation of the T1Gd abnormality. This is usually not the case (inconsistency between clinicians), and we chose in this paper to use the raw clinical segmentation as input. However, this has actually very little impact on the delivered dose because the necrotic core is surrounded by regions targeted with an important dose. In the sequel, we use 100 random samples from the posterior to compute the empirical mean and the standard deviation of the tumor cell density for the three scenarios.

Figure 8 shows the personalized tumor cell density profile extracted along two different lines at the second time point in the case of scenario 3, along with the boundaries of the segmentation samples. This allows to visualize the uncertainty in the computed tumor cell density. Figure 5 shows the axial, coronal, and sagittal views of the 3D empirical mean and standard deviation of the tumor cell density at the second time point for scenario 3. Those two figures highlights the two sources of uncertainty for scenario 3: the segmentation, and the tumor infiltration.

The integral dose constraint is set equal to the total dose a patient would receive during a treatment following standard guidelines [11]. For that, we simulate a clinical target volume (CTV) by expanding the T1Gd abnormality visible on the second time point with a 2 cm margin. To respect the boundaries of the tumor progression - much like a clinician would do - we define the CTV as the isoline of the average tumor cell density using only the second time point which is totally included in a 2 cm expansion of the T1Gd abnormality. The clinical radiotherapy planning prescribes 60 Gy to the CTV, and 0 elsewhere (see Figure 9 Top). The corresponding dose distribution resulting from IMRT planning is shown in Figure 10 (Top). We then set $d^{\text{int}} = 4.4e + 07\text{Gy.mm}^3$ which corresponds to the IMRT dose delivered to the brain tissues (i.e. excluding the skull and cerebrospinal fluid). The radiosensitivity parameter $\bar{\alpha}$ is set to 0.35 1/Gy.

Figure 9 shows the prescription MAP doses in the three scenarios: i) using only the second time point, ii) using the two time points, iii) using the two time points and the segmentation uncertainty. In accordance with the histograms of invisibility index (Figure 6), we can see that the MAP dose using a single time point is more shallow compared to the doses using two time points (see the arrows on the different views of Figure 9). Furthermore, there is almost no difference between scenario 2 and 3, i.e. whether or not segmentation uncertainty is incorporated. Figure 10 shows the corresponding IMRT optimization of the MAP doses. It is apparent that the differences between the doses is largely mitigated by the smoothing effect of IMRT planning. This is confirmed by Figure 13 (Left) which shows the DICE coefficient of the 50 Gy isolines of the different dose distributions before and after the IMRT optimization: the DICE coefficient is greater (on average 95%) after IMRT than before (on average 91%).

Figure 11 shows the probabilistic prescription doses in the three scenarios. In this case, the difference between the scenarios is small (see the DICE coefficients on Figure 13 right). However, we can note that the infiltration of the prescription

dose is greater for the scenarios taking into account the two time points, contrary to what is the case for the MAP doses (see the arrows on the axial view of Figure 11). This is because the larger uncertainty in the invisibility index leads to a smoother falloff of the dose. Moreover, the effect of taking into account the uncertainty in the segmentations can be seen on the coronal view (see the arrows in Figure 11). A part of the tumor near the cerebellum receives higher doses in the third scenario. This is because this tumor is located near boundaries of tumor progression and therefore, the delineation of the segmentation has a big impact. Figure 12 shows that this effect is still present after IMRT optimization (see the arrows on the coronal view of Figure 12).

Figure 14 shows a sagittal view of the corrected dose (Top) for different values of $\delta$, and the corresponding IMRT dose (Bottom). The total amount of dose is the same for the three presented doses since they respect the dose constraint. However, because of the correction factor, the dose prescribed inside the brainstem (outlined in white) is being redistributed to other regions of the brain where the brainstem cell density and the uncertainty in the tumor cell density is lower. The white arrows on Figure 14 and 5 (sagittal view) show how the part of the dose where the tumor cell density is lower and the uncertainty higher (arrow on the right) is reduced whereas the dose where the tumor cell density is higher and the uncertainty lower (arrow on the left) is not redistributed. This translates in reduced delivered dose after the IMRT optimization. This can be more clearly observed by looking at the dose volume histograms on Figure 15. One can see that, with increasing values of $\delta$, the dose delivered to the brainstem is reduced.

The method was applied to a second patient for which the two MRI acquisitions are 31 days apart. Figure 16 shows the histograms of the invisibility index for the three different scenarios. Including the second time point results in a larger uncertainty and lower invisibility index. Including the segmentation uncertainty does not have a noticeable impact on the invisibility index. This results in a prescription MAP dose which is more concentrated to the target volume for the scenarios 2 and 3 compared to scenario 1 (see Figure 17). However - and similarly to the first patient - the prescription probabilistic dose is slightly more spread out for the scenarios 2 and 3 compared to scenario 1 (see Figure 18). This is due to the increase in uncertainty when taking into account the second time point. The figures for the second patient can be found in the supplementary material.

## VII. DISCUSSION

Extensions of the model could be considered in order to include the noise level $\sigma$, the thresholds $\tau_1$, and $\tau_2$, as free parameters of the model. We detail subsequently the motivations behind our choices, and the difficulty which could arise in considering such extensions. The value of $\sigma$ is related to the level of noise in the mean Hausdorff distance. We chose a value of 5mm which is in agreement with the Hausdorff distances corresponding to the Maximum A Posteriori which ranges from 4.8 to 12.4 mm for the presented patients. This level of noise allows to explore the parameter space, and to
<figure><img src="image_7.png" /><figcaption>Figure 7. Posterior density of the joint probability $P(D, \rho|S)$ using only the clinical segmentations (Left) and taking into account the uncertainty in the segmentations (Middle). The colorbar indicates the negative log likelihood of the samples (yellow unlikely, blue very likely). The most probable sample is indicated with the crossing solid red lines, the mean is indicated with the dashed red lines. The histograms of the random variable $Z_i$ are on the right: the bar plot corresponds to the number of time the sample $j$ of the segmentation $i$ has been selected (i.e. $\bar{Z}_{ij} = 1$).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 8. Visualization of the 3D tumor cell density profile extracted along 2 lines (in orange on Figure 4, identified by an orange number). The empirical mean of the tumor cell density is the solid black line, and the shaded area encloses the 10th to 90th percentiles. The colored crosses corresponds to the boundaries of the different segmentations visible on Figure 4 with the same color code.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 10. IMRT MAP doses in Gray for the clinical plan and the three different personalized plans. From top to bottom: clinical plan, using only the second time point, using the two time points, using the two time points and the segmentation uncertainty. From left to right: axial, coronal, and sagittal views. The arrows emphasize the difference of falloff between the different scenarios.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9. Prescription MAP doses in Gray for the clinical plan and the three different personalized plans. From top to bottom: clinical plan, using only the second time point, using the two time points, using the two time points and the segmentation uncertainty. From left to right: axial, coronal, and sagittal views.</figcaption></figure>

focus on a region of interest which is in accordance with the lowest distances found. It was set manually after a few experiments. The noise level $\sigma$ could be considered as a parameter to be tuned during the personalization. However, it raises the question of computing the normalization factor of the likelihood which depends on $\sigma$ in order to properly compute the Metropolis-Hastings acceptance ratio. We followed [30] to set the value of $\tau_1$ and $\tau_2$. The threshold values could be included as a parameter of the model. However, we did not include it in this study for several reasons. First, there is a biological interpretation behind those values (detectable
<figure><img src="image_11.png" /><figcaption>Figure 11. Prescription probabilistic doses in Gray for the clinical plan and the three different personalized plans. From top to bottom: clinical plan, using only the second time point, using the two time points, using the two time points and the segmentation uncertainty. From left to right: axial, coronal, and sagittal views. The arrows emphasize the difference of falloff between the different scenarios.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 13. DICE coefficient of the dose binarized with a 50 Gy threshold for the different scenarios: using only one time point (OTP), using two time points (TTP), using two time points and the segmentation uncertainty (TTPS). The DICE coefficient is presented for the MAP doses (Left) and probabilistic doses (Right), and for the prescription doses (light blue) and the IMRT doses (light pink). One can note that the TTP and TTPS scenarios are the closest, and that the IMRT optimization reduces the differences between the doses.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 14. Prescription (Top) and IMRT (Bottom) doses in Gray for increasing values of $\delta$ (from left to right). We can see that with increasing $\delta$, the dose around the brainstem is re-distributed (see the arrows).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 12. IMRT probabilistic doses in Gray for the clinical plan and the three different personalized plans. From top to bottom: clinical plan, using only the second time point, using the two time points, using the two time points and the segmentation uncertainty. From left to right: axial, coronal, and sagittal views. The arrows emphasize the difference of falloff between the different scenarios.</figcaption></figure>

threshold of tumor cell density), and we think the different possible values cannot be distinguished by the model since, for two sets of thresholds ($\tau_1, \tau_2$) and ($\tau'_1, \tau'_2$), several values of $D$ and $\rho$ will result in similar segmentations. This was investigated in [21] where the models with different values of ($\tau_1, \tau_2$) were personalized. It was noted that changing the thresholds only resulted in adjusted parameters. It was quantified by showing that the statistical evidence of the different models was not noticeably different. Second, this adds complexity and increases the computational cost of the method.

<figure><img src="image_15.png" /><figcaption>Figure 15. Dose volume histogram of the corrected dose for the brainstem (solid line) and the target volume $T$ (dashed line) for different values of $\delta$. The x axis is the dose and the y axis if the percentage of volume targeted with this dose. Increasing the value of $\delta$ reduces the dose delivered to the brainstem while keeping the dose delivered to the target volume $T$ approximately constant.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 16. Normalized histogram of the distribution of the invisibility index $\lambda = \sqrt{D/\rho}$ for the second patient. The distribution using a single time point $t_2$ is more peaked (in black) than using two time points (in red), or two time points and the segmentation uncertainty (in blue).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 18. Prescription probabilistic doses in Gray for the three different personalized plans for the second patient. From top to bottom: clinical plan, using only the second time point, using the two time points, using the two time points and the segmentation uncertainty. From left to right: axial, coronal, and sagittal views.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17. Prescription MAP doses in Gray for the three different personalized plans for the second patient. From top to bottom: clinical plan, using only the second time point, using the two time points, using the two time points and the segmentation uncertainty. From left to right: axial, coronal, and sagittal views.</figcaption></figure>

## VIII. CONCLUSION

We presented the proof of concept for a method combining a computational model of tumor growth with a dose response model in order to optimize radiotherapy planning, which takes into account the uncertainty in the model parameters and the clinical segmentations. We presented and compared three different scenarios. In the first one, we only consider one MRI acquisition before therapy, as it would usually be the case in clinical practice. In the second one, we use two time points in order to personalize the model and plan radiotherapy. In the third one, we include uncertainty in the segmentation process. Based on these scenarios, we proposed three principled approaches to define patient specific dose prescriptions, and discussed the difference between them. The MAP dose minimizes surviving tumor cells after irradiation of the most probable situation, while the probabilistic dose allows one to take into account the uncertainty by minimizing the expected surviving tumor cells. We showed that including a second

time point increased the uncertainty in the invisibility index and resulted in more shallow probabilistic doses. However, the difference between the prescription doses are partly smoothed out by IMRT optimization. We also showed that including the uncertainty in the segmentation did not change the results substantially. Considering the large number of potential sources of uncertainty for radiotherapy planning, we think the most conservative approach is to use the probabilistic dose personalized with time points and segmentation uncertainty. We also proposed an approach to redistribute dose in order to take into account the uncertainty in the tumor cell density, and the presence of neighboring OARs such as the brainstem. We believe that this method could be beneficial in situations when an organ at risk is located in an area that may or may not be infiltrated by tumor cells.

In the future, the inclusion of the fractionation scheme of the delivered dose could be optimized. In this case, the personalization on two different time points before therapy would be crucial in order to estimate the speed of growth of the tumor, and simulate its progression during radiotherapy. The impact of the presented planning should also be further investigated on a larger cohort of patients. To that end, the model should be extended in order to take into account the complex therapy the patient is undergoing. As such, the model could be applied to any presented patient given the therapy schedule. Finally, it should be investigated if more conformal dose delivery techniques such as proton therapy lead to IMRT planning more conformal to the prescribed dose.

## ACKNOWLEDGMENT

Part of this work was funded by the European Research Council through the ERC Advanced Grant MedYMA 2011-291080.

## REFERENCES

[1] J. Murray, *Mathematical biology*. Springer, 2002, vol. 2.
[2] P. Tracqui, G. Cruywagen, D. Woodward, G. Bartoo, J. Murray, and E. Alvord, "A mathematical model of glioma growth: the effect of chemotherapy on spatio-temporal growth," *Cell proliferation*, vol. 28, no. 1, pp. 17-31, 1995.

[3] R. Rockne, E. Alvord Jr, J. Rockhill, and K. Swanson, "A mathematical model for brain tumor response to radiation therapy," *Journal of mathematical biology*, vol. 58, no. 4-5, pp. 561-578, 2009.

[4] K. Swanson, E. Alvord, and J. Murray, "Virtual resection of gliomas: effect of extent of resection on recurrence," *Mathematical and Computer Modelling*, vol. 37, no. 11, pp. 1177-1190, 2003.

[5] E. Stretton, E. Mandonnet, E. Geremia, B. H. Menze, H. Delingette, and N. Ayache, "Predicting the location of glioma recurrence after a resection surgery," in *Proceedings of 2nd International MICCAI Workshop on Spatiotemporal Image Analysis for Longitudinal and Time-Series Image Data (STIA'12)*, ser. LNCS. Nice: Springer, October 2012.

[6] O. Saut, J.-B. Lagaert, T. Colin, and H. M. Fathallah-Shaykh, "A multilayer grow-or-go model for gbm: effects of invasive cells and anti-angiogenesis on growth," *Bulletin of mathematical biology*, vol. 76, no. 9, pp. 2306-2333, 2014.

[7] E. Scribner, O. Saut, P. Province, A. Bag, T. Colin, and H. M. Fathallah-Shaykh, "Effects of anti-angiogenesis on glioblastoma growth and migration: model to clinical predictions," *PloS one*, vol. 9, no. 12, p. e115018, 2014.

[8] F. Raman, E. Scribner, O. Saut, C. Wenger, T. Colin, and H. M. Fathallah-Shaykh, "Computational trials: Unraveling motility phenotypes, progression patterns, and treatment options for glioblastoma multiforme," *PloS one*, vol. 11, no. 1, 2016.

[9] K. Kristiansen, S. Hagen, T. Kollevold, A. Torvik, I. Holme, M. Stat, R. Nesbakken, R. Hatlevoll, M. Lindgren, A. Brun et al., "Combined modality therapy of operated astrocytomas grade III and IV. Confirmation of the value of postoperative irradiation and lack of potentiation of bleomycin on survival time: a prospective multicenter trial of the Scandinavian Glioblastoma Study Group," *Cancer*, vol. 47, no. 4, pp. 649-652, 1981.

[10] I. J. Barani and D. A. Larson, "Radiation therapy of glioblastoma," in *Current Understanding and Treatment of Gliomas*. Springer, 2015, pp. 49-73.

[11] W. Mason, R. Del Maestro, D. Eisenstat, P. Forsyth, D. Fulton, N. Laperrière, D. Macdonald, J. Perry, B. Thiessen, C. G. R. Committee et al., "Canadian recommendations for the treatment of glioblastoma multiforme," *Current Oncology*, vol. 14, no. 3, p. 110, 2007.

[12] D. Corwin, C. Holdsworth, R. C. Rockne, A. D. Trister, M. M. Mrugala, J. K. Rockhill, R. D. Stewart, M. Phillips, and K. R. Swanson, "Toward patient-specific, biologically optimized radiation therapy plans for the treatment of glioblastoma," *PLOS ONE*, 2013.

[13] C. Holdsworth, D. Corwin, R. Stewart, R. Rockne, A. Trister, K. Swanson, and M. Phillips, "Adaptive IMRT using a multiobjective evolutionary algorithm integrated with a diffusion-invasion model of glioblastoma," *Physics in medicine and biology*, vol. 57, no. 24, p. 8271, 2012.

[14] R. Rockne, J. Rockhill, M. Mrugala, A. Spence, I. Kalet, K. Hendrickson, A. Lai, T. Cloughesy, E. Alvord Jr, and K. Swanson, "Predicting the efficacy of radiotherapy in individual glioblastoma patients in vivo: a mathematical modeling approach," *Physics in medicine and biology*, vol. 55, no. 12, p. 3271, 2010.

[15] R. C. Rockne, A. D. Trister, J. Jacobs, A. J. Hawkins-Daarud, M. L. Neal, K. Hendrickson, M. M. Mrugala, J. K. Rockhill, P. Kinahan, K. A. Krohn et al., "A patient-specific computational model of hypoxia-modulated radiation resistance in glioblastoma using 18f-fmiso-pet," *Journal of The Royal Society Interface*, vol. 12, no. 103, p. 20141174, 2015.

[16] J. Unkelbach, B. H. Menze, E. Konukoglu, F. Dittmann, N. Ayache, and H. A. Shih, "Radiotherapy planning for glioblastoma based on a tumor growth model: implications for spatial dose redistribution," *Physics in medicine and biology*, vol. 59, no. 3, p. 771, 2014.

[17] J. Unkelbach, B. H. Menze, E. Konukoglu, F. Dittmann, M. Le, N. Ayache, and H. A. Shih, "Radiotherapy planning for glioblastoma based on a tumor growth model: improving target volume delineation," *Physics in medicine and biology*, vol. 59, no. 3, p. 747, 2014.

[18] M. Lè, J. Unkelbach, N. Ayache, and H. Delingette, "Gpssi: Gaussian process for sampling segmentations of images," in *Medical Image Computing and Computer-Assisted Intervention-MICCAI 2015*. Springer, 2015, pp. 38-46.

[19] B. Menze, A. Jakab, S. Bauer, J. Kalpathy-Cramer, K. Farahani et al., "The multimodal brain tumor image segmentation benchmark (BRATS)," *Medical Imaging, IEEE Transactions on*, vol. 34, no. 10, pp. 1993-2024, 2015.

[20] A. Gooya, K. M. Pohl, M. Bilello, L. Cirillo, G. Biros, E. R. Melhem, and C. Davatzikos, "GLISTR: glioma image segmentation and registration," *IEEE TMI*, vol. 31, no. 10, pp. 1941-1954, 2012.

[21] M. Lè, H. Delingette, J. Kalpathy-Cramer, E. R. Gerstner, T. Batchelor, J. Unkelbach, and N. Ayache, "Bayesian personalization of brain tumor growth model," in *Medical Image Computing and Computer-Assisted Intervention-MICCAI 2015*. Springer, 2015, pp. 424-432.

[22] E. Konukoglu, O. Clatz, P.-Y. Bondiau, H. Delingette, and N. Ayache, "Extrapolating glioma invasion margin in brain magnetic resonance images: Suggesting new irradiation margins," *Medical image analysis*, vol. 14, no. 2, pp. 111-125, 2010.

[23] H. Yoshida and M. Nagaoka, "Multiple-Relaxation-Time LBM for the convection and anisotropic diffusion equation," *Journal of Computational Physics*, vol. 229, no. 20, 2010.

[24] D. Yu, R. Mei, L.-S. Luo, and W. Shyy, "Viscous flow computations with the method of lattice Boltzmann equation," *Progress in Aerospace Sciences*, vol. 39, no. 5, pp. 329-367, 2003.

[25] O. Clatz, M. Sermesant, P.-Y. Bondiau, H. Delingette, S. K. Warfield, G. Malandain, and N. Ayache, "Realistic simulation of the 3-D growth of brain tumors in MR images coupling diffusion with biomechanical deformation," *IEEE TMI*, vol. 24, no. 10, pp. 1334-1346, 2005.

[26] K. Swanson, R. Rostomily, and E. Alvord, "A mathematical modelling tool for predicting survival of individual patients following resection of glioblastoma: a proof of principle," *British journal of cancer*, vol. 98, no. 1, pp. 113-119, 2008.

[27] A. L. Baldock, K. Yagle, D. E. Born, S. Ahn, A. D. Trister, M. Neal, S. K. Johnston, C. A. Bridge, D. Basanta, J. Scott et al., "Invasion and proliferation kinetics in enhancing gliomas predict idh1 mutation status," *Neuro-oncology*, vol. 16, no. 6, pp. 779-786, 2014.

[28] C. E. Rasmussen, "Gaussian processes to speed up hybrid Monte Carlo for expensive Bayesian integrals," vol. 7, 2003, pp. 651-659.

[29] J. O. Deasy, A. I. Blanco, and V. H. Clark, "Cerr: a computational environment for radiotherapy research," *Medical physics*, vol. 30, no. 5, pp. 979-985, 2003.

[30] H. L. Harpold, E. C. Alvord Jr, and K. R. Swanson, "The evolution of mathematical modeling of glioma proliferation and invasion," *Journal of Neuropathology & Experimental Neurology*, vol. 66, no. 1, 2007.