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

^{1}

^{2}

^{3}

^{4}

^{5}

^{*}

^{†}

^{‡}

Department of Nuclear Medicine, Hôpital Erasme, Université Libre de Bruxelles, Route de Lennik 808, 1070 Brussels, Belgium

Center for Microscopy and Molecular Imaging (CMMI), Université Libre de Bruxelles, Rue Adrienne Bolland 8, 6041 Charleroi, Belgium

Laboratory of Image Synthesis and Analysis (LISA), École Polytechnique de Bruxelles, Université Libre de Bruxelles, Avenue Franklin Roosevelt 50, 1050 Brussels, Belgium

Department of Pathology, Hôpital Erasme, Université Libre de Bruxelles, Route de Lennik 808, 1070 Brussels, Belgium

Department of Radiology, Hôpital Erasme, Université Libre de Bruxelles, Route de Lennik 808, 1070 Brussels, Belgium

Author to whom correspondence should be addressed.

These authors contributed equally to this work.

These authors contributed equally to this work.

Academic Editor: Emilio Quaia

Received: 6 September 2021
/
Revised: 22 October 2021
/
Accepted: 25 October 2021
/
Published: 29 October 2021

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.

Gliomas are the most common primary brain tumors. Diffuse gliomas, which include their most aggressive form glioblastoma (GBM), are known to be highly infiltrative [1], with the presence of tumor cells reported as far as 4 $\mathrm{c}$$\mathrm{m}$ 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 whole extent of the tumor cell invasion. Instead, tumor-induced alterations of the microenvironment 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-vascularization induced by VEGF, allowing gadolinium-based contrast agents to diffuse freely into brain tissues [3].

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 [6]. 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 tumor invasion extent would also be beneficial for improved surgery and radiotherapy planning [7]. However, due to their highly invasive nature and their long processing times, the number and frequency of biopsies are restricted, limiting the use of pathological examination as a proper tool for assessing tumor invasion and highlighting the complementarity of radiological examination [8]. In this matter, quantitatively linking glioma cell invasion patterns observed histologically to MR-visible abnormalities would be of great interest. Such information would allow noninvasive assessment of the glioma extent beyond the visible outlines of the tumor 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 1990s [9] are the most widely used, with potential applications for patient follow-up and improved radiotherapy planning [7]. These models rely on a reaction-diffusion equation to capture the spatial-temporal evolution of a tumor cell density function:
where $c(\mathit{r},t)$ is the tumor cell density at position $\mathit{r}$ and time t normalized by the maximum carrying capacity ${c}_{\mathrm{max}}$ of brain tissues (i.e., $c(\mathit{r},t)\in [0,1],\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\forall \mathit{r},t$), $\mathit{D}\left(\mathit{r}\right)$ is the symmetric tumor cell diffusion tensor at position $\mathit{r}$, $\rho $ is the tumor cell proliferation rate, ${c}_{0}\left(\mathit{r}\right)$ is the initial tumor cell density at position $\mathit{r}$, and ${\mathit{n}}_{{\partial}_{\mathsf{\Omega}}}\left(\mathit{r}\right)$ is a unit normal vector pointing outwards the boundary ${\partial}_{\mathsf{\Omega}}$ of the brain domain $\mathsf{\Omega}$ at position $\mathit{r}\in {\partial}_{\mathsf{\Omega}}$. The reaction term appearing in Equation (1) is logistic, which accounts for the (linear) saturation of the tumor cell proliferation due to resource restrictions and ensures a maximum cell density of ${c}_{\mathrm{m}ax}$ (i.e., a maximum normalized cell density c of 1). In the case of a constant and isotropic diffusion coefficient d, Equation (1) is also referred to as the Fisher’s equation, whose asymptotic behavior on the infinite cylinder has been extensively studied previously [10,11]. For the sake of completeness, it should, however, be noted that other reaction terms have been proposed previously for reaction diffusion models, namely, exponential $\rho \phantom{\rule{0.166667em}{0ex}}c(\mathit{r},t)$ (no saturation) and Gompertzian $\rho \phantom{\rule{0.166667em}{0ex}}c(\mathit{r},t)\phantom{\rule{0.166667em}{0ex}}\mathrm{ln}\left(1/c(\mathit{r},t)\right)$ (exponential saturation) [12,13]. 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 ${\partial}_{\mathsf{\Omega}}$.

