Age and Glaucoma-Related Characteristics in Retinal Nerve Fiber Layer and Choroid: Localized Morphometrics and Visualization Using Functional Shapes Registration

Optical coherence tomography provides high-resolution 3D imaging of the posterior segment of the eye. However, quantitative morphological analysis, particularly relevant in retinal degenerative diseases such as glaucoma, has been confined to simple sectorization and averaging with limited spatial sensitivity for detection of clinical markers. In this paper, we present point-wise analysis and visualization of the retinal nerve fiber layer and choroid from cross-sectional data using functional shapes (fshape) registration. The fshape framework matches two retinas, or generates a mean of multiple retinas, by jointly optimizing the surface geometry and functional signals mapped on the surface. We generated group-wise mean retinal nerve fiber layer and choroidal surfaces with the respective layer thickness mapping and showed the difference by age (normal, younger vs. older) and by disease (age-matched older, normal vs. glaucomatous) in the two layers, along with a more conventional sector-based analysis for comparison. The fshape results visualized the detailed spatial patterns of the differences between the age-matched normal and glaucomatous retinal nerve fiber layers, with the glaucomatous layers most significantly thinner in the inferior region close to Bruch's membrane opening. Between the young and older normal cases, choroid was shown to be significantly thinner in the older subjects across all regions, but particularly in the nasal and inferior regions. The results demonstrate a comprehensive and detailed analysis with visualization of morphometric patterns by multiple factors.

Introduction

Volumetric optical coherence tomography (OCT) has emerged as a preferred diagnostic tool in ophthalmology for noninvasive, in-vivo, micrometer-resolution imaging of the eye. Recent progress in OCT imaging has allowed the acquisition of highly detailed 3D images from which morphometric measurements can be derived. Optic nerve head measurements and the thickness of the peri-papillary retinal layers have been used clinically to detect and monitor glaucoma progression (Leung, 2014). While these measurements are useful individually, a lack of clear spatial references for biomarkers limits the spatial and anatomical correspondence across multiple images. For example, in the common sectoral layer thickness analysis of OCT scans, the peri-papillary area is split into quadrants (superior, inferior, nasal, and temporal) which are then further subdivided into circumferential areas. This analysis relies on averaging over the local sectors to reduce noise, and to mitigate potential inconsistency in sectoral placement across individuals. Such a sectoral averaging approach is limited by a minimum size of sectors to achieve comparisons in the same vicinity in the different individuals. Hence, the sectorization approach reduces the spatial sensitivity of the measurements due to averaging over larger areas and could potentially impact clinical assessment. This motivates the need to develop tools that can generate measurements on a point-to-point basis, eliminating the need to average data in local regions.

Most previous studies involving registration of OCT images have averaged multiple serially acquired images for noise reduction and motion correction (Jørgensen et al., 2007; Young et al., 2011), or utilized rigid alignment of time-course images (Niemeijer et al., 2009). Relatively little attention has been given to registration of cross-sectional OCT data. Chen et al. (2014) performed intensity-based non-rigid registration of macular OCT scans by a combination of rigid alignment of foveae, and A-scan-wise affine and non-rigid registration using radial basis functions for refined alignment of the retinal layers (Chen et al., 2014). A more recent work (Anthony et al., 2016) by the same group applies this registration technique to perform voxel based morphometry in macular OCT of healthy controls and multiple sclerosis patients. Our group's approach has focused on retinal surface-based registrations and atlas generation. In (Gibson et al., 2010), 3D optic cup surfaces were registered to a single template surface, first by rigid and nonrigid intensity-based volumetric registration, followed by spherical mapping and spherical demons registration of the surfaces. This work was further expanded upon in Lee et al. (2015) which represented the retinal surfaces utilizing the framework of mathematical currents. Two surfaces were brought into close proximity by minimizing a functional of reproducing kernel Hilbert space norm-based energy and a dissimilarity term, then registered by spherical demons to establish homology. More recently, we introduced the functional shape (fshape), framework (Charlier et al., 2015; Lee et al., 2017). In this framework, the retinal surface (shape or geometry) and any value mapped on the surface, for example, retinal layer thickness (function or signal), are considered together as a single mathematical object called fshape. One fshape can be matched to another or the mean of multiple fshapes can be computed by joint optimization of the energies associated to varifold-based dissimilarity measures of geometry and function. For group analyses, the algorithm can generate population atlases and establish homology across the database, facilitating comparison of morphometric measurements in localized regions. Moreover, fshape mean computation or matching can be performed with flexible constructions of fshapes that can include multiple geometric shapes and function parameters; this allows the building of specific sets of geometry and function features to compare across multiple groups.

In this paper, we demonstrate the use of this algorithm for investigating the effect of age and glaucoma on retinal nerve fiber layer (RNFL) and choroid. Loss of RNFL in the optic nerve head (ONH) region is a well-known hallmark of glaucoma that leads to irreversible vision loss (Medeiros et al., 2005). Currently, the RNFL thickness profile along a circular scan centered at the ONH and the sectoral average thickness maps are used in clinics to assess the disease progression. Studies have shown regional patterns in glaucomatous RNFL thinning, with most significant changes in the inferior peripapillary region (Leung et al., 2010; Mwanza et al., 2011). Aging has been also associated with RNFL loss (Budenz et al., 2007; Parikh et al., 2007). The effect of glaucoma in the choroid has been more debated, with some studies reporting glaucoma-related thickness changes (Song et al., 2016; Li et al., 2017), and others reporting no changes (Ehrlich et al., 2011; Maul et al., 2011). Recent works on OCT angiography (Lee E. J. et al., 2016; Mammo et al., 2016) suggest vascular impairment in glaucoma, and this motivates simultaneous, localized analysis of the two layers to investigate possible connection in the structural changes due to glaucoma.

In this work, we aim to (i) examine the spatial RNFL thickness patterns by age and by presence of glaucoma by comparing the reference group of older healthy eyes with younger healthy eyes and with age-matched glaucoma eyes, and (ii) study whether there is spatial relationship between age-related changes and glaucoma-related changes.

