Investigating changes in axonal density and morphology of glaucomatous optic nerves using fixel-based analysis

PURPOSE To characterize neurodegeneration of glaucomatous optic nerves (ONs) in terms of changes in axonal density and morphology using fixel-based analysis (FBA), a novel framework for analyzing diffusion-weighted MRI (DWI). Furthermore, we aimed to explore the potential of FBA measures as biomarkers of glaucomatous ON degeneration. METHODS DWI scans were obtained from 15 glaucoma patients and 15 controls. ONs were tracked and segmented into their three anatomical segments; intraorbital (IO), intracanalicular (ICAN) and intracranial (ICRAN). For each segment, FBA measures were computed, which included fiber density (FD; a measure of axonal density), fiber-bundle cross-section (FC; an estimate of morphological changes), and fiber density and cross-section (FDC). Peripapillary retinal nerve fiber layer (pRNFL) thickness and visual field mean deviation (VFMD) were assessed for glaucoma patients. ANCOVA was used to compare FBA values between the two groups, and Spearman’s correlation analysis was used to test the correlation between FBA measures and pRNFL thickness and VFMD. RESULTS All glaucomatous ON segments showed a significant loss of FD and FDC compared to the controls, while a loss of FC was found in the IO and ICRAN segments only. FD and FDC values of the IO and ICAN segments of glaucomatous ONs showed significant correlations with pRNFL thickness and VFMD. CONCLUSIONS Glaucomatous ONs exhibit lower FD and FC compared to controls, indicating axonal loss and gross atrophy. The correlation between FBA measures of glaucomatous ONs and established clinical tests of glaucoma demonstrates the potential of FBA measures as biomarkers of glaucomatous ON degeneration. C H AP TE R 3


INTRODUCTION
Glaucoma is a progressive, irreversible optic neuropathy, characterized by the death of retinal ganglion cells (RGCs) and a distinct pattern of visual field loss. 1 Clinical evaluation of glaucoma involves the assessment of both functional loss and structural retinal degeneration.
A hallmark of glaucomatous structural degeneration is the loss of retinal nerve fiber layer (RNFL) thickness. The RNFL is composed of RGCs axons, and therefore exhibits a decrease in thickness following the loss of RGCs due to glaucoma.
Glaucomatous structural loss, however, is not confined to the RNFL, but it rather manifests throughout the entire length of the optic nerve (ON), 2 as the ON is primarily formed of myelinated RGC axons. Therefore, probing the ON in glaucoma patients could present a new approach for detecting glaucomatous structural degeneration. Indeed, diffusionweighted MRI (DWI) studies of glaucoma patients have reported degenerative changes in glaucomatous ONs beyond the retina. [3][4][5][6][7][8][9][10][11][12] Furthermore, DWI measures of glaucomatous ONs were found to correlate with ophthalmic measures of glaucoma, indicating the potential of DWI measures as biomarkers of structural glaucomatous degeneration.
To date, most DWI studies of glaucomatous ONs have relied on diffusion tensor imaging (DTI) for data analysis. DTI uses a tensor model to characterize water diffusion on a voxel level, producing measures of white matter (WM) structural integrity such as fractional anisotropy (FA) and mean diffusivity (MD). 13 However, such DTI measures are nonspecific, as they reflect a wide range of WM structural changes. For example, a change in FA could be a response to changes in membrane permeability, myelination, axonal size or axonal packing density. 14 Furthermore, due to its single-compartmental nature, the tensor model is incapable of distinguishing different tissue types within the same voxel. 15 This can affect the reliability of tensor-derived measures, especially in relatively small structures such as the ON, where partial volume effects from surrounding tissues are inevitable. Therefore, adopting higher-order models which overcome some of the limitations of DTI could improve the prospects of DWI measures as biomarkers of glaucomatous degeneration.
Fixel-based analysis (FBA) is a novel framework for analyzing DWI, which overcomes some of the limitations of DTI by using a multi-compartmental model to analyze distinct WM fiber populations within a voxel (termed "fixels"). 16 By doing so, FBA is able to produce measures that directly correspond to biological changes occurring in WM fibers (Figure 1). These measures are (1) fiber density (FD), (2) fiber-bundle cross-section (FC), and (3) fiber density and cross-section (FDC). FD examines microstructural WM changes by estimating the axonal density of a specific WM fiber population within a voxel. 17 16 ) In this study, we aimed to characterize glaucomatous WM degeneration along the different anatomical segments of the ON using FBA. Furthermore, we aimed to evaluate the potential of FBA measures as biomarkers for glaucomatous ON degeneration by determining their correlation with established ophthalmic tests of glaucoma. Finally, for comparison, we performed the same analyses using the conventional voxel-based DTI approach.

