Animals and experimental design
In this study, male WT C57BL/6 J mice (n = 52) and male transgenic APPKM670/671NL/PS1L166P mice were used (n = 67, referred to as APP/PS1 mice) [11]. Mice were housed in the animal facility of the University of Antwerp during the whole experiment. During the study, mice were kept on a normal 12-h/12-h day-night cycle with ad libitum access to food and water. Additional file 1 shows the weight evolution of the mice in the longitudinal cohort.
APP/PS1 mice start developing Aβ plaques from the age of 6 to 8 weeks and show aggressive amyloidosis in subsequent months [11]. As such, we acquired the first dataset when mice were 2 months old, when only low Aβ plaque deposition is present (n = 20 WT mice and n = 19 APP/PS1 mice). We then longitudinally followed these mice by acquiring a DKI datasets when they were 4 and 6 months of age (intermediate Aβ plaque load) and finally at 8 months of age, which corresponds to an extensive Aβ plaque load. After acquisition of this last DKI dataset at 8 months of age, mice were then sacrificed for histological analysis. In addition, three further cohorts of WT and APP/PS1 mice were scanned once each at 2 months (n = 11 WT mice and n = 16 APP/PS1 mice), 4 months (n = 10 WT mice and n = 16 APP/PS1 mice) and 6 months (n = 11 WT mice and n = 16 APP/PS1 mice) of age and killed thereafter for histological analysis. The complete experimental design is shown in Fig. 1.
MRI scan acquisition
At 2, 4, 6 and 8 months of age, mice were subjected to 1H-MRI scanning using a 7-T PharmaScan MRI scanner with a 16-cm diameter horizontal bore (Bruker, Bremen, Germany). This system is equipped with a standard Bruker cross-coil setup using a quadrature volume coil for excitation and an array mouse surface coil for signal detection. The system was interfaced with a Linux PC running TopSpin 2.0 and Paravision 5.1 software (Bruker BioSpin, Ettlingen, Germany). Anaesthesia was induced using 2% isoflurane (Abbott, Maidenhead, UK) in a gas mixture of 30% O2 and 70% N2 at a flow rate of 600 ml/minute. During MR image acquisition, the isoflurane concentration was initially set at 2% and subsequently lowered when required to maintain a stable respiration rate of (110 ± 10) breaths per minute, which was monitored using a pressure-sensitive pad. In addition, body temperature was monitored via a rectal probe and was held constant between 37.0 °C and 37.3 °C using warm air coupled to a feedback unit (SA Instruments, Stony Brook, NY, USA). PC-sam monitoring software (SA Instruments) was used to measure respiration rate and to measure and control body temperature. Following MR image acquisition, mice were left to recover separately under a heating lamp before being returned to their respective cages.
To ensure uniform slice positioning, we first acquired axial, sagittal and horizontal multi-slice 2D rapid acquisition and relaxation enhancement (RARE) images using the following parameters: repetition time (TR) = 2500 ms, echo time (TE) = 33 ms, matrix size (256 × 256), field of view (FOV) = (20 × 20) mm2, resolution = (0.078 × 0.078) mm2, nine slices, slice thickness = 0.8 mm, and RARE factor = 8. Following correct positioning, we acquired three DKI scans, each with 20 unique diffusion gradient directions and seven b values (400, 800, 1200, 1600, 2000, 2400 and 2800 s/mm2). In addition, seven images without diffusion weighting (b0) were included for each DKI scan. This yielded DKI data comprising 21 b0 images and diffusion weighting for 60 diffusion directions with seven b values. Images were collected with a multi-slice two-shot spin-echo/echo planar imaging sequence with the following parameters: TR = 7000 ms, TE = 23.25 ms, δ = 4 ms, Δ = 12 ms, acquisition matrix = (96 × 96), spatial resolution = (0.214 × 0.214) μm2, 28 horizontal slices, and slice thickness = 0.20 mm. This resulted in a total acquisition time of 1 h, 43 minutes. Next, a high-resolution 3D anatomical image was acquired using a T2-weighted 3D RARE sequence in the same horizontal orientation as the DKI data. The following parameters were used: TR = 3185 ms, TE = 44 ms, and spectral width = 50 kHz, averages = 1, RARE factor = 8, matrix size = (265 × 64 × 50), FOV (20.5 × 13.0 × 10.0) mm3, resolution = (0.080 × 0.203 × 0.200) mm3, and total acquisition time = 21 minutes.
MRI analysis
Prior to the actual image analysis, we performed a visual inspection for quality of the acquired raw data (e.g., ghosting and/or movement) combined with a semi-automated quality control to avoid bias in the diffusion parameter estimation due to MR acquisition artefacts. This semi-automated data quality control consisted of (1) validating the signal decay with increasing b values, (2) validating if signals obtained with high b values did not reach the noise level, (3) validating the magnitude of the signal-to-noise ratio, and (4) validating the parametric estimation error using a chi-square test (see below). Once quality was assured, image pre-processing and analysis were initiated. The data were corrected for motion and eddy current artefacts using the Functional Magnetic Resonance Imaging of the Brain (FMRIB) Software Library (FSL) [13].
The DKI model is shown below:
$$ {S}_{dki}\left(b,g\right)=S(0)\exp \left(-b{\sum}_{i,j=1}^3{g}_i{g}_j{D}_{ij}+\frac{b^2}{6}{\left({\sum}_{i=1}^3\frac{D_{ii}}{3}\right)}^2{\sum}_{i,j,k,l=1}^3{g}_i{g}_j{g}_k{g}_l{W}_{ij kl}\right) $$
(1)
with S(0) being the signal intensity without diffusion weighting, D being the rank 2 diffusion tensor, W being the rank 4 kurtosis tensor, b being the b value, and g being the diffusion gradient direction. The DKI model was voxel-wise fitted to the diffusion-weighted images. The diffusion tensor and kurtosis tensor quantify the Gaussian diffusion profile and the deviation from a Gaussian diffusion distribution, respectively [10, 14]. The diffusion tensor and the diffusion kurtosis tensor were estimated simultaneously using conditional least squares estimators while imposing positivity on the kurtosis coefficients [15]. The conditional least squares estimator explicitly accounts for the Rician MR data distribution, for which the noise level has been estimated [16]. Rotational invariant parameter maps of axial diffusivity, radial diffusivity (RD), mean diffusivity (MD) and fractional anisotropy (FA) were computed from the diffusion tensor, and further refered to as diffusion tensor (DT) metrics, and axial kurtosis (AK), radial kurtosis (RK) and mean kurtosis (MK) maps were computed from the diffusion kurtosis tensor, and further refered to as diffusion kurtosis (DK) metrics, using MATLAB software (Mathworks Inc., Natick, MA, USA) [14, 17].
A study-based atlas was constructed with Advanced Normalization Tools [18] software using 25 randomly selected 3D T2-weighted MRI datasets across both genotypes and all ages. We delineated 23 grey matter and white matter regions of interest (ROIs) on this atlas using AMIRA (version 5.4) (Additional file 2). All individual 3D T2-weighted MRI scans were then normalised to the atlas. The inverse transformation of the normalisation was used to map the atlas ROIs onto the native space of the individual 3D T2-weighted MRI scans. Afterwards, all individual b0 MR images were co-registered to their corresponding 3D T2-weighted MRI scans using FSL software. Using the inverse transformation of the latter, we mapped the individual 3D ROIs to the individual DT and DK metric maps. On the basis of optimal contrast of the FA map to differentiate grey matter and white matter, we finally manually checked and, if necessary, corrected the contours of the grey matter and white matter ROIs to limit a possible partial volume effect. Finally, ROI-averaged DT and DK metrics (MD, axial diffusivity, RD, FA and MK, AK and RK) were extracted from the respective metric maps.
Histological analysis
For histological examination, mice were sacrificed by means of cervical dislocation. The complete brain was dissected and fixed in Fade4 fixative. Fixed brains were sent to HistoGeneX (Antwerp, Belgium), where 5-μm-thick sagittal paraffin-embedded sections were cut from the left hemisphere. Sectioning was started at the middle of the brain, and ten sagittal sections were acquired at 150-μm intervals, covering the whole hemisphere. Immunohistochemical (IHC) staining for Aβ plaques was initiated by depigmenting slides using potassium permanganate for 3 minutes, followed by oxalic acid for 1 minute. Slides were pre-treated in formic acid for 10 minutes to retrieve epitopes and were then incubated for 15 minutes at room temperature with a mouse anti-Aβ (clone 4G8) antibody (1:20,000, SIG-39200; Eurogentec, Angers, France). Next, slides were incubated with a labelled polymer (Dako EnVision + System-HRP Labeled Polymer Anti-Mouse, K4001; Dako). Finally, the substrate was visualised using 3,3′-diaminobenzidine (DAB) chromogen (Dako Liquid DAB+ substrate chromogen system; Dako) for 5 minutes. All steps were performed using the automated Lab Vision Autostainer 480S (Thermo Scientific, Waltham, MA, USA). In every 4G8 staining run, an immunoglobulin G (IgG) control (mouse IgG2b; Dako) was also included. For myelin basic protein (MBP) staining, epitope retrieval was performed in Target Retrieval Solution (Dako) for 30 minutes at 97 °C. After endogenous peroxidase activity was quenched, the slides were incubated for 30 minutes at room temperature with a mouse anti-MBP (clone SMI-94) antibody (1:5000, SMI-94R; Covance Antibody Products, Princeton, NJ, USA). Next, slides were incubated with a labelled polymer (Dako EnVision + System-HRP Labeled Polymer Anti-Mouse, K4001). Finally, the substrate was visualised using DAB chromogen for 5 minutes. All steps were performed on the automated Lab Vision Autostainer 480S. In every MBP staining run, an IgG control (mouse IgG1; Abcam, Cambridge, UK) was included. For microgliosis (ionised calcium-binding adapter molecule 1 [IBA1]) IHC staining, epitope retrieval was performed in citrate buffer (pH 6; Lab Vision) for 30 minutes at 97 °C. After quenching endogenous peroxidase activity, the slides were incubated for 30 minutes at room temperature with a rabbit anti-IBA1 antibody (1:5000, 019-19741; Wako Pure Chemical Industries, Osaka, Japan). Next, slides were incubated with a labelled polymer (Dako EnVision + System-HRP Labelled Polymer Anti-Rabbit, K4003). Finally, the substrate was visualised using the DAB chromogen (Dako Liquid DAB+ substrate chromogen system) for 5 minutes. All steps were performed using the automated Lab Vision Autostainer 480S. In every IBA1 staining run, an IgG control (rabbit IgG; Dako) was included. For astrogliosis (glial fibrillary acidic protein [GFAP]) IHC staining, epitope retrieval was performed using cell conditioning solution (Ventana Medical Systems, Tucson, AZ, USA). The slides were incubated for 28 minutes at 37 °C with a rabbit anti-GFAP antibody (1:7500, Z0334; Dako). Next, slides were incubated with an OmniMap anti-rabbit HRP detection system (Ventana Medical Systems). All steps were performed on the automated Ventana Discovery® XT platform (Ventana Medical Systems). In every GFAP staining run, an IgG control (rabbit IgG; Dako) was also used. For all four staining runs, stained slides were scanned using the MIRAX digital slide scanner (Carl Zeiss Microscopy, Göttingen, Germany).
Quantification of histology
Quantification of the stained slides was performed by DCILabs (Keerbergen, Belgium). In short, the histological images are colour de-convolved to separate the colours on the stained slides [19]. This generates a grey value image of the MBP-, GFAP-, 4G8- and IBA1-stained slides. Hysteresis thresholds were applied on these images to determine the percentage of area positive for the respective staining (percent optical density [%O.D.]) and for each of these ROIs. To quantify the amount of 4G8-positive Aβ plaques, a mask based on an intensity cutoff value to distinguish Aβ plaques from background signal in the grey value images was compiled, which then automatically counted the amount of objects in the mask. To assess the presence of elongated structures in the image, an anisotropic measure was defined as the ratio of the difference and the sum of the eigenvectors of the local structure tensor [20]. All analyses were performed with software developed in-house using OpenSlide C interface (openslide.org) to read the MIRAX images.
3D stacking of histological images and co-registration to MRI
Co-registration of all images of the histological staining was performed, and an overview of the whole co-registration pipeline is shown in Fig. 2. In short, because the MBP staining contains considerable anatomical information, we used the MBP-stained slide of each sectioned interval (150 μm apart) to create a 3D histological reference space. These MBP-stained slides were stacked onto a 3D dataset using an algorithm previously described by Lowe et al. [21] and Ourselin et al. [22]. Next, we co-registered this 3D histological stack to the 3D T2-weighted MRI atlas as described by Ourselin et al. [22] and Modat et al. [23]. Next, a moment-matching algorithm was used to match the 4G8, GFAP and IBA1 staining with the MBP images in order to propagate all histological stains to the corresponding MR images. We then combined all MRI and histological data and created a relational database which contained the following per voxel: the MRI atlas coordinates, the ROI information, all diffusion and diffusion kurtosis metrics, and the corresponding histological values. On the basis of this database, we performed the statistical analyses described below.
Statistical analysis
ROI-based MRI analysis
To evaluate significant differences for each ROI between the different genotypes (WT and APP/PS1 mice), as well as over the different ages of the mice that were followed longitudinally, we used a marginal model with age- and genotype-specific fixed effects. The latter models the evolution of all DT and DK metrics over time while taking into account the association between these DTI and DKI metrics of any given subject at any given age [24,25,26]. A Bonferroni correction was applied for the different ages.
Linear discriminant analysis
Given a set of MRI parameters (MK, AK, RK, MD, axial diffusivity, RD and FA), the linear discriminant analysis (LDA) algorithm searches for a linear combination of MRI parameters for which the genotype misclassification error (MCE; the proportion of incorrectly classified mice) is minimised. The LDA was done using only the DT metrics (MD, axial diffusivity, RD and FA), only the DK metrics (MK, AK and RK) or a combination of both DT and DK metrics. Two-fold cross-validation was performed by randomly splitting the data into a training group and a test group. The training group was composed by randomly selecting ten WT mice and ten APP/PS1 mice. The test group contained the remaining ten WT mice and nine APP/PS1 mice. This process was repeated 1000 times, and the test data MCE was computed for each iteration. The MCE denotes the proportion of mice in the test group assigned to the wrong genotype [27]. To test whether the MCEs were significantly different, we performed a pairwise Kolmogorov-Smirnov test, and the Bonferroni correction was applied for the three tests performed (DT vs DK, DK vs DT + DK and DT vs DT+ DK).
Least absolute shrinkage and selection operator analysis
To identify which MRI metric contributed the most to the correct classification of the genotype, we applied least absolute shrinkage and selection operator (LASSO) logistic regression [28]. In short, we randomly divided all mice into ten distinct groups. A classifying model is trained using MRI metrics from nine of these groups, and the resulting model is applied to the remaining group to determine how good the classifying model is at determining the correct genotypes with different weights of L1 regularisation. We repeated this ten times, withholding a different dataset each time, to achieve a ten-fold cross validation. The result of LASSO analysis indicates if an MRI metric was retained as a classifier (1) or not (0). Moreover, because the 10 distinct groups were chosen randomly, we repeated this whole procedure 100 times. For each metric, we then calculated the percentage of the times it was used to classify the genotypes as an indication of the importance of each metric to correctly classify the genotype.
ROI-based correlation between MRI and histology
To evaluate the potential of the DT and DK metrics as possible markers for the histological features of the motor cortex, we calculated the respective Pearson’s correlation coefficients. A Bonferroni multiplicity correction was applied for the different parameters investigated, and a multiplicity-corrected p value < 0.05 was considered significant.
Voxel-based statistical parametric mapping
To assess voxel-level differences between WT and APP/PS1 mice, we co-registered every individual 3D T2-weighted MR image to the 3D T2-weighted MRI atlas which we had previously constructed. We then transformed all DT and DK metric maps to this 3D T2-weighted MRI atlas and spatially smoothed these maps using a Gaussian kernel. Next, statistical analysis was done for each DT and DK metric separately, using statistical parametric mapping [29]. A false discovery rate correction for multiple comparisons was performed [30], and only voxels with p < 0.05 were visualised.
Voxel-based correlation between MRI and histology
We applied a Bayesian multivariate linear regression [31] to estimate the predictive value of the DT metrics, the DK metrics, or a combination of the DT and DK metrics for the histological outcome at a voxel-based level. Evaluation was done using a leave-one-animal-out methodology. Prediction of the histological values based on the DT or DK metrics of a test mouse was done using a model that was trained on the basis of data of all mice except the test mice. The latter was then validated by cross-validating the correlation coefficients, and the average of the correlation coefficients was calculated after applying a Fisher Z-transformation [32].
