Multi-Contrast Magnetic Resonance Imaging of Visual White Matter Pathways in Patients With Glaucoma

Abstract

Purpose: Glaucoma is a disorder that involves visual field loss caused by retinal ganglion cell damage. Previous diffusion magnetic resonance imaging (dMRI) studies have demonstrated that retinal ganglion cell damage affects tissues in the optic tract (OT) and optic radiation (OR). However, because previous studies have used a simple diffusion tensor model to analyze dMRI data, the microstructural interpretation of white matter tissue changes remains uncertain. In this study, we used a multi-contrast MRI approach to further clarify the type of microstructural damage that occurs in patients with glaucoma.

Methods: We collected dMRI data from 17 patients with glaucoma and 30 controls using 3-tesla (3T) MRI. Using the dMRI data, we estimated three types of tissue property metrics: intracellular volume fraction (ICVF), orientation dispersion index (ODI), and isotropic volume fraction (IsoV). Quantitative T1 (qT1) data, which may be relatively specific to myelin, were collected from all subjects.

Results: In the OT, all four metrics showed significant differences between the glaucoma and control groups. In the OR, only the ICVF showed significant between-group differences. ICVF was significantly correlated with qT1 in the OR of the glaucoma group, although qT1 did not show any abnormality at the group level.

Conclusions: Our results suggest that, at the group level, tissue changes in OR caused by glaucoma might be explained by axonal damage, which is reflected in the intracellular diffusion signals, rather than myelin damage. The significant correlation between ICVF and qT1 suggests that myelin damage might also occur in a smaller number of severe cases.

Glaucoma is a chronic progressive optic neuropathy in which retinal ganglion cell damage causes visual field loss and optic nerve damage (Fig. 1A).¹ Because glaucoma is the leading cause of blindness, with a high prevalence rate in elderly populations,² there is an urgent need to understand how glaucoma affects nerve fiber pathways and brain areas.³ 

Figure 1.

 

Clinical features of patients with glaucoma and schematic illustration of the visual system. (A) Glaucomatous fundus photograph of a representative patient with glaucoma (Glc-015; left panel, right eye; right panel, left eye). The small panel at the top right of each main panel depicts associated visual field loss, as measured using the Humphrey Field Analyzer. The superior neural losses correspond to deficits on the lower visual fields of both eyes. (B) A schematic illustration of the early visual pathways from the eye to the primary visual cortex. Retinal ganglion cells receive visual information from photoreceptor cells via two types of intermediate cells. The optic nerve is composed of retinal ganglion cell axons and glial cells, which extend from the optic disc to the optic chiasm and continue as the optic tract to the lateral geniculate nucleus (LGN). From the LGN, fibers of the optic radiation carry visual information to the primary visual cortex in the occipital lobe of the brain.

Recently, diffusion tensor imaging (DTI) has allowed researchers to investigate tissue abnormalities of visual white matter tracts, such as the optic tract (OT) and optic radiation (OR) (Fig. 1B), in living humans.^4,5 DTI studies have previously demonstrated that patients with glaucoma show abnormalities in diffusivity measurements, such as decreased fractional anisotropy (FA)⁶ and increased mean diffusivity (MD) and radial diffusivity (RD) along the OT and OR, with numerous research groups replicating such results^7–18 (for a review on past glaucoma DTI studies, see Nuzzi et al.¹⁹). However, a major limitation is that DTI-based metrics do not directly correlate with specific types of tissue properties, such as axon diameter and myelination.^5,20–24 

