Linking neural and clinical measures of glaucoma with diffusion magnetic resonance imaging (dMRI)

Purpose To link optic nerve (ON) structural properties to clinical markers of glaucoma using advanced, semi-automated diffusion magnetic resonance imaging (dMRI) tractography in human glaucoma patients. Methods We characterized optic neuropathy in patients with unilateral advanced-stage glaucoma (n = 6) using probabilistic dMRI tractography and compared their results to those in healthy controls (n = 6). Results We successfully identified the ONs of glaucoma patients based on dMRI in all patients and confirmed that dMRI measures of the ONs correlated with clinical markers of glaucoma severity. Specifically, we found reduced fractional anisotropy (FA) in the ONs of eyes with advanced, as compared to mild, glaucoma (F(1,10) = 55.474, p < 0.0001, FDR < 0.0005). Furthermore, by comparing the ratios of ON FA in glaucoma patients to those of healthy controls (n = 6), we determined that this difference was beyond that expected from normal anatomical variation (F(1,9) = 20.276, p < 0. 005). Finally, we linked the dMRI measures of ON FA to standard clinical glaucoma measures. ON vertical cup-to-disc ratio (vCD) predicted ON FA (F(1,10) = 11.061, p < 0.01, R2 = 0.66), retinal nerve fiber layer thickness (RNFL) predicted ON FA (F(1,10) = 11.477, p < 0.01, R2 = 0.63) and ON FA predicted perceptual deficits (visual field index [VFI]) (F(1,10) = 15.308, p < 0.005, R2 = 0.52). Conclusion We describe semi-automated methods to detect glaucoma-related structural changes using dMRI and confirm that they correlate with clinical measures of glaucoma.


Methods
We characterized optic neuropathy in patients with unilateral advanced-stage glaucoma (n = 6) using probabilistic dMRI tractography and compared their results to those in healthy controls (n = 6).

