Deep Learning for Reaction-Diffusion Glioma Growth Modeling: Towards a Fully Personalized Model?

Simple Summary

Mathematical tumor growth models have been proposed for decades to capture the growth of gliomas, an aggressive form of brain tumor. However, the estimation of the tumor cell-density distribution at diagnosis and model parameters from partial observations provided by magnetic resonance imaging are ill-posed problems. In this work, we propose a deep learning-based approach to address these problems. 1200 synthetic tumors are first generated using the mathematical model over brain geometries of 6 volunteers. Two deep convolutional neural networks are then trained to (i) reconstruct a whole tumor cell-density distribution and (ii) evaluate the model parameters from partial observations provided in the form of threshold-like imaging contours, with state-of-the-art results. From the estimated cell-density distribution and parameter values, the spatio-temporal evolution of the tumor can ultimately be accurately captured by the mathematical model. Such an approach could be of great interest for glioma characterization and therapy planning.

Abstract

Reaction-diffusion models have been proposed for decades to capture the growth of gliomas, the most common primary brain tumors. However, ill-posedness of the initialization at diagnosis time and parameter estimation of such models have restrained their clinical use as a personalized predictive tool. In this work, we investigate the ability of deep convolutional neural networks (DCNNs) to address commonly encountered pitfalls in the field. Based on 1200 synthetic tumors grown over real brain geometries derived from magnetic resonance (MR) data of six healthy subjects, we demonstrate the ability of DCNNs to reconstruct a whole tumor cell-density distribution from only two imaging contours at a single time point. With an additional imaging contour extracted at a prior time point, we also demonstrate the ability of DCNNs to accurately estimate the individual diffusivity and proliferation parameters of the model. From this knowledge, the spatio-temporal evolution of the tumor cell-density distribution at later time points can ultimately be precisely captured using the model. We finally show the applicability of our approach to MR data of a real glioblastoma patient. This approach may open the perspective of a clinical application of reaction-diffusion growth models for tumor prognosis and treatment planning.

1. Introduction

Gliomas are the most common primary brain tumors, and remain associated with a poor prognosis. Among them, diffuse gliomas are known to be highly infiltrative tumors in which invading tumor cells can be found as far as 4 $\mathrm{c}$$\mathrm{m}$ from the gross tumor [1]. However, the limited sensibility of magnetic resonance imaging (MRI)—the modality of choice for glioma imaging—to changes occurring at the cellular level makes the delineation of the whole tumor extent tedious and often leads to sub-optimal treatment plannings.

Knowledge of the whole tumor cell distribution, within and beyond the outlines of the tumor visible on MRI, would instead allow us to refine surgery or radiation therapy planning. Current standards for the planning of such therapies indeed rely on the addition of a fixed margin to account for tumor infiltration, defining the clinical target volume [2,3]. From an estimated tumor cell distribution, dose deposition could instead be redistributed, with a higher dose delivered to areas more likely to contain tumor cells and, on the other hand, a lower dose delivered to surrounding healthy tissues, while keeping the total dose unchanged [2]. Furthermore, evaluation of the tumor growth dynamics from repeated medical imaging data would also be of great interest to better characterize the tumor, anticipate its growth, and identify its probable migration pathways as well as areas prone to recurrence.

Reaction-diffusion tumor growth models have been studied for decades to circumvent the limitations of current medical imaging techniques and improve treatment planning in gliomas [2,4,5,6,7,8,9,10]. These models rely on partial differential equations (PDEs) to capture the spatio-temporal evolution of a tumor cell-density function over the brain domain, driven by tumor cell migration and proliferation. The most commonly used form involves a logistic reaction term and is referred to as the Fisher equation [11]:

$$\frac{\partial c(\mathit{x},t)}{\partial t}=d{\nabla }^{2}c(\mathit{x},t)+\rho c(\mathit{x},t)\left(1-c(\mathit{x},t)\right)$$

(1)

where $c(\mathit{x},t)$ is the tumor cell-density at location $\mathit{x}$ and time t normalized by the maximum carrying capacity $c_{\max }$ of brain tissues (i.e., $c(\mathit{x},t)\in [0,1],\forall \mathit{x},t$), d is the tumor cell diffusion coefficient, and $\rho$ is the tumor cell proliferation rate. A property of the well-studied Equation (1) is that, under certain conditions and for constant coefficients, it admits a traveling wave solution on the infinite cylinder with propagation speed $v=2\sqrt{d\rho }$, whose profile decays exponentially with decay constant $\lambda =\sqrt{d/\rho }$ as the distance to the origin tends to infinity and for large but finite times t [10,12].

Since their first introduction by Murray and colleagues in the early 1990s [4], reaction-diffusion growth models have been continuously improved to successively integrate (i) a variable tumor cell diffusion rate in white versus gray matter [13] and (ii) an anisotropic diffusion tensor field accounting for the preferred migration of tumor cells along white matter tracts, whose orientation can be assessed by diffusion tensor imaging (DTI) [14]. These improvements led to the formulation that is used throughout this work, presented in Section 2.1. Tumor-induced mass effect [6,7], necrosis, and hypoxia [15,16], as well as the effects of surgery [5], chemotherapy [4,17,18], and radiotherapy [18,19] have also been introduced into reaction-diffusion glioma growth models, but are not considered in this work. For a more detailed overview of reaction-diffusion glioma growth modeling, the reader is referred to [2,4,5,6,7,8,9,10].

Although reaction-diffusion models have shown promising results for patient follow-up and improved radiotherapy planning [2], their clinical application is still prone to severe limitations. Indeed, the estimation of the parameter values and the tumor cell-density distribution at diagnosis time—required to predict the tumor evolution at later times—but also the validation of such models in vivo, implies to extract information on the tumor at the cellular level from medical imaging data. To address this issue, Swanson and colleagues have proposed in [8] to model the imaging function of MRI sequences $I_{\mathrm{seq}}(\mathit{x},t)$—indicating whether a tumor-induced abnormality is visible at location $\mathit{x}$ and time t on the image—as a simple tumor cell-density threshold function:

$$I_{\mathrm{seq}}(\mathit{x},t)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}c(\mathit{x},t)\ge c_{\mathrm{seq}}\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(2)

where $c_{\mathrm{seq}}$ is the tumor cell-density detectability threshold of the sequence. The abnormalities considered in [8] were the hyper-intense enhancing tumor core visible on T1-weighted sequences with injection of gadolinium-based contrast agent (T1Gd) and the peritumor vasogenic edema visible on T2-weighted sequences with or without fluid-attenuated inversion-recovery (T2/T2 FLAIR). Based on these assumptions, the authors suggested that the outlines of these abnormalities would correspond to iso-contours of the tumor cell-density function—i.e., hyper-surfaces along which c has a constant value:

$$\begin{array}{c}\hfill c(\mathit{x},t)=c_{\mathrm{seq}},\forall \mathit{x}\in \partial {\mathsf{\Omega }}_{\mathrm{abn}}\end{array}$$

(3)

where $\partial {\mathsf{\Omega }}_{\mathrm{abn}}$ is the boundary of the visible abnormality.

Building upon this work, Konukoglu and colleagues proposed in [10] a fast marching approach to construct an approximate solution of Equation (1) at imaging time which satisfies Equation (3). This approach has the interesting property of not attempting to dynamically solve the model but seeks to extrapolate the tumor invasion beyond its MR-visible margins within the reaction-diffusion framework, with applications for radiotherapy planning. It has, nevertheless, two major limitations: First, it requires the ability to extract a tumor cell-density iso-contour from the image, from which the whole tumor cell distribution is built. However, we showed in a previous work based on histological data that the outlines of the edema visible on T2 FLAIR MR images do not coincide with a cell-density iso-contour [20]. The proposed explanation is that, due to spatial discontinuities of the tumor cell-density function at interfaces between white and gray matter as well as along the brain domain boundary, Equation (2) does not necessarily imply Equation (3). The second limitation of this approach is that the method still needs to specify the diffusivity and proliferation rate of the tumor, which are unknown at imaging time and need to be adjusted to each tumor.