This study aimed to improve our understanding of the white matter consequences of glaucoma by combining two neuroimaging approaches. One approach is to analyze diffusion magnetic resonance imaging (dMRI) data using biophysical models, which aim to parameterize dMRI signals as a function of biologically meaningful parameters, such as axon density, to provide better explanations for observed abnormalities in white matter. Neurite density and orientation dispersion imaging (NODDI)²⁵ is one such model that assumes three types of microstructural environments (intracellular, extracellular, and free diffusion in cerebrospinal fluid [CSF]). NODDI provides three metrics: intracellular volume fraction (ICVF), orientation dispersion index (ODI), and isotropic volume fraction (IsoV), each of which may reflect intra-axonal volume, orientation dispersion of axons, and volume of the CSF, respectively. Previous studies have provided evidence that NODDI-based metrics correlate with specific properties of axons,^26–28 such that orientation dispersion measured by histology showed stronger correlation with ODI than with FA.^29,30 Although there are debates on model assumptions and interpretive limitations, NODDI may be a promising neuroimaging method to evaluate white matter damage specific to axons in clinical populations.³¹ The other approach uses MRI data acquisition schemes other than the dMRI. For example, quantitative T1 (qT1) mapping aims to quantify T1 relaxation time by combining structural MRI data acquired with multiple parameters.^32–34 It is hypothesized that qT1 is relatively specific to myelin levels in the white matter based on comparisons with histological data.³⁵ In this study, we focused on NODDI metrics (ICVF, ODI, and IsoV) and qT1 to evaluate white matter tissue properties in patients with glaucoma and age-matched controls, which are hypothesized to be correlated with different types of microstructural properties and can be measured by MRI sequences available for clinical studies and publicly available software. In addition to these main metrics of interest, we further evaluated DTI metrics (such as FA) as well as macromolecular tissue volume (MTV), which is derived from normalized proton-density mapping and provides a quantitative metric correlated with lipid volume fraction,^32,36–39 to enable comparisons with previous studies and improve understanding of microstructural changes occurring in the white matter tracts of patients with glaucoma. 

Materials and Methods

This study was approved by the ethics committee of the Jikei University School of Medicine, National Institute of Information and Communications Technology, and Tamagawa University. All subjects provided written informed consent to participate in the study. The study design followed the tenets of the Declaration of Helsinki. The data and code for reproducing figures and statistical analyses are available at https://github.com/htakemur/MultiContrastGlaucoma. 

Subjects

A total of 17 patients with glaucoma (eight females; mean age, 56.6 years; age range, 24–72 years) (Table 1) participated in this study. Thirty volunteers with normal visual function and no visual field defects also participated in this study as control subjects (14 females; mean age, 51.4 years; age range, 36–71 years) (Table 2). With reference to a previous study,⁴⁰ this sample size was predicted to be sufficient to identify a large effect size in group difference (d′ = 1.10) using a two-tailed, two-sample t-test (evaluated by G*Power 3.1.9.6). 

Table 1.

 

Patient Profiles

 

Table 2.

 

 Control Subject Profiles

 

Clinical Features of the Patients With Glaucoma

Patients with glaucoma who participated in this study were diagnosed with open-angle glaucoma (primary and secondary open-angle glaucoma) by experienced ophthalmologists at the Department of Ophthalmology, Jikei University School of Medicine. All patients underwent a comprehensive ophthalmologic examination, including measurement of visual fields using the Humphrey Field Analyzer (HFA) 24-2 or 30-2 Swedish interactive thresholding algorithm standard program by HFA (Carl Zeiss Meditec, Dublin, CA), and measurement of retinal nerve fiber layer thickness using the Cirrus HD-OCT system (Carl Zeiss Meditec, Dublin, CA) or Spectralis (Version 4; Heidelberg Engineering, Heidelberg, Germany). We note that optical coherence tomography (OCT) data were not used in subsequent statistical analyses because the data were acquired using the OCT system of two vendors, and it is known that retinal nerve fiber layer thickness measurements have vendor dependency.⁴¹ The diagnosis of glaucoma was based on glaucomatous optic neuropathy and visual field defects consistent with optic changes.⁴² Figure 1A shows the fundus photography and visual field loss of a representative patient with glaucoma. The results of HFA in all patients are described in Table 1 as a unit of mean deviation from healthy populations. We did not have a clinical record of disease onset, as most patients did not have subjective symptoms in the early phase of glaucoma. 

Several different subtypes of patients with glaucoma participated in this study, such as primary open-angle glaucoma, normal tension glaucoma (NTG), pseudoexfoliation glaucoma, and secondary open-angle glaucoma (Table 1). Except for NTG, all subtypes cause retinal ganglion cell damage due to high intraocular pressure (IOP). Even in NTG, lowering IOP is effective in preventing disease progression, and IOP can be a major cause of the disease.⁴³ Moreover, there is no evidence that glaucoma subtype is an independent risk factor for the progression of visual field loss.⁴⁴ In this study, we did not separately analyze these subtypes (1) to ensure high statistical power and (2) because these subtypes have similar causes of diseases. 

MRI Data Acquisition

MRI images were acquired using a MAGNETOM Trio, A Tim System (with a 32-channel head coil; Siemens, Erlangen, Germany) at the Tamagawa University Brain Science Institute, Machida, Japan. 

Structural MRI Data Acquisition

