Simulating radiotherapy effect in high-grade glioma by using diffusive modeling and brain atlases.

Applying diffusive models for simulating the spatiotemporal change of concentration of tumour cells is a modern application of predictive oncology. Diffusive models are used for modelling glioblastoma, the most aggressive type of glioma. This paper presents the results of applying a linear quadratic model for simulating the effects of radiotherapy on an advanced diffusive glioma model. This diffusive model takes into consideration the heterogeneous velocity of glioma in gray and white matter and the anisotropic migration of tumor cells, which is facilitated along white fibers. This work uses normal brain atlases for extracting the proportions of white and gray matter and the diffusion tensors used for anisotropy. The paper also presents the results of applying this glioma model on real clinical datasets.

1. Introduction

Cancer causes more than 13% of all deaths worldwide with an estimation of 12 million deaths in 2030 [1, 2]. In United States, 2.5% of cancer deaths are caused by brain tumors [3], with grade IV astrocytoma, namely, glioblastoma multiforme (GBM), accounting for 23.1% of them [4-6]. Despite extensive research on GBM during the last decades, mortality has not changed significantly over the last years, with average life expectancy ranging 12 months after diagnosis and only 4% of treated patients being alive after 5 years [7-12].

Unfortunately, the detection rates of the exact boundaries of GBM with common imaging techniques, such as magnetic resonance imaging (MRI), X-ray computed tomography (CT), and positron emission tomography (PET), are still poor [12-14]. Unlike solid tumors, for which simple exponential or geometric expansion could represent expansion of tumor volume, the GBM growth rate cannot be determined as the classical doubling rate [15], because GBM does not form a solid and compact mass with cells invading the surrounding lesion.

Clearly, new mathematical formulations are necessary when studying this specific glioma case. Since early 90s, there has been a vast amount of research towards simulating and formulating the mechanisms of GBM development, both in macroscopic and microscopic levels. Microscopic models [16-20] study the intracellular biological interactions in cell level, while macroscopic models [21-27] study the tumor behavior, velocity, and mass deformation with using anatomical information derived from medical images.

Gompertzian and volumetric macroscopic models [28-30] did not produce realistic clinical cancer growth representations as they lack taking invasion of tumor cells into account. The class of macroscopic models that has been widely used for simulating this diffusive behavior is diffusive modeling [31]. These models make use of equations based on diffusion reaction schemes of glioma cell concentration, in order to couple the migration of GBM cells (diffusion term) and the proliferation of cells (reaction/source term).

The first diffusive GBM model was introduced in 1995 by Tracqui [21] who used the diffusion reaction equation (DRE) for simulating the spatiotemporal change of glioma cell concentration. The DRE writes as follows: where c(x, t) is the tumor concentration in voxel x = (x, y, z) at time t, D is the diffusion coefficient, ∇ and div are the gradient and divergence operators respectively, and f(c) is the net cell proliferation rate. Later, Swanson et al. [22] altered this equation by taking into account the high velocity of GBM in myelin sheaths by setting different diffusion coefficients in white and gray matter (WM, GM). The new equation was changed to where D(x) is the local diffusion coefficient, being either D or D for x being in GM and WM respectively. In that model, D is assumed to be five times larger than D for high grade glioma, according to the clinical observations of Giese et al. [32, 33].

The next change of the DRE equation was proposed by Price et al. who used T2- and Diffusion Tensor Imaging (DTI) MRIs to identify abnormalities caused by GBM [6, 34]. Therefore, Jbabdi et al. [23] introduced diffusion tensors instead of gradient diffusion coefficients. The new advanced equations for simulating anisotropic growth of GBM in WM writes as follows: where D(x) is the local diffusion tensor which describes the directional tumor cell diffusion, that is, a 3-by-3 symmetric positive definite matrix. Continuing, Clatz et al. [25] used biomechanics to simulate the deformation of brain structures caused by the development of GBM. In 2007, Hogea et al. [27] introduced an advection term in (2), with a two-way coupling with the underlying tissue elastic deformation. Konukoglu et al. [26] adapted the model to successive sessions of MRIs by using Eikonal Equations for coupling the reaction-diffusion dynamics with tumor evolution.

The characteristics of GBM lead to aggressive treatment strategies that most often include surgery, irradiation and chemotherapy. The mathematical approaches that were used for simulating GBM growth inevitably raised the need for incorporating techniques for therapy simulation. Including therapy parameters in GBM models could help the clinicians optimize therapy schemes, as they could predict the response of patients to different therapeutic plans.