Materials and Methods

Participants and Image Acquisition

Thirty-eight eyes from five young healthy participants, five older healthy participants, and twelve older glaucoma patients were included in the study. The participant demographics are listed in Table 1. Before being included in the study, each participant was subject to optic nerve head OCT imaging, dilated stereoscopic examination of the optic nerve, stereo disc photography, intraocular pressure (IOP) measurement, and reproducible Humphrey perimetry at the Eye Care Center at Vancouver General Hospital. Eyes with retinal disease other than primary open-angle glaucoma, uveitis, IOP lower than 10 mmHg or greater than 20 mmHg, or optic neuropathy from causes other than glaucoma were excluded. The mean glaucoma duration at the time of imaging was 3.69 ± 3.80 years. A custom 1,060-nm swept-source OCT system by the Biomedical Optics Research Group at SFU was used for imaging. Each image consisted of 400 B-scans, with 400 A-scans per B-scan and 1,024 pixels per A-scan. The axial voxel resolution was 2.7 μm, the axial coherence length was ~6 μm, and the lateral pixel resolution ranged from 11.9 to 14.5 μm depending on the eye's axial length. The A-scan rate of 100 kHz resulted in ~1.6 s of acquisition time per volume.

TABLE 1

Table 1. Participant demographics.

Preprocessing, Segmentation, and Layer Thickness Measurement

Images with artifacts, such as large lateral motion artifact or the ONH not being at the center of the field of view were excluded from this analysis. The image underwent axial motion correction by B-scan cross-correlation and 3D bounded-variation smoothing for reducing the effect of speckles and enhancing the visibility of anatomical structures, with no additional normalization. An example of OCT B-scan before and after processing is shown in Figures 1a,b. Retinal nerve fiber layer (RNFL) and choroid were segmented automatically by delineating inner limiting membrane (ILM), RNFL-ganglion cell layer boundary, Bruch's membrane (BM), and the choroid-sclera boundary using a 3D graph-cut based algorithm (Li et al., 2006; Lee et al., 2013a,b). All automated segmentation results were checked by a trained rater, and incorrect segmentation was manually corrected using Amira (version 5, FEI). Bruch's membrane opening (BMO) was manually segmented on 80 radial slices extracted from the image volume. An example of the segmentation is shown in Figures 1c,d. A best-fit 3D BMO ellipse was computed using principal component analysis and least square fitting. The segmented RNFL and choroid were cropped at 0.25 mm from the BMO ellipse to account for the ambiguity in the retinal layer boundary close to the optic cup (Lee et al., 2014). Recent studies reported on inconsistencies resulting from using the conventional optic cup as a reference due to its ambiguous definition based on 2D projection fundus images, and showed the BMO was an viable alternative reference as a robust anatomical structure defined in 3D space (Chauhan and Burgoyne, 2013; Gardiner et al., 2014). The layer thickness was measured at each point as the closest 3D Euclidean distance between the posterior and anterior surfaces of the layer. Prior to the fshape registration step, all corresponding surfaces were rigidly aligned by matching the BMO ellipse centroid.

FIGURE 1

Figure 1. Example OCT B-scan (a) without processing, (b) after BV smoothing, and (c) with the segmentation of inner limiting membrane (purple), posterior boundary of retinal nerve fiber layer (cyan), Bruch's membrane (yellow), Bruch's membrane opening (green), and choroid-sclera boundary (red). 3D surface generation of the retinal layer boundaries is shown in (d).

Fshape Registration

The framework of functional shapes (or fshapes) provides a quantitative measure of inter-subject variability in the RNFL and choroid. In this section, we briefly summarize the algorithm which is detailed in Lee et al. (2017). Let the ith subject's RNFL or choroid thickness be represented by (Xⁱ, fⁱ), where X is the layer surface (geometry) and f is the surface-indexed function (thickness here) mapped on X. Let a template exemplar for this database be denoted by (X_∗, f_∗), consisting of a template surface geometry X_∗ and an associated surface-indexed signal (thickness here) mapping f_∗. Given (X_∗, f_∗) as the template fshape, and (Xⁱ, fⁱ) as the ith target fshape to be registered to the template, the fshape framework will estimate a smooth deformation ϕⁱ of the template geometry X_∗ and a residual ζⁱ to be added to the template function f_∗, such that after transformation of the template fshape with the mapping ϕⁱ to the ith target fshape, the geometry of the transformed template matches the geometry of the ith target $X^i\approx {\varphi }^i(X_{*})$ and the transformed template function plus residual matches the function of the ith target i.e., $f^i\approx \left(f_{*}+{\zeta }^i\right)◦{{\varphi }^i}^{-1}$. Hence, for each of the target fshapes (Xⁱ, fⁱ) for all i = 1..N a pair (ϕⁱ, ζⁱ) consisting of a smooth deformation and a function residual are estimated such that:

$\begin{array}{cc}\begin{array}{l}\left(X^i,f^i\right)\approx \left({\varphi }^i,{\zeta }^i\right)·\left(X_{\ast },f_{\ast }\right)=\left({\varphi }^i\left(X_{\ast }\right),\left(f_{\ast }+{\zeta }^i\right)◦{({\varphi }^i)}^{-1}\right)\\ =\left({\tilde{X}}^i,{\tilde{f}}^i\right)\end{array}& \left(1\right)\end{array}$

The function residual added to the template function estimated by fshape matching effectively becomes the representative of the target function in the geometry of the template. By fixing a template geometry and function, a group of target fshapes can be brought into the coordinates of the template such that the residuals representing different target function values are now indexed on the same template geometry. The choice of the template is an important consideration and the mean of the observations in the database is a standard choice. A database fshape mean is hence estimated, to be used as a template, by an adaptive gradient descent algorithm that jointly optimizes a total fshape dissimilarity measure between each target fshape (Xⁱ, fⁱ) and the transformed template $\left({\tilde{X}}^i,{\tilde{f}}^i\right)$ for all i taken together. The process is summarized in Figure 2. A similar mean template generation procedure was repeated for the choroid. This maps all the target observations' function values (layer thickness) into a common coordinate system of the mean template where statistical analysis can be applied on the residuals (ζⁱ) to compute the point-wise variability of the function values at each point on the template geometry across the database.

