Forty-four patients suffering from de novo brain tumours have been recruited and acquired at the University Hospital of Padova from July 2017 to March 2021. All procedures were in accordance with the ethical standards of the institutional research committee and with the 1964 Helsinki declaration plus later amendments. All participants provided informed, written consent in accordance with the local University Hospital Institutional Review Board.
Figure 1 shows a comprehensive overview of the analyses introduced in this section.
2.1 MRI Acquisition
Data acquisition was performed with a 3T Siemens Biograph mMR-PET/MR scanner equipped with a 16-channels head-neck coil. The multi-shell dMRI protocol comprised a total of 100 diffusion weighted images (DWIs) (TR/TE 5355/104 ms; voxel size 2x2x2 mm3; FOV 220x220 mm2; 68 slices; multiband accelerator factor=2): 10 images at b=0 s/mm2, 30 DWIs at b-value=710 s/mm2 and 60 DWIs at b-value=2855 s/mm2. This diffusion HARDI protocol is the optimized two shells NODDI protocol as described in (Zhang et al. 2012). Each diffusion direction was acquired with reverse phase encoding directions, i.e., anterior-posterior and posterior-anterior directions, for distortion correction purposes.
In addition, the acquisition protocol included anatomical imaging, which comprised a 3D T2-weighted (T2w) Fluid Attenuated Inversion Recovery (FLAIR) image (TR/TE 5000/395 ms; voxel size 1x1x1 mm3; FOV 250x250 mm2), two 3D T1-weighted (T1w) magnetization-prepared rapid acquisition gradient echo (MPRAGE, TR/TE 2400/3.2 ms; voxel size 1x1x1 mm3; FOV 256x256 mm2; 160 slices) acquired both before and after contrast agent injection and a T2w image (TR/TE 3200/536 ms; voxel size 1x1x1 mm3; FOV 256x256 mm2; 160 slices).
2.2 Tumour Segmentation and Structural Pre-processing
The anatomical images of each patient were linearly registered to the patient naïve T1w image with the Advanced Normalization Tools (ANTs (Avants et al. 2011), v. 2.0.1). Employing these images, two masks were manually delineated through the ITK-SNAP software (http://www.itksnap.org/) by an expert neuroradiologist with more than five years of experience. The first mask, the T, included the tumour core (contrast agent enhancing and non-enhancing regions) and the necrosis, where present. The second mask, the T+O, was created by adding the oedema area to the T mask. In addition, each tumour was labelled by the same neuroradiologist as left, right or bilateral according to the location of its core and to the mainly involved hemisphere.
Structural pre-processing was applied to the T1w image of each patient and consisted in bias field correction (N4BiasFieldCorrection (Tustison et al. 2010)), skull-stripping (Multi-Atlas Skull Stripping (Doshi et al. 2013)), tissue segmentation (into GM, WM and cortico-spinal fluid, with the unified segmentation tool (Ashburner and Friston 2005) of SPM12 v. 7771) and diffeomorphic non-linear registration (as implemented in ANTs SyN algorithm) to the symmetric MNI152 atlas. The last step was performed excluding the T+O area from the computation of the cost function as suggested in (Andersen, Rapcsak, and Beeson 2010).
For each patient, the T and T+O masks were then mapped into the MNI152 exploiting the estimated diffeomorphic non-linear transformations.
Finally, since direct and indirect disconnection mapping methods differ in the tracking of potential disconnections within subcortical areas, for each patient, we used the AAL3 atlas (Rolls et al. 2020) to segment the following regions: thalamus, caudates, putamen, pallidum, and hippocampus. Such areas were discarded not to bias the subsequent methods’ comparison.
2.3 Disconnection Maps Computation
2.3.1 Direct Disconnection Map
Diffusion MRI Pre-processing: The acquired diffusion weighted volumes were visually inspected to identify and remove those images affected by interslice instabilities (J. D. Tournier, Mori, and Leemans 2011) which were deemed excessively corrupted for subsequent pre-processing techniques to correct. The rest of the pre-processing was executed in its entirety within the MRtrix3 Software (J.-D. Tournier et al. 2019) and was featured by an initial denoising step based on random matrix theory (Veraart, Fieremans, and Novikov 2016), and a subsequent call to the tools TOPUP (Andersson, Skare, and Ashburner 2003) and eddy (Andersson and Sotiropoulos 2016) from the FMRIB Software library (FSL) for B0 inhomogeneity, eddy current and motion joint correction.
T1w segmentation results (including GM, subcortical parcellation, lesion and tumour masks) were registered to the naïve B0 volume using ANTs, by applying an affine transformation previously estimated on the patient’s naïve T1w image.
Diffusion tractography specifications: The patient structural connectome reconstruction was performed in its entirety within the MRtrix3 software. We firstly performed multi-shell multi-tissue spherical deconvolution (Jeurissen et al. 2014) to recover the orientation distribution functions for each voxel. Subsequently, we computed the structural connectome by employing Anatomically Constrained Tractography (R. E. Smith et al. 2012), tracking individual fibres with a second-order Integration over Fiber Orientation Distributions algorithm (J.-D. Tournier, Calamante, and Connelly 2010). Standard streamline termination criteria values were used. The number of generated streamlines for each patient initially amounted to 100M, which were quantitatively reduced to 10M via the Spherical-deconvolution Informed Filtering of Tractograms framework (R. Smith et al. 2013).
Disconnection Maps Computation: For each patient, following the quantification of the diffusion tractogram, we computed the patterns of structural disconnections by taking the two following steps:
- We computed the subset of the streamlines in the tractogram that featured an overlap with the tumoral lesion (respectively, with T and, after, with T+O masks).
- In the dMRI space (i.e., the naïve space of the tractogram), we computed how many altered streamlines (the subset in step 1) were passing through each voxel of the brain. We labelled direct structural disconnection (dSD) maps these voxel-wise frequency maps.
We used in-house MATLAB (ver. 2020b, The Mathworks, Natick, MA) scripts to perform the creation of the dSD maps.
As further volumetric analyses of dSD maps required their binarization, the definition of a threshold to define significant disconnection was necessary. We defined such threshold with the following procedure:
- Similarly to the dSD map quantification, we computed for each subject in the dMRI space a voxel-wise streamline density map, this time considering the entirety of the tractogram.
- We brought the individual streamline density maps to the MNI152 space via the previously estimated diffeomorphic transformations.
- In the common MNI152 space, we computed the population average of the individual maps, omitting the lesioned ROI on a patient-by patient basis.
- To balance the presence of a higher number of patients with left tumours, we symmetrised the frequency template by flipping (left-right) the obtained map, summing the flipped and non-flipped maps and dividing by 2.
- Using the same transform as in 2, we projected back the population-averaged streamline density map to the T1w space of each subject. We referred this resulting map as AvgDensity.
- Finally, we considered the structural disconnection value of voxel with coordinates (x,y,z) as significant if the following criterion was met:
dSD(x,y,z)/AvgDensity(x,y,z) > 10%
Additional investigations were performed also using more stringent and permissive thresholds (i.e., 5%,15%,20%,25%).
Hence for each patient we obtained a mask of significantly disconnected voxels.
2.3.2 Indirect Disconnection Maps Computation
Indirect structural disconnection (iSD) maps were quantified with the BCB Toolkit v. 4.2 (Foulon et al. 2018; Thiebaut de Schotten, Foulon, and Nachev 2020). In choosing which healthy controls tractography atlas to use within the toolbox, we opted for the extended diffusion dataset provided by the toolkit authors, for which the structural connectome of 180 healthy from the Human Connectome Project’ 7T controls was quantified (full tractography specification in (Thiebaut de Schotten et al. 2017)).
In summary, using the BCB Toolkit, for each patient, the lesion masks (both T and T+O) in the MNI152 space were registered to each control naïve space using affine and diffeomorphic deformations and subsequently used as seed for the tractography in TrackVis (http://trackvis.org). Tractographies from the lesions were transformed in visitation maps, binarized and brought to the MNI152. Finally, the percentage overlap map was computed by summing at each point in the MNI152 space the normalized visitation map of each healthy subject. Hence, in the resulting disconnectome map (i.e., iSD), the value in each voxel considers the tracts’ interindividual variability and indicate a probability of disconnection from 0 to 100% for the given lesion.
As no ad-hoc studies are available regarding the recommended use of iSD maps in tumours, we set the probabilistic threshold to 0.5 as in the software defaults. Additional investigations were performed also using more stringent and permissive thresholds (i.e., 0.7 and 0.3 respectively).
2.4 Metrics of Comparison
Having obtained a total of four different SD maps (i.e., two for each methodology, different in terms of the employed input mask: T and T+O), we wanted to compare them both intra-methodology (i.e., same approach, T versus T+O maps) and inter-methodology (i.e., same lesion mask, direct versus indirect maps).
The framework we employed for their comparison had three simple metrics, useful to quantify different similarity features:
- The difference in volume (ΔVol): This metric allowed us to quantify the difference in the extension of alteration which is detected by two disconnection maps.
- The Dice Similarity index (Dice): The Sørensen-Dice similarity index is a well-known metric of comparison between digital images and is defined by the following formula:
where A, B are the two binary matrices for which the similarity needs to be tested. The Dice index quantifies how similar the shape of structural alterations is between two different approaches.
- The Correlation at the intersection (Corr): We computed the Pearson correlation in a region of interest defined by the intersection of the two matrices in exam. Unlike the Dice index which evaluates only the similarity between shapes, the correlation analysis considers where the hotspots of alteration are in the two SD maps, and measures their spatial agreement.
While the ΔVol and Dice metrics were calculated using the binarized disconnection maps, the Corr index was computed using the thresholded SD maps.
The Dice and the ΔVol indexes were computed both considering the entire disconnection maps and subdividing them into their ipsilateral and contralateral components (bilateral tumours were excluded from this last analysis).
To summarise the obtained indexes, we computed the median value and the 25th / 75th percentiles across all subjects for each metric of comparison.
2.5 Statistical Analysis
To test for statistically significant differences in the ΔVol, both considering intra-method and inter-method results at the whole brain or considering only the ipsilateral or contralateral hemisphere, a Wilcoxon rank sum test (significance level α=0.05) was employed.
To assess whether there was a linear relationship between the comparison metrics and the extension of the input mask, we performed a correlation analysis (Spearman Correlation, significance level α= 0.05) between the three indexes (i.e., ΔVol, Dice and Corr) and the volume of the input mask, separately for T and T+O masks. To test the sensitivity of the comparison metrics to the set of thresholds used for the dSD/iSD maps, we defined their normalized range of variation (nRV) as
where is the median value across the dataset of the comparison metric between two SD maps, given the indirect and direct thresholds. The investigated thresholds values were [5%, 10%, 15%, 20%, 25%] and [0.3, 0.5, 0.7].