ETHICAL APPROVAL
This study was approved by the ethics board of the University Medical Center Groningen (UMCG). The study adhered to the tenets of the Declaration of Helsinki. All participants provided written informed consent prior to participation. 50

PARTICIPANTS
Eighteen glaucoma patients and 18 control subjects were initially recruited for this crosssectional study. Inclusion criteria for glaucoma patients was having been diagnosed with glaucoma in at least one eye. Inclusion criteria for controls included having an intraocular pressure (IOP) of ≤ 21 mmHg, no visual field defects, and a visual acuity of 0.8 or higher (+0.1 logMAR or less) in both eyes. Exclusion criteria for controls and glaucoma patients included having any ophthalmic disease that affects optic nerve structure or function (apart from glaucoma in the glaucoma group), having a neurodegenerative disorder, having an unsuccessful eye exam or MRI scan during data collection, having a contraindication to being inside an MRI scanner, and the detection of any apparent pathology on an acquired MRI scan. Two controls were excluded for having uncorrectable MRI artifacts compromising their DWI scans. Three glaucoma patients were excluded, all for having incomplete eye exams. In total, 15 glaucoma patients and 15 controls were included in this study.

OPHTHALMIC TESTS
All participants underwent visual acuity, tonometry, optical coherence tomography (OCT) and standard automated perimetry exams. Visual acuity was tested using Snellen charts with optimal correction for the viewing distance. Tonometry was performed using a Tonoref non-contact tonometer (Nidek, Hiroishi, Japan). A Canon OCT-HS100 scanner (Canon, Tokyo, Japan) was used to measure the average peripapillary RNFL (pRNFL) thickness. For glaucoma patients, perimetry was performed using a Humphrey Field Analyzer (HFA; Carl Zeiss Meditec, Jena, Germany) using the 30-2 grid and the Swedish Interactive Threshold Algorithm (SITA), and the results were expressed as visual field mean deviation (VFMD). For controls, perimetry was performed using Frequency Doubling Technology (FDT; Carl Zeiss Meditec, Jena, Germany) with the C20-1 screening mode to exclude controls with any visual field defects.

DWI DATA ACQUISITION
Whole-brain DWI scans were acquired on a Siemens MAGNETOM Prisma 3T MRI scanner (Siemens, Erlangen, Germany) with a 64-channel head coil. The following parameters were used for data acquisition: echo-time/repetition time (TE/TR) = 85/5500 ms, 66 axial slices, slice thickness = 2 mm, voxel size = 2 mm isotropic, field of view (FoV) = 210 x 210 mm, two DWI shells (b=1000 s/mm 2 and b=2500 s/mm 2 ), 64 diffusion gradient directions for each shell in two phase-encoding directions (anteroposterior and posteroanterior), and three non-diffusion weighted images (b=0 s/mm 2 ) in each phase-encoding direction (six in total).
During scanning, we ensured the inclusion of the entire length and thickness of both ONs in the FoV of each participant. The participants were asked to close their eyelids and reduce eye movements during the scan as much as possible. 51

DATA PREPROCESSING
The images were first denoised, 18 then corrected for susceptibility, 19 motion, and eddy current induced distortions 20 using FMRIB's Software Library (FSL v5.011, https://fsl.fmrib. ox.ac.uk/fsl), and finally corrected for bias field inhomogeneity. 21 The images were then visually inspected for any uncorrected motion artifacts of the ONs.

FIXEL-BASED ANALYSIS
The complete FBA pipeline has been previously described in detail by Raffelt et al. 16 Here, we briefly give an overview of the main steps of FBA, and specifically highlight any performed analyses which deviate from the default FBA pipeline. All described steps have been performed in MRtrix3 (www.mrtrix.org). 22 A summary of the main steps of this study's pipeline is illustrated in Figure 2.
First, average tissue response functions were produced from the WM, grey matter, and cerebrospinal fluid of all included subjects. 23 Then, all images were upsampled to an isotropic voxel size of 1.3 mm. Next, brain masks were computed for the upsampled images and manually edited to include the ONs. Multi-tissue constrained spherical deconvolution 24 was then performed using the produced set of average tissue response functions to estimate fiber orientation distributions (FODs) for each subject, followed by joint bias field correction and intensity normalization. 17 A group-specific FOD template was produced from all 30 included subjects, 25

and all
FODs were then non-linearly registered to the FOD template space to attain spatial correspondence. 26 Subsequently, the FOD lobes were segmented to identify the number and orientation of fixels within each voxel, and the fixels were assigned FD values based on the amplitude of the FOD lobes. Then, the fixels were reoriented in FOD template space based on the non-linear transformation of their corresponding voxels, and then assigned to their spatially corresponding template fixels. The FC was then computed for each fixel based on the warps produced during registration. Finally, the FDC was calculated based on the computed FD and FC values.

