Initial Condition Assessment for Reaction-Diffusion Glioma Growth Models: A Translational MRI-Histology (In)Validation Study

Reaction-diffusion models have been proposed for decades to capture the growth of gliomas. Nevertheless, these models require an initial condition: the tumor cell density distribution over the whole brain at diagnosis time. Several works have proposed to relate this distribution to abnormalities visible on magnetic resonance imaging (MRI). In this work, we verify these hypotheses by stereotactic histological analysis of a non-operated brain with glioblastoma using a 3D-printed slicer. Cell density maps are computed from histological slides using a deep learning approach. The density maps are then registered to a postmortem MR image and related to an MR-derived geodesic distance map to the tumor core. The relation between the edema outlines visible on T2-FLAIR MRI and the distance to the core is also investigated. Our results suggest that (i) the previously proposed exponential decrease of the tumor cell density with the distance to the core is reasonable but (ii) the edema outlines would not correspond to a cell density iso-contour and (iii) the suggested tumor cell density at these outlines is likely overestimated. These findings highlight the limitations of conventional MRI to derive glioma cell density maps and the need for other initialization methods for reaction-diffusion models to be used in clinical practice.


Introduction
Gliomas are the most common primary brain tumors.Diffuse gliomas, which include its most aggressive form glioblastoma (GBM), are known to be highly infiltrative [1], with the presence of tumor cells reported as far as 4 cm from the gross tumor [2].The early diagnosis and follow-up of gliomas usually rely on magnetic resonance imaging (MRI).However, whereas recent advents in MR technologies have given access to a significantly deeper insight into the tumor biology, none of the routinely acquired MR sequences allows to directly assess the extent of the tumor cell invasion.Instead, tumor-induced alterations of the micro-environment are seen on MR images such as peritumor vasogenic edema visible on T2/T2 FLAIR sequences and the enhancing tumor core visible on T1-weighted sequences with injection of gadolinium-based contrast agent (T1Gd).Peritumor vasogenic edema originates from an increase in the blood-brain-barrier (BBB) permeability induced by the release of vascular endothelial growth factor (VEGF) by tissues under hypoxic stress [3,4] combined with changes in the brain hydrodynamic pressure [5].The formation of an enhancing tumor core results from a breakdown of the BBB subsequent to neo-vascularisation induced by VEGF, allowing gadolinium-based contrast agents to diffuse into brain tissues [3].
The pathological and molecular examination carried out on resected or biopsied tissue samples remains the gold standard method to confirm the diagnosis and determine the histological type and grade.According to the 2016 WHO classification of the central nervous system tumors, the identification of infiltrative patterns is essential in the differential diagnosis of diffuse versus pilocytic astrocytomas, a more-circumscribed neoplasm.An overall assessment of the invasion extent would also be beneficial for surgery and radiotherapy planning.However, due to their high invasiveness and long processing times, the number and frequency of biopsy procedures are restricted, limiting the use of pathological examination as a proper tumor invasion assessment tool and highlighting the complementarity of radiological examination [6].In this matter, quantitatively linking glioma cell invasion patterns observed histologically to MR-visible abnormalities would be of great interest.Such information would allow to non-invasively assess glioma extent and progression, within and beyond the abnormality outlines, while providing a better interpretation of the observed abnormalities.
Mathematical glioma growth modeling has addressed the problem of estimating glioma cell distribution within brain tissues and predicting its temporal evolution.Among the investigated models, reaction-diffusion models first introduced by Murray and colleagues in the early 1990's [7] are probably the most widely used, with potential applications for patient follow-up and improved radiotherapy planning [8].These models rely on a reaction-diffusion equation to capture the spatial-temporal evolution of a tumor cell density function: (2) (3) where c(r, t) is the tumor cell density at position r and time t normalized by the maximum carrying capacity c max of brain tissues (c(r, t) ∈ [0, 1], ∀r, t), D(r) is the symmetric tumor cell diffusion tensor at position r, ρ is the tumor cell proliferation rate, c 0 (r) is the initial tumor cell density at position r, and n ∂Ω (r) is a unit normal vector pointing outwards the boundary ∂ Ω of the brain domain Ω at position r ∈ ∂ Ω .Equation (2) specifies the initial condition of the problem.Equation (3) provides no-flux Neumann boundary conditions reflecting the inability of tumor cells to diffuse across ∂ Ω .
Reaction-diffusion models as the one in Equations ( 1) to (3) are particularly attractive for clinical applications since they only have a few parameters that could be assessed from patient imaging data.For instance, based on prior observations that tumor cells rather migrate along white matter tracts, Jbabdi and colleagues proposed to derive the tumor cell diffusion tensor D(r) from diffusion tensor imaging (DTI) data [9].For a more detailed overview of reaction-diffusion glioma growth modeling and its potential clinical applications, the reader is referred to [7,8,[10][11][12].
One problem arising when attempting to solve Equations ( 1) to (3) from actual imaging data of newly diagnosed glioma patients is to estimate the initial cell density c 0 (r) at every location of the brain domain Ω.Early works on glioma growth modeling proposed to relate MR-visible abnormalities observed at time t to the tumor cell density function c(r, t).In [11], the authors suggested to model the MR imaging process as a simple cell density threshold function, that is: where I T1Gd (r, t) and I T2 (r, t) are respectively the imaging functions of the T1Gd and T2/T2 FLAIR MR sequences indicating whether the abnormality is visible at location r and time t on the sequence, and c enhancing and c edema are the corresponding tumor cell density detection thresholds.Based on these assumptions, the authors suggested that the outlines of the tumor enhancing core in T1Gd images and of the vasogenic edema in T2/T2 FLAIR images would correspond to iso-contours of the tumor cell density function: c(r, t) = c enhancing for r ∈ ∂Ω enhancing c edema for r ∈ ∂Ω edema (6) where ∂ enhancing and ∂ edema are respectively the enhancing core and edema outlines.The authors also suggested hypothetical values for c enhancing and c edema of 0.80 and 0.16 respectively [11], although no rationale was provided for these values.Building upon this work, Konukoglu and colleagues proposed a fast-marching approach to construct an approximate solution of Equations ( 1) to (3) at imaging time satisfying Equation ( 6) [12].More recently, the same group suggested that, for a spatially constant and isotropic diffusion coefficient d white and from a certain distance to the tumor core, the tumor cell density in white matter would approximately decrease exponentially with the distance d to the core [8]: λ white (7) where λ white is the infiltration length of tumor cells in white matter given by d white /ρ.Provided two iso-cell density contours and a distance map to the tumor core, the value of λ white can theoretically be assessed using Equation (7).A similar reasoning had been previously applied in [7] for glioma growth modeling in computed tomography (CT) images, leading to the same expression as Equation (7) for the initial condition c 0 (r).
Nevertheless, these attempts to derive a tumor cell density distribution from MR images rely on the existence of cell density iso-contours in MR images and are based on unverified assumptions in Equations ( 4) to (7).However, as will be further discussed, the extent of vasogenic edema is known to be impacted by the administration of corticosteroids and anti-angiogenic treatments [3].Furthermore, as previously observed by our group in [13], the reaction-diffusion model is highly sensitive to the provided initial tumor cell density distribution c 0 (r).The validation of the aforementioned hypotheses is thus crucial for the model to be usable in clinical routine.In this work, we propose to verify these assumptions, as well as the value of 0.16 for c edema suggested in [11], through a translational MRI/histology study conducted in a case of non-operated GBM.To this end, stereotactic histological analyses are performed using a 3D-printed slicer designed from antemortem MRI data.Cell density maps are computed automatically from the scanned histological slides using a deep convolutional neural network and related to an MR-derived geodesic distance map to the tumor core.The relation between the edema outlines and the geodesic distance to the core is also investigated.Our results highlight the limitations of using routine MRI to derive glioma cell density maps and point out the need for other validated initialization methods to make reaction-diffusion growth models usable in clinical practice.

Clinical case
For the needs of this work, the case of a deceased 89-year-old female patient with GBM was studied retrospectively.
The patient underwent an MRI examination in November 2017 in the context of a clinical frontal syndrome, which revealed a massive right-frontal expanding lesion with intense heterogeneous enhancement post-injection of gadolinium contrast agent and a necrotic core, surrounded by a large area of perilesional edema.Considering the patient's age, a consensual decision was taken not to perform surgery and a diagnosis of GBM was made exclusively based on MRI.
A corticotherapy (methylprednisolone) and a palliative chemotherapy-based treatment (temozolomide) were initiated right after diagnosis.In December 2017, after one cycle of temozolomide, the patient died of a septic shock caused by bowel perforation.An autopsy was carried out within 12 hours after death as part of which the patient's brain was collected and fixed in formalin for 24 months.Upon microscopic examination, an infiltrating glioma with high cellularity, astrocytic phenotype cells and nuclear atypia was observed along with glomeruloid vascular proliferation and areas of pseudo-palissading necrosis.All procedures performed in this study were approved by the Hospital-Faculty Ethics Committee of Hôpital Erasme (accreditation 021/406) under reference P2019/245 and are in accordance with the 1964 Helsinki declaration and its later amendments or comparable ethical standards.

MR image acquisitions
The T1 (TR = 8 ms, TE = 2.9 ms, TI = 950 ms, FA = 8°) and T2 FLAIR (TR = 4800 ms, TE = 320 ms, TI = 1650 ms) MR images routinely acquired at diagnosis time on a 3T Achieva scanner (Philips Healthcare, The Netherlands)referred to as the 'in vivo' images hereafter -were retrospectively used in this work for registration guidance and delineation of the vasogenic edema.
Additionally, a T1 BRAVO 'ex vivo' acquisition (TR = 8.264 ms, TE = 3.164 ms, TI = 450 ms, FA = 12°) of the brain placed inside the 3D-printed slicer (see below) was performed on a 3T Signa PET/MR scanner (GE Healthcare, USA) right before slicing.It should however be noted that brain fixation has caused convergence of the white and gray matter T1 values, as reported in [14,15].Consequently, the acquired 'ex vivo' T1 image rather has a proton-density (PD) contrast.Also note that drainage of the extracellular fluid made it impossible to delineate edema regions on postmortem T2 FLAIR images, motivating the use of a registered antemortem T2 FLAIR image for edema delineation hereunder.

Slicer design and tissue sampling
To relate histological observations to the abnormalities in MR images and to the MR-derived distance map to the tumor core (see below), a brain slicer was designed based on the in vivo T2 FLAIR image and 3D-printed.Such a slicer allows the brain to be re-positioned in antemortem imaging orientation and facilitates the cutting of sagittal brain slices.The slicer design procedure is illustrated in Figure A1.A similar slicer design approach was previously adopted in [16].
Ten guides were also designed to ease the collection of sample blocks from brain slices that are compatible with our histological processing chain.These consist of plates with grooves from which the brain slice volume was subtracted.
The brain slicing and samples collection procedure is illustrated in Figure 1.More details on the design steps of the slicer and the cutting guides are available in Appendix A.

Sample processing and analysis
28 tissue samples were selected within the tumor core, as well as within and beyond the peritumoral edematous region based on the in vivo T2 FLAIR image.The samples were formalin-fixed and paraffin-embedded (FFPE, ISO 15189).5 µm slides were cut from each sample and stained with hematoxylin and eosin (HE).The 28 stained slides were scanned in 20× mode (0.46 µm/px) on a calibrated NanoZoomer 2.0-HT digital slide scanner (Hamamatsu Photonics, Japan) for numerical processing.The stained slides were independently examined by an experienced pathologist blinded to MRI for the presence of pseudo-palisading necrosis, tumor cells (in block or infiltrating), glomeruloid vascular proliferation and edema.As will be further discussed, immunohistochemistry (IHC) staining was also investigated but did not provide satisfactory results due to over-fixation of the brain tissues.

Cell density maps
Cell density maps were computed from the scanned HE slides to highlight tumor cell invasion in normal brain tissues.Each scanned slide was first resampled to an isotropic pixel size of 1 µm × 1 µm and divided into adjacent tiles of 100 px × 100 px.Cell nuclei within each tile were automatically counted using a weakly-supervised deep learning approach detailed in Appendix B. The counting result was divided by the actual tissue area within the tile, defined as the number of tissue pixels (see Appendix B) times the pixel area (10 −6 mm 2 ).The computed cell density was finally stored as a 2D image with a pixel size of 0.1 mm × 0.1 mm where each pixel exactly corresponds to one 100 px × 100 px tile of the resampled slide.The cell density map computation procedure is illustrated in Figure 2.
In addition, volume cell density values were extrapolated from the computed surface cell density values, since the former are the actual values of interest for the reaction-diffusion growth model.Under the assumptions that: 1. Cell nuclei are approximately spherical, 2. Cell nuclei distribution is locally isotropic, 3. Tile dimensions are sufficiently large to contain multiple cells, 4. Tile dimensions are sufficiently small for the cell density to be considered homogeneous and isotropic, and denominating l the side length of the square tile, a volume cell density value c volume can be extrapolated for a cube with the same side length l from the surface cell density c surface by:   A1).(d) Whole computed cell density map.