Conclusion
We describe semi-automated methods to detect glaucoma-related structural changes using dMRI and confirm that they correlate with clinical measures of glaucoma. PLOS  Introduction Vision loss is a major cause of disability worldwide that is particularly common among the elderly, conferring a greater risk of injury and diminished quality of life [1][2]. Glaucoma is the leading cause of irreversible vision loss worldwide and is projected to affect nearly 80 million individuals by 2020 [3]. It is clinically defined by characteristic patterns of visual field impairment and optic nerve (ON) damage [4]. There is growing evidence that glaucoma may be a neurodegenerative condition [5][6]. Thus, research seeking to understand the possible underlying neurodegenerative processes associated with glaucoma may help guide the development of more robust treatment paradigms and have applications for improving our understanding of other neurodegenerative diseases. In addition, current glaucoma treatments are limited to the reduction of intraocular pressure (IOP) to prevent progressive visual field loss. However, many patients continue to lose vision despite treatment, and no treatments are available to reverse damage that has already occurred to the visual system [4]. Therefore, the development of robust methods enabling earlier diagnosis and more precise quantification of disease progression are essential to limiting glaucomatous damage and improving patient outcomes. MRI-based in vivo human studies using voxel-based morphometry (VBM) and diffusion tensor imaging (DTI) to explore gray-and white-matter cortical changes associated with glaucoma have shown reduced gray-matter (GM) volume in late stages of the disease, along with significant rarefaction along the optic radiations [7][8][9]. While these results are consistent with animal models and post-mortem pathological studies of human subjects [6], there is a need for more precise evaluation of neurological changes, particularly at the level of the ON. Animal model studies of the early visual system using diffusion MRI (dMRI) demonstrate the ability to detect changes in the neural structure of the ONs from damage occurring within the retina [10][11]. In a meta-analysis of studies examining the ONs of human glaucoma patients using various DTI methods, Li, et al. (2014) noted significant decreases in fractional anisotropy (FA) and increases in mean diffusivity (MD) in the ONs of glaucoma patients compared to controls [12]. A number of studies in this meta-analysis also examined the correlation between various measures of disease severity (including glaucoma stage and optical coherence tomography [OCT] measurements) and structural white-matter changes. Generally, increasing glaucoma disease severity is associated with greater white-matter disruption (i.e. decreasing FA and increasing MD) [13][14][15][16][17].
While these studies have quantified the diffusion properties of glaucomatous ONs and their relationship with various clinical measures of disease severity, they rely on older dMRI methodologies. In particular, these studies employ relatively low-resolution, single phase-encoding direction diffusion scans and sample the ONs using manually placed regions of interest (ROIs) [13][14][15][16][17][18]. By sampling from only small portions of the ONs, these techniques are limited in their ability to measure the full extent of optic neuropathy. Moreover, manual ON segmentation is time-intensive, may introduce operator error, and limits the ability to translate this technique into widespread clinical practice. Diffusion MRI methods have advanced substantially since the publication of these studies, and there is an opportunity to investigate and validate methods that rely on semi-automated techniques to assess the ON using dMRI.
Recent dMRI based probabilistic tractography methods have been developed that can more precisely evaluate white-matter changes in the human visual system. These methods have demonstrated altered white-matter structure (including reduced FA) in patients with amblyopia [19]. In this study, we apply this technique to patients with asymmetric glaucomatous optic neuropathy in each eye to evaluate dMRI methods as a diagnostic tool for visual disorders, linking changes in white-matter structure to retinal and perceptual changes in glaucoma. Recent methodological advances make it possible to reliably identify the microstructural properties of the ON with limited user input [20]. Using a pair of diffusion scans acquired with opposite phase-encoding directions, a low-noise field-corrected volume can be created [21][22], allowing the ONs to be isolated using probabilistic tractography. This provides a unique opportunity to quantify changes across the entire length of the ONs in glaucoma patients. Further, we purposefully selected patients with asymmetric glaucomatous ON damage to allow for within-subjects comparisons of ON properties, quantifying differences in eyes with "advanced" versus "mild" glaucoma as defined by the American Academy of Ophthalmology [23].
We used an advanced dMRI tractography method to identify and analyze the ONs of six asymmetric glaucoma patients and six controls. Using both within-subject analyses and comparison to controls, we evaluated structural changes in the ONs associated with glaucoma. Furthermore, we assessed the relationship between these MRI-based neural measures, clinical measures of ON and retinal structure (e.g. vertical cup-to-disc ratio and average peripapillary retinal nerve fiber layer thickness), and perceptual measures (e.g. visual field index).

Participants
Our study was conducted in accordance with the Code of Ethics of the World Medical Association (Declaration of Helsinki) and was approved by the University of Wisconsin-Madison Institutional Review Board. Informed consent was obtained from all participants (written) and all participants completed MRI screening with consultation and approval obtained from their physicians as needed to ensure they could safely participate.
Glaucoma patients. Six glaucoma patients (4 female) aged 19-66 years (mean 53.3 ± 17.4) were recruited from the Glaucoma Service of the University of Wisconsin Hospitals and Clinics (Table 1). All patients had a diagnosis of either primary open-angle, pigment dispersion, pseudoexfoliation, or chronic angle-closure glaucoma and a history of IOPs greater than 22 mmHg. Selection criteria included a Snellen best-corrected visual acuity of 20/25 or better in the eye with "mild" glaucoma and 20/200 or better in the eye with "advanced" glaucoma. Patients with any history of neurodegenerative diseases, normal-tension glaucoma, diabetic retinopathy, advanced macular degeneration, uveitis, or previous (non-surgical) eye trauma were excluded.
Control subjects. Control subjects were recruited from the University of Wisconsin-Madison. Six gender-matched subjects aged 21-34 years (mean 24 ± 5.3) were included in the analysis. All subjects had Snellen best-corrected visual acuity of 20/20 or better and had no prior medical history of neurologic or ocular pathology other than refractive error. Eye dominance was determined as follows: subjects were instructed to form a small aperture using both hands (right and left hands overlapping so a small opening is formed with the inner sides of the palms and thumbs) and to fixate on a distant object through that opening with both eyes open. Without moving their head or hands, subjects were then instructed to close their left eye and were asked whether or not they could still see the object. This same task was repeated with the right eye closed. The eye with which they could see the fixation target was recorded as the "dominant" eye. Ocular dominance was successfully determined for 6/6 control subjects.