Indeed, simulating therapy can be performed by adjusting the proliferation term f(c) of the DRE, which reflects the rate of proliferating cells (B(c)) minus the dead cells due to treatment (T(c)). In order to simulate radiotherapy and chemotherapy, someone should adjust the treatment term T(t). Swanson et al. [35] predicted survival time after resection by estimating the net rates of proliferation and diffusion, and their ratio from pre-treatment gadolinium-enhanced T1-weighted and T2-weighted MR images. Rockne et al. [36] simulated radiotherapy by introducing the radiobiology parameters α and β into T(t), which are interpreted biologically as repairable single and lethal double-strand breaks to the cell's DNA, respectively [37]. Lastly, some important works have been proposed for non diffusive models of glioma, taking into account the biological mechanisms involved in tumor and normal tissue [38, 39].

The most commonly used formalisms of B(c) are either the geometric law

the logistic law: or the Gompertzian law: where ρ is the geometrical growth parameter and c is the maximum tumor cell concentration parameter.

In 2011, Roniotis et al. [40, 41] used the proportions of white and GM for calculating the diffusion coefficient of the DRE, rather than using discrete diffusion coefficients [42]. The proportions of white and GM, as well as the diffusion tensor were extracted by the open SRI-24 brain atlas, thus no DTI processing was needed. This study extends our prior successes with in silico prediction of tumor growth to incorporate the effects of radiotherapy in real clinical datasets. We use the proportional model of Wahl et al. and apply the linear quadratic model [43] for modeling radiotherapy. Compared to previous publications dealing with tumor growth modeling, our approach includes several improvements:

the use of proportional local diffusion coefficients extracted from brain atlases, instead of discrete diffusion coefficients,

the application of isotropic diffusion in GM and anisotropic diffusion in WM on the proportional diffusive model [40],

the application of the linear quadratic model on our previous model [40] for simulating radiotherapy,

the simulation of radiotherapy in multiple fractions by using the linear quadratic model,

the evaluation of radiotherapy simulation by using Jaccard, Dice, and Volume similarity metrics.

2. Materials and Methods

2.1. Diffusive Modeling Using Brain Matter Proportions

This study makes use of the proportional diffusion model presented in [40], which is an extended version of (3) with where W(x) is a symmetric matrix that describes the anisotropy of cell diffusion along the 3-dimensional directions for each position x and D(x) denotes the proportional diffusion coefficient. In this case, W(x) writes as where w(x) ∈ [0,1], i = x, y, z is the directional diffusion weight, which denotes the contribution of the respective direction to the anisotropic migration of GBM cells in position x. Thus, W denotes the contribution of each axis to the direction of white fibers, while D actually denotes the velocity of cells along this direction. D and ρ can be approximated by using the method of Harpold et al. [44].

In order to make a total use of the proportions P of WM and P of GM that can be provided by brain atlases, such as the STI-24 atlas [45, 46], D(x) does not take two discrete values D and D [32]. Instead, the model uses a continuous pattern of coefficients by computing D(x) as follows: This paper extends this idea by improving (5), so as to turn off anisotropic migration in GM tissue. In order to apply anisotropy only to WM tissues, the identity matrix I3 is used for changing (5) to the following formalism:

2.2. The SRI-24 Atlas

The SRI24 atlas [45, 46] is an MRI-based atlas of normal adult human brain anatomy, generated by registering images of 24 normal brains. This atlas provides the proportions P of WM and P of GM in each position of the brain, which can be registered to the clinical datasets of the patients. The data is provided in 155 slices with 240 × 240 pixels, for both stripped and unstripped skull.

Apart from the matter proportions, the atlas provides the dominant eigenvectors of the diffusion tensor. These eigenvectors have been extracted from DTI, by computing the covariance matrices of the distribution of the 3D Gaussian probability which simulates the diffusion of water. Thus, they can be directly used for simulating the anisotropic migration of glioma cells along white matter fibers, as they represent the directions towards which the water diffusion extends mostly.

By mapping clinical MRI data to a brain matter atlas, it is possible to approximate the required tissue information/composition in the tumor area, since this is not possible to extract from the real MRI. Hence, using local diffusion coefficients on real MRI images becomes possible, even if the lesion “hides” the underlying tissue.

2.3. Simulating Radiotherapy Using the Linear Quadratic Model

