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.


Introduction
Gliomas are the most common primary brain tumours and remain associated with a poor prognosis.Among them, diffuse glioma are known to be highly infiltrative tumours.This, combined with the limited sensibility of magnetic resonance imaging (MRI)-the modality of choice for glioma imaging-, makes the delineation of such tumours tedious and often leads to sub-optimal treatment plannings.
Reaction-diffusion tumour growth models have been studied for decades to circumvent the limitations of current medical imaging techniques and improve treatment planning in gliomas [1][2][3][4][5][6][7][8].These models rely on partial differential equations (PDEs) to capture the spatio-temporal evolution of a tumour cell density function over the brain domain, driven by tumour cell migration and proliferation.The most commonly used form involves a logistic reaction term and is referred to as the Fisher equation [9]: where c(r, t) is the tumour cell density at location r and time t normalised by the maximum carrying capacity c max of brain tissues (i.e.c(r, t) ∈ [0, 1], ∀r, t), d is the tumour cell diffusion coefficient, and ρ is the tumour cell proliferation rate.A property of the well-studied Equation ( 1) is that, under certain conditions and for constant coefficients, it admits a travelling wave solution on the infinite cylinder with propagation speed v = 2 √ d ρ, whose profile decays exponentially with decay constant λ = d/ρ as the distance x to the origin tends to infinity and for large but finite times t [7,10].
Since their first introduction by Murray and colleagues in the early 1990's [1], reaction-diffusion growth models have been continuously improved to successively integrate (i) a variable tumour cell diffusion rate in white versus grey matter [11] and (ii) an anisotropic diffusion tensor field accounting for the preferred migration of tumour cells along white matter tracts, whose orientation can be assessed by diffusion tensor imaging (DTI) [12].These improvements led to the formulation that is used throughout this work, presented Section 2.1.Tumour-induced mass effect [3,4], necrosis, and hypoxia [13,14], as well as the effects of surgery [2], chemotherapy [1,15,16], and radiotherapy [16,17] 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 modelling, the reader is referred to [1][2][3][4][5][6][7][8].
Although reaction-diffusion models have shown promising results for patient follow-up and improved radiotherapy planning [8], their clinical application is still prone to severe limitations.Indeed, the estimation of the parameters and initial condition-i.e. the tumour cell density distribution at diagnosis time-of such models, as well as their validation in vivo, require to extract information on the tumour cell density distribution from medical imaging data.To address this issue, Swanson and colleagues have proposed in [5] to model the imaging function of MRI sequences I seq (r, t)-indicating whether a tumour-induced abnormality is visible at location r and time t on the image-as a simple tumour cell density threshold function: where c seq is the tumour cell density detectability threshold of the sequence.The abnormalities considered in [5] were the hyper-intense enhancing tumour core visible on T1-weighted sequences with injection of gadolinium-based contrast agent (T1Gd) and peritumour vasogenic oedema 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 tumour cell density function: c(r, t) = c seq , ∀r ∈ ∂Ω abn (3) where ∂Ω abn is the boundary of the visible abnormality.
Building upon this work, Konukoglu and colleagues proposed in [7] 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 tumour 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 to be able to extract a tumour cell density iso-contour from the image, from which the whole tumour cell distribution is built.However, we showed in a previous work based on histological data that the outlines of the oedema visible on T2 FLAIR MR images do not coincide with a cell density iso-contour [18].The proposed explanation is that, due to spatial discontinuities of the tumour cell density function at interfaces between white and grey 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 requires to specify the diffusivity and proliferation rate of the tumour, which are unknown at imaging time and need to be adjusted to each tumour.
The assessment of the model parameter values from medical imaging data has also been addressed previously.In [19], the definition of the asymptotic speed of the tumour front v = 2 √ d ρ is used to estimate the tumour cell diffusivity d white and d grey in white and grey matter using a fast marching approach.However, the method does not allow to separate the individual contributions of d and ρ to v, hence ρ is supposed constant for all tumours.Furthermore, this estimation is only valid for large times for which the travelling-wave approximation holds.The approach was then further extended in [6] to take into account the transient speed evolution and the curvature of the tumour front, but still considers a constant ρ value for all tumours.Besides, these fast marching formulations make the assumption that the outlines of the peritumour vasogenic oedema visible on T2 MR images correspond to an iso-contour of the travelling 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 [20], a Bayesian approach is used to estimate both the diffusivity and proliferation rate 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 λ = d/ρ of the tumour but less accurately the tumour front propagation velocity v = 2 √ d ρ, based on synthetic and real glioblastoma (GBM) MRI data.In [16], parameter values of a 2-species reaction-diffusion model incorporating tumour-induced mass effect and response to chemoradiation are estimated based on tumour 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 initialise the model and fit the parameters were built piecewise from the enhancing/non-enhancing tumour 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 tumour cell distribution.
The tumour source localisation is another widely addressed problem in reaction-diffusion glioma growth modelling.In [21], an inverse problem approach is used to estimate the tumour source location from a given tumour cell density distribution with promising results.However, to be applicable in clinical practice, the method still requires the ability to derive a whole tumour cell density distribution from medical imaging data.
Finally, several works have attempted to jointly estimate the tumour source location and model parameters from patient imaging data.In [22], a PDE-constrained optimisation approach is used to assess the source location and parameter values of a reaction-diffusion glioma growth model including an additional advection term.The tumour growth model is coupled to a linear elastic model for the tumour-induced mass effect.Two optimisation 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 estimate 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 tumour 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 [23], the fast marching approach of Konukoglu and colleagues [6,19] is used to assess the diffusivity ratio d white /d grey along with the tumour source location, but a fixed proliferation rate ρ was again considered.More recently, a Bayesian framework has been proposed to simultaneously estimate the tumour source, emergence time, diffusivity, and proliferation rate from a combination of T1Gd, T2 FLAIR, and [ 18 F]Fluoroethyl-L-Tyrosine ([ 18 F]FET) positron emission tomography (PET) images in [24].However, the study reported that these last three parameters cannot be individually assessed from a single imaging time point.However, none of the aforementioned works have jointly addressed the estimation of the initial condition and individual diffusivity and proliferation rate parameters of the model.Besides, most of these works considered a spatially constant diffusivity coefficient in white matter and an identical proliferation rate for all tumours.The introduction of an arbitrary diffusion tensor field D(r) and tumour-specific proliferation rate ρ 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 [25].One interesting property of deep neural networks is their ability to approximate any function under certain conditions [26].This property makes the technique attractive to address the aforementioned limitations.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 assess their parameter values, from partial observations provided in the form of thresholding 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: 1. Reconstructing a whole brain tumour 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 peritumour vasogenic oedema on T1Gd and T2/T2 FLAIR MR images, respectively.
2. 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 oedema 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 1,200 synthetic tumours grown over the brain geometries derived from the MR data of 6 healthy subjects.We also show the applicability of our approach on MR data of a real glioblastoma patient.