FIGURE 2

Figure 2. Computation of fshape mean template.

Sectorization

In order to compare the fshape analysis with more conventional “intrinsic” analysis where an individual-specific coordinate system is placed on each geometry separately such as using sectoral analysis, the peri-papillary retinal layers were divided into regional sectors (Lee et al., 2014) as shown in Figure 3. Unlike standard sectors defined by fixed distances from the center of the optic disk on the enface projection image that do not take into account the individually varying sizes of Bruch's membrane opening (BMO), the sectors in this study were defined in 3D in each eye by the distance from the BMO, which is a more reliable anatomical landmark than the optic disk (Lee et al., 2014). This normalizes the sectors for different retinal tilts and BMO/optic disks sizes. The sectors were first delineated by elliptical annuli at constant distances (0.25, 0.75, 1.25 mm) from the BMO ellipse. These were further divided by superior, nasal, inferior, and temporal sectors, and additionally into superior-nasal (SN), inferior-nasal (IN), inferior-temporal (IT), and superior-temporal (ST) sectors. The first four sectors are 60° wide and the latter are 30° wide. For each sector, the thickness measurements for all points in the sector for that individual eye are taken and averaged to create one scalar number representing the average sectoral thickness value.

FIGURE 3

Figure 3. Sectorization of the layer surface with reference to the Bruch's Membrane Opening (BMO) as reference.

Statistical Analysis

Group Analysis

Two analyses were conducted (1) the fshape analysis of computing the mean fshape geometry and function and the residual function for each subject indexed on the common mean template, enabling point-wise comparison across the subjects in the database and (2) sectoral averages within each subject that enable comparison across subjects by the sectors defined over the subject's layer geometry.

The database consisted of members from three groups: young normal (Group A), older normal (Group B), and older glaucomatous (Group C) individuals. These groups enable analyses of two main questions (1) the effect of age on RNFL and choroid layer thicknesses in healthy young and older individuals (Group A and B comparison) and, (2) the effect of glaucoma between age-matched individuals (Group B and Group C comparison). These two questions were analyzed by point-wise (fshape) and sector-wise group mean thickness difference maps and two-sample t-test maps.

Regression Analysis

The effect of age and glaucoma on RNFL and choroidal thicknesses was additionally examined by point-wise and sector-wise linear regression to quantify trend, on Group A and B for the effect of age, and Group B and C for the effect of glaucoma. Each layer thickness (RNFL or choroid) values from Group A and B were fitted to layer thickness = a*Age + b to estimate the rate of change (mm/year) in the cross-sectional healthy subjects. To assess change as a function of visual field mean deviation (VFMD), a measure of glaucomatous loss measured in decibels (dB), the layer thickness (RNFL or choroid) values from Group B and C were fitted to layer thickness = a*VFMD + b to estimate the rate of change per VFMD (mm/dB).

Point-Wise Visualization of RNFL and Choroid Thickness Pattern

To visualize the measurements across the database highlighting the variability and trends, the fshape measures for each subject on the common template were normalized by using a z-score computed by subtracting the mean of a reference group and dividing by the standard deviation of the measures in the reference group. The z-score was calculated point-wise as ${\text{z}}_{\text{k}}=({\text{x}}_{\text{k}}-{\bar{\text{x}}}_{\text{k}}^{\text{h}})/{\sigma }_{\text{i}}^{\text{h}}$, where x_k is the thickness at k^th point for, ${\bar{\text{x}}}_{\text{k}}^{\text{h}}$ is the mean thickness of the reference group at k^th point, and ${\sigma }_k^h$ is the standard deviation of the reference group at k^th point. The reference groups were chosen to be Group A for age comparison and Group B for glaucoma comparison such that the average measures of the young normal (Group A) or older healthy (Group B) subjects are removed to visualize the residual effects of age and glaucoma, respectively. To present a compact visualization, measurement of each subject was unraveled into a column format by subdividing the mean template surface into sectors, and subdividing each sector into smaller sub-sectors, and arranging the z-score values by their sectors in a column format consistent across the database while preserving spatial adjacencies.

Results

This section will present the results of point-wise (using fshape) and sectoral analysis of RNFL and choroid thickness across the three cohorts chosen for this study. All figures are right-eye oriented, with the left side temporal and the right side nasal. Group averages of RNFL thickness are visualized in Figure 4. The top row shows the point-wise average computed using fshapes and the bottom row shows the results using sectoral averaging. The point-wise fshape mean RNFL templates display detailed salient features from multiple, cross-sectional eyes in each group showing the characteristic hourglass pattern of healthy RNFL in both Group A (young normal) and Group B (older normal) whereas Group C (older glaucoma) show marked glaucomatous thinning. There is good correspondence between the fshape mean templates constructed from point-wise registration and sectoral average maps taken from unregistered RNFL thickness averages in each sector in each individual.

FIGURE 4

Figure 4. Retinal nerve fiber layer (RNFL) thickness using fshapes point-wise registration (top row) and sectorization (bottom row) for Group A (Young Normal), Group B (Older Normal), and group C (Older Glaucoma). The RNFL thickness in Group B is similar to Group A, indicating that RNFL is relatively better preserved with age, whereas in Group C, this layer undergoes marked thinning. All images are in the right-eye orientation.

The group averages of choroidal thickness by fshape (top row) and sectoral averaging (bottom row) is shown in Figure 5. All three groups display thicker choroid in the nasal and superior regions and thinner choroid in the inferior region. The choroid is visibly thinner with age as seen in Group B compared to Group A. As with RNFL thickness, there is overall correspondence between the fshape mean templates and sectoral average maps.

FIGURE 5