This paper investigates an extension of the preceding proportional glioma model to include the effects of the radiotherapy using the Linear Quadratic Model. According to the model, the probability of cells surviving S following a single dose of radiation R(x, t) was observed to follow this relationship [47]: where linear α (Gy−1) and quadratic β (Gy−2) are the radiobiology parameters, which are interpreted biologically as repairable single and lethal double-strand breaks to the cell's DNA, respectively [37]. In our case, where the clinical datasets have been extracted from patients with radiation in a number of fractions with the same dose, (9) turns to the following equation [48]: where r is the dose per fraction and R = nr is the total dose for n fractions. In general, fast-proliferating cells, like GBM, have a high tissue response α/β. In our experiment, we use a constant ratio α/β = 10 that is widely used in clinical applications for highly developing cancer and has been extracted from empirical dose response data [36]. Thus, if (4b) is used for approximating cell proliferation, the overall proliferation term f(c) for (3) turns out to be the following:

2.4. Clinical Datasets

The model uses T1- and T2-MRIs (255 × 255 pixel slices) taken from six patients diagnosed with malignant glioblastoma multiforme and cured with radiotherapy. The datasets were provided by Saarland University, Germany, for the needs of the ContraCancrum project [49]. For each dataset there is information provided for the dose and the fractions. Moreover, there are two or more sessions taken on different dates for tracking glioma development and evaluation simulation result.

2.5. The Workflow of Simulation

The collected MRI datasets have different dimensions from the SRI-atlas. Thus, they have to be firstly registered and interpolated to the atlas, before applying the DRE equation. The Optimized Automatic Registration 3D approach by MIPAV is used for registering the clinical datasets to the atlas [50].

After registration, the datasets are delineated by expert radiologists from the University of Saarland, by using DoctorEye. DoctorEye is a platform for annotating, segmenting and visualizing MRI slices [51]. Figure 1(a) presents the 99th slice of one patient dataset after registration and interpolation to the atlas. Figure 1(b) presents the 99th slice of the SRI24 normal brain atlas. Respectively, Figures 1(d), 1(e), and 1(f) present the 99th image of white matter proportion, gray matter proportion and the one dominant eigenvector. Finally, Figure 1(c) presents the 99th slice of the MRI set with the tumor area annotated by the clinician.

Figure 1

(a) the 99th slice of the patient dataset after registration and interpolation, (b) the 99th slice of normal brain atlas, (c) the 99th image after skull stripping and segmentation of tumor, (d) the 99th image of white matter, (e) gray matter proportion, and (f) the one dominant eigenvector images extracted from SRI24 atlas.

In order to solve the DRE, the mathematical framework of Finite Differences that has been presented in [52] is used. More specifically, the Crank Nikolson numerical scheme has been implemented for solving the DRE in three dimensions. Parameters D, ρ, and initial concentration of cells c0, have been extracted from bibliography according to [41, 44]. Thus, we used ρ = 0.012 day−1, D = 0.010 mm2/day, D = 0.002 mm2/day, c0 = 106 cells/mm3. For the solver of finite differences, we used a grid with dimensions ΔX = ΔY = ΔZ = 1 mm and a stable time step of ΔT = 1 h.

3. Results and Discussion

3.1. Radiobiology Parameters

Before applying the diffusive model, it is appropriate to estimate the values of the radiobiology parameters α and β. The method used is that of Rockne et al. [36], by projecting the temporal tumor curves for different values of a on the real datasets. Estimating radius is ideal for spherical tumors however GBM has very complex shape. Thus, instead of tumor radii, the estimated total number of GBM cells is used for identifying the value of a that minimizes the error between simulated and actual target total GBM cells.

Figure 2 presents a set of simulated curves extracted for one specific patient after applying 7 different values of a. The real tumor total cells are marked on the curve for three different MRI provided (the three green points on day 1, 15, and 70). The patient was treated with a total dose of 60 Gy at 2 Gy/fraction and 5 fractions/week for 6 weeks. Treatment starts 21 days after the first MRI and finishes on day 56. Each different curve is extracted by simulating radiotherapy for 0.003 increments of a, starting from zero (no effect) and reaching 0.021 Gy−1. One can see that the optimal value for a, that minimizes the error from the total number of cells is a = 0.006 Gy−1. The same process is followed for approximating the radiobiology parameters for the rest five patients. Table 1 presents the patient specific parameters for all the patients.

Figure 2