Methods
The following sections successively describe the reaction-diffusion model considered, the MR image acquisition and processing steps, the tumour dataset synthesis, the architecture and training procedure of the DCNNs, and the validation methodology adopted for our approach.

The Reaction-Diffusion Model
The reaction-diffusion tumour growth model that is used throughout this work is described by the following set of equations [3,23,27]: where c(r, t) is the tumour cell density at location r and time t normalised by the maximum carrying capacity c max of brain tissues (i.e.c(r, t) ∈ [0, 1], ∀r, t), D(r) is the symmetric tumour cell diffusion tensor at location r, ρ is the tumour cell proliferation rate, c 0 (r) is the initial tumour cell density at location r, and n ∂Ω (r) is a unit normal vector pointing outwards the boundary ∂ Ω of the brain domain Ω at location r ∈ ∂ Ω .Equation ( 5) specifies the initial condition of the problem whereas Equation ( 6) provides no-flux Neumann boundary conditions reflecting the inability of tumour cells to diffuse across ∂ Ω .

MR Data Acquisition
For the needs of this work, 6 healthy volunteers were enrolled for an MRI acquisition comprising a T1 BRAVO (echo time: ∼3 ms, repetition time: ∼8.To validate our approach (see Section 2.6), similar T1 BRAVO, T2 FLAIR, and DTI images as well as an additional T1Gd (echo time: 3.2 ms, repetition time: 8 ms, flip angle: 8°, pixel bandwidth: 255 Hz px −1 , slice thickness/spacing: 1/1 mm, matrix: 232 px×231 px) image acquired on a 3 T Achieva scanner (Philips Medical Systems, The Netherlands) of a 69-year-old male patient with IDH-wildtype GBM were retrospectively analysed.
The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Hospital-Faculty Ethics Committee of Hôpital Erasme (accreditation 021/406, protocol P2018/311, 3 May 2018).Informed consent was obtained from all subjects involved in the study.

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 [28].A water diffusion tensor field D water 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.

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 mm × 1 mm × 1 mm 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 maximisation of their mutual information [29].All resampling and registration steps were performed using an in-house software in C++ based on ITK [30] and VTK [31].

Skull Stripping
The brain volume was then segmented on the registered T2 FLAIR image using the Otsu thresholding [32] followed by a morphological erosion with structuring element radius 2 voxels (vx), a largest component extraction, a morphological dilation of radius 2 vx, a morphological closing of radius 8 vx, and a morphological dilation of radius 1 vx.

Brain Tissue Segmentation
The extracted brain domain was further segmented on the registered T1 BRAVO image into white matter, grey matter, and cerebrospinal fluid using an in-house implementation of the MICO intensity-based clustering algorithm comprising a bias field correction step [33].Manual corrections were further applied to the mis-segmented basal nuclei and falx cerebri.This last step is crucial to prevent the migration of tumour cells between brain hemispheres via other routes than the corpus callosum, as highlighted previously [8].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.

Tumour Segmentation
For validation purpose (see Section 2.6), the enhancing core and peritumour vasogenic oedema were also semiautomatically segmented on the T1Gd and T2 FLAIR images of the GBM patient using combinations of thresholding, connected component extraction, and morphological operations.

Tumour Diffusion Tensor
A tumour diffusion tensor field D(r) was built piecewise from the DTI-derived water diffusion tensor D water (r) and the segmented brain domain Ω as follows: where D white (r) and D grey are respectively the tumour cell diffusion tensor fields within the white and grey matter domains Ω white and Ω grey , and 0 is the null tensor, with: where d grey is the mean diffusivity of tumour cells in grey matter.
The white matter tumour cell tensor field D white (r) was built from the DTI-derived water diffusion tensor D water (r) using the method proposed by Jbabdi and colleagues in [12].This step is motivated by the observation that, for mechanical [34] and molecular [35] reasons, tumour cells preferentially migrate along rather than across brain fibres, similarly to diffusing water molecules.The method consists in increasing the degree of anisotropy of the water diffusion tensor and scaling its mean diffusivity MD = Tr (D water ) /3 to account for the difference in diffusive behaviour between tumour cells and water molecules [12]: where d white is the mean diffusivity of tumour cells in white matter, e i (r) is the i th eigenvector of D water (r), and λi (a) = l i (a)λ i , with: where λ i is the i th eigenvalue of D water (r), a ≥ 1 is a multiplicative factor introduced to increase the anisotropy of the tensor, and c l , c p , and c s are respectively the linear, planar, and spherical anisotropy measures of D water (r).
For reasons that will become clearer in Section 2.4, a unit tumour cell diffusion tensor was built at this stage by fixing d white to 1.0.A value of 0.1 was chosen for the d grey /d white ratio, as proposed previously in [11] to account for the restricted migration of tumour cells in grey compared to white matter.This ratio was supposed constant among all subjects as it is expected to depend exclusively on the structural organisation of healthy white versus grey matter and not on the tumour characteristics.Similarly, the anisotropy factor a was fixed to 10 for all subjects, as suggested in [12].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://lisaserver.ulb.ac.be/owncloud/index.php/s/KwEPG65gh1U7xNM.Further details on these data are available in Appendix B.

Dataset Synthesis
A synthetic tumour dataset was generated from the processed MR data described hereabove.200 fictitious tumours 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) to (6), totalling 1,200 synthetic tumours.For each tumour, the cell density distribution was sampled at four imaging time points t 1 , t 2 , t 3 , and t 4 .
Each tumour seed consisted of a 3 vx × 3 vx × 3 vx neighbourhood chosen among all segmented white matter voxels whose initial (i.e. at time t = t 0 ) normalised tumour cell density c was set to 1.For each simulated tumour, an infiltration depth λ = d white /ρ, a tumour front propagation velocity v = 2 √ d white ρ, and two imaging time intervals ∆t 1 and ∆t 2 were randomly chosen from uniform distributions (floating point for λ and v, integer for ∆t 1,2 ).The value ranges of the uniform distributions are provided in Table 1 and are of the same order of magnitude as in [17].
In this manner, a wide diversity of tumours can be uniformly sampled within a realistic range of infiltration depths and propagation velocities.Independently sampling d white and ρ values from uniform distributions may instead have resulted in too small (i.e. with too small values of d white and ρ) or too large (i.e. with too large values of d white and ρ) tumours.The empirical joint distribution of the sampled (λ, v) values as well as the corresponding joint distribution of (d white , ρ) and marginal distributions of d white and ρ are depicted in Figure 3.   2.
Table 2: Parameter values used for the tumour simulations in Figure 4.The model was solved by a forward Euler finite difference approach with a time step ∆t = 0.5 d using a GPU implementation of the 3D standard stencil referenced in [36] for the computation of the divergence term in Equation ( 4).The processing time was around 1 ms per iteration on a GeForce GTX 1080 GPU (NVIDIA, USA), leading to total simulation times in range 0.72-1.08s per tumour.All simulations were performed using a Python wrapping of our Tumour Growth Simulation ToolKit (TGSTK) written in C++/CUDA language based on VTK [31] and publicly available at https://github.com/cormarte/tgstk.The toolkit documentation can be found at https: //cormarte.github.io/tgstk/html.

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)[37] in Python on a GeForce RTX 3090 GPU (NVIDIA, USA).