Cell density maps to ex vivo T1 registration
Substantial deformations of the brain occurred between the in vivo MR acquisitions and the histological analyses.The ex vivo T1 image space was thus used as the reference space for the analyses and the cell density maps were registered to the ex vivo T1 image as follows.The ex vivo T1 image was first resampled by linear interpolation to an isotropic voxel size of 0.5 mm.The 2D cell density maps were artificially extended to 3D by addition of a thickness of 0.5 mm and resliced to sagittal orientation.The density maps were finally rigidly registered to the corresponding slice of the resampled ex vivo T1 image based on user-defined landmark pairs using an in-house software in C++ based on VTK [17] and ITK [18].
The cell density map registration process is greatly facilitated by the use of a 3D-printed slicer since it allows to impose the brain slicing orientation.The complex histology slide to MR image registration process in 3D is thus reduced to a simple MR slice selection followed by the identification of at least 3 landmark pairs in-plane.Furthermore, the computed cell density maps have the great advantage of providing spatial tissue information at an intermediate scale between histological and radiological images, with a contrast similar to T1-weighted MR images.The cell densities of white and gray matter are indeed substantially different, as is their T1 and PD values, which eased the identification of landmarks pairs.

Edema delineation
To verify the assumptions in Equations ( 5) and ( 6), the edema region has to be delineated in the reference ex vivo T1 image space.However, as previously mentioned, the drainage of the extracellular fluid made it impossible to discern vasogenic edema on the ex vivo MR images.The in vivo T2 FLAIR image was thus registered to the reference ex vivo T1 image.To this end, the in vivo T1 image acquired on the same day was first registered on the ex vivo T1 image using rigid followed by B-spline transforms using the Elastix software [19].The computed transforms were then successively applied to the in vivo T2 FLAIR image.The Elastix parameter files used for registration are available in Appendix C. The edema segmentation was finally performed semi-automatically on the registered T2 FLAIR image using a combination of thresholding and morphological operations.