Clinical measures
Clinical measures of ON structure and function were assessed for each of the six glaucoma patients during clinical ophthalmologic exams by a glaucoma specialist (Y.L.). These measures included Snellen best-corrected visual acuity (VA), vertical cup-to-disc ratio (vCD), average peripapillary retinal nerve fiber layer thickness (RNFL), and visual field index (VFI). The vCD was determined by direct visualization of the ON using slit-lamp biomicroscopy, average peripapillary RNFL thickness was measured using Cirrus Spectral-Domain Optical Coherence Tomography (Carl Zeiss Meditec, Inc., Dublin, CA, USA) with all scans having adequate signal strength (>7/10), and VFI was measured using the Humphrey visual field 24-2 SITA-Standard testing algorithm (Carl Zeiss Meditec, Inc. Dublin, CA, USA) on visual fields with adequate reliability indices (<33% fixation losses, false positives, and false negatives). Clinical measures were not assessed for control subjects.

Magnetic resonance imaging data acquisition
Brain imaging data was obtained at the Waisman Center in Madison, WI using a GE Discovery Medical 3T MRI scanner (GE Healthcare, Inc., Chicago, IL, USA) equipped with a 32-channel head coil. First, a 10-minute structural whole-brain T1-weighted anatomical scan (2.93 ms TE; 6.70 ms TR; 1 mm 3 isotropic voxels) was acquired. Then, a 15-minute diffusion sequence with two 48-direction diffusion-weighted scans (6 b 0 ), collected in the anterior to posterior (AP) and posterior to anterior (PA) directions (76.7 ms TE; 8.1 s TR; 2x2x2 mm 3 isotropic voxels; b = 2000 s/mm 2 ; reconstruction matrix FOV: LR 212 mm x AP 212 mm x FH 144 mm).