The estimation of the model parameter values from medical imaging data has also been addressed previously. In [21], the definition of the asymptotic speed of the tumor front $v=2\sqrt{d\rho }$ is used to estimate the tumor cell diffusivity $d_{\mathrm{white}}$ and $d_{\mathrm{gray}}$ in white and gray matter using a fast marching approach. However, the method does not allow us to separate the individual contributions of d and $\rho$ to v, hence $\rho$ is supposed constant for all tumors. Furthermore, this estimation is only valid for large times for which the traveling-wave approximation holds. The approach was then further extended in [9] to take into account the transient speed evolution and the curvature of the tumor front, but still considers a constant $\rho$ value for all tumors. Besides, these fast marching formulations make the assumption that the outlines of the peritumor vasogenic edema visible on T2 MR images correspond to an iso-contour of the traveling wave arrival time function. However, this hypothesis might not be verified due to discontinuities appearing at the brain boundary voxels, which could have been reached long before the imaging time. In [22], a Bayesian approach is used to estimate both the diffusion and proliferation parameters of the model from two imaging contours obtained by Equation (2) at two different times. The method was found to accurately estimate the infiltration length $\lambda =\sqrt{d/\rho }$ of the tumor, but less accurately the tumor front propagation speed $v=2\sqrt{d\rho }$, based on synthetic and real glioblastoma (GBM) MRI data. In [18], parameter values of a two-species reaction-diffusion model incorporating tumor-induced mass effect and response to chemoradiation are estimated based on tumor cell-density distributions derived from longitudinal T1Gd, T2 FLAIR, and diffusion-weighted (DW) MR data, with promising results. However, the cell-density distributions used to initialize the model and fit the parameters were built piecewise from the enhancing/non-enhancing tumor regions delineated on T1Gd/T2 FLAIR images as well as average diffusion coefficient (ADC) maps derived from DW-MR data, and are therefore not guaranteed to be solution of Equation (1) nor to reflect the actual tumor cell distribution.

Tumor source localization is another widely addressed problem in reaction-diffusion glioma growth modeling. In [23], an inverse problem approach is used to estimate the tumor source location from a given tumor cell-density distribution, with promising results. However, to be applicable in clinical practice, the method still requires the ability to derive a whole tumor cell-density distribution from medical imaging data.

Finally, several works have attempted to jointly estimate the tumor source location and model parameters from patient imaging data. In [24], a PDE-constrained optimization approach is used to assess the source location and parameter values of a reaction-diffusion glioma growth model including an additional advection term. The tumor growth model is coupled to a linear elastic model for the tumor-induced mass effect. Two optimization criteria are used in the study: (i) the squared difference between the true and estimated cell-density fields at given imaging times and (ii) the squared distance between the true and estimated position of manually defined landmarks on staggered scans, that are displaced as the surrounding brain tissues are deformed under mass effect. However, the first criterion requires the knowledge of the true tumor cell-density field, which again cannot be derived directly from imaging data. Promising results were obtained for the landmark-based criterion on a real glioma case but strong assumptions are made on the initial cell-density field—supposedly Gaussian—and no ground truth was available to assess the model parameter estimation. In [25], the fast marching approach of Konukoglu and colleagues [9,21] is used to assess the diffusivity ratio $d_{\mathrm{white}}/d_{\mathrm{gray}}$ along with the tumor source location, but a fixed proliferation rate $\rho$ was again considered. More recently, a Bayesian framework has been proposed to simultaneously estimate the tumor source, emergence time, diffusivity, and proliferation rate from a combination of T1Gd, T2 FLAIR, and [¹⁸F]fluoroethyl-L-tyrosine ([¹⁸F]FET) positron emission tomography (PET) images in [26]. However, the study reported that these last three parameters cannot be individually assessed from a single imaging time point. Finally, in [27] a numerical approach is proposed to solve the ill-posed problem of estimating both the tumor initial location as well as its diffusivity and proliferation rate, and is further applied on 206 GBM cases from the BraTS dataset [28] in [29], for which actual parameter values are however not known. Although the method has shown promising results and only requires a single imaging time point, it still has three major limitations: (i) it only allows for the estimation of unscaled dimensionless diffusion and proliferation parameters as the time between tumor emergence and imaging is unknown, which restricts the absolute comparison of the estimated parameter values between tumors and prevents the use of the model as a personalized prediction tool for the tumor evolution over time, (ii) it makes strong assumptions regarding the initial tumor cell-density distribution—a sparse set of Gaussian distributions with identical standard deviation and a maximum density value of 1 over the set—and (iii) it relies on many user-defined parameters, to some of which the method is reported sensitive [27,29].

In addition to their various limitations highlighted hereabove, none of the aforementioned works have jointly addressed the estimation of the tumor cell-density distribution at diagnosis time and individual diffusion and proliferation parameters of the model, which could however made it possible to anticipate the growth of the tumor at further times using the model and thus adapt treatment strategies. Besides, most of these works considered a spatially constant diffusion coefficient in white matter and/or an identical proliferation rate for all tumors, which is not realistic. The introduction of an arbitrary diffusion tensor field $\mathit{D}\left(\mathit{x}\right)$ and a tumor-specific proliferation rate $\rho$ would make the addressed problems even more challenging.

Over the last five years, the advent of deep-learning techniques—and in particular deep convolutional neural networks (DCNNs)—has opened tremendous opportunities in the field of medical imaging, achieving state-of-the-art performance in many image classification and segmentation challenges [30]. One interesting property of deep neural networks is their ability to approximate any function under certain conditions [31]. This property makes the technique attractive for the problems addressed in this work. DCNNs may indeed be used to approximate solutions of PDEs such as Equation (1) over complex domains, and for spatially variable coefficients, as well as to estimate their parameter values from partial observations provided in the form of threshold-like imaging contours.

In this work, we investigate the ability of DCNNs to address common pitfalls encountered in the clinical application of reaction-diffusion glioma growth models. In particular, we focus on the following two tasks:

• Reconstructing a whole brain tumor cell-density distribution compatible with the reaction-diffusion model from a pair of contours obtained through a threshold-like imaging process as in Equation (2) for two different detectability threshold values at a given imaging time. These contours may for example correspond to the outlines of the enhancing core and peritumor vasogenic edema on T1Gd and T2/T2 FLAIR MR images, respectively.
• Estimating the values of the diffusion and proliferation parameters of the model from three imaging contours: (i) two contours obtained for a first detectability threshold value (e.g., the vasogenic edema outlines) at two different imaging times and (ii) a third contour obtained for a second detectability threshold value (e.g., the enhancing core outlines) at the second imaging time.

We demonstrate the ability of DCNNs to perform these tasks accurately based on 1200 synthetic tumors grown over brain geometries derived from the MR data of six healthy subjects. We also show the applicability of our approach on MR data of a real glioblastoma patient.

2. Materials and Methods

2.1. The Reaction-Diffusion Model

The reaction-diffusion tumor growth model that is used throughout this work is described by Equations (4)–(6) [6,25,32]:

$$\left\{\begin{array}{ccc}\frac{\partial c(\mathit{x},t)}{\partial t}=\nabla ·(\mathit{D}(\mathit{x})\mathbf{▽}c(\mathit{x},t))+\rho c(\mathit{x},t)(1-c(\mathit{x},t))\hfill & \hfill \forall \mathit{x}\in \mathrm{\Omega },\forall t>0& \hfill (4)\\ c(\mathit{x},0)=c_0(\mathit{x})\hfill & \hfill \forall \mathit{x}\in \mathrm{\Omega }& \hfill (5)\\ \mathit{D}(\mathit{x})\nabla c(\mathit{x},t)·{\mathit{n}}_{{\partial }_{\mathsf{\Omega }}}\left(\mathit{x}\right)=0\hfill & \hfill \forall \mathit{x}\in \partial \mathrm{\Omega }& \hfill (6)\end{array}\right.$$

where $c(\mathit{x},t)$ is the tumor cell-density at location $\mathit{x}$ and time t normalized by the maximum carrying capacity $c_{\max }$ of brain tissues (i.e., $c(\mathit{x},t)\in [0,1],\forall \mathit{x},t$), $\mathit{D}\left(\mathit{x}\right)$ is the symmetric tumor cell diffusion tensor at location $\mathit{x}$, $\rho$ is the tumor cell proliferation rate, $c_0\left(\mathit{x}\right)$ is the initial tumor cell-density at location $\mathit{x}$, and ${\mathit{n}}_{{\partial }_{\mathsf{\Omega }}}\left(\mathit{x}\right)$ is a unit normal vector pointing outwards the boundary ${\partial }_{\mathsf{\Omega }}$ of the brain domain $\mathsf{\Omega }$ at location $\mathit{x}\in {\partial }_{\mathsf{\Omega }}$. Equation (5) specifies the initial condition of the problem whereas Equation (6) provides no-flux Neumann boundary conditions reflecting the inability of tumor cells to diffuse across ${\partial }_{\mathsf{\Omega }}$.

2.2. MR Data Acquisition

For the needs of this work, 6 healthy volunteers were enrolled for an MRI acquisition comprising a T1 BRAVO (echo time: $\sim 3 \mathrm{m}\mathrm{s}$, repetition time: $\sim 8.3 \mathrm{m}\mathrm{s}$, inversion time: $450 \mathrm{m}\mathrm{s}$, flip angle: 12${}^{°}$, pixel bandwidth: 244 $\mathrm{Hz}$ ${\mathrm{voxel}}^{-1}$, slice thickness/spacing: $1/1$ $\mathrm{m}$$\mathrm{m}$, matrix: $240 \mathrm{voxel}\times 240 \mathrm{voxel}\times 172 \mathrm{voxel}$, field of view: $250 \mathrm{m}\mathrm{m}\times 175 \mathrm{m}\mathrm{m}\times 174 \mathrm{m}\mathrm{m}$), a T2 FLAIR (echo time: $\sim 119 \mathrm{m}\mathrm{s}$, repetition time: $7.2$ $\mathrm{s}$, inversion time: $\sim 2040 \mathrm{m}\mathrm{s}$, flip angle: 90${}^{°}$, bandwidth: 122 $\mathrm{Hz}$ ${\mathrm{voxel}}^{-1}$, phase/slice acceleration factor: 2/2, slice thickness/spacing: $1.4/0.7$ $\mathrm{m}$$\mathrm{m}$, matrix: $256 \mathrm{voxel}\times 256 \mathrm{voxel}\times 252 \mathrm{voxel}$, field of view: $256 \mathrm{m}\mathrm{m}\times 245.8 \mathrm{m}\mathrm{m}\times 176.4 \mathrm{m}\mathrm{m}$), and an EPI-DTI (echo time: $77.1$ $\mathrm{m}$$\mathrm{s}$, repetition time: 7 $\mathrm{s}$, inversion time: 108 $\mathrm{m}$$\mathrm{s}$, flip angle: 90${}^{°}$, bandwidth: $1953.12$ $\mathrm{Hz}$ ${\mathrm{voxel}}^{-1}$, phase/slice acceleration factor: 2/1, multiband factor: 3, slice thickness/spacing: $2/2$ $\mathrm{m}$$\mathrm{m}$, matrix: $120 \mathrm{voxel}\times 120 \mathrm{voxel}$, slices: 72, field of view: $240 \mathrm{m}\mathrm{m}\times 240 \mathrm{m}\mathrm{m}\times 144 \mathrm{m}\mathrm{m}$, directions: 32, b-value: 1000 $\mathrm{s}$ ${\mathrm{mm}}^{-2}$) sequence. To correct for susceptibility-induced distortions (see Section 2.3.1), a second DTI acquisition with reversed phase-encode polarity and only 6 directions was additionally performed. All acquisitions were performed on a 3 $\mathrm{T}$ Signa PET/MR scanner (GE Healthcare, Chicago, IL, USA) with a Nova 32-channel head coil (Nova Medical, Houston, TX, USA).

To illustrate our approach (see Section 2.6), similar T1 BRAVO, T2 FLAIR, and DTI images as well as an additional T1Gd (echo time: $3.2 \mathrm{m}\mathrm{s}$, repetition time: $8 \mathrm{m}\mathrm{s}$, flip angle: 8${}^{°}$, pixel bandwidth: 255 $\mathrm{Hz}$ ${\mathrm{voxel}}^{-1}$, slice thickness/spacing: $1/1$ $\mathrm{m}$$\mathrm{m}$, matrix: $232 \mathrm{voxel}\times 231 \mathrm{voxel}\times 175 \mathrm{voxel}$, field of view: $230 \mathrm{m}\mathrm{m}\times 190.3 \mathrm{m}\mathrm{m}\times 144 \mathrm{m}\mathrm{m}$) image acquired on a 3 $\mathrm{T}$ Achieva scanner (Philips Medical Systems, Eindhoven, The Netherlands) of a 69-year-old male patient with IDH-wildtype GBM were retrospectively analyzed.

2.3. Processing

2.3.1. DTI Data Processing

The acquired DTI data were first corrected for susceptibility-induced distortion, eddy currents, and patient motion using the topup and eddy tools available as part of FSL [33]. A water diffusion tensor field ${\mathit{D}}_{\mathrm{water}}\left(\mathit{x}\right)$ was then reconstructed from the corrected DTI data by least-squares fitting using FSL’s dtifit tool. The whole FSL script used for DTI data processing is available in Appendix A.

2.3.2. Resampling and Registration

The acquired T1 BRAVO, T1Gd (GBM patient only), and T2 FLAIR images as well as the corrected DTI data and the derived water diffusion tensor field were resampled to an isotropic voxel size of $1\mathrm{m}\mathrm{m}\times 1\mathrm{m}\mathrm{m}\times 1\mathrm{m}\mathrm{m}$ by linear interpolation. To correct for patient motion throughout the acquisition, the T1 BRAVO and T2 FLAIR images were rigidly registered to the $b=0$ DTI image used as reference by maximization of their mutual information [34]. All resampling and registration steps were performed using an in-house software in C++ based on ITK [35] and VTK [36].

2.3.3. Skull Stripping

The brain volume was then segmented on the registered T2 FLAIR image using the Otsu thresholding [37] followed by a morphological erosion with structuring element of radius 2 $\mathrm{voxel}$, a largest component extraction, a morphological dilation of radius 2 $\mathrm{voxel}$, a morphological closing of radius 8 $\mathrm{voxel}$, and a morphological dilation of radius 1 $\mathrm{voxel}$.

2.3.4. Brain Tissue Segmentation

The extracted brain domain was further segmented on the registered T1 BRAVO image into white matter, gray matter, and cerebrospinal fluid using an in-house implementation of the MICO intensity-based clustering algorithm comprising a bias field correction step [38]. Manual corrections were further applied to the mis-segmented basal nuclei and falx cerebri. This last step is crucial to prevent the migration of tumor cells between brain hemispheres via routes other than the corpus callosum, as highlighted previously [2]. The segmentation results were finally merged into a single brain map. An example of T2 FLAIR and T1 BRAVO images with the corresponding brain mask and segmented brain map is depicted in Figure 1.

2.3.5. Tumor Segmentation

For illustration purposes (see Section 2.6), the enhancing core and peritumor vasogenic edema were also semi-automatically segmented on the T1Gd and T2 FLAIR images of the GBM patient using combinations of thresholding, connected component extraction, and morphological operations.

2.3.6. Tumor Diffusion Tensor

A tumor diffusion tensor field $\mathit{D}\left(\mathit{x}\right)$ was built piecewise from the DTI-derived water diffusion tensor ${\mathit{D}}_{\mathrm{water}}\left(\mathit{x}\right)$ and the segmented brain domain $\mathsf{\Omega }$ as follows:

$$\mathit{D}\left(\mathit{x}\right)=\left\{\begin{array}{cc}{\mathit{D}}_{\mathrm{white}}\left(\mathit{x}\right)\hfill & \mathrm{if}\mathit{x}\in {\mathsf{\Omega }}_{\mathrm{white}}\hfill \\ {\mathit{D}}_{\mathrm{gray}}\hfill & \mathrm{if}\mathit{x}\in {\mathsf{\Omega }}_{\mathrm{gray}}\hfill \\ \mathbf{0}\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(7)

where ${\mathit{D}}_{\mathrm{white}}\left(\mathit{x}\right)$ and ${\mathit{D}}_{\mathrm{gray}}$ are the tumor cell diffusion tensor fields within the white and gray matter domains ${\mathsf{\Omega }}_{\mathrm{white}}$ and ${\mathsf{\Omega }}_{\mathrm{gray}}$, respectively, and $\mathbf{0}$ is the null tensor, with:

$${\mathit{D}}_{\mathrm{gray}}=\left[\begin{array}{ccc}{d}_{\mathrm{gray}}& 0& 0\\ 0& d_{\mathrm{gray}}& 0\\ 0& 0& d_{\mathrm{gray}}\end{array}\right]$$

(8)

where $d_{\mathrm{gray}}$ is the mean diffusivity of tumor cells in gray matter.

The white matter tumor cell tensor field ${\mathit{D}}_{\mathrm{white}}\left(\mathit{x}\right)$ was built from the DTI-derived water diffusion tensor ${\mathit{D}}_{\mathrm{water}}\left(\mathit{x}\right)$ using the method proposed by Jbabdi and colleagues in [14]. This step is motivated by the observation that, for mechanical [39] and molecular [40] reasons, tumor cells preferentially migrate along rather than across brain fibers, similarly to diffusing water molecules. The method consists of increasing the degree of anisotropy of the water diffusion tensor and scaling its mean diffusivity $\mathrm{MD}=Tr\left({\mathit{D}}_{\mathrm{water}}\right)/3$ to account for the difference in diffusive behavior between tumor cells and water molecules, $Tr\left(\mathit{A}\right)$ being the trace of matrix $\mathit{A}$ [14]:

$${\mathit{D}}_{\mathrm{white}}\left(\mathit{x}\right)=\frac{3d_{\mathrm{white}}}{{\sum }_{i=1}^3{\tilde{\lambda }}_i\left(a\right)}\sum _{i=1}^3{\tilde{\lambda }}_i\left(a\right){\mathit{e}}_i\left(\mathit{x}\right){\mathit{e}}_i^{\top }\left(\mathit{x}\right)$$

(9)

where $d_{\mathrm{white}}$ is the mean diffusivity of tumor cells in white matter, ${\mathit{e}}_i\left(\mathit{x}\right)$ is the ith eigenvector of ${\mathit{D}}_{\mathrm{water}}\left(\mathit{x}\right)$, and ${\tilde{\lambda }}_i\left(a\right)=l_i\left(a\right){\lambda }_i$, with:

$$\begin{array}{c}\left[\begin{array}{c}{l}_1\left(a\right)\\ l_2\left(a\right)\\ l_3\left(a\right)\end{array}\right]=\left[\begin{array}{ccc}a& a& 1\\ 1& a& 1\\ 1& 1& 1\end{array}\right]\left[\begin{array}{c}{c}_l\\ c_p\\ c_s\end{array}\right],\end{array}$$

(10)

$$c_l=\frac{{\lambda }_1-{\lambda }_2}{{\lambda }_1+{\lambda }_2+{\lambda }_3},c_p=\frac{2({\lambda }_2-{\lambda }_3)}{{\lambda }_1+{\lambda }_2+{\lambda }_3},c_s=\frac{3{\lambda }_3}{{\lambda }_1+{\lambda }_2+{\lambda }_3},$$

(11)

where ${\lambda }_i$ is the ith eigenvalue of ${\mathit{D}}_{\mathrm{water}}\left(\mathit{x}\right)$, $a\ge 1$ is a multiplicative factor introduced to increase the anisotropy of the tensor, and $c_l$, $c_p$, and $c_s$ are the linear, planar, and spherical anisotropy measures of ${\mathit{D}}_{\mathrm{water}}\left(\mathit{x}\right)$, respectively.

For reasons that will become clearer in Section 2.4, a unit tumor cell diffusion tensor was built at this stage by fixing $d_{\mathrm{white}}$ to 1. A value of 0.1 was chosen for the $d_{\mathrm{gray}}/d_{\mathrm{white}}$ ratio, as proposed previously in [13] to account for the restricted migration of tumor cells in gray compared to white matter. This ratio was supposed constant among all subjects as it is expected to depend exclusively on the structural organization of healthy white versus gray matter and not on the tumor characteristics. Similarly, the anisotropy factor a was fixed to 10 for all subjects, as suggested in [14]. An example of processed DTI data is depicted in Figure 2. The processed MR data of the 6 volunteers used in this study are publicly available at https://0-doi-org.brum.beds.ac.uk/10.5281/zenodo.6563613 (accessed on 17 May 2022). Further details on these data are available in Appendix B.

2.4. Dataset Synthesis

A synthetic tumor dataset was generated from the processed MR data described hereabove. 200 fictitious tumors were grown over the segmented brain domain of each of the 6 volunteers from randomly picked seeds and parameter values using the reaction-diffusion model in Equations (4)–(6), totaling 1200 synthetic tumors. For each tumor, the cell-density distribution was sampled at four imaging time points $t_1$, $t_2$, $t_3$, and $t_4$.

Each tumor seed consisted of a $3\mathrm{voxel}\times 3\mathrm{voxel}\times 3\mathrm{voxel}$ neighborhood chosen among all segmented white matter voxels whose initial (i.e., at time $t=t_0$) normalized tumor cell-density c was set to 1. For each simulated tumor, an infiltration depth $\lambda =\sqrt{d_{\mathrm{white}}/\rho }$, a tumor front propagation speed $v=2\sqrt{d_{\mathrm{white}}\rho }$, and two imaging time intervals $\mathsf{\Delta }{t}_1$ and $\mathsf{\Delta }{t}_2$ were randomly chosen from uniform distributions (floating point for $\lambda$ and v, integer for $\mathsf{\Delta }{t}_1$ and $\mathsf{\Delta }{t}_2$). The value ranges of the uniform distributions are provided in

Table 1. Value ranges and units of the uniform distributions used to sample the tumor growth model parameters for the generation of the synthetic tumor dataset.

	Min	Max	Units
$\mathsf{\lambda }$	0.5	2.0	$\mathrm{m}\mathrm{m}$
$\mathit{v}$	28.28	48.99	$\mathrm{m}\mathrm{m}{\mathrm{year}}^{-1}$
$\mathsf{\Delta }{\mathit{t}}_{\mathbf{1}}$	90	180	$\mathrm{d}$
$\mathsf{\Delta }{\mathit{t}}_{\mathbf{2}}$	90	180	$\mathrm{d}$

and are of the same order of magnitude as in [19].

Two additional time intervals, $\mathsf{\Delta }{t}_3$ and $\mathsf{\Delta }{t}_4$, were fixed to 90 d for verification purposes (see Section 2.6). Starting from tumor emergence time $t_0$, the time intervals $\mathsf{\Delta }{t}_{1-4}$ univocally define the four imaging time points $t_i=t_0+{\sum }_{j=1}^i\mathsf{\Delta }{t}_j,i=1,\dots ,4$.

For each sampled couple of $(\lambda ,v)$ values, a white matter diffusion rate value $d_{\mathrm{white}}$ and a proliferation rate value $\rho$ were derived as:

$$\begin{array}{cc}\hfill d_{\mathrm{white}}& =\frac{\lambda v}{2}\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \rho & =\frac{v}{2\lambda }\hfill \end{array}$$

(13)

In this manner, a wide diversity of tumors can be uniformly sampled within a realistic range of infiltration depths $\lambda$ and propagation speeds v. Independently sampling $d_{\mathrm{white}}$ and $\rho$ values from uniform distributions may instead have resulted in tumors that are too small (i.e., with too small $d_{\mathrm{white}}$ and $\rho$ values resulting in an empty ${\mathsf{\Gamma }}_1$ contour, see Section 2.5.1 and Section 2.5.2) or too large (i.e., with too large $d_{\mathrm{white}}$ and $\rho$ values resulting in a brain domain almost filled with tumor cells). The empirical joint distribution of the sampled $(\lambda ,v)$ values as well as the corresponding joint distribution of $(d_{\mathrm{white}},\rho )$ and marginal distributions of $d_{\mathrm{white}}$ and $\rho$ are depicted in Figure 3.

For each synthetic tumor, a tumor cell diffusion tensor field $\mathit{D}\left(\mathit{x}\right)$ was then obtained by multiplying the unit (unscaled) diffusion tensor derived as described in Section 2.3.6 by the derived value of $d_{\mathrm{white}}$. As a reminder, the ratio $d_{\mathrm{gray}}/d_{\mathrm{white}}$ was considered constant among tumors in this work (see Section 2.3.6). A tumor was finally grown from the sampled seed, tumor cell diffusion tensor field $\mathit{D}\left(\mathit{x}\right)$, and proliferation rate $\rho$ using the model and the simulated tumor cell distributions at times $t_{1-4}$ were stored. Examples of synthetic tumors are depicted in Figure 4. The corresponding model parameter values are provided in

Table 2. Parameter values used for the tumor simulations in Figure 4.

	${\mathit{d}}_{\mathbf{white}}$ [${\mathrm{mm}}^2{\mathrm{year}}^{-1}$]	$\mathit{\rho }$ [${\mathrm{year}}^{-1}$]	$\mathit{\lambda }$ [$\mathbf{m}\mathbf{m}$]	v [$\mathrm{mm}{\mathrm{year}}^{-1}$]	${\mathit{t}}_1$ [$\mathrm{d}$]	${\mathit{t}}_2$ [$\mathrm{d}$]	${\mathit{t}}_3$ [$\mathrm{d}$]	${\mathit{t}}_4$ [$\mathrm{d}$]
Tumor 1	10.87	31.77	1.71	37.16	175	328	418	508
Tumor 2	15.07	13.95	0.96	29.00	146	316	406	496
Tumor 3	41.49	11.31	0.52	43.33	137	242	332	422

.

The model was solved by a forward Euler finite difference approach using a GPU implementation of the 3D standard stencil referenced in [41] for the computation of the divergence term in Equation (4). A time step $\mathsf{\Delta }t=0.5 \mathrm{d}$ was chosen, ensuring the numerical stability of the method within the sampled parameter range. The processing time was around 1 $\mathrm{m}$$\mathrm{s}$ per iteration on a GeForce GTX 1080 GPU (NVIDIA, Santa Clara, CA, USA), leading to total simulation times in range 0.72–1.08 s per tumor. All simulations were performed using a Python wrapping of our Tumour Growth Simulation ToolKit (TGSTK) written in C++/CUDA language based on VTK [36] and publicly available at https://github.com/cormarte/tgstk (accessed on 17 May 2022). The toolkit documentation can be found at https://cormarte.github.io/tgstk/html (accessed on 17 May 2022).

2.5. Deep Convolutional Neural Networks

Two deep convolutional neural networks were implemented and trained to respectively address the cell-density estimation and model parameter estimation problems introduced hereabove. All training steps were performed using the TensorFlow framework (version 2.5.0) [42] in Python on a GeForce RTX 3090 GPU (NVIDIA, USA). The Python code used for dataset generation and network training is available at https://github.com/cormarte/DeepLearningGliomaGrowthModeling (accessed on 17 May 2022).

2.5.1. Cell-Density Estimation

The first problem addressed was to reconstruct a tumor cell-density distribution from (i) two imaging contours—${\mathsf{\Gamma }}_1$ and ${\mathsf{\Gamma }}_2$—obtained through an imaging process described by Equation (2) for two detectability threshold values $c_1$ and $c_2$ at a given imaging time and (ii) a unit (unscaled) tumor cell diffusion tensor field derived from DTI data as presented in Section 2.3.6.

This approach is motivated by the asymptotic properties of the traveling wave solution admitted by Equation (1) for constant coefficients and on the infinite cylinder, whose profile decreases exponentially with decay constant $\lambda$ (see Section 1). For such a solution, the value of $\lambda$ can be trivially estimated given the distance between two cell-density iso-contours, and an approximate tumor cell-density distribution can subsequently be reconstructed for a sufficiently large distance to the tumor core. Here, we assess the ability of deep neural networks to build an approximate solution of Equations (4)–(6) in the more general case of a complex domain and variable anisotropic diffusion tensor field, and from 2 imaging contours obtained by Equation (2)—which do not necessarily coincide with cell-density iso-contours as discussed previously [20].

In a clinical setting, the value of the white matter diffusion rate $d_{\mathrm{white}}$ used to scale the tumor cell diffusion tensor in Equation (9) is unknown, and its estimation will be addressed in the next section. Therefore this information is not considered for this problem. In contrast, the preferred migration directions of tumor cells along the white matter tracts can be assessed from clinical DTI acquisitions as described in Section 2.3.6 and may be used for the estimation of the tumor cell-density distribution. This motivates the introduction of the unit unscaled diffusion tensor field as an input of this problem, in addition to the ${\mathsf{\Gamma }}_1$ and ${\mathsf{\Gamma }}_2$ contours.

To address this problem, a 3D DCNN based on the U-Net architecture [43] was implemented, as it has been successfully applied to many medical imaging problems previously [30,44]. The network consists of 4 down-sampling blocks, 4 up-sampling blocks, and 1 output block. Each down-sampling block is made of 2 convolutional layers with kernel size $3\times 3$ and stride 1, followed by a bias-adding layer and a rectified linear unit (ReLU) activation layer. A convolutional layer with kernel size $2\times 2$ and stride 2 is added at the end of the block to reduce the feature map dimensions by a factor 2. The up-sampling blocks are identical to the down-sampling blocks except that the last convolutional layer is replaced by a transposed convolution layer with kernel size $2\times 2$ and stride 2, followed by a bias-adding layer and a ReLU activation layer to expand the feature map dimensions by a factor 2. Skip connections are added between the output of the second ReLU activation layer of each down-sampling block and the input of the corresponding upsampling block with the same spatial dimensions, implemented as a concatenation operation. The output block has the same structure as the down-sampling blocks except that the last convolutional layer is replaced by a convolutional layer with kernel size $1\times 1$ and stride 1 followed by a bias-adding layer but no activation layer to merge the last 32 feature maps into a single tumor cell-density map. The network architecture with its feature map dimensions is depicted in Figure 5.

The network takes as input tensors of shape $192\times 192\times 128\times 8$ (width × height × depth × channels). The first 2 channels are fed with the two binary regions respectively delimited by ${\mathsf{\Gamma }}_1$ and ${\mathsf{\Gamma }}_2$. These regions were obtained by thresholding each generated tumor cell distribution at the second imaging time point $t_2$ with threshold values $c_1$ and $c_2$ of 0.80 and 0.16, as previously proposed in [8] for the enhancing core and edema outlines, respectively. The last 6 channels correspond to the 6 independent components of the unit (unscaled) tumor cell diffusion tensor (see Section 2.3.6).

To evaluate the generalization ability of the model, the dataset was further split into training and test sets in proportion 83–17% on a patient basis (i.e., the 200 tumors generated from the MR data of the last patient are kept for evaluation purposes). The network was trained using the Adam optimizer [45] (learning rate: 10⁻⁴, ${\beta }_1$: 0.9, ${\beta }_2$: 0.999, $ϵ$: 10⁻⁶) and the mean absolute error (MAE) loss over mini-batches of size 1. Data augmentation was performed by applying random shifts in range $\pm 15\mathrm{voxel}$ in the three spatial dimensions to each batch. Rotations were not applied as they would also imply transformation of the tensor components [46], resulting in longer execution times for on-the-fly augmentation or larger dataset size for offline augmentation. The training was stopped early after no improvement in the test loss for 100 epochs. The network parameter values that provided the best test loss value (MAE = 1.24 × 10 ⁻⁴) were kept, which occurred after 876 epochs.

2.5.2. Parameter Estimation

The second problem addressed is to estimate the value of the model parameters $d_{\mathrm{white}}$ and $\rho$ (or equivalently of the derived parameters $\lambda =\sqrt{d_{\mathrm{white}}/\rho }$, and $v=2\sqrt{d_{\mathrm{white}}\rho }$) from (i) three imaging contours: two imaging contours—${\mathsf{\Gamma }}_1$ and ${\mathsf{\Gamma }}_2$—obtained by Equation (2) for two different threshold values $c_1$ and $c_2$ at imaging time $t_2$ and a third imaging contour ${\mathsf{\Gamma }}_3$ obtained for the same $c_2$ threshold value at the earlier imaging time $t_1$ and (ii) the unit (unscaled) tumor cell diffusion tensor field. The time interval $\mathsf{\Delta }{t}_2$ between $t_1$ and $t_2$ is also considered as an input of the problem.

As for the cell-density estimation problem, the motivation for such inputs lies in the properties of the asymptotic traveling wave solution of Equation (1) (see Section 1), whose profile decay constant $\lambda =\sqrt{d/\rho }$ can be assessed from two cell-density iso-contours at a given time point as mentioned hereabove. In addition, the propagation speed of the tumor front $v=2\sqrt{d\rho }$ can similarly be assessed from the distance between a same cell-density iso-contour taken at two different time points given their temporal spacing. The knowledge of $\lambda$ and v can finally be used to assess the individual values of d and $\rho$. Here again, we assess the ability of deep neural networks to generalize these properties in the case of a complex domain and variable anisotropic diffusion tensor field from 3 threshold-like imaging contours obtained by Equation (2). The same remark as for the previous section holds regarding the possibility of deriving a unit tumor diffusion tensor from DTI data, which can provide additional information for the estimation of the steepness and speed of the tumor front.

The DCNN implemented for this task is a convolutional encoder. The network consists of 6 convolutional down-sampling blocks and a fully connected output block. Each down-sampling block is made of 2 convolutional layers with kernel size $3\times 3$ and stride 1, followed by a bias-adding layer and a ReLU activation layer. A convolutional layer with kernel size $2\times 2$ and stride 2 is added at the end of the block to reduce the feature map dimensions by a factor 2. The output block flattens the $3\times 3\times 2\times 8$ output of the last down-sampling block and concatenates a $1\times 1$ (width × channels) tensor to the flattened vector to feed the imaging time interval $\mathsf{\Delta }{t}_2$. A fully connected layer followed by a bias-adding layer but no activation layer is finally used to merge the last 145 components into 2 scalar values for $\lambda$ and v.

The network takes as input tensors of shape $192\times 192\times 128\times 9$ (width × height × depth × channels). The first 3 channels are fed with the binary regions respectively delimited by ${\mathsf{\Gamma }}_1$, ${\mathsf{\Gamma }}_2$, and ${\mathsf{\Gamma }}_3$. These regions were respectively obtained by thresholding each generated tumor cell distribution at time $t_2$ with threshold values of 0.80 (${\mathsf{\Gamma }}_1$) and 0.16 (${\mathsf{\Gamma }}_2$), and the distribution at time $t_1$ with a value threshold value of 0.16 (${\mathsf{\Gamma }}_3$) [8]. The last 6 channels correspond to the 6 independent components of the unit (unscaled) tumor cell diffusion tensor (see Section 2.3.6). To account for their different value range and scale, the target values of $\lambda$ and v were standardized using the theoretical mean and variance of the respective uniform distributions from which they were sampled.

The same training/test splitting as for the tumor cell-density estimation network was applied to the dataset. The network was trained using the Adam optimizer [45] (learning rate: 10⁻⁴, ${\beta }_1$: 0.9, ${\beta }_2$: 0.999, $ϵ$: 10⁻⁶) and the mean squared error (MSE) loss. Data augmentation was performed by applying random shifts in range $\pm 15\mathrm{voxel}$ in the three spatial dimensions to each input batch. Early stopping was applied if no improvement was observed in the test loss for 100 epochs. The network parameter values that provided the best test loss value (MSE = 6.75 × 10 ⁻²) were kept, which occurred after 628 epochs.

2.6. Verification

To verify and illustrate our approach, we conducted the following numerical experiment: Starting from the tumor cell-density distribution estimated at time $t_2$ from ${\mathsf{\Gamma }}_1$ and ${\mathsf{\Gamma }}_2$ using our first network (Figure 5), as well as the values of $d_{\mathrm{white}}$ and $\rho$ estimated from ${\mathsf{\Gamma }}_1$, ${\mathsf{\Gamma }}_2$, and ${\mathsf{\Gamma }}_3$ using our second network (Figure 6) and Equations (12) and (13), we computed a tumor cell-density distribution using the reaction-diffusion model at times $t_3$ and $t_4$, 90 d and 180 d later, respectively. We then compared the estimated distributions to the actual tumor cell-density distributions at times $t_3$ and $t_4$—i.e., those obtained for the true cell-density distribution at time $t_2$ as well as the true values of $d_{\mathrm{white}}$ and $\rho$—using the MAE computed voxelwise within the $c>0.01$ contour. This latter restriction prevents background or weakly invaded voxels to artificially lower the MAE. The Hausdorff distance $d_{\mathrm{H}}$ and the average symmetric surface distance (ASSD) $d_{\mathrm{S}}$ between the imaging contours obtained from the true and estimated tumor cell-density distributions for threshold values of $c_1$ and $c_2$ were also computed for each test tumor and time point, as given by:

$$\begin{array}{cc}\hfill d_{\mathrm{H}}(A,B)& =max\left\{\underset{b\in B}{max}\left\{\underset{a\in A}{min}d(a,b)\right\},\underset{a\in A}{max}\left\{\underset{b\in B}{min}d(a,b)\right\}\right\},\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill d_{\mathrm{S}}(A,B)& =\frac{1}{\left|A\right|+\left|B\right|}\left(\sum _{b\in B}\underset{a\in A}{min}d(a,b)+\sum _{a\in A}\underset{b\in B}{min}d(a,b)\right),\hfill \end{array}$$

(15)

where $d(a,b)$ is the Euclidian distance between elements a and b, and $\left|\mathit{X}\right|$ is the cardinal of set X.

It should be noted that minor post-processing was applied to the estimated tumor cell-density distributions at time $t_2$ provided by the first network prior to the computation of the densities at times $t_3$ and $t_4$. First, the cell-density of non-brain voxels (i.e., cerebrospinal fluid and background voxels) was set to 0. Indeed, small (∼10⁻⁵) but non-zero values were observed for some of these voxels in the predicted tumor cell-density distributions. Second, maximum densities were clipped to 1 as small overshootings were also occasionally observed. Third, voxels located outside the largest connected region with densities above 1 × 10⁻⁶ were also set to 0 since small local maxima (∼10⁻⁵) were sporadically observed far from the tumor core, which gave rise to new tumor foci throughout the simulation. These post-processing steps allow us to correct for inaccuracies in the non-constrained output of our convolutional network and ensure numerical stability of the reaction-diffusion model solution at later times.

Sensitivity analyses of our approach were also performed by evaluating both network outputs on the test set after application of a systematic variation of $\pm 10\%$ on the threshold values $c_1$ and $c_2$ used to generate the ${\mathsf{\Gamma }}_1$, ${\mathsf{\Gamma }}_2$, and ${\mathsf{\Gamma }}_3$ input contours, reflecting the inherent uncertainties in these values.

Finally, to demonstrate the applicability of our approach in a clinical context, a cell-density map was generated from the retrospective MR data of the GBM patient (see Section 2.2 and Section 2.3). To this extent, the segmented enhancing core and edema regions (see Section 2.3.5) were provided to the first network along with the derived unit tumor cell diffusion tensor.

3. Results

The distribution of the mean absolute error computed over the test set between the true and estimated tumor cell-density distributions at time $t_2$ within the $c>0.01$ contour is summarized by a boxplot in Figure 7 (first plot). Boxplots of the Hausdorff distance and ASSD distributions computed over the test set between the true and estimated imaging contours at time $t_2$ for threshold values $c_1=0.80$ and $c_2=0.16$ are provided in Figure 8 (first plots). The corresponding median values are provided in

Table 3. Median values of the voxelwise mean absolute error (MAE) between the true and estimated cell-density distributions within the $c>0.01$ contour as well as of the Hausdorff distance and average symmetric surface distance (ASSD) between the true and estimated imaging contours for threshold values of $c_1=0.80$ and $c_2=0.16$ computed over the test set for times $t_2$, $t_3$, and $t_4$.

	${\mathit{t}}_2$	${\mathit{t}}_3$	${\mathit{t}}_4$
Median MAE [10−2]	0.96	1.38	2.20
Median Hausdorff (${\mathit{c}}_{\mathbf{1}}$) [$\mathrm{m}\mathrm{m}$]	1.00	1.41	2.24
Median Hausdorff (${\mathit{c}}_{\mathbf{2}}$) [$\mathrm{m}\mathrm{m}$]	1.00	2.00	3.00
Median ASSD (${\mathit{c}}_{\mathbf{1}}$) [$\mathrm{m}\mathrm{m}$]	0.03	0.12	0.19
Median ASSD (${\mathit{c}}_{\mathbf{2}}$) [$\mathrm{m}\mathrm{m}$]	0.02	0.12	0.18

. An example of true and estimated tumor cell-density distributions at time $t_2$ from the test set is depicted in Figure 9 (first column), along with the corresponding absolute error map as well as the true and estimated imaging contours for threshold values $c_1=0.80$ and $c_2=0.16$. Additional examples are provided in Appendix C. All predicted tumor cell-density distributions at time $t_2$ used in Figure 7, Figure 8 and Figure 9 were provided by the first network (Figure 5).

The distributions of the relative error on the values of $\lambda$ and v computed at time $t_2$ over the test set as well as on the values of $d_{\mathrm{white}}$ and $\rho$ derived with Equations (12) and (13) are summarized by boxplots in Figure 10. The corresponding median relative errors were 3.41%, 3.30%, 5.86%, and 2.75% for $\lambda$, v, $d_{\mathrm{white}}$, and $\rho$, respectively. The true versus predicted values of $\lambda$ and v as well as of $d_{\mathrm{white}}$ and $\rho$ from the test set are plotted in Figure 11. The corresponding Lin’s concordance correlation coefficients ($\mathrm{CCC}$) [47] were 0.99, 0.95, 0.97, and 0.99 for $\lambda$, v, $d_{\mathrm{white}}$, and $\rho$, respectively.

As for imaging time $t_2$, the distributions of the mean absolute error computed over the test set between the true and estimated tumor cell distributions at times $t_3$ and $t_4$ within the $c>0.01$ contour are summarized by boxplots in Figure 7 (second and third plot, respectively). Boxplots of the Hausdorff distance and ASSD distributions computed over the test set between the true and estimated imaging contours at times $t_3$ and $t_4$ for threshold values $c_1=0.80$ and $c_2=0.16$ are also provided in Figure 8 (second and third plots). The corresponding median values are provided in Table 3. The true and estimated tumor cell-density distributions at times $t_3$ and $t_4$ are depicted in Figure 9 (second and third column, respectively) for the same test case as for time $t_2$, along with the corresponding absolute error maps as well as the true and estimated imaging contours for threshold values $c_1=0.80$ and $c_2=0.16$. Additional examples are provided in Appendix C. A loss of accuracy in the estimated tumor cell-density distributions over simulated time is observed in Figure 7, Figure 8 and Figure 9 and Table 3. The estimated tumor cell-density distributions at times $t_3$ and $t_4$ used in Figure 7, Figure 8 and Figure 9 and Table 3 were computed using the reaction-diffusion model as described in Section 2.6 from (i) the cell-density distribution predicted at time $t_2$ provided the first network (Figure 5) and (ii) the predicted model parameter values provided by the second network (Figure 6).

The results of the sensitivity analyses are summarized in

Table 4. Median relative errors (MRE) and Lin’s concordance correlation coefficients (CCC) [47] between the true and predicted values of $\lambda$ and v provided by the second network (Figure 6) and of $d_{\mathrm{white}}$ and $\rho$ derived using Equations (12) and (13), computed over the test set for all combinations of $\pm 10\%$ perturbations on the $c_1$ and $c_2$ threshold values used to generate the input contours ${\mathsf{\Gamma }}_1$, ${\mathsf{\Gamma }}_2$, and ${\mathsf{\Gamma }}_3$.

	Perturbation on ${\mathit{c}}_1|{\mathit{c}}_2$
	−10% |−10%		−10% |+10%		+10% |−10%		+10% |+10%
	MRE	CCC	MRE	CCC	MRE	CCC	MRE	CCC
$\mathsf{\lambda }$	12.51%	0.92	18.55%	0.83	25.54%	0.82	18.90%	0.88
$\mathit{v}$	3.48%	0.94	3.88%	0.93	4.62%	0.91	4.12%	0.93
${\mathit{d}}_{\mathrm{white}}$	13.43%	0.90	19.93%	0.81	28.25%	0.80	22.47%	0.86
$\mathsf{\rho }$	13.00%	0.94	20.54%	0.89	17.04%	0.86	13.31%	0.91

and

Table 5. Median values of the voxelwise mean absolute error (MAE) between the true and estimated cell-density distributions within the $c>0.01$ contour as well as of the Hausdorff distance and average symmetric surface distance (ASSD) between the true and estimated imaging contours for threshold values of $c_1=0.80$ and $c_2=0.16$ computed over the test set for times $t_2$, $t_3$, and $t_4$. The predicted cell-density estimations and parameter values used for the calculations were obtained for all combinations of $\pm 10\%$ perturbations on the $c_1$ and $c_2$ threshold values used to generate the input contours ${\mathsf{\Gamma }}_1$, ${\mathsf{\Gamma }}_2$, and ${\mathsf{\Gamma }}_3$ of both networks.

	Perturbation on ${\mathit{c}}_1|{\mathit{c}}_2$
	−10% |−10%			−10% |+10%			+10% |−10%			+10% |+10%
	${\mathit{t}}_{\mathbf{2}}$	${\mathit{t}}_{\mathbf{3}}$	${\mathit{t}}_{\mathbf{4}}$	${\mathit{t}}_{\mathbf{2}}$	${\mathit{t}}_{\mathbf{3}}$	${\mathit{t}}_{\mathbf{4}}$	${\mathit{t}}_{\mathbf{2}}$	${\mathit{t}}_{\mathbf{3}}$	${\mathit{t}}_{\mathbf{4}}$	${\mathit{t}}_{\mathbf{2}}$	${\mathit{t}}_{\mathbf{3}}$	${\mathit{t}}_{\mathbf{4}}$
Median MAE [10−2]	3.10	3.14	3.38	2.94	2.95	3.26	3.70	3.41	3.63	3.86	3.46	3.57
Median Hausdorff (${\mathit{c}}_{\mathbf{1}}$) [$\mathrm{m}\mathrm{m}$]	2.24	2.45	3.32	2.24	2.45	3.32	3.0	3.0	3.32	3.0	3.0	3.67
Median Hausdorff (${\mathit{c}}_{\mathbf{2}}$) [$\mathrm{m}\mathrm{m}$]	1.00	2.00	3.00	1.0	2.0	3.0	1.0	2.0	3.16	1.0	2.0	3.16
Median ASSD (${\mathit{c}}_{\mathbf{1}}$) [$\mathrm{m}\mathrm{m}$]	0.63	0.51	0.39	0.63	0.49	0.39	0.90	0.56	0.40	0.90	0.61	0.43
Median ASSD (${\mathit{c}}_{\mathbf{2}}$) [$\mathrm{m}\mathrm{m}$]	0.12	0.16	0.19	0.11	0.14	0.18	0.12	0.17	0.24	0.11	0.15	0.22

. Performance indices are reported for all possible combinations of $\pm 10\%$ perturbations on the $c_1$ and $c_2$ values used to generate the input ${\mathsf{\Gamma }}_1$, ${\mathsf{\Gamma }}_2$, and ${\mathsf{\Gamma }}_3$ imaging contours of both CNNs.

Finally, the estimated tumor cell-density distribution for the studied GBM patient provided by the first network (see Figure 5) is depicted in Figure 12 along with the T1Gd and T2 FLAIR images with superimposed segmented enhancing core and edema contours, respectively.

4. Discussion

Reaction-diffusion models have been studied for decades to capture the growth of gliomas, but the ill-posedness of their initialization at imaging time and estimation of their parameter values has restrained their use as a personalized predictive clinical tool. In this work, we showed the ability of DCNNs to circumvent these limitations, opening a wide range of opportunities in the field. Our approach only requires (i) deriving a unit diffusion tensor field from clinical DTI data as described herein, accounting for the preferential migration of tumor cells along white matter tracts and (ii) extracting three imaging contours obtained through a cell-density threshold-like process described by Equation (2) for two different threshold values and time points.

Regarding the second requirement, the outlines of the peritumor vasogenic edema and enhancing core have been proposed previously [8], visible on T2 FLAIR and T1Gd MR images acquired in routine for glioma follow-up, respectively. Nevertheless, it is worth noticing that peritumor vasogenic edema does not strictly speaking correspond to a region of tumor cell invasion but results from an accumulation of extracellular fluid originating from tumor-induced alterations of the blood–brain barrier [48,49] and changes in hydrodynamic pressure [50]. Consequently, the T2 FLAIR imaging process might not be accurately described by Equation (2), as also supported by our previous histological analysis in [20]. Furthermore, anti-angiogenic drugs are known to dramatically reduce vasogenic edema without however stopping tumor progression [48]. Therefore, other MR sequences or modalities could be better suited for the estimation of the tumor cell-density distribution and parameters of reaction-diffusion glioma growth models. For instance, ADC maps derived from DW-MRI data could more accurately reflect tumor cell invasion, as proposed in [18,51]. PET imaging with radio-labeled amino acids could also provide additional information to this extent, as suggested in [26,52].

Once the aforementioned prerequisites are met, our approach makes it possible to (i) extrapolate a whole brain-tumor cell-density distribution within and beyond the visible outlines of the tumor that is compatible with the reaction-diffusion model in Equations (4)–(6) and (ii) individually assess the value of the diffusion and proliferation parameters of the model. Extrapolating tumor invasion is of utmost interest for radiotherapy planning since it would allow us to define personalized margins which more accurately target the tumor while avoiding irradiation of the healthy tissues, as previously discussed in [2,10]. The independent assessment of the diffusivity and proliferation parameters of the model is for its part of great interest to better characterize the tumor [22]. The combination of both gives access to a fully personalized tool, initialized from clinical imaging data and allowing us to anticipate the spatial–temporal growth of gliomas. Such a tool could, for example, be of considerable interest for dose fractionation optimization in radiotherapy using a reinforcement learning approach, as used in [53]. Furthermore, as it only depends on post-processed data (binary segmentations and a DTI-derived water diffusion tensor) rather than raw MR data, the proposed approach may be robustly extended to other scanners and centers. In addition, the method is by design robust to variations in the time interval between the two required MR acquisitions since the interval is provided as an input of the second network for the estimation of the model parameters, which makes it well-adapted to the clinical reality.

The proposed method was found to provide accurate estimations of the three-dimensional tumor cell distribution from only two imaging contours at a single time point, with a median voxelwise MAE below 10⁻² within the $c>0.01$ contour—as evaluated on 200 synthetic tumors grown over the real brain domain of a test subject not used for network training. Our method also provided accurate estimates of the individual diffusion and proliferation parameters of the model from three imaging contours extracted from two time points for the same test tumors, with median relative errors of 5.86% and 2.75%, respectively (see Figure 10), and strong concordance ($\mathrm{CCC}\ge 0.95)$ with the true parameter values (see Figure 11). Furthermore, we showed that the spatio-temporal evolution of the tumor cell-density distribution at later time points (90 d and 180 d later) can be accurately captured from the estimated distribution at imaging time and parameter values using the reaction-diffusion model. The ASSD between the true and estimated imaging contours obtained for threshold values of $c_1=0.80$ and $c_2=0.16$ were indeed found to be lower than or equal to the pixel spacing ($1\mathrm{m}\mathrm{m}\times 1\mathrm{m}\mathrm{m}\times 1\mathrm{m}\mathrm{m}$) in most cases (see Figure 8). Nevertheless, a loss of accuracy in the estimated tumor cell-density over simulated time was observed (see Figure 7, Figure 8 and Figure 9 and Table 3), imputed to the amplification of errors originating from uncertainties in the estimated model parameter values and tumor cell-density distribution at imaging time. In particular, artefactual local maxima in the tumor cell-density distributions predicted by the CNN were found to give rise to new tumor foci over time. Post-processing steps were introduced to circumvent these effects (see Section 2.6), but residual artefacts were still observed, resulting in a large Hausdorff distance though small ASSD values for a few isolated cases (see outliers in Figure 8a,b). Our approach was also found to be robust to uncertainties in the tumor cell-density threshold values defining the input imaging contours of both CNNs. Indeed, all combinations of $\pm 10\%$ perturbations on both threshold values used to generate the contours resulted in an increase in median relative error within reasonable ranges of 14.82–22.39% and 10.25–17.79% for the diffusivity and proliferation rate, respectively. Finally, we also demonstrated the applicability of our proposed method to actual MR data of a GBM patient, for which we were able to reconstruct a tumor cell-density distribution compatible with the imaging data. Nevertheless, the lack of biopsy samples combined with the multiple treatments undergone by the patient prevented the validation of the estimated distribution, which was left for a future prospective study.

Compared to the current state-of-the-art approach in [27], our method appears to perform at least as well or better in most cases on the model parameter estimation problem, by comparison of the reported relative errors on the parameter values in [27]. Nevertheless, both methods were not evaluated on the same synthetic tumor dataset and no open-source code was available in [27] for further comparison. Besides, once trained, the CNNs used in this work no longer depend on any arbitrary parameters, as opposed to the method in [27] for which adequate sparsity level and observation operator weighting need to be selected, in addition to the many parameters involved in the numerical scheme (Gaussian standard deviation, coarsening levels, initial solution guesses, tolerance thresholds, termination criteria, …) [27,29]. Our method also performs significantly faster: parameter estimation on a single case is performed within around $0.1$ $\mathrm{s}$ versus 1000 $\mathrm{s}$ for [27], which however remains largely sufficient for the application in both cases. Furthermore, the proposed approach allows for absolute parameter estimation, as opposed to [27] which only provides non-dimensionalized estimates that could not be scaled as the time between tumor emergence and imaging is actually not known. Consequently, comparison of the estimated parameter values between tumors or their application to predict tumor evolution over time using the model are prevented. Nevertheless, contrary to [27], our method requires two imaging time points to compensate for the ill-posedness of the problem and to allow for dimensionalized parameter estimation—but in return makes no explicit assumption on the initial tumor cell distribution. This latter requirement of our approach implies that the tumor diffusivity and proliferation rate remain constant between the scans—which is in any case also implicitly assumed between the tumor emergence and the single scan time for the forward problem described in [27]. As a consequence, our method is expected to be sensitive to any treatment administered between the scans that would significantly impact the tumor model parameters (chemotherapy, radiotherapy) or solving domain (surgery). Finally, as opposed to [27], our method considers a spatially variable tumor cell diffusion tensor accounting for the preferential migration of tumor cells along white matter tracts, but is restricted to monofocal tumors in the present form.

As a future work, tumor-induced mass effect should be further integrated into the reaction-diffusion model since it is known to cause substantial deformations of the brain parenchyma and distortions of the white matter tracts as the tumor grows, which should also be taken into account for accurate treatment planning. Such effects have been previously considered [6,7], but would introduce additional parameters to be estimated. In addition, transient brain deformations would hardly be integrated into a regular grid-based approach such as the finite difference method used in this work without loss of precision. A finite element formulation over an unstructured mesh could be used instead but would be much more computationally expensive—hence less suited for the generation of large high-resolution datasets such as the one described herein. Tumor-induced destruction of the white matter tracts should also be further considered, as an accurate capture of the original orientation of the brain fibers within the tumor region is required for the evaluation of the cell-density distribution and model parameter values using both DCNNs. A solution to this problem has been previously proposed in [25]—though subject to limitations—in which symmetry of the brain is exploited to artificially reconstruct the missing brain fibers. More advanced methods should be investigated in this sense. Necrosis could also be integrated into the model as proposed in [15,16], which would have avoided the counter-intuitive correspondence between the hyper-dense ($c\sim 1$) region of the estimated tumor cell-density distribution and the necrotic area visible on MRI in Figure 12. Furthermore, the deep neural networks presented herein remain little flexible as they would need to be retrained if different imaging threshold values were considered, although transfer learning could be used to benefit from the lower-level features learned herein and avoid retraining the networks from scratch [54]. Ultimately, the threshold values could be fed to the networks along with the binary contours, but this would make the problem even more complex and would therefore require an even larger training dataset. Although real medical imaging data were used in this work, the verification of our approach still relied on healthy subject data. Therefore, the underlying hypothesis was made that the reaction-diffusion model defined by Equation (4)–(6) and used for tumor synthesis is indeed able to accurately capture the growth of real gliomas, which has never been extensively demonstrated so far to the best of our knowledge. Validation of our approach on actual glioma patient data should be further performed, but longitudinal imaging data with stereotactic biopsies of untreated glioma patients remain scarce. Including the effects of treatments into reaction-diffusion models has also been proposed previously [4,5,17,18,19], but again introduces additional parameters, increasing the complexity of the problem. Alternately, the method could be applied to large publicly available datasets such as the BRaTS dataset [28]. Whereas ground-truth cell-density distributions and model parameter values are unknown for such datasets, indirect validation by investigation of the predictive performance of the estimated model parameters for tumor bio-markers could still be performed, as attempted in [29].

This work also highlights the added value of DCNNs for the resolution of ill-posed problems that are hardly solved by classical optimization methods, and provides encouraging results towards the full personalization of reaction-diffusion glioma growth models from medical imaging data, which has remained unsolved for decades.

5. Conclusions

We proposed a deep learning-based approach to simultaneously address the problems of estimating the tumor cell-density distribution at diagnosis and parameter values of a reaction-diffusion glioma growth model from patient magnetic resonance imaging data. We demonstrated the accuracy of our approach on synthetic tumors grown over actual brain domains of healthy volunteers. We also showed the applicability of our method on MR data of a real glioblastoma patient. Our promising results which point towards the full personalization of glioma reaction-diffusion models may open up tremendous possibilities in the field.