Cell Density Estimation
The first problem addressed was to reconstruct a tumour cell density distribution from (i) two imaging contours-Γ 1 and Γ 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) tumour cell diffusion tensor field derived from DTI data as presented in Section 2.3.6.
This approach is motivated by the properties of the asymptotic travelling wave solution admitted by Equation ( 1) for constant coefficients and on the infinite cylinder, whose profile decreases exponentially with decay constant λ (see Section 1).For such a solution, the value of λ can be trivially estimated given the distance between two cell density iso-contours, and an approximate tumour cell density distribution can subsequently be reconstructed for sufficiently large x.Here, we assess the ability of deep neural networks to build an approximate solution of Equations ( 4) to (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 [18].
In a clinical setting, the value of the white matter diffusion rate d white used to scale the tumour cell diffusion tensor in Equation ( 9) is unknown-and its assessment will be addressed in the next section.Therefore this information is not considered for this problem.In contrast, the preferred migration directions of tumour 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 tumour cell density distribution.This motivates the introduction of the unit unscaled diffusion tensor field as an input of this problem, in addition to the Γ 1 and Γ 2 contours.
To address this problem, a 3D DCNN based on the U-Net architecture [38] was implemented, as it has been successfully applied to many medical imaging problems previously [25,39].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 × 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 × 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 × 2 and stride 2, followed by a bias-adding layer and a ReLU activation layer to expand the feature maps 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 × 1 and stride 1 followed by a bias-adding layer but no activation layer to merge the last 32 feature maps into a single tumour 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 × 192 × 128 × 8 (width × height × depth × channels).The first 2 channels are fed with the two binary regions respectively delimited by Γ 1 and Γ 2 .These regions were obtained by thresholding each generated tumour 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 [5] for the enhancing core and oedema outlines, respectively.The last 6 channels correspond to the 6 independent components of the unit (unscaled) tumour cell diffusion tensor (see Section 2.3.6).
To evaluate the generalisation 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 tumours generated from the MR data of the last patient are kept for evaluation purpose).The network was trained using the Adam optimiser [40] (learning rate: 10 −4 , β 1 : 0.9, β 2 : 0.999, : 10 −6 ) and the mean absolute error (MAE) loss over mini-batches of size 1.Data augmentation was performed by applying random shifts in range ±15 vx in the three spatial dimensions to each batch.Rotations were not applied as they would imply transformation of the tensor components, 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 −4 ) were kept, which occurred after 876 epochs.

