HAL Id: hal-01403847
https://hal.inria.fr/hal-01403847
Submitted on 27 Nov 2016
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entic 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 diusion de documents
scientiques 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ê, 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 Batche-
lor, et al.. Personalized Radiotherapy Planning Based on a Computational Tumor Growth Model.
IEEE Transactions on Medical Imaging, Institute of Electrical and Electronics Engineers, 2016, pp.11.
�10.1109/TMI.2016.2626443�. �hal-01403847�
JOURNAL SUBMISSION 1
Personalized Radiotherapy Planning
Based on a Computational Tumor Growth Model
Matthieu L
ˆ
e
1
, Herv
´
e Delingette
1
, Jayashree Kalpathy-Cramer
2
, Elizabeth R. Gerstner
3
,
Tracy Batchelor
3
, Jan Unkelbach
4
, Nicholas Ayache
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 tu-
mor growth model, personalization, uncertainty, segmentation,
glioblastoma
I. INTRODUCTION
H
IGH 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 tis-
sue [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
1
Asclepios Project, Inria Sophia Antipolis, France.
2
Martinos Center for Biomedical Imaging, Harvard-MIT Division of Health
Sciences and Technology, Charlestown, MA, USA.
3
Department of Neurology, Massachusetts General Hospital, Boston, MA,
USA.
4
Department of Radiation Oncology, Massachusetts General Hospital and
Harvard Medical School, Boston, MA, USA.
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.
equation in order to model the impact of chemo or radio-
therapy [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.
JOURNAL SUBMISSION 2
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.
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
τ
1
(resp. τ
2
) corresponding to the T1Gd (resp. T2-FLAIR) abnormality.
(Bottom) Comparison between the clinician segmentation and the contours
from the model.
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].
JOURNAL SUBMISSION 3
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 prin-
cipled 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 abnor-
malities at the first and second time points are generated. Let
S
0
i
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
k
i
k=1,...,K
denote sets of K plausible
segmentations per modality and time point, where each S
k
i
is
a plausible sample from S
0
i
, 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 segmen-
tation approaches which have been proposed for glioblastoma
[20].
III. TUMOR GROWTH MODEL
The tumor growth model is based on the reaction-diffusion
equation,
∂u
∂t
= ∇(D.∇u)
| {z }
Diffusion
+ ρu(1 − u)
| {z }
Logistic Proliferation
(1)
D∇u.
−→
n
∂Ω
= 0 (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 ρ. Equation (2) enforces Neumann boundary conditions
on the brain domain Ω. Following [21], we define the diffusion
tensor as D = d
w
I in the white matter, and D = d
w
/10 I in
the gray matter, where 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 τ
1
the value of the tumor cell density u
corresponding to the frontier of the T1Gd abnormality, and
τ
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
JOURNAL SUBMISSION 4
index λ =
p
D/ρ . By construction of the initialization, the
T1Gd abnormality falls exactly on the threshold τ
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 approxi-
mately 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 θ = (D, ρ), knowing the clinical
segmentations S
0
3
on the T1Gd and S
0
4
on the T2-FLAIR at
the second time point. To cast the problem in a probabilis-
tic framework, we follow the Bayes rule: P (θ|S
0
3
, S
0
4
) ∝
P (S
0
3
, S
0
4
|θ) P (θ). The likelihood is modeled as
P (S
0
3
, S
0
4
|θ) ∝ exp
−
H(D, ρ, S
0
3
, S
0
4
)
2
σ
2
(3)
where H(D, ρ, S
0
3
, S
0
4
) is the 95th percentile of the symmet-
ric Hausdorff distance between the border of the segmentation
S
0
4
, and the isoline at τ
2
of the simulated tumor cell density
u using (D, ρ), and initialized with the segmentation S
0
3
.
We further model the prior as log-uniform and independent
between the parameters,
P (θ) = P (D)P (ρ) (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 λ =
p
D/ρ. 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 θ = (D, ρ), knowing the clinical
segmentations S
0
i
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
0
i
}
i=1,2,3,4
|θ) ∝ exp
−
1
σ
2
P
4
i=2
H
i
(D, ρ, S
0
1
, S
0
i
)
3
!
2
(5)
where H
i
(D, ρ, S
0
1
, S
0
i
) is the 95th percentile of the sym-
metric Hausdorff distance between the border of the segmen-
tation S
0
i
for i = 2, 3, 4, and the isoline of the simulated
tumor cell density u using (D, ρ), and initialized with the
segmentation S
0
1
. We model the prior as described in Section
IV-A.
We sample from the posterior distribution using the Gaus-
sian Process Hamiltonian Monte Carlo (GPHMC) algorithm
first described by [28], and used for tumor growth personal-
ization 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) ∝ P (S
j
i
), and Z
il
= 0 for l 6= j when Z
ij
= 1.
The random variable Z
i
is a measure of the plausibility
of the samples: P (Z
i
) =
Q
K
i=1
P (Z
ij
= 1)
Z
ij
. We are
interested in the posterior probability of the model parameter
θ = (D, ρ, Z
1
, Z
2
, Z
3
, Z
4
), knowing the observations S. We
model the likelihood as
P (S|θ) ∝ exp
−
1
σ
2
P
4
i=2
H
i
(D, ρ, Z
1
, Z
i
)
3
!
2
(6)
where H
i
(D, ρ, Z
1
, Z
i
) is the 95th percentile of the sym-
metric Hausdorff distance between the border of the segmen-
tation indexed by Z
i
, and the isolines of the simulated tumor
cell density u using (D, ρ), and initialized with the contour
selected with Z
1
. We model the prior independent between
the parameters, log-uniform for D and ρ, and uniform for Z
i
(i.e. P (Z
ij
= 1) = 1/K),
P (θ) = P (D)P (ρ)
4
Y
i=1
P (Z
i
) (7)