Data processing
Pre-processing. To improve the quality of our tractography and increase the signal-tonoise ratio in the nasal cavity, the two reverse-encoded (AP and PA) diffusion scans were combined into a single corrected volume using the FSL software (University of Oxford, Oxford, England) [21][22]. Subsequent processing was completed using the mrVista software package (Stanford University, Stanford, California), based on previously published methods [20]. A mean b = 0 image was calculated from the corrected DTI volume and underwent eddy current correction. This corrected b 0 image was co-registered to the AC-PC-aligned T1 image and diffusion tensors were fit to the volume using a least-squares estimate bootstrapped 500 times [24].
ROI placement. We manually identified three ROIs along the brain's visual pathway. The T1 image was used to place the left and right ONs and the optic chiasm (OC) by gross  [25][26][27][28][29][30][31][32][33]. Constrained spherical deconvolution (CSD) estimates were used to generate fibers between two ROI pairs, representing the left and right ONs (ON » OC). Whole-brain tractography was completed using an L max of 6 with 500,000 seeds and a maximum of 5,000,000 fibers. A modified white-matter mask generated using mrVista was used to constrain fibers to the brain while still allowing CSDs to be fit within the nasal cavity, enabling detection of the ONs. Final pathways were restricted to fibers passing between the specified ROIs, omitting any spurious results.
Fiber cleaning. Fiber groups were cleaned using the Automated Fiber Quantification (AFQ) toolkit (Stanford University) [34]. Fibers were removed when they were more than 2.6 standard deviations in distance away from the pathway's fiber core or were more than three standard deviations longer than the pathway's average fiber length. After automated cleaning was complete, all optic nerve pathways were subject to quality assessment, with additional manual cleaning as necessary. Fibers were overlaid on the anatomical T1 volume and any fibers that were found to be anatomically implausible were manually removed. Pathways from all 12 participants were processed using the same AFQ cleaning parameters and were manually refined by the same operator (N.M.) (Fig 1). Additional fiber cleaning details are provided in Supplementary S1 Table,  Voxel-wise tensor properties were extracted from the volumetric region defined by each tractography-generated pathway. The main diffusion properties included in our analysis were mean diffusivity (MD, μm 2 /s) and fractional anisotropy (FA) (Fig 2). MD provides an average measure of pathway diffusivity and is a useful approximation of white-matter density, where large values indicate a diffuse ("weak") pathway, and small values indicate a denser ("strong") pathway [35]. FA provides a measure of diffusion directionality and is indicative of microstructural changes across various visual and non-visual pathologies, where small values indicate multiple intersecting, degenerated, or demyelinated pathways [35][36]. Studies examining white-matter changes in other visual disorders like amblyopia, macular degeneration, and retinitis pigmentosa have shown decreased FA in the optic radiations, among other visual pathways [19][20][37][38][39]. While highly indicative of microstructural changes, FA is nonspecific and provides little information about the nature of the underlying tissue properties (e.g. demyelination, axonal degeneration, crossing fibers, etc.). To more precisely characterize the MD and FA measures, we also assessed the component measures, radial diffusivity (RD, μm 2 /s) and axial diffusivity (AD, μm 2 /s). RD has been demonstrated to be more indicative of changes in white-matter myelination, while AD is more indicative of axonal degeneration [40]. Typically, large RD values are evidence of possible demyelination, while small AD values are evidence of possible axonal degeneration.
Analysis. A combination of within-and between-groups analyses were conducted to evaluate the structural white-matter changes associated with glaucoma-related vision damage. All six asymmetric glaucoma patients had one eye with no or early glaucomatous visual field defects ("mild") and one eye with moderate or advanced glaucomatous visual field defects ("advanced") [23]. Within-subjects, diffusion-tensor properties were compared between the "advanced" and "mild" eyes. Between-subjects, the ratios of "advanced" / "mild" ON FA in glaucoma patients were compared to non-dominant/dominant ratios in controls. This design minimized the possible effect of global changes in white-matter properties as a result of age.
To facilitate these comparisons and normalize pathway lengths, 100 samples were taken along the length of each pathway such that 100 average MD, FA, RD, and AD values were available for each pathway in each of the two groups (glaucoma patients and control subjects). These values were generated for each cross-section using a Gaussian-weighted average, where the calculated "core" of each pathway was selectively weighted over the outlying fibers. From  Linking neural and clinical measures of glaucoma with diffusion magnetic resonance imaging (dMRI) these 100 samples, the central 80 were retained for further analysis to reduce the risk of including measures contaminated by retinal cell bodies or contralateral fiber tracts [19][20]38]. The central 80% of each sample was further subdivided into 10% bins to more precisely quantify differences along the pathway length.
Statistical analysis. Linear mixed effect (LME) models were used to compare the MD, FA, RD, and AD values of the advanced and mild ONs across the middle 80% of samples and at each of the eight 10% bins in glaucoma patients. This model factored glaucoma severity (Group; advanced or mild) as a fixed effect and subject as a random effect. These models followed the general structure: MD/FA/RD/AD~Group + (1|Subject). Reported p-values are from ANOVAs of the fixed "glaucoma severity" effects. Similar LME models were used in the glaucoma and control ratio comparisons, as well as in all correlational data (factoring clinical and neurological measures as fixed effects and subjects as random effects). For all correlations, reported R 2 values are adjusted to the number of predictors included in the LME model.
To compensate for multiple comparisons in the advanced vs. mild tract profile analyses, we implemented the Benjamini and Hochberg False Discovery Rate (FDR) test [41]. In this test, the fixed effect (advanced vs. mild) p-values from each LME model were pooled and an FDR value was calculated for each test. A total of 36 p-values were pooled (eight individual bin tests and one whole-pathway test for each of the four diffusion measures) and tests for which the FDR value was less than our target alpha (α = 0.05) were considered significant.