Distance map
To verify assumption in Equation ( 7), a 3D geodesic distance map to the tumor core across white matter was computed from the ex vivo T1 image.White matter was first segmented using an in-house gradient-based anisotropic diffusion algorithm followed by manual corrections to ensure that no physically incompatible bypass exists between white matter regions.The tumor core was then segmented on the same image using a combination of thresholding and morphological operations.A distance map to the tumor core across the segmented white matter region was finally computed using an adapted implementation of the anisotropic fast marching (AFM) algorithm presented in [20].Note that since no DTI images were available for the patient, the anisotropy of glioma cell diffusion mentioned hereabove could not be taken into account in this work.A unit isotropic metric tensor field was thus provided to the AFM algorithm, hence the abusive use of the term 'distance map' to designate the 'traveling time map' returned by the algorithm.
The relation between the edema extent and the geodesic distance map was also investigated using the Hausdorff distance and the average symmetric surface distance (ASSD), computed between the edema outlines and the contours of the binary region obtained by thresholding the distance map.The Hausdorff distance d Hausdorff and the ASSD d ASSD between two sets A and B are respectively given by [21]: where d(x, y) is the Euclidian distance between elements x and y, and |X| is the cardinal of set X.

Cell density model
As mentioned, the over-fixation of brain tissues prevented any IHC staining.Therefore, HE staining had to be used instead, making it difficult to distinguish infiltrating tumor cells from healthy brain cells on the scanned slides.Consequently, nuclei-based cell density maps were computed which reflect the total cell density, with no distinction between tumor and healthy cells.To address this problem, we propose to verify the following equation instead of Equation ( 7): where c total , c tumor , and c white are respectively the total, tumor, and healthy cell density in white matter, and c core is the tumor cell density iso-value along the tumor core boundary.In this formulation, the cellularity of healthy white matter is supposed to be approximately constant and the invading tumor cells are assumed to be superimposed to the white matter baseline cellularity c white .
To relate the total cell density and the geodesic distance to the tumor core, the registered cell density maps were resampled to the same voxel size as the geodesic distance map (0.5 mm × 0.5 mm × 0.5 mm) and all available pairs of density/distance values among the segmented white matter voxels were extracted.c core , λ white and c white in Equation ( 12) were finally least-squares fitted to the available experimental density/distance pairs using SciPy's 'optimize' module in Python [22].
3 Results     A1) (1 st and 3 rd columns) and corresponding slices of the geodesic distance map to the tumor core (2 nd and 4 th columns) superimposed to the ex vivo T1 image.
The available pairs of density/distance values among all white matter voxels are plotted in Figure 5 for the surface density data (Figure 5(a)) and the volume density data extrapolated using Equation (8) (Figure 5(b)).The fitted model curve given by Equation ( 12) is superimposed in red for each plot and the corresponding parameter values are respectively provided in Table 1 and Table 2.
Table 1: Least-squares fitted values of the cell density model parameters in Equation ( 12) for the surface cell density data plotted in Figure 5(a).(blue).The distribution is rather continuous, suggesting that the edema boundary does not correspond to an iso-distance contour in contrast to the expected step-like distribution superimposed in red in Figure 6.Five percent of the edema boundary voxels are located at a distance greater or equal to 49.4 mm and only 1 % are located at a distance greater than 60.1 mm from the tumor core.The threshold geodesic distance values that provided the smallest Hausdorff distance (24.65 mm) and ASSD (1.97 mm) between the edema outlines and the contour of the corresponding thresholded region in the distance map were 43.5 mm and 35.5 mm, respectively.The corresponding volumes are depicted in Figure 7.
The results of the blinded pathological examination and numerical tile processing are summarized in Table A1, reporting the presence of pseudo-palisading necrosis, tumor cells (tumor block versus infiltrative cells) and glomeruloid vascular proliferation, along with the minimum, maximum and mean cell density and distance values within the corresponding slide.The furthest distance at which suspected infiltrating tumor cells were identified was around 46 mm (slide 7), whereas edema was detected as far as 61.7 mm on the same slide.