Parameter Estimation
The second problem addressed is to estimate the value of the model parameters d white and ρ (or equivalently of the derived parameters λ = d white /ρ and v = 2 √ d white ρ) from (i) three imaging contours: two imaging contours-Γ 1 and Γ 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 Γ 3 obtained for the same c 2 threshold value at the earlier imaging time t 1 and (ii) the unit (unscaled) tumour cell diffusion tensor field.The time interval ∆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 travelling wave solution of Equation (1) (see Section 1), whose profile decay constant λ = d/ρ can be assessed from two cell density iso-contours at a given time point as mentioned hereabove.In addition, the propagation velocity of the tumour front v = 2 √ d ρ 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 λ and v can finally be used to assess the individual values of d and ρ.Here again, we assess the ability of deep neural networks to generalise 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 tumour diffusion tensor from DTI data, which can provide additional information for the estimation of the steepness and velocity of the tumour front.
The DCNN implemented for this task is a convolutional encoder.The network consists of 6 convolutional downsampling blocks and a fully connected output block.Each down-sampling block is made of 2 convolutional layers with kernel size 3 × 3 and stride 1, followed by a bias-adding layer and a ReLU activation layer.A convolutional layer with kernel size 2 × 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 × 3 × 2 × 8 output of the last down-sampling block and concatenates a 1 × 1 (width × channels) tensor to the flattened vector to feed the imaging time interval ∆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 λ and v.The network takes as input tensors of shape 192 × 192 × 128 × 9 (width × height × depth × channels).The first 3 channels are fed with the binary regions respectively delimited by Γ 1 , Γ 2 , and Γ 3 .These regions were respectively obtained by thresholding each generated tumour cell distribution at time t 2 with threshold values of 0.80 (Γ 1 ) and 0.16 (Γ 2 ) and the distribution at time t 1 with a value threshold value of 0.16 (Γ 3 ) [5].The last 6 channels correspond to the 6 independent components of the unit (unscaled) tumour cell diffusion tensor (see Section 2.3.6).To account for their different value range and scale, the target values of λ and v were standardised using the theoretical mean and variance of the respective uniform distributions from which they were sampled.
The same training/test splitting as for the tumour cell density estimation network was applied to the dataset.The network was trained using the Adam optimiser [40] (learning rate: 10 −4 , β 1 : 0.9, β 2 : 0.999, : 10 −6 ) and the mean squared error (MSE) loss.Data augmentation was performed by applying random shifts in range ±15 vx 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 −2 ) were kept, which occurred after 628 epochs.