T1-weighted magnetization prepared rapid acquisition with gradient echo (MPRAGE) images were collected from all subjects (1-mm isotropic voxels; repetition time [TR], 2000 ms; echo time [TE], 1.98 ms). Acquisition of T1-weighted MRI data took 9 minutes 18 seconds per subject. T2-weighted sampling perfection with application optimized contrast using different flip angle evolution (SPACE) images were also acquired from 42 subjects (all 17 patients with glaucoma and 25 control subjects; 1-mm isotropic voxels; TR, 3000 ms; TE, 478 ms) to obtain T1-weighted/T2-weighted ratio maps (T1w/T2w).⁴⁵ Acquisition of T2-weighted MRI data took approximately 11 minutes 14 seconds per subject. 

dMRI Data Acquisition

dMRI data were collected from all subjects using single-shot spin-echo, echo planar imaging (EPI; 32 directions with b = 700 s/mm²; 64 directions with b = 2000 s/mm²; 1.7-mm isotropic voxels; TR, 4500 ms; TE, 94 ms; in-plane acceleration, 2; multiband factor, 3; phase partial Fourier, 6/8; diffusion scheme, monopolar) implemented in a multiband accelerated EPI pulse sequence provided by the Center for Magnetic Resonance Research, Department of Radiology, University of Minnesota; https://www.cmrr.umn.edu/multiband/).⁴⁶ Twelve non-diffusion-weighted (b = 0) measurements were also acquired. dMRI data were acquired with a single diffusion encoding design, as it is a conventional method and can be acquired by using widely available sequences. To minimize the impact of EPI distortion in subsequent analyses, two image sets were acquired with reversed phase-encoding directions (A–P and P–A). Acquisition of the dMRI data took approximately 9 minutes 11 seconds per subject. 

qT1 Data Acquisition

For all subjects, qT1 was measured following the protocols described in previous publications.^32,40 Four fast low-angle shot (FLASH) images with flip angles of 4°, 10°, 20°, and 30° (TR, 12 ms; TE, 2.41 ms) and isotropic 2-mm voxels were acquired. Five additional spin-echo inversion-recovery (SEIR) scans with an EPI readout (TR, 3000 ms; TE, 49 ms; 2× acceleration) were acquired to remove field inhomogeneities. The inversion times were 50, 200, 400, 1200, and 2400 ms. The in-plane resolution and slice thickness of the additional scans were 2 × 2 mm and 4 mm, respectively. Acquisition of qT1 data took approximately 13 minutes 30 seconds per subject. 

MRI Data Analysis

Structural MRI Data Preprocessing

T1-weighted images of individual subjects were aligned to the anterior commissure–posterior commissure (AC-PC) space. This aligned image on the AC-PC coordinate was then used as a common coordinate frame across the dMRI and qT1 datasets. T1w/T2w maps were also calculated by aligning the T2-weighted image with the T1-weighted image in the AC-PC space and calculating the ratio between the image intensities of the T1-weighted image and T2-weighted image for all 42 subjects who participated in the T2-weighted image acquisition. 

dMRI Data Preprocessing