Figure 5. Choroid thickness using fshapes (top row) and sectorization (bottom row) for Group A (Young Normal), Group B (Older Normal) and Group C (Older Glaucoma). The choroid in the older normal subjects is markedly thinner than the choroid in the younger normal subjects revealing a thinning process that seems to be based on chronological age. All images are in the right-eye orientation.

The effect of age and glaucoma in RNFL thickness is shown by difference of group averages in Figure 6. The top row shows the effect of age in RNFL thickness in the healthy subjects comparing the mean RNFL thickness in young normal vis-a-vis the older normal individuals. The bottom row shows the effect of glaucoma in RNFL thickness by comparing age-matched older subjects with and without glaucoma. The left panel shows the point- and sector-wise difference of RNFL thickness from Group B by Group A (top, young vs. older) and from Group C by Group B (bottom, healthy vs. glaucoma) at each corresponding points / sectors. The right panel shows point and sector-wise t-test results indicating where the group difference is significant. The RNFL thickness is found to not change significantly over age across individuals (top row), whereas the difference due to glaucoma is apparent (bottom row). The loss of RNFL thickness is observed to be higher in regions where normal RNFL thickness is higher and suggests some regional correspondence between the degree of glaucomatous RNFL loss and the original RNFL thickness. The t-test significance map between the healthy and glaucoma subjects shows the highest statistical significance in the inferior region.

FIGURE 6

Figure 6. Effect of age (top row) and glaucoma (bottom row) on RNFL thickness. The point- and sector-wise subtraction of mean RNFL thickness between older normal (Group B) and young normal (Group A) is shown in the top left panel, and the p-values from point- and sector-wise t-test is shown in the top right panel. The point- and sector-wise subtraction of RNFL mean between older glaucoma (Group C) and older normal (Group B) is shown in the bottom left panel, and the p-values from point- and sector-wise t-test is shown in the bottom right panel. The RNFL layer is found not to change significantly with age, whereas it changes significantly with glaucoma. The fshapes point-wise comparison shows the pattern of change in greater detail than the sectorization, revealing the localized pattern of glaucomatous RNFL thinning. All images are in the right-eye orientation.

Similar group difference visualization for choroid is shown in Figure 7. Compared to RNFL, choroidal thickness is more different by age (top row, Group B − Group A) than by glaucoma (Group C − Group B). Choroidal thickness of Group B is consistently lower than that of Group A across all regions, but in particular in the nasal and superior regions. The locations of statistical significance, as shown in the t-test map, of the age-related group difference is found particularly in the nasal region. The choroidal thickness difference by glaucoma was not as strong as that due to age and although the choroidal thickness in Group B was overall larger than that of group C, the point- and sector-wise t-test show limited statistical significance.

FIGURE 7

Figure 7. Effect of age (top row) and glaucoma (bottom row) on choroidal thickness. The point- and sector-wise subtraction of mean choroidal thickness across older normal (Group B) and young normal (Group A) is shown in the top left panel, and the p-values from point- and sector-wise t-test is shown in the top right panel. The point-and sector-wise subtraction of mean choroidal thickness between older glaucoma (Group C) and older normal (Group B) is shown in the bottom left panel, and the p-values from point- and sector-wise t-test is shown in the bottom right panel. The choroid is found to thin significantly with age, whereas changes are relatively less with glaucoma. The fshapes point-wise comparison shows the pattern of change in greater detail, revealing the localized pattern of age-related choroidal thinning. All images are in the right-eye orientation.

The relationship of RNFL thickness to aging and severity of glaucoma was examined by point- and sector-wise linear regression in Figure 8. The estimated slope or the rate of change associated with aging (mm/year) in the cross-sectional healthy subjects was plotted in the top row along with the goodness of fit by r². The estimated slope or the rate of change associated with visual field loss (mm/dB) in the cross-sectional age-matched subjects was plotted in the bottom row along with the goodness of fit by r². Aging did not show consistent, significant trend with RNFL thickness, whereas increasing glaucoma severity (more negative VFMD) was correlated to decreasing RNFL thickness. Similar spatial patterns are observed as seen in the group difference maps shown in Figure 6. In the most severely affected regions of superior-temporal and inferior-temporal regions, RNFL thickness change per MD unit decrease exceeded 5 μm.

FIGURE 8

Figure 8. Linear regression plots of RNFL thickness v. age (top row) and v. visual field mean deviation (VFMD). The slope in RNFL thickness = a*Age + b is plotted point- and sector-wise in the top left panel, along with the goodness of fit (r²) in the top right panel. The slope of RNFL thickness = a*VFMD + b is plotted point- and sector-wise in the bottom left panel, along with the goodness of fit (r²) on the bottom right panel. As with the group difference, RNFL thickness is negatively correlated with VFMD. The point-wise maps reveal that the rate of change has distinct spatial pattern with greater thinning along the regions where RNFL is generally thicker in healthy subjects. All images are in the right-eye orientation.

The relationship of choroidal thickness to aging and severity of glaucoma was also examined by point- and sector-wise linear regression in Figure 9. The estimated slope or the rate of change associated with aging (mm/year) in the cross-sectional healthy subjects was plotted in the top row along with the goodness of fit by r². The estimated slope or the rate of change associated with visual field loss (mm/dB) in the cross-sectional age-matched subjects was plotted in the bottom row along with the goodness of fit by r². Again, compared to the RNFL thickness, change in the choroidal thickness was associated more with aging than glaucoma. Aging was correlated to globally decreasing choroidal thickness, most significantly in the nasal and inferior regions. In the most severely affected regions, the average choroidal thickness change per year was ~3–4 μm.

FIGURE 9

Figure 9. Linear regression plots of choroidal thickness v. age (top row) and v. visual field mean deviation (VMFD). The slope in Choroidal thickness = a*Age + b is plotted point- and sector-wise in the top left panel, along with the goodness of fit (r²) in the top right panel. The slope of Choroidal thickness = a*VFMD + b is plotted point- and sector-wise in the bottom left panel, along with the goodness of fit (r²) on the bottom right panel. As with the group difference, choroidal thickness is negatively correlated with age. The point-wise maps show the steepest change in the nasal and inferior regions. All images are in the right-eye orientation.