Discussion
The exponentially decreasing glioma cell density profile with distance to the tumor core suggested in [7] and [8] is compatible with our experimental data, as observed in Figure 5 for both surface (a) and extrapolated volume (b) cell densities.The fitted value of λ white (8.46 mm) for the volume cell density data is in the same order of magnitude as the one used in [8] for simulation (4.2 mm).The baseline volume cell density value of 0.59 cell mm −3 in white matter is also in accordance with the literature, although large variations are to be noted between the reported values [23].In contrast, the assumption of iso-cell density edema outlines in Equation ( 6) was not verified by our experimental data.Indeed, considering a monotonically (exponentially) decreasing relation between the tumor cell density and the distance to the tumor core, iso-density contours should coincide with iso-distance contours.However, our results suggest that edema outlines do not coincide with an iso-distance contour (see Figure 6) and would therefore not correspond to a cell density iso-contour either.This apparent incompatibility of assumptions in Equation ( 5) and Equation ( 6) can be explained by the thresholding behavior of the imaging function suggested in Equation ( 5) from which assumption in Equation ( 6) was deduced in [11].In fact, thresholding a spatial function may give rise to iso-value contours only if the function is sufficiently smooth and continuous.In contrast, the tumor cell density function discretized over the voxel grid has discontinuities at interfaces between white and gray matter and along the brain domain boundary.Indeed, the difference in tumor cell diffusivity between white and gray matter [8,10,12] gives rise to steep tumor cell gradients at the white/gray matter interfaces, resulting in substantial discontinuities of the cell density function at the voxel precision.
Along the brain domain boundary, discontinuities are even more pronounced since no tumor cell are allowed to diffuse outside the brain domain (see Equation ( 3)), resulting in an accumulation of cells along boundaries.Consequently, edema outlines may not correspond to cell density iso-contours even if the thresholding behavior of the edema imaging function in Equation ( 5) turned out to be verified.As an illustration, the edema outlines in Figure 3(b) were split into blue and red segments, respectively corresponding to parts of the edema boundary where tumor cell diffusion is free (blue) and where diffusion is restricted due to a local decrease in tumor cell diffusivity (red).From Figure 6, it can be reasonably assumed -taking a margin of 1 % for registration errors -that the edema extended up to 60.1 mm from the tumor core, which would correspond to the actual 'free' edema extent.Finally, it should be noted that due to the restricted number of available cell density map voxels also belonging to the free-to-diffuse part of the edema outlines, assumption in Equation ( 6) could not be directly verified, which motivates the indirect distance-based reasoning herein.
For reasons mentioned hereabove, the thresholding behavior of the edema imaging function in Equation ( 5) may still be valid even after invalidation of Equation ( 6).Nevertheless, the thresholded distance map regions whose contour respectively minimizes the Hausdorff distance and the ASSD to the edema outlines were not found to accurately coincide with the edema region, as depicted in Figure 7.Both thresholded distance map regions (blue) in Figure 7 were indeed found to extend further in the contralateral hemisphere via the corpus callosum and less far towards the right posterior region, compared to the edema region (red).Still under the hypothesis of a monotonically decreasing relation between the tumor cell density and the distance to the tumor core, the threshold-like imaging function in Equation ( 6) is therefore not compatible with our results.It should however be noted that this apparent inadequacy could result from the use of an isotropic metric tensor for the geodesic map computation instead of a DTI-derived anisotropic metric tensor, which would have allowed to account for the preferential migration of glioma cells along white matter tracts [24].
The iso-density value of 16 % of the maximum cell carrying capacity suggested for the edema contours in [11] was not supported by our experimental data either.Indeed, assuming that Equation ( 6) would still be verified on the free-to-diffuse part of the edema boundary, the proposed cell density model in Equation (12) with the parameter values fitted to volume cell density data in Table 2 suggests an over-cellularity of only 8.6 × 10 1 cell mm −3 with regard to the white matter cellularity baseline c white at a distance of 60.1 mm to the tumor core corresponding to the free edema extent.This over-cellularity corresponds to only 0.08 % of c core ≤ c max , which strongly invalidates the commonly accepted value of 16 %.In addition, since the c core value was likely underestimated as mentioned hereabove, an even lower actual value of c edema is to be expected.These results are confirmed by blinded pathological examination, which did not reveal any noticeable invasion of the brain parenchyma -even within the edematous region -beyond a distance of 46 mm to the tumor core.This overall analysis does however not exclude the possible presence of isolated infiltrating tumor cells at the edema boundary and beyond as observed in [2,25,26].
Although MRI provides high contrast in soft tissues and is currently the standard of care for radiological examination of gliomas, abnormalities visible on conventional MR sequences are not trivially related to the tumor cell invasion extent.A striking illustration of this limitation is the administration of corticosteroid or anti-angiogenic therapies to reduce edema-related symptoms in glioma patients, which does however not stop tumor progression.Consequently, a decoupling arises between tumor progression and its visible effects on the surrounding environment, potentially leading to a misclassification of the disease as responding.The impact of such therapies on the MR-based follow-up of gliomas has been extensively studied through numerical simulations in [3].In this work, we invalidated two commonly made assumptions relating the outlines of visible abnormalities on MRI to the tumor cell density function: (i) Assuming a threshold-like edema imaging function (see Equation ( 5)), the edema contour may not correspond to an iso-contour of the cell density function as soon as the migration of tumor cells is locally restricted or prevented and (ii) at a distance corresponding to the maximum extension of the vasogenic edema, the over-cellularity was found to be negligible in our studied case, as opposed to the previously hypothesized value of 16 %.These results raise the question of the applicability of both previously proposed methods to assess the glioma cell density distribution from routine MR images presented in the introduction.Indeed, whereas the method proposed in [12] allows to compute an accurate approximate solution of Equation ( 1), it still relies on the assumption that iso-density contours can be derived from MR data.In the case of a more simple exponentially decreasing model as in Equation ( 7) [7,8], iso-density contours would still be required to assess the infiltration length parameter λ white .Since a high sensibility of the reaction-diffusion tumor growth models to their initial condition was previously reported by our group [13], deriving an initial spatial cell density distribution that is as reliable as possible is crucial for the model to be applied in clinical practice.To this end, the use of other imaging sequences or modalities such as average diffusion coefficient (ADC) maps [27] or amino-acid positron emission tomography should be investigated for the in vivo assessment of the tumor cell density distribution.
This study was however prone to several limitations.First, due to the scarcity of the human body material analyzeda non-operated brain with GBM -this study was based on a single case and should be further conducted on a larger diffuse glioma cohort of various grades.In addition, the use of murine glioma models for conducting such studies at a larger scale would be of interest but is restricted due to reported substantial differences in cortical [28], glial [29], and endothelial [30] cells between human and mouse brain, leading to only partial capture of human glioma features by such models [31].Second, IHC staining could not be performed on the autopsied material, which has prevented the specific identification of tumor cells.Instead, HE staining was used in this work and the overall cellularity was analyzed.
The assumption was made that the cellularity baseline is approximately constant across healthy white matter and that the over-cellularity observed locally is exclusively attributed to tumor cell invasion since tumor-induced recruitment of inflammatory cells is limited in brain tissues.As a consequence, the identification of isolated infiltrating tumor cells on pathological examination may have been prevented.Third, the substantial deformation of the brain between in vivo and ex vivo imaging -with a volume decrease estimated to 16 % ex vivo based on MRI -may have resulted in partial distortion of the ex vivo MR-derived distance map compared to in vivo, not fully compensated by deformable registration of the in vivo edema outlines.
We would finally like to emphasize that the automated cell density map computation and histological slide to MR image registration procedures described in this work are not limited to the problem addressed herein and could be applied to various histological stainings, imaging modalities, and organs, such as prostate.Besides, the use of tailor-made 3D-printed slicer and cutting guides makes it possible to precisely analyze whole organ slices at low cost even for centers that do not have access to whole organ slice microscopy, opening tremendous possibilities for translational microscopic/macroscopic imaging research [32].