Validation
To validate and illustrate our approach, we conducted the following numerical experiment: Starting from the tumour cell density distribution estimated at time t 2 from Γ 1 and Γ 2 using our first network (Figure 5), as well as the values of d white and ρ estimated from Γ 1 , Γ 2 , and Γ 3 using our second network (Figure 6) and Equations ( 12) and ( 13), we computed a tumour 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 tumour cell density distributions at times t 3 and t 4 -i.e.these obtained for the true cell density distribution at time t 2 as well as the true values of d white and ρ-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 Hausdorff and the average symmetric surface distance (ASSD) d ASSD between the imaging contours obtained from the true and estimated tumour cell density distributions for threshold values of c 1 and c 2 were also computed for each test tumour and time point, as given by: where d(x, y) is the Euclidian distance between elements x and y, and |X| is the cardinal of a set X.
It should be noted that minor post-processing was applied to the estimated tumour cell 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 as small (∼ 10 −5 ) but non-zero values were observed for some of these voxels in the predicted tumour cell distribution.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 −6 were also set to 0 since small local maxima (∼ 10 −5 ) were sporadically observed far from the tumour core, which gave rise to new tumour foci throughout the simulation.These post-processing steps allow 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.
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 Sections 2.2 and 2.3).To this extend, the segmented enhancing core and oedema regions (see Section 2.3.5) were provided to the first network along with the derived unit tumour cell diffusion tensor.