The point-wise nature of fshape metrics can be utilized by simultaneous visualization of all data points from multiple subjects. Figure 9 displays the RNFL fshape thickness in the order of VFMD, which measures the glaucomatous functional loss, for the age-matched normal and glaucomatous eyes of Group B and C. In the top panel, each column represents an eye's point-wise RNFL fshape thickness, which approximates the true RNFL thickness as the sum of the RNFL fshape mean template thickness and the residual (t_i ≈ f_∗ + ζ_i) at each point. Horizontally, the eyes are ordered by VFMD, plotted below the thickness plot in grayscale. Vertically, the points are ordered by sectors from Nasal (N) and counter-clockwise to Inferior Nasal (IN). Within each sector, the points are ordered by the distance from BMO center, from the closest to the farthest. This visualization allows one to compare all eyes at each spatial location. There is observed an overall group-wise difference between Healthy older subjects, Early Glaucoma, and Moderate to Severe Glaucoma subjects. Comparing vertically from the top to the bottom, the healthy eyes show the thickness pattern that follows the mean template for Group B in Figure 4 where the superior and inferior regions are the thickest, and within each region, RNFL is thicker toward BMO and thins farther from BMO.

The bottom panel of Figure 10 visualizes the f-shape thickness in the top panel using z-scores normalized by the group mean and standard deviation of the RNFL thickness of the young healthy group (Group A) at each point, with increasing magnitude indicating greater deviation from the reference group. The differences between the groups are more apparent in z-score plot, with clear regional characteristics. Regionally, inferior RNFL, and to a lesser degree, inferior-temporal RNFL, are consistently thinner in all glaucoma eyes. In other sectors, the magnitude of z-score is greater for the moderate to severe glaucoma group than early glaucoma. The plot also shows glaucomatous RNFL thickness change is greater nearer BMO by the vertical gradations within individual regions in the glaucomatous eyes.

FIGURE 10

Figure 10. Multi-subjects RNFL thickness by visual field mean deviation (VFMD). By mapping to the common template space, all subjects have estimated thickness values (t_i ≈ f_∗ + ζ_i) at the same corresponding spatial points, visualized in the top panel. Vertically, the points are ordered by regions, from Nasal (N) to Inferior Nasal (IN). Within each region, the points are ordered by their distance from BMO center, from the closest to the farthest. Horizontally, the eyes are ordered by VFMD, plotted in grayscale below the thickness plot. In the bottom panel, the same data is visualized z-scores, which normalizes the values by the mean and standard deviation of the reference group (Group A, young healthy) and highlights the trend across regions and increasing VFMD magnitude.

Figure 11 displays the choroidal fshape thickness in the order of age for the normal eyes to show the thinning of the choroid observed with age. In the top panel, each column represents an eye's point-wise choroidal fshape thickness, which approximates the true choroidal thickness at each point. Horizontally, the eyes are ordered by age, plotted below the thickness plot in grayscale. As in Figure 10, the points are vertically ordered by sectors from Nasal (N) and counter-clockwise to Inferior Nasal (IN), and within each sector, the points are ordered by the distance from BMO center, from the closest to the farthest. In the bottom panel, the same data is visualized in z-scores normalized by the group mean and standard deviation of the choroidal thickness of the young healthy group. As shown in Figures 5, 7, there is a marked difference between the young healthy and older healthy eyes, and the choroid appears generally thicker in the superior half than the inferior half. Within the young eyes, the z-score is generally lower for the older eyes. As seen in Figure 7, the nasal region of the older eyes shows the highest magnitudes of z-score, suggesting the age-related choroidal thinning may be the most significant in the region.

FIGURE 11

Figure 11. Multi-subjects choroidal thickness by age. By mapping to the common template space, all subjects have estimated thickness values (t_i ≈ f_∗ + ζ_i) at the same corresponding spatial points, visualized in the top panel. Vertically, the points are ordered by regions, from Nasal (N) to Inferior Nasal (IN). Within each region, the points are ordered by their distance from BMO center, from the closest to the farthest. Horizontally, the eyes are ordered by age, plotted in grayscale below the thickness plot. In the bottom panel, the same data is visualized in z-scores, which normalizes the values by the mean and standard deviation of the reference group (Group A, young healthy), and highlights the trend across regions and increasing age.

Discussion

The 3D OCT images reveal the structure of the ocular posterior segment in great detail so as to enable visualization and quantification of the retinal layer morphometry. Individual measurements of layer thicknesses can be pooled into population-wide assessments of normative layer thicknesses and any changes that may occur as a function of age and disease. These allow insights into whether age and disease have a stereotypical pattern of influence in the retina, with common shape features and localizations, as well as variability across individuals and deviation of a particular subject from a reference population. In this paper, we presented the effect of age and glaucoma on retina nerve fiber layer (RNFL) and choroid using our novel f-shapes approach, which enables a point-wise assessment of retinal morphometrics across individuals via a registration approach. The fshape registration estimates a residual function that is added to the template thickness such that the template geometry after transformation matches the subject geometry, and the template thickness plus the residual function after transformation matches the subject thickness. This maps an individual's layer thicknesses onto the common coordinate system of the template geometry via the specific residual function estimated for that individual. At each location on the template surface, subsequent statistical analysis across the database can reveal trends and features that are common across individuals as a function of age and disease. A more conventional approach utilizes sectorization of the layer surface by calculating the average of all thickness measurements within each sector for each individual eye. Assuming that the sectors are defined with consistent anatomical and spatial correspondence across individuals, within-sector average provides a single scalar summary measure for a given sector that can be statistically analyzed for cross-sectional data taken from the same sector across the individual eyes. However, such approach is limited in spatial sensitivity due to averaging of values in a region. We examined the effect of age and glaucoma in RNFL and choroidal thickness in both of these approaches with four quantitative visualizations: (i) group rages (Figures 4, 5), (ii) group-wise difference and t-test maps (Figures 6, 7), (iii) linear regression with age and visual field mean deviation (VFMD) as predictors (Figures 8, 9), and (iv) multi-subject fshape thickness and z-scores plots (Figures 10, 11). The computation time including automated segmentation and mean template generation was ~40 min on a high-performance GPU cluster.