Conclusion
Through a translational radiological/histological analysis performed on a case of non-operated glioblastoma, we invalidated two commonly made assumptions relating the outlines of the visible abnormalities in magnetic resonance images to the tumor cell density function in the context of reaction-diffusion glioma growth modeling.We showed that, due to local restrictions of the tumor cell migration at brain tissue interfaces and along the brain boundary, the outlines of vasogenic edema in T2 FLAIR images do not generally coincide with a cell density iso-contour, as opposed to what was previously suggested.We also found that the commonly adopted tumor cell density iso-value at the edema outlines is likely overestimated since the over-cellularity at the maximum edema extent was found to be negligible in our studied case.This, however, does not exclude the possible presence of isolated tumor cells migrating beyond edema outlines, as previously reported.This work highlights the limitations of using routine magnetic resonance images to derive cell density maps for reaction-diffusion tumor growth models and points out the need of validating other methods to accurately initialize such models and make them usable for clinical applications.

A Slicer design
The in vivo T2 FLAIR image was first resampled to an isotropic voxel size of 0.5 mm × 0.5 mm × 0.5 mm.A brain mask was generated using the Otsu thresholding [33], followed by a morphological opening of radius 1, a largest component extraction and a morphological closing of radius 15.The brain mask bounding box was then drawn on a binary image with the same spacing and was evenly extended on both sides in each spatial dimension, with final dimensions of 160.5 mm × 190.5 mm × 172 mm.The brain mask volume was subtracted from the box volume.Ten slits of 12 mm × 190.5 mm × 152 mm spaced by 5.5 mm were carved into the box along the sagittal plane for brain slicing and the box was split in half along the axial plane at the level of the largest brain mask outline for brain insertion.Four cylindric fixations of radius 5.5 mm and height 9.5 mm were added to each corner of the upper part of the slicer and subtracted from its lower part.A 3D surface mesh was finally generated from the binary image using the marching cubes algorithm and exported in .stlformat for printing.The brain slicer design is depicted in Figure A1.
Ten cutting guides were also designed to ease the cutting of the brain slices into blocks.Those consisted of 190.5 mm × 149.5 mm × 11 mm plates from which the slice volume with a thickness of 5.5 mm was subtracted.Grooves were carved into the plates to guide the scalpel-cutting of 26 mm × 18 mm × 5.5 mm tissue blocks whose dimensions are compatible with the histological processing chain at our institution.All image processing and design steps were performed using an in-house software in C++ based on VTK [17] and ITK [18].
The brain slicer was printed using 1.75 mm PLA thread on a HICTOP 3DD-17-ATL-FM printer (HIC Technology, China) with a nozzle size of 0.4 mm (layer thickness 0.2 mm, honeycomb filling 5 %) for a total printing time of 96 hours.The printing of the ten cutting guides was parallelized on 5 Prusa i3 MK2 printers (Prusa Research, Czech Republic) with equivalent settings for a printing time of around 10 hours per guide.