dMRI data were preprocessed using the TOPUP and EDDY tools in FSL (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki) to correct for susceptibility-induced distortion and eddy-current artifacts.^47,48 The dMRI data were then aligned to the T1-weighted images in the AC-PC space using a 14-parameter constrained nonlinear coregistration algorithm. The tensor model was then fitted to the dMRI data using a least-squares algorithm to estimate FA, MD, axial diffusivity (AD), and RD. Furthermore, NODDI was also fitted to the dMRI data to obtain ICVF, ODI, and IsoV (Fig. 2B), using the NODDI MATLAB toolbox (http://mig.cs.ucl.ac.uk/index.php?n=Tutorial.NODDImatlab) using MATLAB (MathWorks, Natick, MA, USA). 

Figure 2.

 

(A) Visual white matter tracts identified by probabilistic tractography in representative subjects in each group (left, Ctl-008; right, Glc-004; magenta, OT; green, OR), overlaid on an axial slice of a T1-weighted image, located inferior to the tracts. (B) Tissue measurement maps (qT1, ICVF, ODI, and IsoV) in a representative healthy subject (Ctl-008).

qT1 Data Preprocessing

Both the FLASH and SEIR scans were processed using the mrQ software package (https://github.com/mezera/mrQ) to produce qT1 and MTV maps (Fig. 2B),³² and the qT1 maps were registered to the T1-weighted images using the FSL FLIRT tool.⁴⁹ 

Tract Identification and Evaluation

Tract Identification

We identified the OT and OR by analyzing the dMRI data using the same method as that used in previous studies (Fig. 2A).^40,50 Briefly, tracking was performed between regions of interest, which were defined manually (lateral geniculate nucleus [LGN]; see Supplementary Materials and Methods)⁴⁰ or by FreeSurfer segmentation (optic chiasm and primary visual cortex [V1]).^51–53 Further details of the tract identification methods are provided in the Supplementary Materials and Methods. 

Tissue Property Evaluation

We evaluated the tissue properties of OT and OR using Automated Fiber Quantification (AFQ) (https://github.com/yeatmanlab/AFQ).⁵⁴ Briefly, each streamline was resampled to 100 equidistant nodes, and tissue properties (qT1, ICVF, ODI, IsoV, FA, MD, AD, RD, and MTV) were calculated at each node of each streamline. The properties at each node were summarized by taking a weighted average of the tissue measurements on each streamline within that node. The weight of each streamline was based on the Mahalanobis distance from the tract core to minimize the impact of the partial volume effect with neighboring tissues and individual variabilities of the tract volume.⁵⁴ Data from the left and right hemispheres were averaged. We excluded the first and last 10 nodes because they are susceptible to crossing with U-fibers and partial voluming with gray matter. Although the results of the remaining 80 nodes are plotted as tract profiles (Figs. 3, 4), we averaged the data of 80 nodes to obtain a single-number summary of each metric per subject for statistical comparisons. For the dMRI metrics, we report the average results of the two runs. 

Figure 3.

 

Tissue properties along the OT. The four panels show the qT1 (left upper panel), ICVF (left lower panel), ODI (right upper panel), and IsoV (right lower panel) measurements along the OT. Profiles of individual patients with glaucoma are depicted as thin blue curves. Thick curves show the mean of each group (glaucoma, blue; healthy control, black). The lighter gray shades indicate a range of ±2 SD from the control mean, and the darker gray band shows ±1 SD from the control mean. The horizontal axis describes the normalized position along the tract (left, anterior; right, posterior).

Figure 4.

 

Tissue properties along the OR (left upper panel, qT1; left lower panel, ICVF; right upper panel, ODI; right lower panel, IsoV). Other conventions are identical to those of Figure 3.

Tissue Properties Controlled for Age

Linear regression was performed to predict the tissue properties (qT1, ICVF, ODI, and IsoV) along the OT and OR according to the ages of the subjects, which were pooled across both glaucoma and control groups. We then calculated the residual, which is the difference between the measured tissue properties and tissue properties predicted by age.⁵⁵ We used the residual as the tissue property controlled for age to evaluate the impact of age variability on the main analysis. 

Statistical Analyses

Group Comparison

We evaluated the effect sizes of the group differences in each MRI metric (qT1, ICVF, ODI, IsoV, FA, MD, AD, RD, and MTV) by calculating Cohen's d′. We evaluated the statistical significance of the group differences by performing a two-tailed, two-sample t-test for major metrics of interest (qT1, ICVF, ODI, and IsoV). We defined the significance level (α) as 0.05. 

Correlation Between Measurements in Patients With Glaucoma

We also evaluated between-metric correlations in the extent of the patient's tissue abnormality in each metric (qT1, ICVF, ODI, and IsoV) in the OT and OR. To this end, we first calculated the degree of deviation from the control mean in each patient with glaucoma in units of standard deviation of the control subjects. We then calculated the inter-patient Pearson correlation of each metric in the same tract and evaluated its statistical significance. We defined α as 0.008, which is equivalent to P = 0.05, with Bonferroni correction for six comparisons. 

Correlation Between MRI Measurements and Visual Field Test Results

We also assessed the correlation between tissue properties of the OT and the OR and the visual field test scores (mean deviation value of HFA) in patients with glaucoma. To this end, we normalized the MRI measurements of the patients with glaucoma by calculating the degree of deviation from the control mean, using units of standard deviation in the control subjects. The visual field test scores of each patient with glaucoma (Table 1) were averaged across the two eyes. We then calculated the Pearson correlation between the MRI measurements and the visual field test results in patients with glaucoma. In addition, we performed a supplementary analysis dividing the visual field test scores (total deviation of HFA) in each data point into left and right visual fields. We then calculated the correlation between the mean total deviation in each hemifield and the tissue properties of the OT and OR in the contralateral hemisphere. 

Results

We successfully identified the OT and OR in all hemispheres (for representative subjects, see Fig. 2A). We found that the estimated volume of the OT was significantly larger in the glaucoma group compared with that in the control group (mean estimated OT volume in glaucoma, 621.65 mm³; mean estimated OT volume in control, 368.87 mm³; d′ = −1.97; t₄₅ = −6.84; 95% confidence interval [CI], −327.20 to −178.37; P = 0.00000002). In contrast, estimated OR volume was not significantly different across groups (mean estimated OR volume in glaucoma, 7256.79 mm³; mean estimated OR volume in control, 7024.93 mm³; d′ = −0.07; t₄₅ = −6.85; 95% CI, −8737.53 to 8273.81; P = 0.96). Group differences in the estimated OT size are difficult to interpret and most likely occur because generating coherent OT streamlines in glaucoma group data is difficult due to reduced diffusion anisotropy, rather than true anatomical differences in the tract size. To minimize the impact of apparent tract volume difference, in subsequent analyses we focused on the tract profiles of four tissue measurements (ICVF, ODI, IsoV, and qT1) along the OT and OR, which were calculated by weighted average based on the distance from the tract core (see Materials and Methods). 

Glaucoma Affected All Types of Tissue Measurements in the OT

Figure 3 depicts the profile of OT in the control (black) and glaucoma (blue) subjects. The gray-shaded regions in Figure 3 show the standard deviation from the control mean for qT1 (left upper panel), ICVF (left lower panel), ODI (right upper panel), and IsoV (right lower panel), along the length of the OT. The individual blue curves depict data from individual patients with glaucoma. We found significant group differences in all measurements; the patients with glaucoma showed significantly higher qT1 (two-tailed, two-sample t-test d′ = −0.93; t₄₅ = −3.05; 95% CI, −0.09 to −0.02; P = 0.004), significantly lower ICVF (d′ = 1.49; t₄₅ = 5.14; 95% CI, 0.06–0.13; P = 0.000006), significantly higher ODI (d′ = −2.01; t₄₅ = −6.85; 95% CI, −0.08 to −0.04; P = 0.00000002), and significantly lower IsoV (d′ = 1.55; t₄₅ = 5.39; 95% CI, 0.05–0.12; P = 0.000002). These differences were preserved after controlling for age (Supplementary Fig. S1). 

Glaucoma Affected ICVF in the OR, But Not Other Measurements

Figure 4 shows the profile of the OR in the control and glaucoma subjects. We observed significantly lower ICVF in patients with glaucoma (d′ = 1.13; t₄₅ = 3.80; 95% CI, 0.02–0.06; P = 0.0004) (Fig. 4). By contrast, we did not find significant differences in qT1 (d′ = −0.12, t₄₅ = −0.40; 95% CI, −0.03 to 0.02; P = 0.69), ODI (d′ = −0.42; t₄₅ = −1.38; 95% CI, −0.02 to 0.00; P = 0.18), or IsoV (d′ = 0.19; t₄₅ = 0.66; 95% CI, −0.01 to 0.01; P = 0.52). Similar results were obtained when the tissue properties were controlled for age (Supplementary Fig. S2). 

Do Multiple MRI Measurements Detect Similar Types of Tissue Abnormalities Within the Same Tract?

For each MRI metric along the OT, we quantified the deviation from the control mean for each patient with glaucoma and calculated the correlation between them. We found a significant correlation between the extent of abnormalities in ICVF and IsoV (R = 0.77; P = 0.0003) (Fig. 5A). The correlations between other pairs did not reach statistical significance (α = 0.008; qT1-ICVF, R = −0.49, P = 0.05; qT1-ODI, R = 0.51, P = 0.04; qT1-IsoV, R = −0.60; P = 0.01; ICVF-ODI, R = −0.12, P = 0.65; ODI-IsoV, R = −0.57, P = 0.02), suggesting that these metrics may reflect multiple underlying sources of tissue abnormalities in the OT (see Supplementary Table S1 for the correlations among MRI metrics in healthy controls). 

Figure 5.

 

Correlations between multiple MRI measurements in patients with glaucoma. The two-dimensional scatterplots depict the extent of deviation of each individual patient with glaucoma (N = 17) from the control mean, with units of 1 SD of the control group. Individual dots are data points from individual patients with glaucoma. Black thick lines depict regression lines, and cyan curves indicate the 95% confidence intervals of linear regression estimated by the bootstrapping method. (A) Correlation between ICVF and IsoV in the OT. (B) Correlation between ICVF and qT1 in the OR.

We also investigated the correlation between the extent of abnormalities in each MRI metric along the OR in patients with glaucoma. We found a significant correlation between qT1 and ICVF (R = −0.64; P = 0.005) (Fig. 5B), whereas other pairings were not significantly correlated (qT1-ODI, qT1-IsoV, ICVF-ODI, ICVF-IsoV, and ODI-IsoV: R = −0.33, 0.46, 0.39, −0.11, and 0.00, respectively; P = 0.19, 0.06, 0.13, 0.69 and 1.00, respectively). A similar negative correlation was found between qT1 and ICVF in healthy controls (Supplementary Table S2). 

Relationship Between White Matter Tissue Measurement and Visual Field Loss

We assessed the extent to which this inter-patient variability in MRI measurements could be explained by the severity of visual field loss by evaluating the correlations between MRI measurements and visual field test data (Table 1) among patients with glaucoma. In the OT, the ODI was significantly correlated with the visual field test (R = −0.61; P = 0.009) (Supplementary Fig. S3A). This negative correlation remained after controlling the ODI for differences in age (R = −0.51) (Supplementary Fig. S3B). However, qT1, ICVF, and IsoV along the OT did not show significant correlations with the visual field test (R = −0.30, 0.30, and 0.39; P = 0.25, 0.24, and 0.13, respectively, for qT1, ICVF, and IsoV). None of the MRI measurements along the OR was significantly correlated with the visual field test (R = −0.14, 0.31, −0.03, and 0.01; P = 0.59, 0.23, 0.92, and 0.97, respectively, for qT1, ICVF, ODI, and IsoV). 

We also performed a supplementary analysis to investigate whether tissue abnormalities in the ODI along the OT could be related to the pattern of visual field loss. To this end, we compared the visual field test results in the left and right visual fields with the ODI along the contralateral OT. However, we found that, after separating visual fields and hemispheres, correlation between the visual field test and ODI along the OT became smaller (R = −0.31 and −0.50 for left OT/right visual field and right OT/left visual field, respectively). We speculate that this result may be due to limitations in the signal-to-noise ratio of the dMRI measurements along the OT or spatial precision of the visual field test performed in a clinical setting. 

Comparisons Across Multiple MRI-Based Tissue Measurement Metrics

Finally, we also evaluated other MRI metrics, such as DTI-based metrics (FA, MD, RD, and AD) and MTV, to enable comparisons with previous studies using these metrics (for profiles of the OT and OR in these metrics, see Supplementary Figs. S4 and S5). Figure 6 shows how each MRI metric in individual patients with glaucoma deviated from the control mean. In this plot, the vertical axis indicates the effect size of the difference (d′) between individual patients with glaucoma and the control mean for each metric. We found that FA was lower in patients with glaucoma in both the OT and OR (effect size of group difference, d′ = 2.39 and 1.50 in OT and OR, respectively). We also found higher MD in patients with glaucoma, although the difference in the OT was not consistent across all patients (d′ = −0.39 and −1.03 in OT and OR, respectively). In both OT and OR, we found much higher RD in the patients with glaucoma (d′ = −1.63 and −1.44, respectively), although AD showed a smaller effect in OT and inconsistent results across patients in OR (d′ = 1.09 and −0.23, respectively). Therefore, profound differences in FA can be mostly explained by higher RD rather than lower AD. 

Figure 6.

 

Evaluation of different MRI metrics derived from the diffusion tensor model, quantitative MRI, and NODDI (left panel, OT; right panel, OR). The vertical axis of each plot represents the extent to which individual patients with glaucoma (color bars; Glc001-017) deviate from the control mean (N = 30). The vertical axis unit is the effect size (d′) of the difference between individual patients with glaucoma (each colored bar) and the control mean.

In the OT, the patients with glaucoma consistently showed lower MTV (d′ = 1.28), which is another type of quantitative MRI method based on calibrated proton density measurements.³² In contrast, in the OR, we did not find consistently lower MTV among patients with glaucoma (d′ = 0.02), similar to the observations for qT1 (d′ = −0.12). Therefore, similar to our previous study on Leber's hereditary optic neuropathy (LHON),⁴⁰ qT1 and MTV showed abnormalities in the OT and a lack of abnormalities in the OR. 

We also evaluated the extent to which DTI parameters can provide information similar to NODDI parameters. To this end, we calculated the correlation of the extent of deviation of patients with glaucoma from control mean between NODDI and DTI parameters (Supplementary Tables S3 and S4). Notably, ICVF showed negative correlation with RD (R = −0.81 and −0.88 in OT and OR, respectively) whereas ODI showed a negative correlation with AD (R = −0.90 and −0.85 in OT and OR, respectively). 

We also evaluated the T1w/T2w ratio, which has been used in recent neuroimaging studies (Supplementary Fig. S6).^45,56,57 In the OT, patients with glaucoma showed slightly lower T1w/T2w ratios than the controls (effect size of group difference, d′ = 0.63); however, in the OR, a group-level difference was almost absent (d′ = 0.04). In both tracts, two patients with glaucoma showed higher T1w/T2w ratios than the controls (Supplementary Fig. S5; Glc-013 and Glc-014). However, because these two patients did not show clear opposite trends in dMRI and quantitative MRI metrics (Fig. 6), it is unclear whether these trends reflect a true tissue difference or measurement bias. 

Discussion

Previous anatomical studies on patients with glaucoma and non-human primate models revealed that glaucoma causes tissue changes in the LGN^58–61 and V1,⁵⁸ indicating that glaucoma affects tissues of the subsequent visual areas receiving inputs from retinal ganglion cells. More recently, neuroimaging studies have demonstrated abnormalities in diffusivity measurements (such as FA) along the OT and OR in patients with glaucoma.^{7,11,12,18,19,62–77} However, the types of microstructural changes that occur as a consequence of glaucoma are not yet understood, as DTI-based metrics are not specific to myelin damage or axonal loss.^21–23 

We combined NODDI and qT1 to overcome the limitations of DTI. We found that all four MRI metrics (qT1, ICVF, ODI, and IsoV) (Fig. 3) indicated abnormalities in the OT of patients with glaucoma. This suggests that OT tissue abnormalities in glaucoma may be caused by multiple factors. In contrast, only ICVF indicated abnormalities in the OR of patients with glaucoma (Fig. 4), which suggests that tissue abnormalities in the OR are associated with more specific neurobiological mechanisms. Because it is natural to interpret ICVF in white matter voxels as a measurement of apparent axonal density,⁷⁸ our finding may indicate that glaucoma causes axonal damage in the OR rather than myelin damage. This interpretation is also in line with a previous study demonstrating that diffusivity measurements are sensitive to Wallerian degeneration, in which axonal loss has been considered as one of the main determinants of tissue damage.⁷⁹ Because the OT is a direct continuation of the optic nerve and the OR is one synapse away from retinal ganglion cells, distinct results between OT and OR suggest that direct damage and transsynaptic effects cause different types of microstructural abnormalities. However, these interpretations remain speculative because histological validation of the ICVF along the OR remains to be performed. In the future, it is important to establish a theory on how different types of degeneration (anterograde, retrograde, transsynaptic, and Wallerian degeneration) affect each MR metric along the visual pathways. 

We also found that ICVF was correlated with qT1 in the OR of patients with glaucoma (Fig. 5B), whereas, on average, qT1 did not show any abnormalities (Fig. 4). This result indicates that patients with glaucoma with severe ICVF abnormalities tend to exhibit some degree of abnormality in qT1. Although speculative, these results suggest the possibility that, in the OR of patients with glaucoma, axonal damage (quantified by ICVF) occurs earlier than myelin damage (quantified by qT1). Myelin damage may occur at a later stage of disease progression. In contrast, You et al.¹⁸ proposed that myelin damage in the OR occurs in an early stage of glaucoma based on a correlation with latency of visual evoked potential (VEP), although it is not clear whether VEP latency can be a specific marker solely for assessing myelin levels, as conduction velocity is also dependent on axonal properties.⁸⁰ It is necessary to perform future longitudinal studies or preclinical MRI studies on animal models,⁶⁰ in which control of the duration since disease onset is practical, to clarify the time course of axonal and myelin damage in the OR of patients with glaucoma. 

Previous studies have shown that tissue abnormalities have been observed not only in glaucoma but also in LHON, which also damages retinal ganglion cells.^40,81 Abnormalities in the DTI metric along the OT and OR have also been reported in disorders causing photoreceptor damage (retinitis pigmentosa⁸² and macular degeneration^81,83–85) and amblyopia.^86–88 Therefore, OT and OR tissue changes are not specific to disorders causing retinal ganglion cell damage but rather are a generalized phenomenon among visual disorders. However, it is not clear whether the present findings based on qT1 and NODDI can be generalized to photoreceptor diseases and amblyopia, as these previous studies used a simplistic DTI approach rather than a multi-contrast approach. The extension of the approach used in this study to populations with other types of disease will clarify the similarity and differences in microstructural damages occurring in the OT and OR across different types of disorders. 

Takemura et al.⁴⁰ used a similar combination of MRI methods (DTI and qT1) to evaluate white matter tissue properties in patients with LHON. Although their study did not use NODDI, we can compare our DTI results obtained from patients with glaucoma (Fig. 6) with those obtained from patients with LHON. This comparison revealed consistent findings: (1) both diffusivity and qT1 demonstrated abnormalities in the patients’ OT, (2) abnormal diffusivity measurements were made from the patients’ OR, and (3) qT1 measurements in the patients’ OR showed no evidence of abnormality. These results suggest that main findings may be generalizable to LHON and glaucoma, both of which cause retinal ganglion cell damage. However, the data from patients with glaucoma showed a notable difference from the data from patients with the LHON.⁴⁰ Specifically, although the LHON study did not show a large abnormality in RD along the OT, the glaucoma data acquired in this study did (Fig. 6), suggesting that glaucoma and LHON may cause different types of tissue changes. We speculate that these differences may be related to differences in pathology between glaucoma and LHON. Specifically, although visual field damage is specific to the fovea in LHON,⁸⁹ patients with glaucoma often exhibit damage within a wide range of visual fields. Therefore, such differences in visual field damage may explain the differences in diffusivity along the OT. In addition, the differences in disease progression may affect diffusivity along the OT; although dendrite and soma damage precedes axonal damage in glaucoma,^90–92 it is unclear whether a similar order of progression occurs in LHON. 

This study has several limitations. First, although the sample size used in this study was sufficient to identify large group differences between glaucoma and control subjects, statistical power was still limited when evaluating correlations between MRI measurements and severity of visual field loss in patients with glaucoma. Because visual field loss in glaucoma varies across patients, it is also difficult to perform more direct comparisons between the type of visual field loss and tissue properties of subcomponents of the OR^40,85 unless we obtain a very large sample size and divide patients into different subgroups. Second, because this study focused entirely on open-angle glaucoma, the results may not be fully generalizable to other types of glaucoma populations.⁹³ Third, it remains unclear how the time from the onset of glaucoma is related to the degree of abnormality in white matter tissue properties. Given that most of the patients did not have subjective symptoms of visual field loss in the early phase of glaucoma, an accurate estimation of disease onset was very difficult. To address this question, preclinical MRI studies on animal models⁶⁰ are essential. Fourth, although the NODDI metrics and qT1 are hypothesized to be relatively specific markers of microstructural properties,^{26,27,30,33–35} they do not have a complete one-to-one relationship with the properties of myelin or axons.^26,94 Specifically, NODDI has limitations regarding uncertainty in parameter estimations⁹⁵ and dependency on the selection of assumed parameters.⁹⁶ In addition, although qT1 and other types of quantitative MRI maps are sensitive to myelin content, they are not fully specific markers for myelin levels,⁹⁷ most likely because other types of microstructural properties also impact qT1.⁹⁴ A possible way to overcome these limitations is to acquire dMRI data with advanced diffusion encoding^98,99 or to build a model to predict the fractions of certain types of lipids from multiple quantitative MRI metrics.¹⁰⁰ Fifth, it will be important to extend this work by further testing MRI measurements in cases of congenital glaucoma and glaucoma with full visual field loss to test the extent to which the findings of this study can be generalized to these cases. Finally, it is important to understand the extent to which MRI-quantifiable white matter tissue changes are reversible. This requires a longitudinal MRI study to evaluate the effect of specific treatment strategies on MRI-based metrics. Although such a study is not trivial to perform, we hope that the extension of a quantitative multi-contrast MRI approach will continue to improve our understanding of the white matter consequences of glaucoma. 

Acknowledgments

The authors thank Yusuke Sakai for supporting the data analysis. 

Supported by Japan Society for the Promotion of Science KAKENHI (JP20K18396 to SO; JP17H04684 and JP21H03789 to HT; JP18K16939 and JP21K09729 to HH; JP20K22273 and JP21K12617 to AM; JP19K09982 to YM). 

Disclosure: S. Ogawa, None; H. Takemura, None; H. Horiguchi, None; A. Miyazaki, None; K. Matsumoto, None; Y. Masuda, None; K. Yoshikawa, None; T. Nakano, None