OPTIC NERVE TRACKING AND SEGMENTATION
In order to analyze the ONs, deterministic tractography was used to track the ONs in the  ii) The result of deterministic tracking between the ROIs, color-coded by directionality. iii) ON segmentation. Each optic nerve is segmented into 14 sub-segments, and each sub-segment is assigned to an anatomical ON segment. Alternating shades of the same color represent different subsegments of the same anatomical segment. Green: intraorbital segment; red: intracanalicular segment; blue: intracranial segment. 53 Following ON tracking, approximately 3 mm of the distal ends of the ONs were excluded from the analysis to avoid signal contamination from the posterior poles of the eye globes. 12,27 The ONs were then segmented into 14 sub-segments of equal coronal thickness, each segment measuring 2.6 mm (or the equivalent of two slices) in thickness. The resulting ON subsegments were then numbered in ascending order, starting at the sub-segments closest to the eye globes and ending at the optic chiasm. Then, each sub-segment was assigned to were assigned to the IO segment. The ICRAN was identified as the segment emerging from the proximal end of the optic canal to enter the cranial cavity and ending at the optic chiasm.
Finally, the FD, FC and FDC values were computed for each of the 14 sub-segments. The FD was computed using the Apparent Fiber Density framework described by Raffelt et al.,17 where the amplitude of an FOD along a WM fiber population is assumed to represent the within-voxel volume occupied by the axons of that specific fiber population. The FC was calculated using the non-linear warps required to co-register the fixels of each subject to the population template, where a positive FC indicates a cross-sectional area greater than that of the population template and a negative FC represents a cross-sectional area smaller than that of the population template. The FDC was computed as described by Raffelt et al., 16 to incorporate both FD and FC, providing a measure which takes into account both micro-and macrostructural changes.

DTI ANALYSIS
For comparison, a conventional voxel-based DTI analysis was also performed. First, diffusion tensors were produced from the b=0 s/mm 2 and b=1000 s/mm 2 shells of the preprocessed diffusion data. Then, FA and MD maps were derived from the tensors. The maps were then warped to the FOD template space using the same non-linear subject-to-template warps produced during FOD registration. Voxel-based masks for each of the 14 ON sub-segments were also created. Finally, the FA and MD of each sub-segment were computed.

STATISTICAL ANALYSIS
As some of the glaucoma patients were affected unilaterally and others bilaterally, a single ON was included from each patient for statistical analysis. For unilaterally affected patients, the glaucomatous ON was included; for bilaterally affected patients, a single ON was randomly selected. For the control group, a single random ON was included for each 54 subject, ensuring to keep the same right-to-left side ratio as the glaucoma group. In total, nine right-sided and six left-sided ONs were included from both groups.
To analyze the patients' demographics and clinical characteristics of the included ONs, independent-samples t-test was used for parametric continuous variables, Mann-Whitney U test was used for non-parametric continuous variables, and χ 2 test was used for categorical variables. Statistical significance was reported at P < 0.05.
To compare FBA and DTI measures of each ON sub-segment between the two groups, analysis of covariance (ANCOVA) was used, while controlling for age and sex as nuisance co-variates. To study the three ON anatomical segments, FBA measures were first averaged for each segment, then ANCOVA was also used to compare the measures of each anatomical segment between the two groups. False discovery rate (FDR) was used to correct for multiple comparisons, and statistical significance was reported at P FDR < 0.05.
A within-subject comparison between the average FBA measures of each ON segment was performed for the controls and glaucoma patients using within-subjects ANOVA, and statistical significance was reported at a post-hoc Bonferroni corrected P < 0.05. Spearman's rank correlation analysis was used to study the correlation between the average FBA and DTI measures of each anatomical segment of the glaucomatous ONs and the average pRNFL thickness and VFMD. A Spearman's correlation coefficient (ρ) ≤ 0.3 signified a very weak correlation, 0.3 < ρ ≤ 0.5 signified a weak correlation, 0.5 < ρ ≤ 0.7 signified a moderate correlation, and ρ > 0.7 signified a strong correlation. Statistical significance was reported at P < 0.05.

DEMOGRAPHICS AND CLINICAL CHARACTERISTICS
A summary of participants' demographics and the clinical characteristics of the eyes corresponding to the included ONs can be found in Table 1

FIXEL-WISE ANALYSIS
All sub-segments of the glaucomatous ONs showed a significant decrease (P FDR < 0.05) in  (Table 2).

DTI ANALYSIS
Sub-segments 2 to 14 of the glaucomatous ONs showed a significant decrease