The curves depict the total number of GBM cells for one patient for eight values of a ranging from 0 (no effect of radiotherapy) to 0.021 Gy−1 after applying the diffusive model. The same plot shows the three different points of real cell numbers, estimated after segmenting the tumors on three MRI datasets on the respective days. The optimal value of a lies around 0.006 Gy−1.

Table 1

Summary of model parameters for all six patients with glioblastoma multiforme included in the study.

Patient	Age	a(Gy−1)	R(Gy)	r(Gy)	n
1	47	0.114	60	2	30
2	64	0.003	48	3	16
3	52	0.006	60	2	30
4	41	0.211	58	2	29
5	73	0.000	42	3	14
6	56	0.013	60	2	30

3.2. Simulation Result

After defining all simulation parameters, the diffusive model has been applied on the clinical datasets of MRIs. Figure 3 presents the simulated results for the 4th patient 21 days after the end of radiotherapy (100 days after diagnosis). The hot areas depict the resulting simulated tumor projected on the real MRI dataset. The real dataset has been extracted on the same day of simulation.

Figure 3

The simulated cell concentration (hot areas) projected on the real patient data on the 100th day after diagnosis and 21 days after the last fraction of radiotherapy.

For this specific patient, radiotherapy caused shrinkage of the tumor, which cannot be clearly seen on the simulation results of Figure 3. In order to have a better understanding of the shrinkage effect of radiotherapy, we use a line intersecting the point with the maximum concentration. Then, concentration is estimated along this line and presented on a graph. This is applied both on the real initial MRIs (day 1) and the simulated final result. Figure 4(a) shows the initial MRI with the intersection line projected, while Figure 4(b) presents the graph of the concentration of GBM along these line. The red curve shows the initial estimated concentration along this line, while the blue curve shows the simulated concentration of radiotherapy along the same line. In order to have a better visualization of the shrinkage effect, we also provide Figure 6 that depicts the change of the total number of cells in time (for 100 days). It is obvious that the tumor has shrunk considerably.

Figure 4

The graph presents the tumor cell concentration for the initial day (diagnosis day) and the 100th day after diagnosis (21th day after the last fraction of radiotherapy). These values have been calculated along the yellow line shown on the left MRI slice, which comes across the point with the maximum concentration value.

Figure 6

The temporal change of the number of glioma cells for the 4th patient for 100 days. The last radiotherapy session was carried on day 79.

3.3. Evaluation of Model

We have made simulations on a series of six datasets in order to compare the simulated tumor with the final actual tumor, as annotated by the radiologists. For evaluation, we used a scheme with solid metrics and objective comparison. Therefore, the annotated final tumor serves as a golden ground truth and the Jaccard (JC), Dice (DS) and Volume Similarity (VS) metrics are adopted [53].

The evaluation metrics are defined as: where TP (true positive) is the number of tumor voxels belonging to both ground truth and simulated result, FP (false positive) is the number of tumor voxels belonging to simulated result and not belonging to ground truth, and FN (false negative) is the number of tumor voxels belonging to ground truth but not belonging to simulated tumor.

Figure 5 depicts the scatter plots of the three metrics for each patient. The average resulting JC, DS, and VS values for all patients are 95.21%, 98.09% and 99.62%. These values are 2.1%, 4.01%, and 1.89% higher than the respective results for discrete diffusion coefficients and 6.45%, 8.49%, and 4.07% than the uniform tumor growth model.

Figure 5

Scatter plot of JC, DS, and VS (%) for six cases where the proportional, the discrete, and the uniform growth models have been applied.

4. Conclusions

Despite the limited number of datasets, the evaluation results indicate that there might be a slight improvement in using the proposed model on glioblastoma multiforme. This is expected due to the several improvements introduced.

By using the proportional diffusion coefficients, all available information is fully exploited for characterizing tissues, without needing to truncate tissue areas to either white or gray matter.

Moreover, the use of atlases contributes to defining the matter proportions on the areas underlying the tumor area. This could not be true without the use of atlas, as the tumor hides the tissue matter needed for defining the diffusion coefficients at the tumor area. Also, by using anisotropy only on the portions of white matter, instead of all type of tissue, we better approximate the glioma migration along white fibers.

Radiotherapy is usually applied in fractions and the linear quadratic model has been adjusted to this, contrary to previous diffusive models which use the quadratic model for a single dose. By using the formalism of (10), the diffusive model incorporates this idea.