With age, RNFL showed relatively little difference between the young and older healthy subjects, with regression estimating no strong relationship between age and RNFL thickness. In previous studies using OCT measurements, RNFL thickness has been negatively associated with age (Budenz et al., 2007; Parikh et al., 2007; Bendschneider et al., 2010; Sarunic et al., 2010; More et al., 2011). In this study, the mean age of the young healthy subjects was thirty, and that of the older healthy subjects was fifty-seven. The age difference between the two groups may be too small for any marked difference, especially with the small sample size in the study. Older age was, however, associated with markedly thinner choroid, and the point-wise registration showed the difference was more significant in the nasal and inferior regions. With the recent advancement in OCT devices enabling the posterior boundary of the choroid to be imaged, age-related choroid thinning has been reported by multiple groups (Manjunath et al., 2010; Maul et al., 2011; Barteselli et al., 2012). Our result suggests overall thinning of peripapillary choroid with age, but with regional differences. That the older healthy subjects had comparable RNFL but thinner choroid compared to the young healthy subjects may indicate the age-related choroidal thinning does not directly and concurrently impact RNFL thickness. Recent work using speckle-variance OCT (SV-OCT) (Mammo et al., 2016) has shown degradation of RNFL microcapillaries in glaucoma, and it has been suggested the glaucomatous tissue loss may be driven by changes to the microvasculature. Although the choroid is a vascular layer, it mainly supplies the outer layers of the retina that are unaffected in glaucoma, and may therefore be separate from the factors that drive the RNFL tissue and capillary loss in glaucoma.

Glaucoma, as expected, was observed to be associated with significant thinning of RNFL. The pattern of loss visualized in the point-wise maps forms the hourglass crescent-shaped ridge, particularly in the inferior arm, highest toward the middle of the ridge, decreasing further away from the ridge center. The same hourglass-like pattern is observed in young normal subjects. The results of our study are consistent with known patterns of RNFL loss in glaucoma, but more importantly, it shows that the glaucomatous RNFL loss occur in a specific, uneven pattern that follows the initial RNFL thickness, suggesting that the time of onset and significance of the RNFL loss may be proportional to the initial RNFL thickness in the region. The loss significance was also higher in the region closer to BMO. The fshape RNFL loss map elucidates the results of previous studies that found the highest diagnostic ability of the RNFL loss in the inferior and temporal-inferior sectors (Sehi et al., 2009; Mwanza et al., 2011). Although the choroidal thickness in the older glaucoma subjects seemed to be somewhat less than in the age-matched healthy subjects, it was still statistically comparable, indicating that glaucoma pathology may not have significant, direct impact on the choroid as it does on the RNFL. The role of choroid in glaucoma has been debated, and it is likely complex and multifaceted. Disturbed autoregulation of the choroid has been suggested as part of the disease pathology (Hayreh, 1969; Ulrich et al., 1996). Multiple studies using OCT images (Ehrlich et al., 2011; Maul et al., 2011; Li et al., 2013) consistently reported no changes in the peripapillary choroidal thickness in primary open angle glaucoma (POAG); however, Li et al. (2017) reported thicker temporal peripapillary choroidal area in POAG patients using enhanced depth imaging OCT, and Song et al. (2016) reported global and all 12 clock-hour peripapillary choroidal thickness thinner in OAG patients using swept-source OCT.

Figures 10, 11 presented a large-data visualization with each subject's point-wise thickness values color mapped in a single column, and multiple subjects' data columns displayed concurrently, arranged in the order of visual field mean deviation (VFMD, Figure 10) and age (Figure 11). This method allows for presentation of all data points from multiple subjects in a way that highlights the trend and discrepancies in the data. Normalizing the data by the mean and standard deviation of the young, healthy group as the reference removes the baseline in the data and further brings out the differences with respect to controls.

The patterns of change shown in the point-wise and sector-wise approaches were overall consistent. The point-wise registration was able to show localized features in higher resolution compared to the sectorization, revealing detailed regional patterns and potentially furthering our understanding of the disease mechanism. This approach may be useful for characterizing the focal localized patterns that are often seen in glaucoma, both in RNFL and in other subsurface structures such as lamina cribrosa and peripapillary tissues. In this work, we also examined multiple factors (age, glaucoma) in multiple layers (RNFL, choroid) concurrently for a more complete picture in understanding our data. The results showed glaucomatous thinning of RNFL, and age-related thinning of choroid and how the spatial patterns of the tissue loss in the two layers were localized and distinct. Limitations of this work include relatively small sample size, inclusion of fellow eyes, and light beam angle-related uncertainty in retinal layer thickness measurement, although given the relatively small size of the field of view in the images in this study, this effect is likely limited. Our future work will include expanding the analyses presented here to other retinal layers and to the macular region to build a comprehensive presentation of the retinal morphometrics, and incorporating retinal vessels and capillary density measures from SV-OCT as metrics along with layer thickness.

Ethics Statement

The study followed the tenets of the Declaration of Helsinki, and informed and written consents were obtained from the participants. Ethics review for the study was approved by Simon Fraser University (SFU) and University of British Columbia (UBC).

Author Contributions

Designed research: SL, KP, NC, BC, AT, PJM, MVS, and MFB. Analyzed and interpreted data: SL and MFB. Wrote manuscript: SL, MLH, MVS, and MFB. Reviewed and approved manuscript: KP, NC, BC, AT, PJM, MVS, and MFB.

Funding

This work was supported by Canadian Institutes of Health Research, Natural Sciences and Engineering Research Council of Canada, Michael Smith Foundation for Health Research, Pacific Alzheimer Research Foundation, and Brain Canada.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work is based on our conference presentation (Lee S. et al., 2016). It contains extended experimental study with comparison with sectorization and regression analyses.