WITHIN-SUBJECT COMPARISON OF ON SEGMENTS
For the glaucoma group, FD was significantly higher in the IO and ICAN segments compared to the ICRAN segment (P = 0.044 and P = 0.001, respectively) and FDC was higher in the ICAN segment compared to the ICRAN segment (P = 0.010). For the control group, FD was significantly higher in the IO segment compared to the ICAN and ICRAN segments (P = 0.023 and P = 0.001, respectively) and FDC was higher in the IO segment compared 58 to the ICAN and ICRAN segments (P = 0.016 and P = 0.006, respectively). FC showed no difference between segments in both groups.

CORRELATION WITH OPHTHALMIC TESTS
Average FBA and DTI measures of the anatomical segments of glaucomatous ONs showed varying degrees of correlation with the average pRNFL thickness and the VFMD ( Table 3).
The strongest correlation was found between the FDC of the IO segment and VFMD (ρ 2 = 0.83, P < 0.001), while the FC of all segments showed no correlation with either of the ophthalmic tests.  While all anatomical segments of the glaucomatous ONs showed a significant decrease in FD, the effect size of this decrease appears to be considerably higher at the IO segment.
From a histopathological perspective, we expect axonal loss to be uniform throughout the length of the ON. Therefore, this variation in effect size along the ON is most likely due to technical limitations of DWI. The skull base and air-filled paranasal sinuses adjacent to the ONs create susceptibility-induced magnetic field inhomogeneities, which can cause geometric distortions and signal loss. 27,31 These magnetic susceptibility artifacts are more likely to affect the ICAN and ICRAN segments, as they are in closer proximity to the magnetic field distortions compared to the IO segment, which explains the lower effect sizes found in these segments. This, in addition to the within-subject FD differences we found between the ON segments, highlights the importance of assessing the entire length of the ON to overcome some of the limitations of DWI of the ON and to gain a more complete picture of glaucomatous ON degeneration.

SEGMENT-SPECIFIC ATROPHY OF GLAUCOMATOUS ONs
Glaucomatous ONs showed atrophic changes (as indicated by a decrease in FC) at the IO  Figure 3B). This finding further demonstrates the importance of assessing the entire length of glaucomatous ONs.  (Figure 4). This could be due to the fact that our DTI analysis utilized only the low b-value shell for tensor fitting, while FBA used the high b-value shell as well, which is associated with lower signal-to-noise ratio and greater geometric distortions, 37 to which the ICAN and ICRAN segments are susceptible due to their anatomical proximity to the skull base and paranasal sinuses.

FBA MEASURES AS BIOMARKERS OF GLAUCOMATOUS ON DEGENERATION
A moderate to strong correlation was found between the FD of the glaucomatous ONs and the average pRNFL thickness and VFMD. This correlation is expected, as a decrease in FD represents a loss of RGC axons, which reflects a loss of structure and function of the RGCs at the level of the retina. A primate study found a similar correlation between the RNFL thickness and nerve fiber density in ONs with experimentally induced glaucoma. 29 FC of the glaucomatous ONs, however, showed no significant correlation with either of the ophthalmic tests. This is most likely due to the difference in scale between these tests, as the FC is a coarse measurement of gross morphological changes while the ophthalmic tests operate on a much finer scale. The FDC showed the strongest correlation with the ophthalmic tests, specifically at the IO segment. As FDC is a combination of both FD and FC, it is supposed to represent the overall information carrying capacity of the ONs, which explains its strong correlation with the ophthalmic tests. This finding suggests that FDC of the IO segment offers the most accurate representation of glaucomatous ON degeneration and hence shows the most potential for being used as a biomarker of glaucomatous ON degeneration.

COMPARISON TO PREVIOUS FBA STUDIES
In a previous study, we utilized FBA to investigate post-chiasmatic visual pathway changes in POAG patients to better understand the underlying pathophysiology of the central degeneration found in these patients. 38 We found a significant loss of FD, FC and FDC in the optic tracts (OTs) of POAG patients, which is in line with our current findings as the OTs

STUDY LIMITATIONS AND FUTURE DIRECTIONS
The present study has a number of limitations. Firstly, the number of participants is moderate.
However, our group sizes are still in line with previous DWI studies of glaucoma. 3,12,[43][44][45][46] Future studies could benefit from a larger sample size and a separate analysis of different stages of glaucoma. Furthermore, a normative database including representative samples of normal ONs of different age groups, genders and ethnicities would be needed to make FBA measures of ONs clinically applicable.
Another limitation is the relatively low resolution of the DWI scans in comparison to ON size. However, this limitation is not specific to our study, but is rather a general limitation of imaging small structures like the ON using DWI. Optimizing a DWI protocol for ON imaging and using a dedicated surface coil could increase image resolution and reduce scanning time considerably, making FBA a more feasible clinical tool for assessing ON structural changes in glaucoma. 62

CONCLUSIONS
Using FBA, we were able to detect axonal loss (as indicated by FD loss) and gross atrophy