B Deep learning-based nuclei counting
For each of the 28 stained slides resampled to a pixel size of 1 µm × 1 µm, 50 tiles of 100 px × 100 px were randomly chosen after exclusion of background tiles.Background tiles were defined as tiles containing less than 10 % of tissue pixels, empirically chosen as pixels whose the lowest RGB value is above 220 for our calibrated scanner.The 1600-tile dataset was then weakly supervised using an in-house annotation program in Python allowing the user to locate each cell nucleus within a tile by a simple mouse click.For each annotated tile, a nuclei presence probability map was generated by the superposition of 2D Gaussian blobs centered on user-pointed nucleus locations with full width at half maximum (FWHM) of 3 px (i.e. standard deviation σ ≈ 1.27 px), truncated to a radius of 4 σ.The supervised dataset was split into training and validation sets in proportion 80 %-20 %.
A U-Net [34] deep convolutional neural network was implemented for cell nuclei detection, consisting of 2 downsampling blocks, 3 up-sampling blocks and an output block.Each down-sampling block is made of two 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 max pool layer with kernel size 2 × 2 and stride 2 is added at the end of the block to reduce the feature maps dimensions by a factor 2. The up-sampling blocks are identical to the down-sampling blocks except that the max pool layer is replaced by a deconvolution 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. The output block has the same structure as the down-sampling blocks except that the max pool layer is replaced by a convolutional layer with kernel size 1 × 1 and stride 1 followed by a bias-adding layer and a sigmoid activation layer to merge the last 128 feature maps into a single presence probability map.The network architecture is depicted in Figure A2 and was implemented in Python using the TensorFlow framework (version 1.13.1)[35].
The network was trained on the weakly-supervised training set using the cross-entropy loss and the Adam optimizer [36] (learning rate: 10 −4 , β 1 : 0.9, β 2 : 0.999, : 10 −6 ).An evaluation was performed on the validation set at the end of each epoch and the training was stopped early after no improvement of the validation metric in 100 epochs.The network parameter values that gave the best validation metric value (0.041) were kept, which was reached at epoch 96.A similar weakly supervised deep-learning approach was used in [37] for cell nuclei detection in histological slides, but Euclidian distance to user-pointed nuclei locations and the mean squared error loss were used instead of Gaussian nuclei presence probability maps and the cross-entropy loss.
After training, each resampled scanned slide was divided into adjacent tiles of 100 px × 100 px from which nuclei presence probability maps were predicted by the trained network.A nuclei count was finally derived for each predicted probability map by detecting its local maxima with a value greater or equal to 0.1 and spaced by at least 3 px.C In vivo to ex vivo registration