References

Anthony, B. J., Chen, M., Carass, A., Jedynak, B. M., Al-Louzi, O., Solomon, S., et al. (2016). “Voxel based morphometry in optical coherence tomography: validation & core findings,” in Proceedings of SPIE 9788. San Diego, CA.

Barteselli, G., Chhablani, J., El-Emam, S., Wang, H., Chuang, J., Kozak, I., et al. (2012). Choroidal volume variations with age, axial length, and sex in healthy subjects: a three-dimensional analysis. Ophthalmology 119, 2572–2578. doi: 10.1016/j.ophtha.2012.06.065

PubMed Abstract | CrossRef Full Text | Google Scholar

Bendschneider, D., Tornow, R. P., Horn, F. K., Laemmer, R., Roessler, C. W., Juenemann, A. G., et al. (2010). Retinal nerve fiber layer thickness in normals measured by spectral domain, O. C. T. J. Glaucoma 19, 475–482. doi: 10.1097/IJG.0b013e3181c4b0c7

CrossRef Full Text | Google Scholar

Budenz, D. L., Anders, D. R., Varma, R., Schuman, J., Cantor, L., Savell, J., et al. (2007). Determinants of normal retinal nerve fiber layer thickness measured by stratus, O. C. T. Ophthalmology 114, 1046–1052. doi: 10.1016/j.ophtha.2006.08.046

CrossRef Full Text | Google Scholar

Charlier, B., Charon, N., and Trouvé, A. (2015). The fshape framework for the variability analysis of functional shapes. Found. Comput. Math. 17:287. doi: 10.1007/s10208-015-9288-2

CrossRef Full Text | Google Scholar

Chauhan, B. C., and Burgoyne, C. F. (2013). From clinical examination of the optic disc to clinical assessment of the optic nerve head: a paradigm change. Am. J. Ophthalmol. 156, 218–227. doi: 10.1016/j.ajo.2013.04.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, M., Lang, A., Ying, H. S., Calabresi, P. A., Prince, J. L., and Carass, A. (2014). Analysis of macular OCT images using deformable registration. Biomed. Opt. Express 5, 2196–2214. doi: 10.1364/BOE.5.002196

PubMed Abstract | CrossRef Full Text | Google Scholar

Ehrlich, J. R., Peterson, J., Parlitsis, G., Kay, K. Y., Kiss, S., and Radcliffe, N. M. (2011). Peripapillary choroidal thickness in glaucoma measured with optical coherence tomography. Exp. Eye Res. 92, 189–194. doi: 10.1016/j.exer.2011.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Gardiner, S. K., Ren, R., Yang, H., Fortune, H., Burgoyne, C. F., and Demirel, S. (2014). A method to estimate the amount of neuroretinal rim tissue in glaucoma: comparison with current methods for measuring rim area. Am. J. Ophthalmol. 157, 540–549. doi: 10.1016/j.ajo.2013.11.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Gibson, E., Young, M., Sarunic, M. V., and Beg, M. F. (2010). Optic nerve head registration via hemispherical surface and volume registration. IEEE Trans. Biomed. Eng. 57, 2592–2595. doi: 10.1109/TBME.2010.2060337

PubMed Abstract | CrossRef Full Text | Google Scholar

Hayreh, S. S. (1969). Blood supply of the optic nerve head and its role in optic atrophy, glaucoma, and oedema of the optic disc. Br. J. Ophthalmol. 53, 721–748. doi: 10.1136/bjo.53.11.721

PubMed Abstract | CrossRef Full Text | Google Scholar

Jørgensen, T. M., Thomadsen, J., Christensen, U., Soliman, W., and Sander, B. (2007). Enhancing the signal-to-noise ratio in ophthalmic optical coherence tomography by image registration–method and clinical examples. J. Biomed. Opt. 12:041208. doi: 10.1117/1.2772879

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, E. J., Lee, K. M., Lee, S. H., and Kim, T. W. (2016). OCT angiography of the peripapillary retina in primary open-angle glaucoma. Invest. Ophthalmol. Vis. Sci. 57, 6265–6270. doi: 10.1167/iovs.16-20287

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S., Beg, M. F., and Sarunic, M. V. (2013a). “Segmentation of the macular choroid in OCT images acquired at 830nm and 1060nm,” in Proceedings of SPIE 8802. Munich. doi: 10.1117/12.2032603

CrossRef Full Text | Google Scholar

Lee, S., Fallah, N., Forooghian, F., Ko, A., Pakzad-Vaezi, K., Merkur, A. B., et al. (2013b). Comparative analysis of repeatability of manual and automated choroidal thickness measurements in nonneovascular age-related macular degeneration. Invest. Ophthalmol. Vis. Sci. 54, 2864–2871. doi: 10.1167/iovs.12-11521

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S., Charon, N., Charlier, B., Popuri, K., Lebed, E., Sarunic, M. V., et al. (2017). Atlas-based shape analysis and classification of retinal optical coherence tomography images using the functional shape (fshape) framework. Med. Image Anal. 35, 570–581. doi: 10.1016/j.media.2016.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S., Han, S. X., Young, M., Beg, M. F., Sarunic, M. V., and Mackenzie, P. J. (2014). Optic nerve head and peripapillary morphometrics in myopic glaucoma. Invest Ophthalmol. Vis. Sci. 55, 4378–4393. doi: 10.1167/iovs.14-14227

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S., Lebed, E., Sarunic, M. V., and Beg, M. F. (2015). Exact surface reistration of retinal surfaces from 3D optical coherence tomography images. IEEE Trans. Biomed. Eng. 62, 609–617. doi: 10.1109/TBME.2014.2361778

CrossRef Full Text | Google Scholar

Lee, S., Popuri, K., Matsubara, J. A., Eadie, B., Merkur, A., Mackenzie, P., et al. (2016). Novel fshape mapping of peripapillary RNFL and choroidal thickness for group-wise comparison of glaucoma and healthy aging. Invest. Ophthalmol. Vis. Sci. 57, 3389.