Diffusion magnetic resonance imaging
ON white-matter pathways were successfully identified and refined in 6/6 glaucoma patients and 6/6 controls. All pathways appeared to be anatomically plausible after cleaning and were amenable to within-subjects and group-wise comparisons. A significant difference in mean FA between the advanced and mild ONs of glaucoma patients was noted in 3/8 bins (F 1 (1,10) = 55.442, p 1 = 2.20e-5, FDR 1 = 4.0e-4; F 2 (1,10) = 18.382, p 2 = 0.0016, FDR 2 = 0.019; F 3 (1,10) = 11.322, p 3 = 0.0072, FDR 3 = 0.048), along with a significant difference across the middle 80% of samples (F(1,10) = 55.474, p = 2.19e-5, FDR = 4.0e-4) (Fig 4). In the same pathway, a significant difference in average MD was noted in 1/8 bins (F(1,10) = 10.885, p = 0.0080, FDR = 0.048). For RD, we found a significant difference in 1/8 bins (F(1,10) = 15.87, p = 0.0026, FDR = 0.023). No significant difference in AD was found (all FDR > 0.05). Our subsequent analyses focus on ON FA because of the magnitude and reliability of the effect (compared to MD and RD). To more precisely characterize the nature of these within-subject effects, we compared the FA ratios of advanced/mild ONs in glaucoma patients to non-dominant/dominant ONs in controls. As in earlier analyses, these ratios were computed from the central 80% of samples, omitting samples from the first and last 10% of the pathway estimates. We found selective FA reductions in the "advanced" ONs of glaucoma patients. As expected, there were no significant differences between non-dominant and dominant ON FA values in control subjects. We found a significant difference between the ON FA ratios of glaucoma patients compared to controls in a LME model including group and age as fixed factors (F(1,9) = 20.276, p = 0.0015) (Fig 5).

Relating dMRI to clinical measures in glaucoma
Reductions in ON FA were found to correlate with clinical measures of glaucoma (vCD, RNFL, and VFI). We found that vCD predicted ON FA (F(1,10) = 11.061, p = 0.0077, R 2 = 0.66), Optic nerve fractional anisotropy ratios in glaucoma patients (n = 6) and controls (n = 6). Comparison of optic nerve fractional anisotropy ratios (%) in glaucoma patients (green) and controls (blue). Glaucoma patient ratios were calculated for "advanced" / "mild" optic nerve fractional anisotropy, and control subject ratios were calculated for non-dominant/dominant optic nerve fractional anisotropy. Mean ratios are indicated by the bold lines, with standard error denoted by the surrounding box. The dashed gray line at 100% represents baseline, i.e. the same fractional anisotropy measures in the two optic nerves. We found a significant difference between fractional anisotropy ratios in glaucoma patients versus controls (p = 0.0015). https://doi.org/10.1371/journal.pone.0217011.g005 Linking neural and clinical measures of glaucoma with diffusion magnetic resonance imaging (dMRI) RNFL predicted ON FA (F(1,10) = 11.477, p = 0.0069, R 2 = 0.63), and ON FA in turn predicted VFI (F(1,10) = 15.308, p = 0.0029, R 2 = 0.52) (Fig 6). Thus, dMRI measures of FA in the optic nerve reliably linked the retinal and perceptual deficits observed in our patient sample.