Figure 1 :
Figure 1: Brain slicing and sample collection procedure.(a) The brain is placed inside the 3D-printed slicer.(b) Sagittal slices are cut carefully.(c) Each brain slice is placed inside its cutting guide.(d) Sample blocks are cut with a scalpel along the grooves and placed into standard cassettes.

Figure 2 :
Figure 2: Cell density map computation procedure.(a) 3 × 3 adjacent tiles (dotted squares) with dimensions 100 px × 100 px and pixel size 1 µm × 1 µm extracted from the resampled slide in panel (c).Cell nuclei detected by the deep convolutional neural network are indicated with cyan dots.(b) Corresponding 3 × 3 pixels (dotted squares) of the cell density map with pixel size 0.1 mm × 0.1 mm whose value is given by the corresponding tile cell count divided by the true tissue area.(c) Whole hematoxylin and eosin stained slide (slide 13, see TableA1).(d) Whole computed cell density map.

Figure 3 Figure 3 :
Figure 3 depicts an example of brain slice in its cutting guide (Figure 3(a)) with the corresponding registered in vivo T2 FLAIR image slice (Figure 3(b)), registered cell density maps (Figure 3(c)) and geodesic distance map (Figure3(d)).An over-cellularity front is visible in Figure3(c), progressing from the frontal necrotic tumor core but rapidly decreasing to reach a normal cellularity of around 1450 cell/mm 2 beyond a geodesic distance of 20 mm (Figure3(d)).The edema outlines, on the other hand, extend to over 50 mm on the depicted image slice (see red and blue delineations in Figure3(b)).The distinction between the red and blue segments of the edema outlines in Figure3(b) will be used in the discussion.
Figure 3: Cell density profile analysis.(a) Brain slice inside its 3D-printed cutting guide.(b) Corresponding slice of the registered in vivo T2 FLAIR image with segmented edema outlines.The blue and red segments of the outline respectively correspond to free and non-free to diffuse parts of the edema boundary (see Section 4).(c) Corresponding slice of the ex vivo T1 image (grayscale) and superimposed registered cell density maps (colored) with their slide number (see Table A1).(d) Corresponding slice of the geodesic distance map to the tumor core across white matter.More examples of registered cell density maps with the corresponding geodesic distance map slice are depicted in Figure 4.The decreasing behavior of the tumor cell density with the distance to the tumor core was observed among all these examples.

