Initial condition assessment for reaction-diffusion glioma growth models: A translational MRI/histology (in)validation study
Corentin Martens, Laetitia Lebrun, Christine Decaestecker, Thomas Vandamme, Yves-Rémi Van Eycke, Antonin Rovai, Thierry Metens, Olivier Debeir, Serge Goldman, Isabelle Salmon, Gaetan Van Simaeys
II NITIAL CONDITION ASSESSMENT FOR REACTION - DIFFUSIONGLIOMA GROWTH MODELS : A
TRANSLATIONAL
MRI/
HISTOLOGY ( IN ) VALIDATION STUDY
A P
REPRINT
Corentin Martens
Department of Nuclear Medicine – Hôpital ErasmeUniversité libre de BruxellesBrussels, Belgium [email protected]
Laetitia Lebrun
Department of Pathology – Hôpital ErasmeUniversité libre de BruxellesBrussels, Belgium
Christine Decaestecker
Laboratory of Image Synthesis and AnalysisUniversité libre de BruxellesBrussels, Belgium
Thomas Vandamme
Laboratory of Image Synthesis and AnalysisUniversité libre de BruxellesBrussels, Belgium
Yves-Rémi Van Eycke
DIAPath – Center for Microscopy and Molecular ImagingUniversité libre de BruxellesCharleroi (Gosselies), Belgium
Antonin Rovai
Department of Nuclear Medicine – Hôpital ErasmeUniversité libre de BruxellesBrussels, Belgium
Thierry Metens
Department of Radiology – Hôpital ErasmeUniversité libre de BruxellesBrussels, Belgium
Olivier Debeir
Laboratory of Image Synthesis and AnalysisUniversité libre de BruxellesBrussels, Belgium
Serge Goldman
Department of Nuclear Medicine – Hôpital ErasmeUniversité libre de BruxellesBrussels, Belgium
Isabelle Salmon
Department of Pathology – Hôpital ErasmeUniversité libre de BruxellesBrussels, Belgium
Gaetan Van Simaeys
Department of Nuclear Medicine – Hôpital ErasmeUniversité libre de BruxellesBrussels, BelgiumFebruary 4, 2021 A BSTRACT
Diffuse gliomas are highly infiltrative tumors whose early diagnosis and follow-up usually rely onmagnetic resonance imaging (MRI). However, the limited sensitivity of this technique makes itimpossible to directly assess the extent of the glioma cell invasion, leading to sub-optimal treatmentplaning. Reaction-diffusion growth models have been proposed for decades to extrapolate gliomacell infiltration beyond margins visible on MRI and predict its spatial-temporal evolution. Thesemodels nevertheless require an initial condition, that is the tumor cell density values at every location a r X i v : . [ phy s i c s . m e d - ph ] F e b PREPRINT - F
EBRUARY
4, 2021of the brain at diagnosis time. Several works have proposed to relate the tumor cell density functionto abnormality outlines visible on MRI but the underlying assumptions have never been verifiedso far. In this work we propose to verify these assumptions by stereotactic histological analysisof a non-operated brain with glioblastoma using a tailored 3D-printed slicer. Cell density mapsare computed from histological slides using a deep learning approach. The density maps are thenregistered to a postmortem
MR image and related to an MR-derived geodesic distance map to thetumor core. The relation between the edema outlines visible on T2 FLAIR MRI and the distanceto the core is also investigated. Our results suggest that (i) the previously suggested exponentialdecrease of the tumor cell density with the distance to the tumor core is not unreasonable but (ii) theedema outlines may in general not correspond to a cell density iso-contour and (iii) the commonlyadopted tumor cell density value at these outlines is likely overestimated. These findings highlightthe limitations of using conventional MRI to derive glioma cell density maps and point out the needof validating other methods to initialize reaction-diffusion growth models and make them usable inclinical practice. K eywords Cellularity · Glioma · Histology · Magnetic Resonance Imaging · Reaction-Diffusion Equation · TumorGrowth Modeling · Gliomas are the most common primary brain tumors. Diffuse gliomas, which include its most aggressive formglioblastoma (GBM), are known to be highly infiltrative [1], with the presence of tumor cells reported as far as 4 cmfrom 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 intothe tumor biology, none of the routinely acquired MR sequences allows to directly assess the extent of the tumorcell invasion. Instead, tumor-induced alterations of the micro-environment are seen on MR images such as peritumorvasogenic edema visible on T2/T2 FLAIR sequences and the enhancing tumor core visible on T1-weighted sequenceswith injection of gadolinium-based contrast agent (T1Gd). Peritumor vasogenic edema originates from an increase inthe blood-brain-barrier (BBB) permeability induced by the release of vascular endothelial growth factor (VEGF) bytissues under hypoxic stress [3, 4] combined with changes in the brain hydrodynamic pressure [5]. The formation ofan 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 goldstandard method to confirm the diagnosis and determine the histological type and grade. According to the 2016WHO classification of the central nervous system tumors, the identification of infiltrative patterns is essential in thedifferential diagnosis of diffuse versus pilocytic astrocytomas, a more-circumscribed neoplasm. An overall assessmentof the invasion extent would also be beneficial for surgery and radiotherapy planning. However, due to their highinvasiveness and long processing times, the number and frequency of biopsy procedures are restricted, limiting theuse of pathological examination as a proper tumor invasion assessment tool and highlighting the complementarity ofradiological examination [6]. In this matter, quantitatively linking glioma cell invasion patterns observed histologicallyto MR-visible abnormalities would be of great interest. Such information would allow to non-invasively assess gliomaextent and progression, within and beyond the abnormality outlines, while providing a better interpretation of theobserved abnormalities.Mathematical glioma growth modeling has addressed the problem of estimating glioma cell distribution within braintissues and predicting its temporal evolution. Among the investigated models, reaction-diffusion models first introducedby Murray and colleagues in the early 1990’s [7] are probably the most widely used, with potential applications forpatient follow-up and improved radiotherapy planning [8]. These models rely on a reaction-diffusion equation to capturethe spatial-temporal evolution of a tumor cell density function: ∂c ( r , t ) ∂t = ∇ · ( D ( r ) ∇ c ( r , t )) + ρ c ( r , t ) (1 − c ( r , t )) ∀ r ∈ Ω , ∀ t > c ( r ,
0) = c ( r ) ∀ r ∈ Ω D ( r ) ∇ c ( r , t ) · n ∂ Ω ( r ) = 0 ∀ r ∈ ∂ Ω (1)(2)(3)2 PREPRINT - F
EBRUARY
4, 2021where c ( r , t ) is the tumor cell density at position r and time t normalized by the maximum carrying capacity c max ofbrain tissues ( c ( r , t ) ∈ [0 , , ∀ r , t ), D ( r ) is the symmetric tumor cell diffusion tensor at position r , ρ is the tumorcell proliferation rate, c ( r ) is the initial tumor cell density at position r , and n ∂ Ω ( r ) is a unit normal vector pointingoutwards the boundary ∂ Ω of the brain domain Ω at position r ∈ ∂ Ω . Equation (2) specifies the initial condition of theproblem. Equation (3) provides no-flux Neumann boundary conditions reflecting the inability of tumor cells to diffuseacross ∂ Ω .Reaction-diffusion models as the one in Equations (1) to (3) are particularly attractive for clinical applications since theyonly have a few parameters that could be assessed from patient imaging data. For instance, based on prior observationsthat tumor cells rather migrate along white matter tracts, Jbabdi and colleagues proposed to derive the tumor celldiffusion tensor D ( r ) from diffusion tensor imaging (DTI) data [9]. For a more detailed overview of reaction-diffusionglioma growth modeling and its potential clinical applications, the reader is referred to [7, 8, 10–12].One problem arising when attempting to solve Equations (1) to (3) from actual imaging data of newly diagnosed gliomapatients is to estimate the initial cell density c ( r ) at every location of the brain domain Ω . Early works on gliomagrowth 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, thatis: I T1Gd ( r , t ) = (cid:26) if c ( r , t ) ≥ c enhancing otherwise (4) I T2 ( r , t ) = (cid:26) if c ( r , t ) ≥ c edema otherwise (5)where I T1Gd ( r , t ) and I T2 ( r , t ) are respectively the imaging functions of the T1Gd and T2/T2 FLAIR MR sequencesindicating whether the abnormality is visible at location r and time t on the sequence, and c enhancing and c edema are thecorresponding tumor cell density detection thresholds. Based on these assumptions, the authors suggested that theoutlines of the tumor enhancing core in T1Gd images and of the vasogenic edema in T2/T2 FLAIR images wouldcorrespond to iso-contours of the tumor cell density function: c ( r , t ) = (cid:26) 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 suggestedhypothetical values for c enhancing and c edema of 0.80 and 0.16 respectively [11], although no rationale was provided forthese values. Building upon this work, Konukoglu and colleagues proposed a fast-marching approach to construct anapproximate solution of Equations (1) to (3) at imaging time satisfying Equation (6) [12]. More recently, the samegroup suggested that, for a spatially constant and isotropic diffusion coefficient d white and from a certain distance to thetumor core, the tumor cell density in white matter would approximately decrease exponentially with the distance d tothe core [8]: c ( r , t ) ∝ exp (cid:18) − d ( r ) λ white (cid:19) (7)where λ white is the infiltration length of tumor cells in white matter given by (cid:112) d white /ρ . Provided two iso-cell densitycontours and a distance map to the tumor core, the value of λ white can theoretically be assessed using Equation (7). Asimilar 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 ( r ) .Nevertheless, these attempts to derive a tumor cell density distribution from MR images rely on the existence of celldensity iso-contours in MR images and are based on unverified assumptions in Equations (4) to (7). However, as will befurther discussed, the extent of vasogenic edema is known to be impacted by the administration of corticosteroids andanti-angiogenic treatments [3]. Furthermore, as previously observed by our group in [13], the reaction-diffusion modelis highly sensitive to the provided initial tumor cell density distribution c ( r ) . The validation of the aforementioned3 PREPRINT - F
EBRUARY
4, 2021hypotheses is thus crucial for the model to be usable in clinical routine. In this work, we propose to verify theseassumptions, as well as the value of 0.16 for c edema suggested in [11], through a translational MRI/histology studyconducted in a case of non-operated GBM. To this end, stereotactic histological analyses are performed using a3D-printed slicer designed from antemortem MRI data. Cell density maps are computed automatically from the scannedhistological slides using a deep convolutional neural network and related to an MR-derived geodesic distance map tothe tumor core. The relation between the edema outlines and the geodesic distance to the core is also investigated. Ourresults highlight the limitations of using routine MRI to derive glioma cell density maps and point out the need for othervalidated 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, whichrevealed a massive right-frontal expanding lesion with intense heterogeneous enhancement post-injection of gadoliniumcontrast agent and a necrotic core, surrounded by a large area of perilesional edema. Considering the patient’s age, aconsensual 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 initiatedright after diagnosis. In December 2017, after one cycle of temozolomide, the patient died of a septic shock causedby bowel perforation. An autopsy was carried out within 12 hours after death as part of which the patient’s brainwas collected and fixed in formalin for 24 months. Upon microscopic examination, an infiltrating glioma with highcellularity, astrocytic phenotype cells and nuclear atypia was observed along with glomeruloid vascular proliferationand areas of pseudo-palissading necrosis. All procedures performed in this study were approved by the Hospital-FacultyEthics Committee of Hôpital Erasme (accreditation 021/406) under reference P2019/245 and are in accordance with the1964 Helsinki declaration and its later amendments or comparable ethical standards.
The T1 (TR = , TE = . , TI =
950 ms , FA = °) and T2 FLAIR (TR = , TE =
320 ms , TI = )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 anddelineation of the vasogenic edema.Additionally, a T1 BRAVO ‘ ex vivo ’ acquisition (TR = 8.264 ms, TE = .
164 ms , TI =
450 ms , FA = °) of the brainplaced 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 matterT1 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.
To relate histological observations to the abnormalities in MR images and to the MR-derived distance map to the tumorcore (see below), a brain slicer was designed based on the in vivo
T2 FLAIR image and 3D-printed. Such a slicer allowsthe brain to be re-positioned in antemortem imaging orientation and facilitates the cutting of sagittal brain slices. Theslicer 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 ourhistological 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 theslicer and the cutting guides are available in Appendix A.
28 tissue samples were selected within the tumor core, as well as within and beyond the peritumoral edematous regionbased on the in vivo
T2 FLAIR image. The samples were formalin-fixed and paraffin-embedded (FFPE, ISO 15189). µ m slides were cut from each sample and stained with hematoxylin and eosin (HE). The 28 stained slides werescanned in × mode ( . µ 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 blinded4 PREPRINT - F
EBRUARY
4, 2021 (a) (b)(c) (d)
Figure 1: Brain slicing and sample collection procedure. (a) The brain is placed inside the 3D-printed slicer. (b) Sagittalslices are cut carefully. (c) Each brain slice is placed inside its cutting guide. (d) Sample blocks are cut with a scalpelalong the grooves and placed into standard cassettes.to MRI for the presence of pseudo-palisading necrosis, tumor cells (in block or infiltrating), glomeruloid vascularproliferation and edema. As will be further discussed, immunohistochemistry (IHC) staining was also investigated butdid 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 µ m × µ m and divided into adjacent tiles of
100 px ×
100 px . Cell nuclei within each tile were automatically counted using a weakly-supervised deep learningapproach detailed in Appendix B. The counting result was divided by the actual tissue area within the tile, defined as thenumber of tissue pixels (see Appendix B) times the pixel area ( − mm ). The computed cell density was finally storedas a 2D image with a pixel size of . × . where each pixel exactly corresponds to one
100 px ×
100 px tileof 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 theformer 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 cubewith the same side length l from the surface cell density c surface by: c volume = √ c surface (8)5 PREPRINT - F
EBRUARY
4, 2021 (a) (b) (c) (d)
Figure 2: Cell density map computation procedure. (a) × adjacent tiles (dotted squares) with dimensions
100 px ×
100 px and pixel size µ m × µ m extracted from the resampled slide in panel (c). Cell nuclei detected by the deepconvolutional neural network are indicated with cyan dots. (b) Corresponding × pixels (dotted squares) of the celldensity map with pixel size . × . whose value is given by the corresponding tile cell count divided bythe true tissue area. (c) Whole hematoxylin and eosin stained slide (slide 13, see Table A1). (d) Whole computed celldensity map. 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 registeredto the ex vivo
T1 image as follows. The ex vivo
T1 image was first resampled by linear interpolation to an isotropicvoxel size of . . The 2D cell density maps were artificially extended to 3D by addition of a thickness of . and resliced to sagittal orientation. The density maps were finally rigidly registered to the corresponding slice ofthe resampled ex vivo T1 image based on user-defined landmark pairs using an in-house software in C++ based onVTK [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 imposethe brain slicing orientation. The complex histology slide to MR image registration process in 3D is thus reducedto a simple MR slice selection followed by the identification of at least 3 landmark pairs in-plane. Furthermore, thecomputed cell density maps have the great advantage of providing spatial tissue information at an intermediate scalebetween histological and radiological images, with a contrast similar to T1-weighted MR images. The cell densities ofwhite and gray matter are indeed substantially different, as is their T1 and PD values, which eased the identification oflandmarks pairs. 6
PREPRINT - F
EBRUARY
4, 2021
To verify the assumptions in Equations (5) and (6), the edema region has to be delineated in the reference ex vivo
T1image space. However, as previously mentioned, the drainage of the extracellular fluid made it impossible to discernvasogenic edema on the ex vivo
MR images. The in vivo
T2 FLAIR image was thus registered to the reference exvivo
T1 image. To this end, the in vivo
T1 image acquired on the same day was first registered on the ex vivo
T1image using rigid followed by B-spline transforms using the Elastix software [19]. The computed transforms were thensuccessively applied to the in vivo
T2 FLAIR image. The Elastix parameter files used for registration are available inAppendix C. The edema segmentation was finally performed semi-automatically on the registered T2 FLAIR imageusing a combination of thresholding and morphological operations.
To verify assumption in Equation (7), a 3D geodesic distance map to the tumor core across white matter was computedfrom the ex vivo
T1 image. White matter was first segmented using an in-house gradient-based anisotropic diffusionalgorithm followed by manual corrections to ensure that no physically incompatible bypass exists between white matterregions. The tumor core was then segmented on the same image using a combination of thresholding and morphologicaloperations. A distance map to the tumor core across the segmented white matter region was finally computed usingan adapted implementation of the anisotropic fast marching (AFM) algorithm presented in [20]. Note that since noDTI images were available for the patient, the anisotropy of glioma cell diffusion mentioned hereabove could not betaken into account in this work. A unit isotropic metric tensor field was thus provided to the AFM algorithm, hence theabusive 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 distanceand the average symmetric surface distance (ASSD), computed between the edema outlines and the contours of thebinary region obtained by thresholding the distance map. The Hausdorff distance d Hausdorff and the ASSD d ASSD betweentwo sets A and B are respectively given by [21]: d Hausdorff ( A, B ) = max (cid:26) max b ∈ B (cid:26) min a ∈ A d ( a, b ) (cid:27) , max a ∈ A (cid:26) min b ∈ B d ( a, b ) (cid:27)(cid:27) (9) d ASSD ( A, B ) = 1 | A | + | B | (cid:32)(cid:88) b ∈ B min a ∈ A d ( a, b ) + (cid:88) a ∈ A min b ∈ B d ( a, b ) (cid:33) (10)where d ( x, y ) is the Euclidian distance between elements x and y , and | X | is the cardinal of set X . As mentioned, the over-fixation of brain tissues prevented any IHC staining. Therefore, HE staining had to beused 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 distinctionbetween tumor and healthy cells. To address this problem, we propose to verify the following equation instead ofEquation (7): c total ( r ) = c tumor ( r ) + c white (11) = c core exp (cid:18) − d ( r ) λ white (cid:19) + c white (12)where c total , c tumor , and c white are respectively the total, tumor, and healthy cell density in white matter, and c core is thetumor cell density iso-value along the tumor core boundary. In this formulation, the cellularity of healthy white matteris supposed to be approximately constant and the invading tumor cells are assumed to be superimposed to the whitematter baseline cellularity c white . 7 PREPRINT - F
EBRUARY
4, 2021To relate the total cell density and the geodesic distance to the tumor core, the registered cell density maps wereresampled to the same voxel size as the geodesic distance map ( . × . × . ) and all available pairs ofdensity/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 inPython [22]. Figure 3 depicts an example of brain slice in its cutting guide (Figure 3(a)) with the corresponding registered in vivo
T2FLAIR image slice (Figure 3(b)), registered cell density maps (Figure 3(c)) and geodesic distance map (Figure 3(d)).An over-cellularity front is visible in Figure 3(c), progressing from the frontal necrotic tumor core but rapidly decreasingto reach a normal cellularity of around / mm beyond a geodesic distance of
20 mm (Figure 3(d)). Theedema outlines, on the other hand, extend to over
50 mm on the depicted image slice (see red and blue delineations inFigure 3(b)). The distinction between the red and blue segments of the edema outlines in Figure 3(b) will be used in thediscussion. (a) (b)
19 2118 2223 2420 (c) (d)
Figure 3: Cell density profile analysis. (a) Brain slice inside its 3D-printed cutting guide. (b) Corresponding sliceof the registered in vivo
T2 FLAIR image with segmented edema outlines. The blue and red segments of the outlinerespectively correspond to free and non-free to diffuse parts of the edema boundary (see Section 4). (c) Correspondingslice of the ex vivo
T1 image (grayscale) and superimposed registered cell density maps (colored) with their slidenumber (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 inFigure 4. The decreasing behavior of the tumor cell density with the distance to the tumor core was observed among allthese examples. 8
PREPRINT - F
EBRUARY
4, 2021
13 17
11 15
27 28
Figure 4: Example of registered cell density maps with their slide numbrt (see Table A1) (1 st and 3 rd columns) andcorresponding 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 surfacedensity data (Figure 5(a)) and the volume density data extrapolated using Equation (8) (Figure 5(b)). The fittedmodel curve given by Equation (12) is superimposed in red for each plot and the corresponding parameter values arerespectively 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 densitydata plotted in Figure 5(a). c core [ cell mm − ] λ white [ mm ] c white [ cell mm − ] c core [ cell mm − ] λ white [ mm ] c white [ cell mm − ] PREPRINT - F
EBRUARY
4, 2021
Distance ( mm ) C e ll d e n s it y ( ce ll mm − ) (a) Distance ( mm ) C e ll d e n s it y ( ce ll mm − ) (b) Figure 5: Scatter plot of the surface cell density (a) and the extrapolated volume cell density (b) versus distance foreach available value pairs among white matter voxels (blue dots) with superimposed fitted model curves (red curves).(blue). The distribution is rather continuous, suggesting that the edema boundary does not correspond to an iso-distancecontour in contrast to the expected step-like distribution superimposed in red in Figure 6. Five percent of the edemaboundary voxels are located at a distance greater or equal to . and only are located at a distance greaterthan . from the tumor core. . . . . Distance ( mm ) F r e qu e n c y Figure 6: Inverse cumulative distribution of the geodesic distance values along the edema outlines. The expecteddistribution under the hypothesis of iso-distance edema outlines is plotted in red.The threshold geodesic distance values that provided the smallest Hausdorff distance ( .
65 mm ) and ASSD ( .
97 mm )between the edema outlines and the contour of the corresponding thresholded region in the distance map were . and . , 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, reportingthe presence of pseudo-palisading necrosis, tumor cells (tumor block versus infiltrative cells) and glomeruloid vascularproliferation, along with the minimum, maximum and mean cell density and distance values within the correspondingslide. The furthest distance at which suspected infiltrating tumor cells were identified was around
46 mm (slide 7),whereas edema was detected as far as . on the same slide. The exponentially decreasing glioma cell density profile with distance to the tumor core suggested in [7] and [8] iscompatible with our experimental data, as observed in Figure 5 for both surface (a) and extrapolated volume (b) celldensities. The fitted value of λ white ( .
46 mm ) for the volume cell density data is in the same order of magnitude as theone used in [8] for simulation ( . ). The baseline volume cell density value of .
59 cell mm − in white matteris also in accordance with the literature, although large variations are to be noted between the reported values [23].10 PREPRINT - F
EBRUARY
4, 2021 (a) (b) (c)(d) (e) (f)
Figure 7: Edema region (red) with superimposed thresholded region of the distance map whose contour minimizes theHausdorff 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.However, high variance was observed in our experimental data (see blue points in Figure 5(a) and (b)), resulting inuncertainties on the fitted parameters. The fitted values of c core were likely underestimated, since cell density values ofup to . × cell mm − and . × cell mm − were respectively observed in our surface and volume cell densitydata (see Figure 5(a) and (b)), and consequently the corresponding fitted values of λ white were probably 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 distanceto the tumor core, iso-density contours should coincide with iso-distance contours. However, our results suggest thatedema outlines do not coincide with an iso-distance contour (see Figure 6) and would therefore not correspond to acell density iso-contour either. This apparent incompatibility of assumptions in Equation (5) and Equation (6) can beexplained by the thresholding behavior of the imaging function suggested in Equation (5) from which assumption inEquation (6) was deduced in [11]. In fact, thresholding a spatial function may give rise to iso-value contours only if thefunction is sufficiently smooth and continuous. In contrast, the tumor cell density function discretized over the voxelgrid has discontinuities at interfaces between white and gray matter and along the brain domain boundary. Indeed, thedifference in tumor cell diffusivity between white and gray matter [8, 10, 12] gives rise to steep tumor cell gradients atthe 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 diffuseoutside 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 imagingfunction in Equation (5) turned out to be verified. As an illustration, the edema outlines in Figure 3(b) were split intoblue 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 bereasonably assumed – taking a margin of for registration errors – that the edema extended up to . fromthe tumor core, which would correspond to the actual ‘free’ edema extent. Finally, it should be noted that due to therestricted 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 stillbe valid even after invalidation of Equation (6). Nevertheless, the thresholded distance map regions whose contour11 PREPRINT - F
EBRUARY
4, 2021respectively minimizes the Hausdorff distance and the ASSD to the edema outlines were not found to accuratelycoincide with the edema region, as depicted in Figure 7. Both thresholded distance map regions (blue) in Figure 7 wereindeed found to extend further in the contralateral hemisphere via the corpus callosum and less far towards the rightposterior region, compared to the edema region (red). Still under the hypothesis of a monotonically decreasing relationbetween the tumor cell density and the distance to the tumor core, the threshold-like imaging function in Equation (6) istherefore not compatible with our results. It should however be noted that this apparent inadequacy could result from theuse 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] wasnot supported by our experimental data either. Indeed, assuming that Equation (6) would still be verified on thefree-to-diffuse part of the edema boundary, the proposed cell density model in Equation (12) with the parameter valuesfitted to volume cell density data in Table 2 suggests an over-cellularity of only . × cell mm − with regard tothe white matter cellularity baseline c white at a distance of . to the tumor core corresponding to the free edemaextent. This over-cellularity corresponds to only .
08 % of c core ≤ c max , which strongly invalidates the commonlyaccepted value of
16 % . In addition, since the c core value was likely underestimated as mentioned hereabove, an evenlower actual value of c edema is to be expected. These results are confirmed by blinded pathological examination, whichdid not reveal any noticeable invasion of the brain parenchyma – even within the edematous region – beyond a distanceof
46 mm to the tumor core. This overall analysis does however not exclude the possible presence of isolated infiltratingtumor 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 examinationof gliomas, abnormalities visible on conventional MR sequences are not trivially related to the tumor cell invasionextent. A striking illustration of this limitation is the administration of corticosteroid or anti-angiogenic therapies toreduce edema-related symptoms in glioma patients, which does however not stop tumor progression. Consequently, adecoupling arises between tumor progression and its visible effects on the surrounding environment, potentially leadingto a misclassification of the disease as responding. The impact of such therapies on the MR-based follow-up of gliomashas been extensively studied through numerical simulations in [3]. In this work, we invalidated two commonly madeassumptions relating the outlines of visible abnormalities on MRI to the tumor cell density function: (i) Assuming athreshold-like edema imaging function (see Equation (5)), the edema contour may not correspond to an iso-contour ofthe cell density function as soon as the migration of tumor cells is locally restricted or prevented and (ii) at a distancecorresponding to the maximum extension of the vasogenic edema, the over-cellularity was found to be negligible inour studied case, as opposed to the previously hypothesized value of
16 % . These results raise the question of theapplicability of both previously proposed methods to assess the glioma cell density distribution from routine MR imagespresented in the introduction. Indeed, whereas the method proposed in [12] allows to compute an accurate approximatesolution of Equation (1), it still relies on the assumption that iso-density contours can be derived from MR data. Inthe case of a more simple exponentially decreasing model as in Equation (7) [7, 8], iso-density contours would stillbe required to assess the infiltration length parameter λ white . Since a high sensibility of the reaction-diffusion tumorgrowth models to their initial condition was previously reported by our group [13], deriving an initial spatial cell densitydistribution that is as reliable as possible is crucial for the model to be applied in clinical practice. To this end, the useof other imaging sequences or modalities such as average diffusion coefficient (ADC) maps [27] or amino-acid positronemission 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 analyzed –a non-operated brain with GBM – this study was based on a single case and should be further conducted on a largerdiffuse glioma cohort of various grades. In addition, the use of murine glioma models for conducting such studies at alarger scale would be of interest but is restricted due to reported substantial differences in cortical [28], glial [29], andendothelial [30] cells between human and mouse brain, leading to only partial capture of human glioma features bysuch models [31]. Second, IHC staining could not be performed on the autopsied material, which has prevented thespecific 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 thatthe over-cellularity observed locally is exclusively attributed to tumor cell invasion since tumor-induced recruitment ofinflammatory cells is limited in brain tissues. As a consequence, the identification of isolated infiltrating tumor cellson pathological examination may have been prevented. Third, the substantial deformation of the brain between invivo and ex vivo imaging – with a volume decrease estimated to
16 % ex vivo based on MRI – may have resulted inpartial distortion of the ex vivo
MR-derived distance map compared to in vivo , not fully compensated by deformableregistration of the in vivo edema outlines.We would finally like to emphasize that the automated cell density map computation and histological slide to MR imageregistration procedures described in this work are not limited to the problem addressed herein and could be appliedto various histological stainings, imaging modalities, and organs, such as prostate. Besides, the use of tailor-made12
PREPRINT - F
EBRUARY
4, 20213D-printed slicer and cutting guides makes it possible to precisely analyze whole organ slices at low cost even forcenters that do not have access to whole organ slice microscopy, opening tremendous possibilities for translationalmicroscopic/macroscopic imaging research [32].
Through a translational radiological/histological analysis performed on a case of non-operated glioblastoma, weinvalidated two commonly made assumptions relating the outlines of the visible abnormalities in magnetic resonanceimages to the tumor cell density function in the context of reaction-diffusion glioma growth modeling. We showedthat, due to local restrictions of the tumor cell migration at brain tissue interfaces and along the brain boundary, theoutlines of vasogenic edema in T2 FLAIR images do not generally coincide with a cell density iso-contour, as opposedto what was previously suggested. We also found that the commonly adopted tumor cell density iso-value at the edemaoutlines is likely overestimated since the over-cellularity at the maximum edema extent was found to be negligible inour studied case. This, however, does not exclude the possible presence of isolated tumor cells migrating beyond edemaoutlines, as previously reported. This work highlights the limitations of using routine magnetic resonance images toderive cell density maps for reaction-diffusion tumor growth models and points out the need of validating other methodsto accurately initialize such models and make them usable for clinical applications.
Acknowledgments
C. Martens is funded by the FRIA grant no.5120417F (F.R.S.-FNRS – Belgian National Fund for Scientific Research).L. Lebrun is funded by the Fonds Erasme. C. Decaestecker is senior research associate with the F.R.S.-FNRS. TheDepartment of Nuclear Medicine at Hôpital Erasme is supported by the Association Vinçotte Nuclear (AVN), theFonds Erasme and the Walloon Region (Biowin). The Department of Pathology at Hôpital Erasme is supported by theFonds Erasme and the Fonds Yvonne Boël. The CMMI is supported by the European Regional Development Fund andthe Walloon Region. The authors would like to thank the technical, medical, and scientific staff in charge of imageacquisition and tissue sample processing for the needs of this work.
References [1] Q. T. Ostrom, G. Cioffi, H. Gittleman, N. Patil, K. Waite, C. Kruchko, and J. S. Barnholtz-Sloan. CBTRUSstatistical report: Primary brain and other central nervous system tumors diagnosed in the United States in2012–2016.
Neuro-Oncol. , 21(S5):1–100, Nov 2019.[2] D. L. Silbergeld and M. L Chicoine. Isolation and characterization of human malignant glioma cells fromhistologically normal brain.
J. Neurosurg. , 86(3):525–31, Mar 1997.[3] A. Hawkins-Daarud, R. C. Rockne, A. R. A. Anderson, and K. R. Swanson. Modeling tumor-associated edema ingliomas during anti-angiogenic therapy and its impact on imageable tumor.
Front. Oncol. , 3(66), Apr 2013.[4] Z.-X. Lin. Glioma-related edema: New insight into molecular mechanisms and their clinical implications.
Chin. J.Cancer , 32(1):49–52, Jan 2013.[5] S. Lu, D. Ahn, G. Johnson, M. Law, D. Zagzag, and R. I. Grossman. Diffusion-tensor MR imaging of intracranialneoplasia and associated peritumoral edema: Introduction of the tumor infiltration index.
Radiology , 232(1):221–8,Jul 2004.[6] P. Wesseling, J. M. Kros, and J. W. M. Jeuken. The pathological diagnosis of diffuse gliomas: Towards asmart synthesis of microscopic and molecular information in a multidisciplinary context.
Diagn. Histopathol. ,17(11):486–94, Nov 2011.[7] P. Tracqui, G. C. Cruywagen, D. E. Woodward, G. T. Bartoo, J. D. Murray, and E. C. Alvord. A mathematicalmodel of glioma growth: The effect of chemotherapy on spatio-temporal growth.
Cell Prolif. , 28(1):17–31, Jan1995.[8] J. Unkelbach, B. H. Menze, E. Konukoglu, F. Dittmann, M. Le, N. Ayache, and H. A. Shih. Radiotherapyplanning for glioblastoma based on a tumor growth model: Improving target volume delineation.
Phys. Med. Biol. ,59(3):747–70, Feb 2014.[9] S. Jbabdi, E. Mandonnet, H. Duffau, L. Capelle, K. Rae Swanson, M. Pélégrini-Issac, R. Guillevin, and H. Benali.Simulation of anisotropic growth of low-grade gliomas using diffusion tensor imaging.
Magn. Reson. Med. ,54(3):616–24, Sep 2005. 13
PREPRINT - F
EBRUARY
4, 2021[10] O. Clatz, M. Sermesant, P.-Y. Bondiau, H. Delingette, S. K. Warfield, G. Malandain, and N. Ayache. Realisticsimulation of the 3-D growth of brain tumors in MR images coupling diffusion with biomechanical deformation.
IEEE Trans. Med. Imag. , 24(10):1334–46, Oct 2005.[11] K. R. Swanson, R. C. Rostomily, and E. C. Alvord. A mathematical modelling tool for predicting survival ofindividual patients following resection of glioblastoma: A proof of principle.
Br. J. Cancer , 98(1):113–9, Jan2008.[12] E. Konukoglu, O. Clatz, P.-Y. Bondiau, H. Delingette, and N. Ayache. Extrapolating glioma invasion margin inbrain magnetic resonance images: Suggesting new irradiation margins.
Med. Image Anal. , 14(2):111–25, Apr2010.[13] C. Martens, T. Metens, O. Debeir, S. Goldman, and G. Van Simaeys. Initial condition assessment from patientMRI data for reaction-diffusion glioma growth models. In
Proc. Intl. Soc. Mag. Reson. Med. 27 , 2019. 2852.[14] M. Tovi and A. Ericsson. Measurement of T1 and T2 over time in formalin-fixed human whole-brain specimens.
Acta Radiol. , 33(5):400–4, Sep 1992.[15] M. R. Raman, Y. Shu, T. G. Lesnick, C. R. Jack, and K. Kantarci. Regional T1 relaxation time constants in exvivo human brain: Longitudinal effects of formalin exposure.
Magn. Reson. Med. , 77(2):774–8, Feb 2017.[16] M. Absinta, G. Nair, M. Filippi, A. Ray-Chaudhury, M. I. Reyes-Mantilla, C. A. Pardo, and D. S. Reich.Postmortem magnetic resonance imaging to guide the pathologic cut: Individualized, 3-dimensionally printedcutting boxes for fixed brains.
J. Neuropathol. Exp. Neurol. , 73(8):780–8, Aug 2014.[17] W. Schroeder, K. Martin, and B. Lorensen.
The Visualization Toolkit . Kitware, Clifton Park, NY, USA, 4th edition,2010.[18] T. Yoo, M. Ackerman, W. Lorensen, W. Schroeder, V. Chalana, S. Aylward, D. Metaxas, and R. Whitaker.Engineering and algorithm design for an image processing API: A technical report on ITK – the Insight Toolkit.
Stud. Health Technol. Inform. , 85:586–92, Feb 2002.[19] S. Klein, M. Staring, K. Murphy, M. A. Viergever, and J. Pluim. elastix: A toolbox for intensity-based medicalimage registration.
IEEE Trans. Med. Imaging , 29(1):196–205, Jan 2010.[20] E. Konukoglu, M. Sermesant, O. Clatz, J.-M. Peyrat, H. Delingette, and N. Ayache. A recursive anisotropic fastmarching approach to reaction diffusion equation: Application to tumor growth modeling. In
Inf. Process. Med.Imaging , pages 687–99, Berlin, Heidelberg, Germany, 2007. Springer.[21] V. Yeghiazaryan and I. Voiculescu. Family of boundary overlap metrics for the evaluation of medical imagesegmentation.
J. Med. Imaging , 5(1):015006, Jan 2018.[22] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson,W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. Jarrod Millman, N. Mayorov, A. R. J. Nelson,E. Jones, R. Kern, E. Larson, C. J. Carey, ˙I. Polat, Y. Feng, E. W. Moore, J. Vanderplas, D. Laxalde, J. Perktold,R. Cimrman, I. Henriksen, E. A. Quintero, C. R Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, and P. vanMulbregt. SciPy 1.0: Fundamental algorithms for scientific computing in Python.
Nat. Methods , 17(S1):261–72,Feb 2020.[23] C. S. von Bartheld, J. Bahney, and S. Herculano-Houzel. The search for true numbers of neurons and glial cells inthe human brain: A review of 150 years of cell counting.
J. Comp. Neurol. , 524(18):3865–95, Dec 2016.[24] S. Jbabdi, P. Bellec, R. Toro, J. Daunizeau, M. Pélégrini-Issac, and H. Benali. Accurate anisotropic fast marchingfor diffusion-based geodesic tractography.
Int. J. Biomed. Imaging , 2008:320195, Dec 2007.[25] P. J. Kelly, C. Daumas-Duport, D. B. Kispert, B. A. Kall, B. W. Scheithauer, and J. J. Illig. Imaging-basedstereotaxic serial biopsies in untreated intracranial glial neoplasms.
J. Neurosurg. , 66(6):865–74, Jun 1987.[26] F. Sahm, D. Capper, A. Jeibmann, A. Habel, W. Paulus, D. Troost, and A. von Deimling. Addressing diffuseglioma as a systemic brain disease with single-cell analysis.
Arch. Neurol. , 69(4):523–6, Apr 2012.[27] N. C. Atuegwu, D. C. Colvin, M. E. Loveless, L. Xu, J. C. Gore, and T. E. Yankeelov. Incorporation of diffusion-weighted magnetic resonance imaging data into a simple mathematical model of tumor growth.
Phys. Med. Biol. ,57(1):225–240, Jan 2012.[28] R. D. Hodge, T. E. Bakken, J. A. Miller, K. A. Smith, E. R. Barkan, L. T. Graybuck, J. L. Close, B. Long,N. Johansen, O. Penn, and Z. Yao. Conserved cell types with divergent features in human versus mouse cortex.
Nature , 573:61–8, Sep 2019.[29] Y. Zhang, S. A. Sloan, L. E. Clarke, C. Caneda, C. A. Plaza, P. D. Blumenthal, H. Vogel, G. K. Steinberg, M. S. B.Edwards, G. Li, and J. A. III Duncan. Purification and characterization of progenitor and mature human astrocytesreveals transcriptional and functional differences with mouse.
Neuron , 89(1):37–53, Jan 2016.14
PREPRINT - F
EBRUARY
4, 2021[30] H. W. Song, K. L. Foreman, B. D. Gastfriend, J. S. Kuo, S. P. Palecek, and E. V. Shusta. Transcriptomiccomparison of human and mouse brain microvessels.
Sci. Rep. , 10(1), Dec 2020.[31] S. de Vleeschouwer, editor.
Glioblastoma , pages 131–40. Codon Publications, Brisbane, Australia, 2017.[32] D. Baldi, M. Aiello, A. Duggento, M. Salvatore, and C. Cavaliere. MR imaging-histology correlation by tailored3D-printed slicer in oncological assessment.
Contrast Media Mol. Imaging , 2019:1071453, May 2019.[33] N. Otsu. A threshold selection method from gray-level histograms.
IEEE Trans. Syst. Man Cybern. Syst. , 9(1):62–6,Jan 1979.[34] O. Ronneberger, P. Fischer, and T. Brox. U-Net: Convolutional networks for biomedical image segmentation. In
Med. Image Comput. Comput. Assist. Interv. arXiv:1412.6980 , Dec 2014.[37] F. Xing, T. C. Cornish, T. Bennett, D. Ghosh, and L. Yang. Pixel-to-pixel learning with weak supervision forsingle-stage nucleus recognition in Ki67 images.
IEEE Trans. Biomed. Eng. , 66(11):3088–97, Nov 2019.
A Slicer design
The in vivo
T2 FLAIR image was first resampled to an isotropic voxel size of . × . × . . A brainmask was generated using the Otsu thresholding [33], followed by a morphological opening of radius 1, a largestcomponent extraction and a morphological closing of radius 15. The brain mask bounding box was then drawn ona binary image with the same spacing and was evenly extended on both sides in each spatial dimension, with finaldimensions of . × . ×
172 mm . The brain mask volume was subtracted from the box volume. Tenslits of
12 mm × . ×
152 mm spaced by . were carved into the box along the sagittal plane for brainslicing 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 . and height . were added to each corner of the upper part of the slicerand subtracted from its lower part. A 3D surface mesh was finally generated from the binary image using the marchingcubes 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 . × . ×
11 mm plates from which the slice volume with a thickness of . was subtracted. Grooves werecarved into the plates to guide the scalpel-cutting of
26 mm ×
18 mm × . tissue blocks whose dimensionsare compatible with the histological processing chain at our institution. All image processing and design steps wereperformed using an in-house software in C++ based on VTK [17] and ITK [18].The brain slicer was printed using .
75 mm
PLA thread on a HICTOP 3DD-17-ATL-FM printer (HIC Technology,China) with a nozzle size of . (layer thickness . , honeycomb filling ) for a total printing time of 96hours. The printing of the ten cutting guides was parallelized on 5 Prusa i3 MK2 printers (Prusa Research, CzechRepublic) 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 µ m × µ m , 50 tiles of
100 px ×
100 px were randomlychosen after exclusion of background tiles. Background tiles were defined as tiles containing less than
10 % of tissuepixels, empirically chosen as pixels whose the lowest RGB value is above 220 for our calibrated scanner. The 1600-tiledataset was then weakly supervised using an in-house annotation program in Python allowing the user to locate eachcell nucleus within a tile by a simple mouse click. For each annotated tile, a nuclei presence probability map wasgenerated by the superposition of 2D Gaussian blobs centered on user-pointed nucleus locations with full width at halfmaximum (FWHM) of (i.e. standard deviation σ ≈ .
27 px ), truncated to a radius of σ . The supervised datasetwas 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 down-sampling blocks, 3 up-sampling blocks and an output block. Each down-sampling block is made of two convolutionallayers with kernel size × and stride 1 followed by a bias-adding layer and a rectified linear unit (ReLU) activationlayer. A max pool layer with kernel size × and stride 2 is added at the end of the block to reduce the feature maps15 PREPRINT - F
EBRUARY
4, 2021 (a) (b)(c) (d)
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.dimensions by a factor 2. The up-sampling blocks are identical to the down-sampling blocks except that the max poollayer is replaced by a deconvolution layer with kernel size × and stride 2 followed by a bias-adding layer and aReLU activation layer to expand the feature maps dimensions by a factor 2. The output block has the same structure asthe down-sampling blocks except that the max pool layer is replaced by a convolutional layer with kernel size × andstride 1 followed by a bias-adding layer and a sigmoid activation layer to merge the last 128 feature maps into a singlepresence probability map. The network architecture is depicted in Figure A2 and was implemented in Python using theTensorFlow 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: − , β : . , β : . , (cid:15) : − ). An evaluation was performed on the validation set at the end of eachepoch and the training was stopped early after no improvement of the validation metric in 100 epochs. The networkparameter values that gave the best validation metric value ( . ) were kept, which was reached at epoch 96. A similarweakly supervised deep-learning approach was used in [37] for cell nuclei detection in histological slides, but Euclidiandistance to user-pointed nuclei locations and the mean squared error loss were used instead of Gaussian nuclei presenceprobability maps and the cross-entropy loss.After training, each resampled scanned slide was divided into adjacent tiles of
100 px ×
100 px from which nucleipresence probability maps were predicted by the trained network. A nuclei count was finally derived for each predictedprobability map by detecting its local maxima with a value greater or equal to . nd spaced by at least .16 PREPRINT - F
EBRUARY
4, 2021
Conv 3x3 stride 1 + bias + ReLUMax pool 2x2 stride 2Deconv 2x2 stride 2 + bias +ReLUConv 1x1 stride 1 + bias + sigmoidConcatenation25x25x512 100x100x128 100x100x1Presenceprobability map100x100x128100x100x3RGB tile 50x50x256 50x50x256Down-samplingblock Up-samplingblock
Figure A2: U-Net architecture with its feature map dimensions used for nuclei detection. The network takes
100 px ×
100 px
RGB tiles with pixel size µ m × µ m as input and predicts a nuclei presence probability map with the samespatial dimensions. C In vivo to ex vivo registration Registration of the in vivo
T2-FLAIR image to the ex vivo
T1 image was performed using the Elastix software (version4.801) [19] as described in Section 2.7. The Elastix parameter files for rigid and B-spline registration are respectivelygiven in Listings A1 and A2Listing A1: Elastix parameter file for rigid in vivo
T1 to ex vivo
T1 registration. // ElastixMain ( FixedImageDimension 3) ( MovingImageDimension 3) ( FixedInternalImagePixelType " double ") ( MovingInternalImagePixelType " double ") // ElastixTemplate ( WriteTransformParametersEachIteration " false ") ( WriteTransformParametersEachResolution " false ") ( UseDirectionCosines " true ") // Registration ( Registration " MultiResolutionRegistration ") ( NumberOfResolutions 1) ( ErodeMask " false ") ( ErodeFixedMask " false ") ( ErodeMovingMask " false ") // Pyramid ( FixedImagePyramid " FixedSmoothingImagePyramid ") ( MovingImagePyramid " MovingSmoothingImagePyramid ") ( ImagePyramidSchedule 1.0 1.0 1.0) ( WritePyramidImagesAfterEachResolution " false ") // Transform ( Transform " EulerTransform ") PREPRINT - F
EBRUARY
4, 2021 ( AutomaticScalesEstimation " true ") ( AutomaticTransformInitialization " true ") ( AutomaticTransformInitializationMethod " GeometricalCenter ") // Metric ( Metric " AdvancedMattesMutualInformation ") ( NumberOfHistogramBins 50) ( NumberOfFixedHistogramBins 50) ( NumberOfMovingHistogramBins 50) ( FixedKernelBSplineOrder 3) ( MovingKernelBSplineOrder 3) ( FixedLimitRangeRatio 0.01) ( MovingLimitRangeRatio 0.01) ( FiniteDifferenceDerivative " false ") ( UseFastAndLowMemoryVersion " true ") ( UseJacobianPreconditioning " false ") ( ShowExactMetricValue " true ") ( ExactMetricSampleGridSpacing 1 1 1) ( CheckNumberOfSamples " false ") ( RequiredRatioOfValidSamples 0.25) // Optimizer ( Optimizer " RegularStepGradientDescent ") ( MaximumNumberOfIterations 500) ( MinimumStepLength 0.00001) ( MaximumStepLength 8.0) ( MinimumGradientMagnitude 0.00001) ( RelaxationFactor 0.8) ( NewSamplesEveryIteration " false ") // ImageSampler ( ImageSampler " Full ") // Interpolator ( Interpolator " LinearInterpolator ") // ResampleInterpolator ( ResampleInterpolator " FinalLinearInterpolator ") // Resampler ( Resampler " DefaultResampler ") ( WriteResultImage " true ") ( WriteResultImageAfterEachResolution " false ") ( WriteResultImageAfterEachIteration " false ") ( ResultImageFormat " mha ") ( ResultImagePixelType " double ") ( CompressResultImage " false ") // Misc ( UseMultiThreadingForMetrics " true ") ( DefaultPixelValue 0.0) Listing A2: Elastix parameter file for B-Spline in vivo
T1 to ex vivo
T1 registration. // ElastixMain ( FixedImageDimension 3) ( MovingImageDimension 3) ( FixedInternalImagePixelType " double ") ( MovingInternalImagePixelType " double ") // ElastixTemplate ( WriteTransformParametersEachIteration " false ") ( WriteTransformParametersEachResolution " false ") ( UseDirectionCosines " true ") PREPRINT - F
EBRUARY
4, 2021 // Registration ( Registration " MultiResolutionRegistration ") ( NumberOfResolutions 5) ( ErodeMask " false ") ( ErodeFixedMask " false ") ( ErodeMovingMask " false ") // Pyramid ( FixedImagePyramid " FixedSmoothingImagePyramid ") ( MovingImagePyramid " MovingSmoothingImagePyramid ") ( ImagePyramidSchedule 16.0 16.0 16.0 8.0 8.0 8.0 4.0 4.0 4.0 2.0 2.0 2.0 1.0 1.01.0) ( WritePyramidImagesAfterEachResolution " false ") // Transform ( Transform " BSplineTransform ") ( BSplineTransformSplineOrder 3) ( FinalGridSpacingInPhysicalUnits 4.0 4.0 4.0) ( GridSpacingSchedule 16.0 16.0 16.0 8.0 8.0 8.0 4.0 4.0 4.0 2.0 2.0 2.0 1.0 1.01.0) ( PassiveEdgeWidth 0) ( UseCyclicTransform " false ") ( HowToCombineTransforms " Compose ") // Metric ( Metric " AdvancedMattesMutualInformation ") ( NumberOfHistogramBins 32) ( NumberOfFixedHistogramBins 32) ( NumberOfMovingHistogramBins 32) ( FixedKernelBSplineOrder 3) ( MovingKernelBSplineOrder 3) ( FixedLimitRangeRatio 0.001) ( MovingLimitRangeRatio 0.001) ( FiniteDifferenceDerivative " false ") ( UseFastAndLowMemoryVersion " true ") ( UseJacobianPreconditioning " false ") ( ShowExactMetricValue " false ") ( ExactMetricSampleGridSpacing 1) ( CheckNumberOfSamples " true ") ( RequiredRatioOfValidSamples 0.25) // Optimizer ( Optimizer " StandardGradientDescent ") ( MaximumNumberOfIterations 12000) ( MaximumNumberOfSamplingAttempts 10) ( SP_a 2500) ( SP_A 1200.0) ( SP_alpha 0.6) ( NewSamplesEveryIteration " true ") // ImageSampler ( ImageSampler " RandomCoordinate ") ( NumberOfSpatialSamples 2000) ( UseRandomSampleRegion " true ") ( SampleRegionSize 100.0 100.0 100.0) ( FixedImageBSplineInterpolationOrder 1) // Interpolator ( Interpolator " BSplineInterpolator ") ( BSplineInterpolationOrder 1) // ResampleInterpolator ( ResampleInterpolator " FinalBSplineInterpolator ") ( FinalBSplineInterpolationOrder 3) PREPRINT - F
EBRUARY
4, 2021 // Resampler ( Resampler " DefaultResampler ") ( WriteResultImage " true ") ( WriteResultImageAfterEachResolution " false ") ( WriteResultImageAfterEachIteration " false ") ( ResultImageFormat " mha ") ( ResultImagePixelType " double ") ( CompressResultImage " false ") // Misc ( UseMultiThreadingForMetrics " true ") ( DefaultPixelValue 0.0) D Pathology results
The results of the pathological examination and numerical tile processing are summarized in Table A1.Table A1: Results of the pathological examination and numerical tile processing. PPN: pseudo-palisading necrosis,GVP: glomeruloid vascular proliferation, susp.: suspected.
Cell density Distance [ mm ][ cell mm − ]Block PPN Tumor cells GVP Edema Min Max Mean Min Max Mean1 No No No No 1.08 20.46 11.77 33.40 40.40 34.952 No No No No 2.25 23.88 13.80 35.05 46.76 40.123 No Infiltrative (susp.) No Yes 1.52 25.25 14.96 7.24 21.18 14.004 No No No No 1.36 20.64 12.50 29.36 41.66 36.845 No No No Yes 1.19 24.85 15.40 11.74 30.27 21.576 No Infiltrative (susp.) No No 3.72 28.44 18.70 32.92 40.35 36.527 No Infiltrative (susp.) No Yes 0.89 38.95 21.44 38.59 61.70 49.048 Yes Block Yes No 3.37 37.29 23.65 0.50 5.17 2.719 Yes Infiltrative Yes Yes 1.02 44.37 31.71 0.50 4.72 1.9110 No Infiltrative (susp.) No Yes 5.52 37.15 25.53 0.50 3.71 1.0611 No No No No 1.73 25.68 15.21 19.85 38.74 30.7412 No Infiltrative Yes Yes 0.89 36.08 19.70 0.50 31.46 14.2413 Yes Block Yes No 1.40 52.46 21.69 0.50 22.85 9.7214 No No No No 1.10 19.25 11.78 17.73 32.20 25.5115 No No No Yes 2.08 22.90 11.41 28.63 54.74 40.1416 No No No No 0.99 26.47 13.52 20.37 48.79 35.9117 No No No No 1.42 23.69 13.10 28.30 56.12 42.8818 No No Yes Yes 2.46 39.29 20.71 7.94 27.01 16.8319 Yes Block Yes Yes 0.99 44.12 19.82 0.50 18.50 6.9220 Yes Block Yes No 6.25 61.05 28.40 0.50 7.37 2.0621 No Infiltrative Yes Yes 4.75 34.20 23.05 0.50 14.62 7.7022 No No No Yes 0.86 30.37 10.79 2.99 36.87 21.3923 No No No No 0.94 25.61 13.23 21.14 49.61 34.4824 No No No No 1.37 26.72 13.72 57.27 90.93 71.0825 No Infiltrative (susp.) No Yes 0.87 39.07 21.03 1.71 21.41 10.8726 Yes Block Yes Yes 6.21 54.24 22.79 0.50 8.81 2.9627 No No No Yes 0.93 37.96 21.33 12.74 30.35 18.4028 No No No No 0.95 34.79 20.76 14.35 35.85 22.57