Results
The distribution of the mean absolute error computed over the test set between the true and estimated tumour cell distributions at time t 2 within the c > 0.01 contour is summarised by a boxplot in Figure 7 (1 st plot).The corresponding median value was 9.58 × 10 −3 .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 (1 st plots).An example of true and estimated tumour cell density distributions at time t 2 from the test set is depicted in Figure 9 (1 st 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 tumour cell distributions at time t 2 used in Figures 7 to 9 were provided by the first network (Figure 5).
The distributions of the relative error on the values of λ and v computed at time t 2 over the test set as well as on the values of d white and ρ derived with Equations ( 12) and ( 13) are summarised by boxplots in Figure 10.The corresponding median relative errors were 3.41 %, 3.30 %, 5.86 %, and 2.75 % for λ, v, d white , and ρ, respectively.The predicted versus estimated values of λ and v as well as of d white and ρ from the test set are plotted in Figure 11.The corresponding Lin's concordance correlation coefficients (CCC) [41] were 0.99, 0.95, 0.97, and 0.99 for λ, v, d white , and ρ, respectively.
As for imaging time t 2 , the distributions of the mean absolute error computed over the test set between the true and estimated tumour cell distributions at times t 3 and t 4 within the c > 0.01 contour are summarised by boxplots in Figure 7 (2 nd and 3 rd plot, respectively).The corresponding median values were 1.38 × 10 −2 and 2.20 × 10 −2 for t 3 and t 4 , 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 (2 nd and 3 rd plots).The true and estimated tumour cell density distributions at times t 3 and t 4 are depicted in Figure 9 (2 nd and 3 rd 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 tumour cell density distributions over simulated time is observed in Figures 7 to 9. The estimated tumour cell density distributions at times t 3 and t 4 used in Figures 7 to 9 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 parameters values provided by the second network (Figure 6).
Finally, the estimated tumour cell 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 oedema contours, respectively.superimposed to the T1 and T2 FLAIR image are depicted in the 4 th and 5 th rows, respectively.The blue, red, and green segments respectively correspond to the target, prediction, and overlapping contour voxels.The predicted cell distribution at time t 2 was provided by the first network (Figure 5).The estimated cell distribution at times t 3 and t 4 were computed using the reaction-diffusion model from the cell distribution predicted at time t 2 and the predicted model parameters values provided by the second network (Figure 6).12) and (13).Horizontal line: median, box: interquartile range, whiskers: ±1.5 interquartile range, asterisks: outliers.12) and (13).For each plot, the identity function is superimposed in red and the corresponding Lin's concordance correlation coefficient CCC is provided.