$$\begin{array}{cc}\hfill \{& \begin{array}{}\hfill \mathrm{(1)}& \frac{\partial c(\mathit{r},\text{}t)}{\partial t}=\u25bd\xb7(\mathit{D}(\mathit{r})\u25bdc(\mathit{r},\text{}t))+\rho \text{}c(\mathit{r},\text{}t)(1-c(\mathit{r},\text{}t))\hfill & \hfill \forall \mathit{r}\in \mathsf{\Omega},\forall t0,\hfill \mathrm{(2)}& c(\mathit{r},\text{}0)={c}_{0}\left(\mathit{r}\right)\hfill & \hfill \forall \mathit{r}\in \mathsf{\Omega},\hfill \mathrm{(3)}& \mathit{D}\left(\mathit{r}\right)\u25bdc(\mathit{r},\text{}t)\xb7{\mathit{n}}_{\partial \mathsf{\Omega}}\left(\mathit{r}\right)=0\hfill & \hfill \forall \mathit{r}\in \partial \mathsf{\Omega},\end{array}\hfill \end{array}$$

Reaction-diffusion models as the one in Equations (1)–(3) are particularly attractive for clinical applications as 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 than across white matter tracts, Jbabdi and colleagues proposed a method to derive a tumor cell diffusion tensor $\mathit{D}\left(\mathit{r}\right)$ from diffusion tensor imaging (DTI) data [14]. For a more detailed overview of reaction-diffusion glioma growth modeling and its potential clinical applications, the reader is referred to the works in [7,9,11,12,15].

One problem arising when attempting to solve Equations (1)–(3) from actual imaging data of newly diagnosed glioma patients is to estimate the initial cell density ${c}_{0}\left(\mathit{r}\right)$ at every location of the brain domain $\mathsf{\Omega}$. To address this issue, early works on glioma growth modeling proposed to relate the MR-visible abnormalities observed at time t to the tumor cell density function $c(\mathit{r},t)$. For example, Swanson and colleagues suggested to model the MR imaging process as a simple cell density threshold function [15]:
where ${I}_{\mathrm{T}1\mathrm{Gd}}(\mathit{r},t)$ and ${I}_{\mathrm{T}2}(\mathit{r},t)$ are, respectively, the imaging functions of the T1Gd and T2/T2-FLAIR MR sequences indicating whether the abnormality is visible at location $\mathit{r}$ and time t on the sequence, and ${c}_{\mathrm{enhancing}}$ and ${c}_{\mathrm{edema}}$ are the corresponding tumor cell density detection thresholds. Based on these assumptions, it was 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:
where ${\partial}_{\mathrm{enhancing}}$ and ${\partial}_{\mathrm{edema}}$ are, respectively, the enhancing core and edema outlines. The authors also suggested hypothetical values for ${c}_{\mathrm{enhancing}}$ and ${c}_{\mathrm{edema}}$ of 0.80 and 0.16, respectively [15], 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)–(3) at imaging time satisfying Equation (6) [11]. More recently, the same group suggested that—for a spatially constant and isotropic diffusion coefficient ${d}_{\mathrm{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 [7]:
where ${\lambda}_{\mathrm{white}}$ is the infiltration length of tumor cells in white matter given by $\sqrt{{d}_{\mathrm{white}}/\rho}$. Provided two iso-cell density contours and a distance map to the tumor core, the value of ${\lambda}_{\mathrm{white}}$ could therefore theoretically be assessed using Equation (7). A similar reasoning had been previously applied in [9] for glioma growth modeling in computed tomography images, leading to the same expression as Equation (7) for the initial condition ${c}_{0}\left(\mathit{r}\right)$.

$$\begin{array}{cc}\hfill {I}_{\mathrm{T}1\mathrm{Gd}}(\mathit{r},t)& =\left\{\begin{array}{cc}1\phantom{\rule{1.em}{0ex}}\hfill & \mathrm{if}\phantom{\rule{4.pt}{0ex}}c(\mathit{r},t)\ge {c}_{\mathrm{enhancing}}\hfill \\ 0\phantom{\rule{1.em}{0ex}}\hfill & \mathrm{otherwise}\hfill \end{array}\right.,\hfill \end{array}$$

$$\begin{array}{cc}\hfill {I}_{\mathrm{T}2}(\mathit{r},t)& =\left\{\begin{array}{cc}1\phantom{\rule{1.em}{0ex}}\hfill & \mathrm{if}\phantom{\rule{4.pt}{0ex}}c(\mathit{r},t)\ge {c}_{\mathrm{edema}}\hfill \\ 0\phantom{\rule{1.em}{0ex}}\hfill & \mathrm{otherwise}\hfill \end{array}\right.,\hfill \end{array}$$

$$c(\mathit{r},t)=\left\{\begin{array}{cc}{c}_{\mathrm{enhancing}}\phantom{\rule{1.em}{0ex}}\hfill & \mathrm{for}\phantom{\rule{4.pt}{0ex}}\mathit{r}\in \partial {\mathsf{\Omega}}_{\mathrm{enhancing}}\hfill \\ {c}_{\mathrm{edema}}\phantom{\rule{1.em}{0ex}}\hfill & \mathrm{for}\phantom{\rule{4.pt}{0ex}}\mathit{r}\in \partial {\mathsf{\Omega}}_{\mathrm{edema}}\hfill \end{array}\right.,$$

$$c(\mathit{r},t)\propto exp\left(-\frac{d\left(\mathit{r}\right)}{{\lambda}_{\mathrm{white}}}\right),$$

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)–(7). However, the extent of vasogenic edema is, for example, known to be impacted by the administration of corticosteroids and anti-angiogenic treatments, independent of the tumor progression [3]. As reaction-diffusion models are highly sensitive to the provided initial tumor cell density distribution ${c}_{0}\left(\mathit{r}\right)$ [16], the validation of the aforementioned hypotheses is crucial for these models to be usable in clinical routine. In this work, we propose to verify these assumptions, as well as the value of 0.16 suggested for ${c}_{\mathrm{edema}}$ in [15], 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.

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, TMZ) were initiated right after diagnosis. In December 2017, after one cycle of TMZ, the patient died of a septic shock caused by bowel perforation. An autopsy was carried out within 12 h 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.

The T1 (TR = 8 $\mathrm{m}$$\mathrm{s}$, TE = $2.9$ $\mathrm{m}$$\mathrm{s}$, TI = 950 $\mathrm{m}$$\mathrm{s}$, FA = 8°) and T2-FLAIR (TR = 4800 $\mathrm{m}$$\mathrm{s}$, TE = 320 $\mathrm{m}$$\mathrm{s}$, TI = 1650 $\mathrm{m}$$\mathrm{s}$) MR images routinely acquired at diagnosis time on a 3T Achieva scanner (Philips Healthcare, Best, 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$ $\mathrm{m}$$\mathrm{s}$, TE = $3.164$ $\mathrm{m}$$\mathrm{s}$, TI = 450 $\mathrm{m}$$\mathrm{s}$, FA = 12°) of the brain placed inside the 3D-printed slicer (see below) was performed on a 3T Signa PET/MR scanner (GE Healthcare, Chicago, IL, USA) right before slicing. It should, however, be noted that brain fixation has caused convergence of the T1 values of white and gray matter, as reported in [17,18]. Consequently, the acquired ex vivo T1 image rather has a proton-density (PD) contrast. Furthermore, 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 hereafter.

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, then 3D printed. Such a slicer allows the brain to be repositioned 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 [19]. 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 of the cutting guides are available in Appendix A.

Twenty-eight 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 (ISO 15189). 5 $\mathsf{\mu}$$\mathrm{m}$ slides were cut from each sample and stained with hematoxylin and eosin (HE). The 28 stained slides were scanned in $20\times $ mode (per-mode=symbol $0.46$ $\mathsf{\mu}$$\mathrm{m}$/$\mathrm{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 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\text{}\mathsf{\mu}\mathrm{m}\times 1\text{}\mathsf{\mu}\mathrm{m}$ and divided into adjacent tiles of $100\text{}\mathrm{px}\times 100\text{}\mathrm{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}$ $\mathrm{m}$${\mathrm{m}}^{2}$). The computed cell density values were finally stored as a 2D image with a pixel size of $0.1\text{}\mathrm{m}\mathrm{m}\times 0.1\text{}\mathrm{m}\mathrm{m}$ where each pixel exactly corresponds to one of the $100\text{}\mathrm{px}\times 100\text{}\mathrm{px}$ tiles 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, as the former are the actual values of interest for the reaction-diffusion growth model in Equations (1)–(3). Under the assumptions that:

- cell nuclei are approximately spherical;
- cell nuclei distribution is locally isotropic;
- tile dimensions are sufficiently large to contain multiple cells;
- 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}_{\mathrm{volume}}$ can be extrapolated for a cube with the same side length l from the surface cell density ${c}_{\mathrm{surface}}$ by

$${c}_{\mathrm{volume}}={{c}_{\mathrm{surface}}}^{\frac{3}{2}}.$$

Substantial deformations of the brain occurred between the antemortem MR acquisitions and the postmortem 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$ $\mathrm{m}$$\mathrm{m}$. The 2D cell density maps were artificially extended to 3D images with a thickness of $0.5$ $\mathrm{m}$$\mathrm{m}$ 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 [20] and ITK [21].

The cell density map registration process was greatly facilitated by the use of a 3D-printed slicer as it allowed to impose the brain slicing orientation. The complex histology slide to MR image registration process in 3D was thus reduced to a simpler 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.

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 available as part of the Elastix software [22]. 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.

To verify the 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 [23]. Note that as no DTI images were available for the patient, the anisotropy of glioma cell diffusion mentioned above 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}_{\mathrm{Hausdorff}}$ and the ASSD ${d}_{\mathrm{ASSD}}$ between two sets A and B are, respectively, given by [24]
where $d(x,y)$ is the Euclidian distance between elements x and y, and $\left|X\right|$ is the cardinal of set X.

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

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

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 without distinction between tumor and healthy cells. To address this problem, we propose to verify the following equation instead of Equation (7):
where ${c}_{\mathrm{total}}$, ${c}_{\mathrm{tumor}}$, and ${c}_{\mathrm{white}}$ are, respectively, the total, tumor, and healthy cell densities in white matter, and ${c}_{\mathrm{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 a white matter baseline cellularity ${c}_{\mathrm{white}}$.

$$\begin{array}{cc}\hfill {c}_{\mathrm{total}}\left(\mathit{r}\right)& ={c}_{\mathrm{tumor}}\left(\mathit{r}\right)+{c}_{\mathrm{white}},\hfill \end{array}$$

$$\begin{array}{cc}\hfill \phantom{\rule{1.em}{0ex}}& ={c}_{\mathrm{core}}exp\left(-\frac{d\left(\mathit{r}\right)}{{\lambda}_{\mathrm{white}}}\right)+{c}_{\mathrm{white}},\hfill \end{array}$$

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\text{}\mathrm{mm}\times 0.5\text{}\mathrm{mm}\times 0.5\text{}\mathrm{mm}$) and all available pairs of density/distance values among the segmented white matter voxels were extracted. The values of ${c}_{\mathrm{core}}$, ${\lambda}_{\mathrm{white}}$, and ${c}_{\mathrm{white}}$ in Equation (12) were finally least-squares fitted to the available experimental density/distance pairs using SciPy’s ‘optimize’ module in Python [25].

Figure 3 depicts an example of brain slice in its cutting guide (Figure 3a) with the corresponding registered in vivo T2-FLAIR image slice (Figure 3b), registered cell density maps (Figure 3c), and geodesic distance map (Figure 3d). An over-cellularity front is visible in Figure 3c, progressing from the frontal necrotic tumor core but rapidly decreasing to reach an apparently normal cellularity of around 1450 $\mathrm{cell}$/$\mathrm{m}$${\mathrm{m}}^{2}$ beyond a geodesic distance of 20 $\mathrm{m}$${\mathrm{m}}^{2}$ (Figure 3d). The edema outlines, on the other hand, extend to over 50 $\mathrm{m}$$\mathrm{m}$ on the depicted image slice (see red and blue delineations in Figure 3b). The distinction between the red and blue segments of the edema outlines in Figure 3b will be used in the discussion.

More examples of registered cell density maps with the corresponding slice of the geodesic distance map 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.

The available pairs of density/distance values among all white matter voxels are plotted in Figure 5 for the surface density data (Figure 5a) and the volume density data extrapolated using Equation (8) (Figure 5b). 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.

The inverse cumulative distribution of the geodesic distance values along the edema outlines—i.e., the fraction of edema boundary voxels located at a geodesic distance greater or equal to a given value on the x-axis—is depicted in Figure 6 (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$ $\mathrm{m}$$\mathrm{m}$ and only 1% are located at a distance greater than $60.1$ $\mathrm{m}$$\mathrm{m}$ from the tumor core.

The distance threshold values which provided the smallest Hausdorff distance ($24.65$ $\mathrm{m}$$\mathrm{m}$) and ASSD ($1.97$ $\mathrm{m}$$\mathrm{m}$) between the edema outlines and the contour of the thresholded distance map region were $43.5$ $\mathrm{m}$$\mathrm{m}$ and $35.5$ $\mathrm{m}$$\mathrm{m}$, respectively. The corresponding circumscribed 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 $\mathrm{m}$$\mathrm{m}$ (slide 7), whereas edema was detected as far as $61.7$ $\mathrm{m}$$\mathrm{m}$ on the same slide.

The exponentially decreasing glioma cell density profile with distance to the tumor core suggested in [7,9] is compatible with our experimental data, as observed in Figure 5 for both surface (a) and extrapolated volume (b) cell densities. Besides, the fitted value of ${\lambda}_{\mathrm{white}}$ ($8.46$ $\mathrm{m}$$\mathrm{m}$) for the volume cell density data is in the same order of magnitude as the one used in [7] for simulations ($4.2$ $\mathrm{m}$$\mathrm{m}$). The baseline volume cell density value of $0.59\times {10}^{5}$ $\mathrm{cell}$ $\mathrm{m}$${\mathrm{m}}^{-3}$ in white matter is also in accordance with the literature, although large variations are observed between the reported values [26]. However, high variance was observed in our experimental data (see the blue point distribution in Figure 5a,b), imputed to errors accumulated throughout the numerous steps of our data processing pipeline. As a result, uncertainties on the fitted parameters are to be expected. In particular, the reported values of ${c}_{\mathrm{core}}$ in Table 1 and Table 2 were likely underestimated as cell density values of up to $6.1\times {10}^{3}$ $\mathrm{cell}$ $\mathrm{m}$${\mathrm{m}}^{-2}$ and $4.8\times {10}^{5}$ $\mathrm{cell}$ $\mathrm{m}$${\mathrm{m}}^{-3}$ were, respectively, observed in our surface and volume cell density data (see Figure 5a,b). Consequently, the corresponding fitted values of ${\lambda}_{\mathrm{white}}$ were likely overestimated.

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 Equations (5) and (6) can be explained by the thresholding behavior of the imaging function proposed in Equation (5), from which assumption in Equation (6) was deduced in [15]. In fact, thresholding a spatial function gives rise to iso-value contours only if the function is sufficiently smooth and continuous. In contrast, the tumor cell density function, discretized over the image voxel grid, has many discontinuities at interfaces between white and gray matter and along the brain domain boundaries. Indeed, the difference in tumor cell diffusivity between white and gray matter [7,11,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 level. Along the brain domain boundary, discontinuities are even more pronounced since no tumor cells are allowed to diffuse outside the brain domain (see Equation (3)), resulting in an accumulation of tumor cells. Consequently, edema outlines may not correspond to a cell density iso-contour even if the proposed thresholding behavior of the edema imaging function in Equation (5) turned out to be valid. As an illustration, the edema outlines in Figure 3b 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 or the presence of brain boundary (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$ $\mathrm{m}$$\mathrm{m}$ from the tumor core, which would correspond to the actual maximum edema extent. Finally, note that, due to the limited number of available cell density map voxels located on the edema outlines, the assumption in Equation (6) could not be directly verified, which motivated the indirect distance-based reasoning herein.

For reasons mentioned above, the thresholding behavior of the edema imaging function proposed 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 farther 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 (5) is therefore incompatible with our results. It should, however, be noted that this apparent inadequacy could still result from residual registration errors between in vivo and ex vivo images. Besides, the use of an isotropic metric tensor for the geodesic map computation instead of a DTI-derived anisotropic metric tensor may also have impacted the presented results. A DTI-derived anisotropic metric tensor would indeed have allowed to account for the well-established preferential migration of glioma cells along white matter tracts [27], originating from both mechanical [28] and molecular [29] factors. Registration of a DTI atlas may have been used to circumvent the lack of DTI data in this study. However, the substantial deformations of the brain observed ex vivo would have required the use of a deformable registration method. Whereas rigid or affine transforms of tensor fields can be trivially performed using matrix products, B-spline transforms of such fields are substantially more challenging as the coordinates system in which the tensor components are expressed is changed locally [30]. One solution could have been to approximate the B-spline transform by an affine transform at every voxel, as suggested in [30], but such developments were out of the scope of this study. We nevertheless investigated the potential impact of an anisotropic metric tensor on the computed distance—or rather ‘arrival time’—maps based on DTI acquisitions of a healthy volunteer. The results presented in Appendix E suggest even shorter arrival times in the contralateral hemisphere due to the high anisotropy in the corpus callosum region, which is not in favor of a coincidence between the edema outlines and a thresholded arrival time map contour (see Figure 7).

The iso-density value of 16% of the maximum cell carrying capacity suggested for the edema contours in [15] 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 (blue contour in Figure 3b), 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\times {10}^{1}$ $\mathrm{cell}$ $\mathrm{m}$${\mathrm{m}}^{-3}$ with regard to the white matter cellularity baseline ${c}_{\mathrm{white}}$ at a distance of $60.1$ $\mathrm{m}$$\mathrm{m}$ from the tumor core—corresponding to the maximum edema extent. This over-cellularity corresponds to only 0.08% of ${c}_{\mathrm{core}}\le {c}_{\mathrm{max}}$, which is far below the commonly accepted value of 16%. In addition, as the value of ${c}_{\mathrm{core}}$ in Table 2 was likely underestimated as mentioned above, an even lower relative value of ${c}_{\mathrm{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 $\mathrm{m}$$\mathrm{m}$ 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,31,32].

Although MRI provides high contrast in soft tissues and is the current 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 to treatment. 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 imaging function for the vasogenic edema in T2/T2-FLAIR images (see Equation (5)), the edema outlines 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% of the maximum carrying capacity. These results raise the question of the applicability of the previously proposed methods for deriving glioma cell density distributions from routine MR images referenced in the introduction. Indeed, whereas the method proposed in [11] 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,9], iso-density contours would still be required to assess the infiltration length parameter ${\lambda}_{\mathrm{white}}$. As a high sensibility of the reaction-diffusion tumor growth models to their initial condition was previously reported by our group [16], deriving an initial spatial cell density distribution from medical imaging data that is as reliable as possible is crucial for the model to be applied in clinical practice.

Beyond the field of reaction-diffusion tumor growth modeling, assessing the tumor cell density distribution in gliomas is also of utmost importance for treatment planning and response assessment [33,34,35,36]. To this extend, advanced MR acquisition and processing methods have been proposed previously. In [37], a linear relation was derived voxel-wise between the cell density and the average diffusion coefficient (ADC) value assessed by diffusion-weighted (DW) MRI. ADC maps were then used to assess the parameters and validate a simple logistic tumor growth model in rats with inoculated GBM cells. This approach was more recently used in patients to personalize a reaction-diffusion glioma growth model including the effects of chemoradiation [38]. Furthermore, microstructure imaging using advanced DW-MRI models could provide finer information on the tumor cell density [34,39] but also on the spatial distribution of tumor cells within a voxel [40]. However, this technique would require higher signal-to-noise ratios than achievable on clinical scanners to be used in practice [39]. In [33], glioma invasion is assessed using choline/N-acetylaspartate ratio maps derived from MR spectroscopy and showed promising relations with histology. In [36], anatomic, DW, and dynamic contrast-enhanced MR sequences are used to predict cell density values assessed by biopsies using a random forest classifier, with moderate-to-strong correlations reported between the observed and predicted density values. A similar approach was adopted in [35], with a reported root mean squared error of 1015 $\mathrm{cell}$ $\mathrm{m}$${\mathrm{m}}^{-2}$ on test data. Aside from MRI, the use of positron emission tomography (PET) imaging to assess glioma cell density has also been investigated previously. In [41], a linear relation was shown between the tumor cell density and the uptake of [^{18}F]fluoroethyl-L-tyrosine and used in [42] along with conventional MRI to personalize a reaction-diffusion tumor growth model. Finally, combined [^{11}C]methionine and [^{18}F]fluorodeoxyglucose PET imaging was used in [43] to improve the detection of glioma cell infiltration. As can be noticed, many imaging techniques have been proposed to assess glioma cell density distributions in vivo but no consensual clinical procedure has seemed to emerge so far, leaving the question open.

This study was however prone to several limitations. First, due to the scarcity of the human body material analyzed—a non-operated brain with GBM, this study was based on a single case and should be further carried out 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 also be of interest but is restricted due to reported substantial differences in cortical [44], glial [45], and endothelial [46] cells between human and mouse brain. Consequently, murine glioma models are known to only partially reflect the characteristics of human gliomas [47]. 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 was approximately constant across healthy white matter and that the over-cellularity observed locally was exclusively attributed to tumor cell invasion—which seems reasonable as 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 entirely compensated by deformable registration of the in vivo edema outlines. Fourth, the effect of the palliative TMZ treatment administered to the patient over the month between the in vivo MRI acquisition and the patient’s death could not be properly assessed as no further medical imaging acquisition was performed after the treatment start. Nevertheless, only limited effects would be expected as the patient underwent a single cycle of TMZ, whereas at least six cycles are classically recommended for GBMs [6,48]. Indeed, although a case of complete response after only one cycle of TMZ has been reported in [49], this case was still presented by the authors as unusual. Besides, the response to TMZ is highly dependent on the MGMT methylation status of the tumor [6,48,49], which could not be determined on this autopsied material. This latter limitation well reflects the difficulties encountered to clinically validate tumor growth models in typically multi-treated glioma patients.

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 [50].

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, contrary to what was suggested previously. We also found that the commonly adopted tumor cell density iso-value at the edema outlines is likely overestimated as the over-cellularity assessed 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 conventional magnetic resonance images to derive tumor cell density maps and points out the need of validating other methods to accurately initialize reaction-diffusion tumor growth models for clinical applications.

Conceptualization, C.M., L.L., C.D., S.G. and G.V.S.; methodology, C.M., L.L., T.V., Y.-R.V.E., A.R. and T.M.; software, C.M. and Y.-R.V.E.; validation, C.M., L.L., C.D., T.M., O.D., S.G., I.S. and G.V.S.; formal analysis, C.M., L.L. and C.D.; investigation, C.M., L.L., Y.-R.V.E., A.R. and T.M.; resources, T.V. and Y.-R.V.E.; data curation, C.M., L.L. and Y.-R.V.E.; writing—original draft preparation, C.M. and L.L.; writing—review and editing, C.M., L.L., C.D., T.V., Y.-R.V.E., A.R., T.M., O.D., S.G., I.S. and G.V.S.; visualization, C.M.; supervision, C.D., T.M., O.D., S.G., I.S. and G.V.S.; project administration, C.M. and L.L.; funding acquisition, C.D., O.D., S.G., I.S. and G.V.S. All authors have read and agreed to the published version of the manuscript.

C.M. is funded by FRIA grant number 5120417F (F.R.S.-FNRS) and the Walloon Region (Belgium, PROTHER-WAL). L.L. is funded by Fonds Erasme. 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). The CMMI is supported by the European Regional Development Fund and the Walloon Region. The Department of Pathology at Hôpital Erasme is supported by Fonds Erasme and Fonds Yvonne Boël.

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 P2019/245, approved on 3 May 2019).

Patient consent was waived by our institutional review board for this retrospective study on residual human body and radiological material in the context of a previously authorized autopsy.

The data presented in this study are available on reasonable request from the corresponding author.

The authors would like to thank the technical, medical, and scientific staff in charge of image acquisition and tissue sample processing for the needs of this work.

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

The following abbreviations are used in this manuscript:

ADC | Average diffusion coefficient |

AFM | Anisotropic fast marching |

ASSD | Average symmetric surface distance |

BBB | Blood–brain barrier |

DTI | Diffusion tensor imaging |

DW | Diffusion-weighted |

FA | Flip angle |

GBM | Glioblastoma |

GVP | Glomeruloid vascular proliferation |

HE | Hematoxylin and eosin |

IHC | Immunohistochemistry |

MR(I) | Magnetic resonance (imaging) |

PD | Proton density |

PET | Positron emission tomography |

PPN | Pseudo-palisading necrosis |

ReLU | Rectified linear unit |

Susp. | Suspected |

T1Gd | T1-weighted sequence with injection of gadolinium-based contrast agent |

T2-FLAIR | T2-weighted sequence with fluid-attenuated inversion recovery |

TE | Echo time |

TI | Inversion time |

TMZ | Temozolomide |

TR | Repetition time |

VEGF | Vascular endothelial growth factor |

WHO | World Health Organization |

The in vivo T2-FLAIR image was first resampled to an isotropic voxel size of $0.5\text{}\mathrm{m}\mathrm{m}\times 0.5\text{}\mathrm{m}\mathrm{m}\times 0.5\text{}\mathrm{m}\mathrm{m}$. A brain mask was generated using the Otsu thresholding [51], followed by a morphological opening of radius 1, a largest connected 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\text{}\mathrm{m}\mathrm{m}\times 190.5\text{}\mathrm{m}\mathrm{m}\times 172\text{}\mathrm{m}\mathrm{m}$. The brain mask volume was subtracted from the box volume. Ten slits of $12\mathrm{m}\mathrm{m}\times 190.5\mathrm{m}\mathrm{m}\times 152\mathrm{m}\mathrm{m}$ spaced by $5.5$ $\mathrm{m}$$\mathrm{m}$ 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 with radius $5.5$ $\mathrm{m}$$\mathrm{m}$ and height $9.5$ $\mathrm{m}$$\mathrm{m}$ 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 stl format 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\text{}\mathrm{m}\mathrm{m}\times 149.5\text{}\mathrm{m}\mathrm{m}\times 11\text{}\mathrm{m}\mathrm{m}$ plates from which the slice volume with a thickness of $5.5$ $\mathrm{m}$$\mathrm{m}$ was subtracted. Grooves were carved into the plates to guide the scalpel-cutting of $26\text{}\mathrm{m}\mathrm{m}\times 18\text{}\mathrm{m}\mathrm{m}\times 5.5\text{}\mathrm{m}\mathrm{m}$ 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 [20] and ITK [21].

The brain slicer was printed using $1.75$ $\mathrm{m}$$\mathrm{m}$ PLA thread on a HICTOP 3DD-17-ATL-FM printer (HIC Technology, Shenzhen, China) with a nozzle size of $0.4$ $\mathrm{m}$$\mathrm{m}$ (layer thickness $0.2$ $\mathrm{m}$$\mathrm{m}$, honeycomb filling 5%) for a total printing time of 96 h. The printing of the ten cutting guides was parallelized on 5 Prusa i3 MK2 printers (Prusa Research, Prague, Czech Republic) with equivalent settings for a printing time of around 10 h per guide.

For each of the 28 stained slides resampled to a pixel size of $1\text{}\mathsf{\mu}\mathrm{m}\times 1\text{}\mathsf{\mu}\mathrm{m}$, 50 tiles of $100\text{}\mathrm{px}\times 100\text{}\mathrm{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 of 3 $\mathrm{px}$ (i.e., standard deviation $\sigma \approx 1.27\mathrm{px}$), truncated to a radius of $4\phantom{\rule{0.166667em}{0ex}}\sigma $. The supervised dataset was split into training and validation sets in proportion 80–20%.

A U-Net [52] deep convolutional neural network was implemented for cell nuclei detection, consisting of two downsampling blocks, three upsampling blocks and an output block. Each downsampling block is made of two convolutional layers with kernel size $3\times 3$ and stride 1 followed by a bias-adding layer and a rectified linear unit (ReLU) activation layer. A max pooling layer with kernel size $2\times 2$ and stride 2 is added at the end of the block to reduce the feature maps dimensions by a factor 2. The upsampling blocks are identical to the downsampling blocks except that the max pooling layer is replaced by a deconvolution layer with kernel size $2\times 2$ and stride 2 followed by a bias-adding layer and a ReLU activation layer to expand the feature maps dimensions by a factor 2. Skip connections are added between the output of the second convolutional layer of each downsampling 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 downsampling blocks except that the max pooling layer is replaced by a convolutional layer with kernel size $1\times 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) [53].

The network was trained on the weakly-supervised training set using the cross-entropy loss and the Adam optimizer [54] (learning rate: 10^{−4}, ${\beta}_{1}$: 0.9, ${\beta}_{2}$: 0.999, $\u03f5$: 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 [55] 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\text{}\mathrm{px}\times 100\text{}\mathrm{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 $\mathrm{px}$.

Registration of the in vivo T2-FLAIR image to the ex vivo T1 image was performed using the Elastix software (version 4.801) [22] as described in Section 2.7. The Elastix parameter files for rigid and B-spline registration are, respectively, given in Listing A1 and Listing A2.

The results of the pathological examination and numerical tile processing are summarized in Table A1.

Cell Density [${10}^{3}$ $\mathbf{cell}$ $\mathbf{m}$${\mathbf{m}}^{-2}$] | Distance [$\mathbf{m}\mathbf{m}$] | |||||||||
---|---|---|---|---|---|---|---|---|---|---|

Block | PPN | Tumor Cells | GVP | Edema | Min | Max | Mean | Min | Max | Mean |

1 | No | No | No | No | 1.08 | 20.46 | 11.77 | 33.40 | 40.40 | 34.95 |

2 | No | No | No | No | 2.25 | 23.88 | 13.80 | 35.05 | 46.76 | 40.12 |

3 | No | Infiltrative (susp.) | No | Yes | 1.52 | 25.25 | 14.96 | 7.24 | 21.18 | 14.00 |

4 | No | No | No | No | 1.36 | 20.64 | 12.50 | 29.36 | 41.66 | 36.84 |

5 | No | No | No | Yes | 1.19 | 24.85 | 15.40 | 11.74 | 30.27 | 21.57 |

6 | No | Infiltrative (susp.) | No | No | 3.72 | 28.44 | 18.70 | 32.92 | 40.35 | 36.52 |

7 | No | Infiltrative (susp.) | No | Yes | 0.89 | 38.95 | 21.44 | 38.59 | 61.70 | 49.04 |

8 | Yes | Block | Yes | No | 3.37 | 37.29 | 23.65 | 0.50 | 5.17 | 2.71 |

9 | Yes | Infiltrative | Yes | Yes | 1.02 | 44.37 | 31.71 | 0.50 | 4.72 | 1.91 |

10 | No | Infiltrative (susp.) | No | Yes | 5.52 | 37.15 | 25.53 | 0.50 | 3.71 | 1.06 |

11 | No | No | No | No | 1.73 | 25.68 | 15.21 | 19.85 | 38.74 | 30.74 |

12 | No | Infiltrative | Yes | Yes | 0.89 | 36.08 | 19.70 | 0.50 | 31.46 | 14.24 |

13 | Yes | Block | Yes | No | 1.40 | 52.46 | 21.69 | 0.50 | 22.85 | 9.72 |

14 | No | No | No | No | 1.10 | 19.25 | 11.78 | 17.73 | 32.20 | 25.51 |

15 | No | No | No | Yes | 2.08 | 22.90 | 11.41 | 28.63 | 54.74 | 40.14 |

16 | No | No | No | No | 0.99 | 26.47 | 13.52 | 20.37 | 48.79 | 35.91 |

17 | No | No | No | No | 1.42 | 23.69 | 13.10 | 28.30 | 56.12 | 42.88 |

18 | No | No | Yes | Yes | 2.46 | 39.29 | 20.71 | 7.94 | 27.01 | 16.83 |

19 | Yes | Block | Yes | Yes | 0.99 | 44.12 | 19.82 | 0.50 | 18.50 | 6.92 |

20 | Yes | Block | Yes | No | 6.25 | 61.05 | 28.40 | 0.50 | 7.37 | 2.06 |

21 | No | Infiltrative | Yes | Yes | 4.75 | 34.20 | 23.05 | 0.50 | 14.62 | 7.70 |

22 | No | No | No | Yes | 0.86 | 30.37 | 10.79 | 2.99 | 36.87 | 21.39 |

23 | No | No | No | No | 0.94 | 25.61 | 13.23 | 21.14 | 49.61 | 34.48 |

24 | No | No | No | No | 1.37 | 26.72 | 13.72 | 57.27 | 90.93 | 71.08 |

25 | No | Infiltrative (susp.) | No | Yes | 0.87 | 39.07 | 21.03 | 1.71 | 21.41 | 10.87 |

26 | Yes | Block | Yes | Yes | 6.21 | 54.24 | 22.79 | 0.50 | 8.81 | 2.96 |

27 | No | No | No | Yes | 0.93 | 37.96 | 21.33 | 12.74 | 30.35 | 18.40 |

28 | No | No | No | No | 0.95 | 34.79 | 20.76 | 14.35 | 35.85 | 22.57 |

To assess the impact of an anisotropic metric tensor on the computed arrival time maps, T1 BRAVO (TR = $8.332$ $\mathrm{m}$$\mathrm{s}$, TE = $3.088$ $\mathrm{m}$$\mathrm{s}$, TI = 450 $\mathrm{m}$$\mathrm{s}$, FA = 12${}^{\xb0}$, pixel bandwidth = 244 $\mathrm{Hz}$ ${\mathrm{px}}^{-1}$, slice thickness/spacing = 1/1 $\mathrm{m}\mathrm{m}$, matrix = $240\text{}\mathrm{px}\times 240\text{}\mathrm{px}$) and EPI DTI (TR = 7000 $\mathrm{m}$$\mathrm{s}$, TE = $77.1$ $\mathrm{m}$$\mathrm{s}$, TI = 108 $\mathrm{m}$$\mathrm{s}$, FA = 90${}^{\xb0}$, pixel bandwidth = $1953.12$ $\mathrm{Hz}$ ${\mathrm{px}}^{-1}$, phase/slice acceleration factor = 2/1, multiband factor = 3, slice thickness/spacing = 2/2 $\mathrm{m}\mathrm{m}$, matrix = $120\text{}\mathrm{px}\times 120\text{}\mathrm{px}$, directions = 32, b-value = 1000 $\mathrm{s}$ $\mathrm{m}$${\mathrm{m}}^{-2}$) MR acquisitions of a healthy volunteer were performed on a 3T Signa PET/MR scanner (GE Healthcare, Chicago, IL, USA). The brain domain was extracted on the T1 image resampled to $1\text{}\mathrm{m}\mathrm{m}\times 1\text{}\mathrm{m}\mathrm{m}\times 1\text{}\mathrm{m}\mathrm{m}$, then segmented into white matter, gray matter, and cerebrospinal fluid using the MICO algorithm [56]. The DTI data were corrected for motion, eddy currents, and susceptibility artifacts using the ‘eddy’ and ‘topup’ tools and a water diffusion tensor field was reconstructed using the ‘dtifit’ tool, all these tools being available as part of the FSL software [57]. A unit tumor diffusion tensor field was derived from the segmented brain map and the DTI-derived water diffusion tensor field using the method proposed by Jbabdi and colleagues in [14] with a white matter/gray matter diffusivity ratio of 10 and an anisotropy factor of 10, then resampled to $1\text{}\mathrm{m}\mathrm{m}\times 1\text{}\mathrm{m}\mathrm{m}\times 1\text{}\mathrm{m}\mathrm{m}$. By unit tensor, it is meant that the tensor mean diffusivity $\mathrm{MD}=\mathrm{Tr}\left(\mathit{D}\right)/3$ is equal to 1 at every white matter voxel. Arrival time maps from a sphere with radius $2.5$ $\mathrm{m}$$\mathrm{m}$ located in the right frontal lobe—similarly to the GBM case studied—were finally computed using an adapted implementation of the AFM algorithm proposed in [23] with the inverse tumor diffusion tensor used as metric [27]. The results are depicted in Figure A3.

The results suggest even shorter arrival times in the contralateral hemisphere in the anisotropic case due to the high anisotropy in the corpus callosum region, which is not in favor of a coincidence between the edema outlines and a thresholded arrival time map contour (see Figure 7).

- Ostrom, Q.T.; Cioffi, G.; Gittleman, H.; Patil, N.; Waite, K.; Kruchko, C.; Barnholtz-Sloan, J.S. CBTRUS Statistical Report: Primary Brain and Other Central Nervous System Tumors Diagnosed in the United States in 2012–2016. Neuro-Oncol.
**2019**, 21, 1–100. [Google Scholar] [CrossRef] [PubMed] - Silbergeld, D.L.; Chicoine, M.L. Isolation and Characterization of Human Malignant Glioma Cells from Histologically Normal Brain. J. Neurosurg.
**1997**, 86, 525–531. [Google Scholar] [CrossRef] [PubMed] - Hawkins-Daarud, A.; Rockne, R.C.; Anderson, A.R.A.; Swanson, K.R. Modeling Tumor-Associated Edema in Gliomas during Anti-Angiogenic Therapy and Its Impact on Imageable Tumor. Front. Oncol.
**2013**, 3, 66. [Google Scholar] [CrossRef] - Lin, Z.X. Glioma-Related Edema: New Insight into Molecular Mechanisms and Their Clinical Implications. Chin. J. Cancer
**2013**, 32, 49–52. [Google Scholar] [CrossRef] - Lu, S.; Ahn, D.; Johnson, G.; Law, M.; Zagzag, D.; Grossman, R.I. Diffusion-Tensor MR Imaging of Intracranial Neoplasia and Associated Peritumoral Edema: Introduction of the Tumor Infiltration Index. Radiology
**2004**, 232, 221–228. [Google Scholar] [CrossRef] [PubMed] - Weller, M.; van den Bent, M.; Preusser, M.; Le Rhun, E.; Tonn, J.C.; Minniti, G.; Bendszus, M.; Balana, C.; Chinot, O.; Dirven, L.; et al. EANO Guidelines on the Diagnosis and Treatment of Diffuse Gliomas of Adulthood. Nat. Rev. Clin. Oncol.
**2021**, 18, 170–186. [Google Scholar] [CrossRef] [PubMed] - Unkelbach, J.; Menze, B.H.; Konukoglu, E.; Dittmann, F.; Le, M.; Ayache, N.; Shih, H.A. Radiotherapy Planning for Glioblastoma Based on a Tumor Growth Model: Improving Target Volume Delineation. Phys. Med. Biol.
**2014**, 59, 747–770. [Google Scholar] [CrossRef] - Wesseling, P.; Kros, J.M.; Jeuken, J.W.M. The Pathological Diagnosis of Diffuse Gliomas: Towards a Smart Synthesis of Microscopic and Molecular Information in a Multidisciplinary Context. Diagn. Histopathol.
**2011**, 17, 486–494. [Google Scholar] [CrossRef] - Tracqui, P.; Cruywagen, G.C.; Woodward, D.E.; Bartoo, G.T.; Murray, J.D.; Alvord, E.C. A Mathematical Model of Glioma Growth: The Effect of Chemotherapy on Spatio-Temporal Growth. Cell Prolif.
**1995**, 28, 17–31. [Google Scholar] [CrossRef] - Cencini, M.; Lopez, C.; Vergni, D. Reaction-Diffusion Systems: Front Propagation and Spatial Structures. In The Kolmogorov Legacy in Physics; Livi, R., Vulpiani, A., Eds.; Lecture Notes in Physics, Springer: Berlin/Heidelberg, Germany, 2003; pp. 328–341. [Google Scholar]
- Konukoglu, E.; Clatz, O.; Bondiau, P.Y.; Delingette, H.; Ayache, N. Extrapolating Glioma Invasion Margin in Brain Magnetic Resonance Images: Suggesting New Irradiation Margins. Med. Image Anal.
**2010**, 14, 111–125. [Google Scholar] [CrossRef] - Clatz, O.; Sermesant, M.; Bondiau, P.Y.; Delingette, H.; Warfield, S.K.; Malandain, G.; Ayache, N. Realistic Simulation of the 3-D Growth of Brain Tumors in MR Images Coupling Diffusion with Biomechanical Deformation. IEEE Trans. Med. Imag.
**2005**, 24, 1334–1346. [Google Scholar] [CrossRef] [PubMed] - Kuang, Y.; Nagy, J.D.; Eikenberry, S.E. Chapter 2—Introduction to Cancer Modeling. In Introduction to Mathematical Oncology; CRC Press: Abingdon, UK, 2016; pp. 45–92. [Google Scholar]
- Jbabdi, S.; Mandonnet, E.; Duffau, H.; Capelle, L.; Swanson, K.R.; Pélégrini-Issac, M.; Guillevin, R.; Benali, H. Simulation of Anisotropic Growth of Low-Grade Gliomas Using Diffusion Tensor Imaging. Magn. Reson. Med.
**2005**, 54, 616–624. [Google Scholar] [CrossRef] - Swanson, K.R.; Rostomily, R.C.; Alvord, E.C. A Mathematical Modelling Tool for Predicting Survival of Individual Patients Following Resection of Glioblastoma: A Proof of Principle. Br. J. Cancer
**2008**, 98, 113–119. [Google Scholar] [CrossRef] [PubMed] - Martens, C.; Metens, T.; Debeir, O.; Goldman, S.; Van Simaeys, G. Initial Condition Assessment from Patient MRI Data for Reaction-Diffusion Glioma Growth Models. Proc. Intl. Soc. Mag. Reson. Med.
**2019**, 27, 2852. [Google Scholar] - Tovi, M.; Ericsson, A. Measurement of T1 and T2 Over Time in Formalin-Fixed Human Whole-Brain Specimens. Acta Radiol.
**1992**, 33, 400–404. [Google Scholar] [CrossRef] [PubMed] - Raman, M.R.; Shu, Y.; Lesnick, T.G.; Jack, C.R.; Kantarci, K. Regional T1 Relaxation Time Constants in Ex Vivo Human Brain: Longitudinal Effects of Formalin Exposure. Magn. Reson. Med.
**2017**, 77, 774–778. [Google Scholar] [CrossRef] - Absinta, M.; Nair, G.; Filippi, M.; Ray-Chaudhury, A.; Reyes-Mantilla, M.I.; Pardo, C.A.; Reich, D.S. Postmortem Magnetic Resonance Imaging to Guide the Pathologic Cut: Individualized, 3-Dimensionally Printed Cutting Boxes for Fixed Brains. J. Neuropathol. Exp. Neurol.
**2014**, 73, 780–788. [Google Scholar] [CrossRef] - Schroeder, W.; Martin, K.; Lorensen, B. The Visualization Toolkit, 4th ed.; Kitware: Clifton Park, NY, USA, 2010. [Google Scholar]
- Yoo, T.; Ackerman, M.; Lorensen, W.; Schroeder, W.; Chalana, V.; Aylward, S.; Metaxas, D.; Whitaker, R. Engineering and Algorithm Design for an Image Processing API: A Technical Report on ITK—The Insight Toolkit. Stud. Health Technol. Inform.
**2002**, 85, 586–592. [Google Scholar] - Klein, S.; Staring, M.; Murphy, K.; Viergever, M.A.; Pluim, J. elastix: A Toolbox for Intensity-Based Medical Image Registration. IEEE Trans. Med. Imaging
**2010**, 29, 196–205. [Google Scholar] [CrossRef] - Konukoglu, E.; Sermesant, M.; Clatz, O.; Peyrat, J.M.; Delingette, H.; Ayache, N. A Recursive Anisotropic Fast Marching Approach to Reaction Diffusion Equation: Application to Tumor Growth Modeling. In Information Processing in Medical Imaging; Springer: Berlin/Heidelberg, Germany, 2007; pp. 687–699. [Google Scholar]
- Yeghiazaryan, V.; Voiculescu, I. Family of Boundary Overlap Metrics for the Evaluation of Medical Image Segmentation. J. Med. Imaging
**2018**, 5, 015006. [Google Scholar] [CrossRef] - Virtanen, P.; Gommers, R.; Oliphant, T.E.; Haberland, M.; Reddy, T.; Cournapeau, D.; Burovski, E.; Peterson, P.; Weckesser, W.; Bright, J.; et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nat. Methods
**2020**, 17, 261–272. [Google Scholar] [CrossRef] [PubMed] - Von Bartheld, C.S.; Bahney, J.; Herculano-Houzel, S. The Search for True Numbers of Neurons and Glial Cells in the Human Brain: A Review of 150 Years of Cell Counting. J. Comp. Neurol.
**2016**, 524, 3865–3895. [Google Scholar] [CrossRef] - Jbabdi, S.; Bellec, P.; Toro, R.; Daunizeau, J.; Pélégrini-Issac, M.; Benali, H. Accurate Anisotropic Fast Marching for Diffusion-Based Geodesic Tractography. Int. J. Biomed. Imaging
**2007**, 15, 320195. [Google Scholar] [CrossRef] - Johnson, J.; Nowicki, M.O.; Lee, C.H.; Chiocca, E.A.; Viapiano, M.S.; Lawler, S.E.; Lannutti, J.J. Quantitative Analysis of Complex Glioma Cell Migration on Electrospun Polycaprolactone Using Time-Lapse Microscopy. Tissue Eng. Part C Methods
**2007**, 2009, 531–540. [Google Scholar] [CrossRef] [PubMed] - Armento, A.; Ehlers, J.; Schötterl, S.; Naumann, U. Molecular Mechanisms of Glioma Cell Motility. In Glioblastoma; De Vleeschouwer, S., Ed.; Codon Publications: Brisbane, Australia, 2017; pp. 73–94. [Google Scholar]
- Gee, J.C.; Alexander, D.C. Diffusion-Tensor Image Registration. In Visualization and Processing of Tensor Fields; Weickert, J., Hagen, H., Eds.; Mathematics and Visualization; Springer: Berlin/Heidelberg, Germany, 2006; pp. 328–341. [Google Scholar]
- Kelly, P.J.; Daumas-Duport, C.; Kispert, D.B.; Kall, B.A.; Scheithauer, B.W.; Illig, J.J. Imaging-Based Stereotaxic Serial Biopsies in Untreated Intracranial Glial Neoplasms. J. Neurosurg.
**1987**, 66, 865–874. [Google Scholar] [CrossRef] - Sahm, F.; Capper, D.; Jeibmann, A.; Habel, A.; Paulus, W.; Troost, D.; von Deimling, A. Addressing Diffuse Glioma as a Systemic Brain Disease with Single-Cell Analysis. Arch. Neurol.
**2012**, 69, 523–526. [Google Scholar] [PubMed] - Ganslandt, O.; Stadlbauer, A.; Fahlbusch, R.; Kamada, K.; Buslei, R.; Blumcke, I.; Moser, E.; Nimsky, C. Proton Magnetic Resonance Spectroscopic Imaging Integrated into Image-guided Surgery: Correlation to Standard Magnetic Resonance Imaging and Tumor Cell Density. Neurosurgery
**2005**, 56, 291–298. [Google Scholar] [CrossRef] [PubMed] - Roberts, T.A.; Hyare, H.; Agliardi, G.; Hipwell, B.; d’Esposito, A.; Ianus, A.; Breen-Norris, J.O.; Ramasawmy, R.; Taylor, V.; Atkinson, D.; et al. Noninvasive Diffusion Magnetic Resonance Imaging of Brain Tumour Cell Size for the Early Detection of Therapeutic Response. Sci. Rep.
**2020**, 10, 9223. [Google Scholar] [CrossRef] - Bobholz, S.A.; Lowman, A.K.; Brehler, M.; Kyereme, F.; Duenweg, S.R.; Sherman, J.; McGarry, S.; Cochran, E.J.; Connelly, J.; Mueller, W.M.; et al. Radio-Pathomic Maps of Cell Density Identify Glioma Invasion Beyond Traditional MR Imaging Defined Margins. bioRxiv
**2021**. [Google Scholar] [CrossRef] - Gates, E.D.H.; Weinberg, J.S.; Prabhu, S.S.; Lin, J.S.; Hamilton, J.; Hazle, J.D.; Fuller, G.N.; Baladandayuthapani, V.; Fuentes, D.T.; Schellingerhout, D. Estimating Local Cellular Density in Glioma Using MR Imaging Data. AJNR Am. J. Neuroradiol.
**2021**, 42, 102–108. [Google Scholar] [CrossRef] - Atuegwu, N.C.; Colvin, D.C.; Loveless, M.E.; Xu, L.; Gore, J.C.; Yankeelov, T.E. Incorporation of Diffusion-Weighted Magnetic Resonance Imaging Data into a Simple Mathematical Model of Tumor Growth. Phys. Med. Biol.
**2012**, 57, 225–240. [Google Scholar] [CrossRef] [PubMed] - Hormuth, D.A.; Al Feghali, K.A.; Elliott, A.M.; Yankeelov, T.E.; Chung, C. Image-Based Personalization of Computational Models for Predicting Response of High-Grade Glioma to Chemoradiation. Sci. Rep.
**2021**, 11, 8520. [Google Scholar] [CrossRef] [PubMed] - McHugh, D.J.; Hubbard Cristinacce, P.L.; Naish, J.H.; Parker, G.J.M. Towards a ‘Resolution Limit’ for DW-MRI Tumor Microstructural Models: A Simulation Study Investigating the Feasibility of Distinguishing between Microstructural Changes. Magn. Reson. Med.
**2019**, 81, 2288–2301. [Google Scholar] [CrossRef] [PubMed] - Szczepankiewicz, F.; Lasič, S.; van Westen, D.; Sundgren, P.C.; Englund, E.; Westin, C.; Ståhlberg, F.; Lätt, J.; Topgaard, D.; Nilsson, M. Quantification of Microscopic Diffusion Anisotropy Disentangles Effects of Orientation Dispersion from Microstructure: Applications in Healthy Volunteers and in Brain Tumors. NeuroImage
**2015**, 104, 241–252. [Google Scholar] [CrossRef] [PubMed] - Stockhammer, F.; Plotkin, M.; Amthauer, H.; van Landeghem, F.K.H.; Woiciechowsky, C. Correlation of F-18-Fluoro-Ethyl-Tyrosin Uptake with Vascular and Cell Density in Non-Contrast-Enhancing Gliomas. J. Neuro-Oncol.
**2008**, 88, 205–210. [Google Scholar] [CrossRef] - Lipková, J.; Angelikopoulos, P.; Wu, S.; Alberts, E.; Wiestler, B.; Diehl, C.; Preibisch, C.; Pyka, T.; Combs, S.; Hadjidoukas, P.; et al. Personalized Radiotherapy Design for Glioblastoma: Integrating Mathematical Tumor Models, Multimodal Scans and Bayesian Inference. IEEE Trans. Med. Imaging
**2019**, 38, 1875–1884. [Google Scholar] [CrossRef] - Kinoshita, M.; Arita, H.; Goto, T.; Okita, Y.; Isohashi, K.; Watabe, T.; Kagawa, N.; Fujimoto, Y.; Kishima, H.; Shimosegawa, E.; et al. A Novel PET Index, 18F-FDG–11C-Methionine Uptake Decoupling Score, Reflects Glioma Cell Infiltration. J. Nucl. Med.
**2012**, 53, 1701–1708. [Google Scholar] [CrossRef] - Hodge, R.D.; Bakken, T.E.; Miller, J.A.; Smith, K.A.; Barkan, E.R.; Graybuck, L.T.; Close, J.L.; Long, B.; Johansen, N.; Penn, O.; et al. Conserved Cell Types with Divergent Features in Human Versus Mouse Cortex. Nature
**2019**, 573, 61–68. [Google Scholar] [CrossRef] - Zhang, Y.; Sloan, S.A.; Clarke, L.E.; Caneda, C.; Plaza, C.A.; Blumenthal, P.D.; Vogel, H.; Steinberg, G.K.; Edwards, M.S.B.; Li, G.; et al. Purification and Characterization of Progenitor and Mature Human Astrocytes Reveals Transcriptional and Functional Differences with Mouse. Neuron
**2016**, 89, 37–53. [Google Scholar] [CrossRef] - Song, H.W.; Foreman, K.L.; Gastfriend, B.D.; Kuo, J.S.; Palecek, S.P.; Shusta, E.V. Transcriptomic Comparison of Human and Mouse Brain Microvessels. Sci. Rep.
**2020**, 10, 12358. [Google Scholar] [CrossRef] - Kijima, N.; Kanemura, Y. Mouse Models of Glioblastoma. In Glioblastoma; De Vleeschouwer, S., Ed.; Codon Publications: Brisbane, Australia, 2017; pp. 131–140. [Google Scholar]
- Blumenthal, D.T.; Gorlia, T.; Gilbert, M.R.; Kim, M.M.; Nabors, L.B.; Mason, W.P.; Hegi, M.E.; Zhang, P.; Golfinopoulos, V.; Perry, J.R.; et al. Is More Better? The Impact of Extended Adjuvant Temozolomide in Newly Diagnosed Glioblastoma: A Secondary Analysis of EORTC and NRG Oncology/RTOG. Neuro-Oncol.
**2017**, 19, 1119–1126. [Google Scholar] [CrossRef] - Ducray, F.; Benouaich-Amiel, A.; Idbaih, A.; Rousseau, A.; Laigle-Donadey, F.; Delattre, J.; Sanson, M. Complete Response After One Cycle of Temozolomide in an Elderly Patient with Glioblastoma and Poor Performance Status. J. Neuro-Oncol.
**2008**, 88, 185–188. [Google Scholar] [CrossRef] - Baldi, D.; Aiello, M.; Duggento, A.; Salvatore, M.; Cavaliere, C. MR Imaging-Histology Correlation by Tailored 3D-Printed Slicer in Oncological Assessment. Contrast Media Mol. Imaging
**2019**, 2019, 1071453. [Google Scholar] [CrossRef] [PubMed] - Otsu, N. A Threshold Selection Method from Gray-Level Histograms. IEEE Trans. Syst. Man Cybern. Syst.
**1979**, 9, 62–66. [Google Scholar] [CrossRef] - Ronneberger, O.; Fischer, P.; Brox, T. U-Net: Convolutional Networks for Biomedical Image Segmentation. In Medical Image Computing and Computer-Assisted Intervention; Springer: Cham, Switzerland, 2015; pp. 234–241. [Google Scholar]
- Abadi, M.; Agarwal, A.; Barham, P.; Brevdo, E.; Chen, Z.; Citro, C.; Corrado, G.S.; Davis, A.; Dean, J.; Devin, M.; et al. TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems. Available online: https://www.tensorflow.org/ (accessed on 6 September 2021).
- Kingma, D.; Ba, J. ADAM: A Method for Stochastic Optimization. arXiv
**2014**, arXiv:1412.6980. [Google Scholar] - Xing, F.; Cornish, T.C.; Bennett, T.; Ghosh, D.; Yang, L. Pixel-to-Pixel Learning with Weak Supervision for Single-Stage Nucleus Recognition in Ki67 Images. IEEE Trans. Biomed. Eng.
**2019**, 66, 3088–3097. [Google Scholar] [CrossRef] - Li, C.; Gore, J.C.; Davatzikos, C. Multiplicative Intrinsic Component Optimization (MICO) for MRI Bias Field Estimation and Tissue Segmentation. Magn. Reson. Imaging
**2014**, 32, 913–923. [Google Scholar] [CrossRef] - Jenkinson, M.; Beckmann, C.F.; Behrens, T.E.; Woolrich, M.W.; Smith, S.M. FSL. NeuroImage
**2012**, 62, 782–790. [Google Scholar] [CrossRef]

${\mathit{c}}_{\mathbf{core}}$ [${10}^{3}$ $\mathbf{cell}$ $\mathbf{m}$${\mathbf{m}}^{-2}$] | ${\mathit{\lambda}}_{\mathbf{white}}$ [$\mathbf{m}\mathbf{m}$] | ${\mathit{c}}_{\mathbf{white}}$ [${10}^{3}$ $\mathbf{cell}$ $\mathbf{m}$${\mathbf{m}}^{-2}$] |
---|---|---|

1.47 | 10.55 | 1.43 |

${\mathit{c}}_{\mathbf{core}}$ [${10}^{5}$ $\mathbf{cell}$ $\mathbf{m}$${\mathbf{m}}^{-3}$] | ${\mathit{\lambda}}_{\mathbf{white}}$ [$\mathbf{m}\mathbf{m}$] | ${\mathit{c}}_{\mathbf{white}}$ [${10}^{5}$ $\mathbf{cell}$ $\mathbf{m}$${\mathbf{m}}^{-3}$] |
---|---|---|

1.05 | 8.46 | 0.59 |

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. |

© 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).