Figure 4 :
Figure 4: Example of registered cell density maps with their slide numbrt (see TableA1) (1 st and 3 rd columns) and corresponding slices of the geodesic distance map to the tumor core (2 nd and 4 th columns) superimposed to the ex vivo T1 image.

c core [10 3 Table 2 : 6 0Figure 5 :
Figure 5: Scatter plot of the surface cell density (a) and the extrapolated volume cell density (b) versus distance for each available value pairs among white matter voxels (blue dots) with superimposed fitted model curves (red curves).

FrequencyFigure 6 :
Figure 6: Inverse cumulative of the geodesic distance values along the edema outlines.The expected distribution under the hypothesis of iso-distance edema outlines is plotted in red.

Figure 7 :
Figure 7: Edema region (red) with superimposed thresholded region of the distance map whose contour minimizes the Hausdorff distance (blue, 1 st row) and average symmetric surface distance (blue, 2 nd row) to the edema contour in axial (1 st column), sagittal (2 nd column) and coronal (3 rd column) planes.

Figure A1 :
Figure A1: Brain slicer design superimposed to the in vivo T2 FLAIR image used as template in axial (a), sagittal (b), coronal (c) and 3D (d) views.

Figure
FigureA2: U-Net architecture with its feature map dimensions used for nuclei detection.The network takes 100 px × 100 px RGB tiles with pixel size 1 µm × 1 µm as input and predicts a nuclei presence probability map with the same spatial dimensions.