Discussion
Reaction-diffusion models have been studied for decades to capture the growth of gliomas but severe limitations regarding the estimation of their initial conditions and parameter values have restrained their use as a proper personalised 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 to (i) derive a unit tensor field from clinical DTI data as described herein, accounting for the preferential migration of tumour cell along white matter tracts and (ii) extract three contours obtained through a cell density threshold-like imaging process described by Equation ( 2) for two different threshold values and time points.
Regarding the second requirement, the outlines of the peritumour vasogenic oedema and enhancing core have been proposed previously [5], respectively visible on T2 FLAIR and T1Gd MR images acquired in routine for glioma follow-up.Nevertheless, it is worth noticing that peritumour vasogenic oedema does not strictly speaking correspond to a region of tumour cell invasion but results from an accumulation of extracellular fluid originating from tumourinduced alterations of the blood-brain barrier [42,43] and changes in hydrodynamic pressure [44].Consequently, the T2 FLAIR imaging process might not be accurately described by Equation (2), as also supported by our previous experimental results in [18].Furthermore, anti-angiogenic drugs are known to dramatically reduce vasogenic oedema without stopping tumour progression [42].Therefore, other MR sequences or modalities could be better suited for the initialisation and parameter estimation of reaction-diffusion glioma growth models.For instance, ADC maps derived from DW-MRI data could more accurately reflect tumour cell invasion as proposed in [16,45].PET imaging with radio-labelled amino acids could also provide additional information to this extent as suggested in [24,46].
Once the aforementioned prerequisites are met, our approach makes it possible to (i) extrapolate a whole brain tumour cell density distribution between and beyond the visible outlines of the tumour that is compatible with the reactiondiffusion model in Equations ( 4) to ( 6) and (ii) individually assess the values of the diffusivity and proliferation parameters of the model.Extrapolating tumour invasion is of utmost interest for radiotherapy planning since it would allow to define personalised margins which more accurately target the tumour while avoiding irradiation of the healthy tissues, as previously discussed in [7,8].The independent assessment of the diffusivity and proliferation parameters of the model is for its part of great interest to better characterise the tumour [20].The combination of both gives access to a fully personalised tool, initialised from clinical imaging data and allowing to anticipate the spatial-temporal growth of gliomas.Such a tool could for example be of considerable interest for dose fractionation optimisation in radiotherapy using a reinforcement learning approach as used in [47].Furthermore, as it only depends on post-processed data (binary segmentations and a DTI-derived water diffusion tensor) and not on raw MR data, the proposed approach may be robustly extended to other scanners and centres.Besides, the method is also robust to variations in the imaging timing as variable imaging time intervals were considered for synthetic dataset generation, making it well adapted to the clinical reality.
The proposed method provided accurate estimations of the tumour cell distribution from only two imaging contours at a single time point, with a mean absolute error below 10 −2 within the c > 0.01 contour for the same tumours, as evaluated on 200 synthetic tumours grown over the brain domain of an unseen test subject.Our method also provided accurate estimates of the individual diffusivity and proliferation rate parameters of the model from three imaging contours extracted from two time points, with median relative errors of 5.86 % and 2.75 %, respectively (see Figure 10), and strong concordance (CCC ≥ 0.95) with the true parameter values (see Figure 11).Furthermore, we showed that the spatio-temporal evolution of the tumour cell density distribution at later time points (90 d and 180 d later) can be accurately captured from the estimated initial distribution 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 mm × 1 mm × 1 mm) in most cases (see Figure 8).Nevertheless, a loss of accuracy in the estimated tumour cell density over simulated time was observed (see Figures 7  to 9), imputed to the accumulation of errors originating from uncertainties on the estimated model parameters values and initial condition.In particular, artefactual local maxima in the tumour cell densities distributions predicted by the CNN were found to give rise to new tumour 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 8(a) and (b)).We also demonstrate the applicability of our proposed method to actual MR data of a GBM patient, for which we were able to reconstruct a tumour 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.
As a future work, tumour-induced mass effect should be further integrated into the reaction-diffusion model since it is known to cause substantial deformations of the brain parenchyma as the tumour grows, which should also be taken into account for accurate treatment planning.Such effects have been previously considered [3,4] but would introduce additional parameters to be assessed.Besides, 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 as the one described herein.Necrosis should also be taken into account by the model as proposed in [13,14], which would have avoided the counter-intuitive correspondence between the hyper-dense (c ∼ 1) region of the estimated tumour 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 [48].
Ultimately, the threshold values could be fed to the networks along with the binary contours, but 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 validation of our approach still relied on healthy subject data.Therefore, the underlying hypothesis was made that the reaction-diffusion model defined by Equations ( 4) to (6) and used for tumour 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 have also been proposed previously [1,2,[15][16][17] but again introduce additional parameters, increasing the complexity of the problem.Finally, the assumption was made throughout this work that the diffusivity and proliferation rate parameters of the model can be considered constant over time in untreated patients.The imaging time interval over which this assumption can reasonably hold should also be determined based on real glioma patient data.
This work further demonstrates the ability of DCNNs to accurately approximate functions over complex domains and provides encouraging results towards the full personalisation of reaction-diffusion glioma growth models, which has remained an unsolved problem for decades.

Conclusion
We proposed a deep learning-based approach to address together the problems of estimating the initial condition 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 tumours 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 towards the full personalisation of glioma reaction-diffusion models may open up tremendous possibilities in the field.

Figure 1 :
Figure 1: Example of processed MR data.(a) Axial slice of the T2 FLAIR image with superimposed segmented brain mask (red).(b) Corresponding slice of the T1 BRAVO image.(c) Segmented brain map obtained with the MICO algorithm followed by manual corrections.

Figure 2 :
Figure 2: Example of processed DTI data.(a) DTI-derived water diffusion tensor field after susceptibility-induced distortion, eddy currents, and patient motion correction using FSL.(b) Tumour diffusion tensor field with increased anisotropy in white matter (a = 10) and scaled diffusivity (d grey /d white = 0.1) built from the water diffusion tensor field in panel (a) and the brain map in Figure 1(c) as described in Section 2.3.6.The subpanel located at row i and column j of panels (a) and (b) corresponds to the tensor component d i,j .

dFigure 3 :
Figure 3: Sampling of the model parameters.(a) Empirical joint distribution of the (λ, v) values sampled from uniform distributions (blue marks) with superimposed sampling domain boundary (red segments).(b) Corresponding joint distribution of the derived (d white , ρ) values using Equations (12) and (13) (blue marks) with superimposed sampling domain boundary (red curves).(c) Empirical marginal distribution of the derived d white values (blue bars) with superimposed theoretical distribution (red curves).(d) Empirical marginal distribution of the derived ρ values (blue bars) with superimposed theoretical distribution (red curves).

ForFigure 4 :
Figure 4: Examples of simulated tumour cell density distributions at times t 1−4 (1 st to 4 th columns, axial slices) from the MR data of the same subject as in Figures 1 and 2. The corresponding model parameter values are provided in Table2.

2 Figure 5 :
Figure 5: 3D U-Net architecture[38] with its feature map dimensions used for cell density estimation.The network takes as input volumes of dimensions 192 × 192 × 128 with 8 channels corresponding to the 2 contours Γ 1 and Γ 2 and the 6 independent components of the unit (unscaled) tumour cell diffusion tensor field and outputs a cell density map with the same spatial dimensions.

2 Figure 6 :
Figure6: 3D convolutional regressor architecture with its feature map dimensions used for parameter estimation.The network takes as input volumes of dimensions 192 × 192 × 128 with 9 channels corresponding to the 3 contours Γ 1 , Γ 2 , and Γ 3 and the 6 independent components of the unit (unscaled) tumour cell diffusion tensor field as well as the time interval δt 2 between Γ 3 and Γ 2 , and outputs estimated values of λ and v.

Figure 10 :
Figure 10: Boxplots of the relative error on the predicted model parameter values evaluated on the test set.(a) Relative errors on the estimated values of λ and v provided by the second network (Figure 6).(b) Corresponding relative errors on the derived values of d white and ρ using Equations (12) and(13).Horizontal line: median, box: interquartile range, whiskers: ±1.5 interquartile range, asterisks: outliers.

d 2 yr − 1 ]Figure 11 :
Figure 11: Scatterplots of the true versus predicted values of the model parameters from the test set.(a)-(b) True versus predicted values of λ and v provided by the second network (Figure 6).(c)-(d) True versus estimated values of d white and ρ derived from the predicted values of λ and v using Equations (12) and(13).For each plot, the identity function is superimposed in red and the corresponding Lin's concordance correlation coefficient CCC is provided.