Discussion
We assessed the utility of probabilistic dMRI tractography in correlating neural and clinical measures of ON damage in patients with asymmetric glaucoma and normal controls. We isolated the ON white-matter pathway in 6/6 glaucoma patients and 6/6 controls. We combined AP and PA phase-encoded dMRI volumes to recover from imaging distortions caused by the adjacent nasal cavities. To our knowledge, this is the first probabilistic tractography study to successfully isolate the ONs in glaucoma patients. Previous work relied largely on manual ON segmentation or ROI-based analyses, while our methodology allows the entire ON to be identified with minimal manual operator input. This technique increases the efficiency of data collection and reduces the risk of operator error. We correlated clinical measures of ON structure and function (i.e. vCD, RNFL, and VFI) with dMRI measures of ON structure (i.e. FA, MD, RD, and AD). Our methods sample diffusion measures along the entire length of the ON (rather than small targeted regions) and provides a more comprehensive account of dMRI measures of disease.
We found significant differences in average FA, MD, and RD of ONs with "advanced" versus "mild" glaucoma. These trends (smaller FA and greater MD and RD in "advanced" ONs) indicate disruption to the visual system white-matter consistent with clinical measures of glaucoma severity. Our results were robust to False Discovery Rate testing to compensate for multiple comparisons and our small sample size. Our findings are mostly consistent with earlier studies, which showed general trends of decreasing FA and increasing MD with increasing glaucoma severity [13][14][15][16][17][18]. However, there was a significant difference in radial, but not axial diffusivity, which would suggest that glaucoma predominantly impacted myelination rather than axonal degeneration of ONs with advanced glaucoma. However, increased RD was noted in only 1/8 bins along the length of the ON, so we are limited in our ability to comment on the precise contributions of axonal degeneration and/or demyelination to the observed changes in FA. Additionally, this result differs from previous work, where changes were noted in both RD Linking neural and clinical measures of glaucoma with diffusion magnetic resonance imaging (dMRI) and AD [15]. This discrepancy may be the result of differences in diffusion sequence parameters but is most likely the result of our small sample size.
To minimize the impact of the difference in mean age between our glaucoma patients and controls in this study, we primarily relied on within-subjects comparisons of the advanced and mildly glaucomatous ONs. Through comparison of ON FA ratios in glaucoma patients and controls, we determined that these structural changes were not global changes in the ONs of glaucoma patients, but rather unilateral changes in the "advanced" glaucomatous eyes. By using a ratio comparison, we compensated for the difference in mean age between our glaucoma patients and controls. Regardless of age, within-subject ratios of the left and right optic nerves of healthy control subjects would be expected to be~100%. Any global changes to white-matter structure as a result of age (and absent any additional ocular or neurological disease) would be expected to affect the left and right ONs equally. Thus, by restricting betweensubject analyses to a ratio comparison, we minimized the potential confound of the differences in age between our patients and controls.
Lastly, we examined the correlation between clinical glaucoma measures (vCD, RNFL, and VFI) and dMRI neural measures (ON FA). We found that vCD and RNFL were correlated with ON FA and that ON FA was predictive of VFI. Thus, our analysis confirms and quantifies correlations between measures of glaucoma-related damage between neural (ON FA), structural (vCD), retinal (RNFL), and functional (VFI) measures.
In summary, we correlated neural dMRI measures and clinical glaucoma measurements in patients with asymmetric glaucoma damage. Our results using current dMRI methods agree with previous studies using older techniques and validate probabilistic tractography methods that assess deficits in the visual pathways of glaucoma patients. Future larger, prospective studies may evaluate dMRI as a possible diagnostic tool for glaucoma evaluation. In addition, studies quantifying the relationship between neural and clinical measures longitudinally may determine whether these methods may be useful for monitoring glaucoma progression and treatment efficacy.  Table. Optic nerve fiber cleaning details for glaucoma patients and control subjects. Detailed description of the specific fiber cleaning results for six glaucoma patients (denoted G1-G6) and six control subjects (denoted C1-C6). The total number of fibers removed, percentage of fibers retained, and the normalized pathway volumes of the left and right optic nerves of all glaucoma patients and control subjects are provided. Values corresponding to the "advanced" glaucomatous eyes are marked in bold. Pathway volume estimates were generated using AFQ and are based on the number of unique fiber coordinates at each sample point. The number of fibers removed, percentage retained, and the normalized pathway volumes varied between subjects based on the number of initial fiber streams isolated. There was not a significant difference between average normalized volumes of glaucomatous and healthy control optic nerves (t(1,22) = -0.4185, p = 0.68). Optic nerve volumes across all subjects were found to be normally distributed in the Shapiro-Wilk test (p > 0.05) [42] and no significant outlier volumes were found with Iglewicz and Hoaglin's outlier test (all modified Z scores < 3.5) [43]. (DOCX)