Leung, C. K. (2014). Diagnosing glaucoma progression with optical coherence tomography. Curr. Opin. Ophthalmol. 25, 104–111. doi: 10.1097/ICU.0000000000000024

PubMed Abstract | CrossRef Full Text | Google Scholar

Leung, C. K. S., Choi, N., Weinreb, R. N., Liu, S., Ye, C., Liu, L., et al. (2010). Retinal nerve fiber layer imaging with spectral-domain optical coherence tomography: pattern of RNFL defects in glaucoma. Ophthalmology 117, 2337–2344. doi: 10.1016/j.ophtha.2010.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, K., Wu, X., Chen, D. Z., and Sonka, M. (2006). Optimal surface segmentation in volumetric images–a graph-theoretic approach. IEEE Trans. Pattern Anal. Mach. Intell. 28, 119–134. doi: 10.1109/TPAMI.2006.19

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L., Bian, A., Zhou, Q., and Mao, J. (2013). Peripapillary choroidal thickness in both eyes of glaucoma patients with unilateral visual field loss. Am. J. Ophthalmol. 156, 1277–1284. doi: 10.1016/j.ajo.2013.07.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, M., Yan, X., Song, Y., Guo, J., and Zhang, H. (2017). Choroidal area assessment in various fundus sectors of patients at different stages of primary open-angle glaucoma by using enhanced depth imaging optical coherence tomography. Medicine 96:e6293. doi: 10.1097/MD.0000000000006293

PubMed Abstract | CrossRef Full Text | Google Scholar

Mammo, Z., Heisler, M., Balaratnasingam, C., Lee, S., Yu, D.-Y., Mackenzie, P., et al. (2016). Quantitative optical coherence tomography angiography of radial peripapillary capillaries in glaucoma, glaucoma suspect, and normal eyes. Am. J. Ophthalmol. 170, 41–49. doi: 10.1016/j.ajo.2016.07.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Manjunath, V., Taha, M., Fujimoto, J. G., and Duker, J. S. (2010). Choroidal thickness in normal eyes measured using cirrus HD optical coherence tomography. Am. J. Ophthalmol. 150, 325–329. doi: 10.1016/j.ajo.2010.04.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Maul, E. A., Friedman, D. S., Chang, D. S., Boland, M. V., Ramulu, P. Y., Jampel, H. D., et al. (2011). Choroidal thickness measured by spectral domain optical coherence tomography: factors affecting thickness in glaucoma patients. Ophthalmology 118, 1571–1579. doi: 10.1016/j.ophtha.2011.01.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Medeiros, F. A., Zangwill, L. M., Bowd, C., Vessani, R. M., Susanna, R., and Weinreb, R. N. (2005). Evaluation of retinal nerve fiber layer, optic nerve head, and macular thickness measurements for glaucoma detection using optical coherence tomography. Am. J. Ophthalmol. 139, 44–55. doi: 10.1016/j.ajo.2004.08.069

PubMed Abstract | CrossRef Full Text | Google Scholar

More, H. L., Chen, J., Gibson, E., Donelan, J. M., and Beg, M. F. (2011). A semi-automated method for identifying and measuring myelinated nerve fibers in scanning electron microscope images. J. Neurosci. Methods 201, 149–158. doi: 10.1016/j.jneumeth.2011.07.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Mwanza, J.-C., Oakley, J. D., Budenz, D. L., and Anderson, D. R. (2011). Ability of cirrus HD-OCT optic nerve head parameters to discriminate normal from glaucomatous eyes. Ophthalmology 118, 241–248. doi: 10.1016/j.ophtha.2010.06.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Niemeijer, M., Garvin, M. K., Lee, K., van Ginneken, B., Abramoff, M. D., and Sonka, M. (2009). “Registration of 3D spectral OCT volumes using 3D SIFT feature point matching,” in Proceedings of SPIE, vol. 7259. Orlando, FL. doi: 10.1117/12.811906

CrossRef Full Text | Google Scholar

Parikh, R. S., Parikh, S. R., Sekhar, G. C., Prabakaran, S., Babu, J. G., and Thomas, R. (2007). Normal age-related decay of retinal nerve fiber layer thickness. Ophthalmology 114, 921–926. doi: 10.1016/j.ophtha.2007.01.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarunic, M. V., Yazdanpanah, A., Gibson, E., Xu, J., Bai, Y., Lee, S., et al. (2010). Longitudinal study of retinal degeneration in a rat using spectral domain optical coherence tomography. Opt. Exp. 18, 23435–23441. doi: 10.1364/OE.18.023435

PubMed Abstract | CrossRef Full Text | Google Scholar

Sehi, M., Grewal, D. S., Sheets, C. W., and Greenfield, D. S. (2009). Diagnostic ability of fourier-domain vs time-domain optical coherence tomography for glaucoma detection. Am. J. Ophthalmol. 148, 597–605. doi: 10.1016/j.ajo.2009.05.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, Y. J., Kim, Y. K., Jeong, J. W., and Park, K. H. (2016). Assessment of open-angle glaucoma peripapillary and macular choroidal thickness using swept-source optical coherence tomography (SS-OCT). PLoS ONE 11:e0157333. doi: 10.1371/journal.pone.0157333

PubMed Abstract | CrossRef Full Text | Google Scholar

Ulrich, A., Ulrich, C., Barth, T., and Ulrich, W. (1996). Detection of disturbed autoregulation of the peripapillary choroid in primary open angle glaucoma. Ophthalmic Surg. Lasers 27, 746–757.

PubMed Abstract | Google Scholar

Young, M., Lebed, E., Jian, Y., Mackenzie, P. J., Beg, M. F., and Sarunic, M. V. (2011). Real-time high-speed volumetric imaging using compressive sampling optical coherence tomography. Biomed. Opt. Express 2, 2690–2697. doi: 10.1364/BOE.2.002690

PubMed Abstract | CrossRef Full Text | Google Scholar