Figure 12 :
Figure 12: T1Gd image (1 st row), T2 FLAIR image (2 nd row), and estimated tumour cell density distribution using the first network (3 rd row) for an IDH-wildtype glioblastoma patient in axial (1 st column), sagittal (2 nd column), and coronal (3 rd column) planes.The contours of the segmented enhancing core and peritumour vasogenic oedema are superimposed in red on the T1Gd (1 st row) and T2 FLAIR (2 nd row) images, respectively.

FundingC
.M. is funded by the Walloon Region, Belgium (PROTHER-WAL).C.D. is a senior research associate at F.R.S.-FNRS.The Department of Nuclear Medicine at Hôpital Erasme is supported by Association Vinçotte Nuclear (AVN), Fonds Erasme, and the Walloon Region (BioWin).

Table 1 :
Value ranges and units of the uniform distributions used to sample the tumour growth model parameters for the generation of the synthetic tumour dataset.Two additional time intervals, ∆t 3 and ∆t 4 , were fixed to 90 d for validation purpose (see Section 2.6).Starting from tumour emergence time t 0 , the time intervals ∆t 1−4 univocally define the four imaging time points t i = t 0 + i j=1 ∆t j , i = 1, . .., 4.For each sampled couple of (λ, v) values, a white matter diffusion rate value d white and a proliferation rate value ρ were derived as: Boxplots of the mean absolute error distribution within the c > 0.01 contour computed voxelwise over the whole test set for times t 2 = ∆t 1 + ∆t 2 ∈ [180, 360] d (see Table 1), t 3 = t 2 + 90 d, and t 4 = t 2 + 180 d.Horizontal line: median, box: interquartile range, whiskers: ±1.5 interquartile range, asterisks: outliers.