Multivariate Model of the Intersubject Variability of the Retinal Nerve Fiber Layer Thickness in Healthy Subjects

Abstract

Purpose: We present and validate a multivariate model that partially compensates for retinal nerve fiber layer (RNFL) intersubject variability.

Methods: A total of 202 healthy volunteers randomly attributed to a training (TS) and a validation (VS) sample underwent complete ophthalmic examination, including Fourier-domain optical coherence tomography (FD-OCT). We acquired FD-OCT data centered at the optic disc (OD) and the macula. Two-dimensional (2D) projection images were computed and registered, to determine the distance between fovea and OD centers (FD) and their respective angle (FA). Retinal vessels were automatically segmented in the projection images and used to calculate the circumpapillary retinal vessel density (RVD) profile. Using the TS, a multivariate model was calculated for each of 256 sectors of the RNFL, including OD ratio, orientation and area, RVD, FD, FA, age, and refractive error. Model selection was based on Akaike Information Criteria. The compensation effect was determined for 12 clock hour sectors, comparing the coefficients of variation (CoV) of measured and model-compensated RNFL thicknesses. The model then was applied to the VS, and CoV was calculated.

Results: The R value for the multivariate model was, on average 0.57 (max = 0.68). Compensation reduced the CoV on average by 18%, both for the TS and VS (up to 23% and 29%), respectively.

Conclusions: We have developed and validated a comprehensive multivariate model that may be used to create a narrower range of normative RNFL data, which could improve diagnostic separation between early glaucoma and healthy subjects. This, however, remains to be demonstrated in future studies.

Glaucoma is an optic neuropathy characterized by the progressive and irreversible loss of retinal ganglion cells and their axons, leading to visual field loss,¹ with no definitive cure established. Despite the efforts to define a set of parameters that can better discriminate between healthy and glaucomatous subjects,² an accurate and early diagnosis still is under development. 

Currently, early diagnosis of glaucoma is partially performed based on imaging techniques, such as optical coherence tomography (OCT).³ This is done mostly by determining the circumpapillary retinal nerve fiber layer (RNFL) thickness and comparing the measurements with a normative database. Despite a relatively high reproducibility,⁴ RNFL thickness measurements present a high intersubject variance,⁵ due to the influence of numerous individual factors. Some of these factors, such as sex, age, and ethnicity,^6–9 are taken into account by the normative databases. However, others are not, including optic disc (OD) and fovea parameters, retinal vessel position, axial length, and refractive error, despite the fact that they also may present an association with circumpapillary RNFL thickness and distribution.^10–19 

Concerning the impact of the retinal vasculature, previous publications focused on the correlation of specific retinal vessels with RNFL thickness.^20–22 To enable a more comprehensive analysis of this factor, we recently developed the parameter of circumpapillary retinal vessel density (RVD)—a function dependent on thickness and distribution of all circumpapillary retinal blood vessels.^23,24 This parameter proved to have a close relationship with RNFL pattern, especially when measured at the RNFL measurement circle of 3.4-mm diameter around the OD. We have shown that interindividual variance may be reduced, on average, by 11% and up to 20% when taking RVD into account.²⁴ 

It is our aim to develop a comprehensive multivariate model to explain a larger part of the interindividual RNFL variability, by including factors that are related to the thickness and/or distribution of circumpapillary RNFL, such as RVD, age, refractive error, OD size and shape, as well as fovea angle (FA) and distance (FD) relative to the OD center. Furthermore, our hypothesis that interindividual variability is reduced by application of our model is validated in an independent sample of healthy volunteers. 

Methods

Subjects

This study was performed in collaboration between the Department of Ophthalmology and Optometry and the Section for Medical Information Management and Imaging of the Center for Medical Statistics Informatics and Intelligent Systems of the Medical University of Vienna. The protocol was approved by the Ethics Committee of the Medical University of Vienna and written informed consent was obtained from all volunteers included. The Guidelines of Good Clinical practice and the Declaration of Helsinki were followed. Both samples (training [TS] and validation [VS]) followed the same criteria. 

Inclusion and Exclusion Criteria

Inclusion and exclusion criteria, as well as description of the examination protocols, have been described previously.²³ In addition to these criteria, only OD examinations acquired with Fourier-domain OCT (FD-OCT, Cirrus; Carl Zeiss Meditec, Inc., Dublin, CA, USA) with quality index higher than 5 were included (ranging from 0–10), due to reports of image quality influence on RNFL thickness measurements.²⁵ Moreover, scans with movement artifacts, in any direction and within the measurement circle, were excluded. 

OCT Acquisition

For each volunteer, the study eye was subjected to pupil dilation using tropicamide drops (Mydriaticum; Agepha Pharma s.r.o., Senec, Slovakia). Two different protocols were acquired using FD-OCT, one centered in the macula (macular cube, 200 × 200) and the other centered in the OD (OD cube, 200 × 200). After each examination, there was a control of the image quality. When the results were unsatisfactory (scan quality lower than 6 and/or movement artifacts within the RNFL measurement circle), the examination was repeated, until the inclusion criteria were met. If a satisfying quality could not be obtained after repeated measurements, the subject was excluded from the analysis. 

Automated Parameter Extraction

The complete set of data was randomly divided into two samples: the TS and the VS. Both samples were subjected to the same methodology for automated parameters extraction. For each subject, 2-dimensional (2D) projection images of both protocols were computed based on the volumetric OCT data. These images were computed by the OPTIMA group, from the Medical University of Vienna, and are based on automatic OCT layer segmentation, using the algorithm developed at the University of Iowa and explained previously.^26–28 Succinctly, based on the automatic layer segmentation of the 3D volumetric data acquired with FD-OCT, each A-scan between the boundary of the myloid and ellipsoid of inner segments (BMEIS) and the outer boundary of the RPE was averaged. This resulted in a 2D representation of the 3D scan interval between the aforementioned retinal structures, with accentuated retinal vasculature contrast. Figure 1 presents an example of one B-scan evidencing the boundaries used to calculate the projection image. 

Figure 1

 

Two-dimensional projection images calculation. The planar representation of the 3D volumetric data is determined by averaging all A-scans within the boundary of the myloid and ellipsoid of inner segments, and the outer boundary of the RPE (solid white lines in the B-scan, left). Black dashed line juxtaposed in the projection image (left) illustrates the scan location.

The parameters referring to OD are: area, orientation (angle between the horizontal axis and the major axis of the OD), and ratio (quotient between major and minor axis). These parameters are extracted using the OD contour information obtained directly from FD-OCT. Optic disc center was determined as the center point of the OD contour. 

An automatic vessel segmentation algorithm, based on morphological operators and on the Frangi method,^29,30 was developed in MATLAB (Version R2012b; Mathworks, Inc., Natick, MA, USA), to obtain the vessel tree from OD and macula centered images. 

Based on the results of the automatic segmentation of the retinal vessel tree, determined in the OD centered image, the thicknesses of all measurable vessels contained within a measurement band extending from 3.28 to 3.64 mm, centered in the OD, were used to determine the RVD. Each vessel thickness was attributed to a sector (in a total of 256), depending on the angle of a line connecting the vessel center point to the OD center and a horizontal line. The thickness values of all vessels within the same sector were summed up. Therefore, for each subject, we obtain a discrete profile of vessel distribution. To generate a semicontinuous profile, similar to RNFL distribution, that expresses vessel density in the vicinity of the RNFL measurement area, each individual discrete profile of vessel distribution was convoluted with a Gaussian window. This operation combines both functions—the discrete profile of retinal vessel distribution and the Gaussian window—by multiplying them and, therefore, expressing the overlap of both functions along the sectors.³¹ The size of the Gaussian window used for this operation was optimized for the maximization of the median of the intersubject correlation of RNFL thickness and RVD. More details on RVD and Gaussian convolution and optimization have been described previously.²³ 

Using both vessel trees of the macula and the OD centered projection images, we automatically registered the two images corresponding to each subject with a rigid transformation, based on a cross-correlation algorithm. From the registered image and considering fovea center as automatically determined in FD-OCT, we obtained the fovea parameters: first, the FD, which corresponds to the distance between OD and fovea centers; second, the FA, which corresponds to the angle between a line connecting fovea and OD centers and a horizontal line passing through the OD center. 

Statistical Analysis

Statistical analysis was performed using SPSS software package (Version 21; SPSS, Inc., Chicago, IL, USA). We developed a multivariate linear model for each sector of the RNFL, which should partially explain the interindividual variance of the RNFL in healthy subjects. The parameters used to develop the multivariate model were: age, spherical refractive error (RefError), RVD, OD area (ODA), orientation (ODO), and ratio (ODR), as well as FA and FD. For each sector of the RNFL we used the RVD value of the same sector. This analysis was performed across all subjects of the TS. 

For each of the 256 sectors of the RNFL, we performed all subset selection over 256 candidate models, which result from considering our 8 putative parameters. Model selection was performed by minimizing the Akaike Information Criteria (AIC), which gives a tradeoff between data fitting and model complexity, by estimating the expected loss in information in choosing a model.³² This process was repeated for all 256 sectors of the RNFL. After obtaining the best model in each sector, we retrieved the regression coefficients of the multiple linear regression of that model for further analyses. 

To evaluate the effect of the final multivariate model on interindividual variation of RNFL thickness distribution, we calculated, for each subject, a compensated RNFL thickness profile by estimating an individual compensation factor. 

For each of the 256 sectors we, thus, defined the compensated RNFL as defined in Equation 1. The term X_i refers to each individual parameter, the term X̄ _i refers to the average value of each parameter, β_i refers to coefficients of the regression model of the relevant sectors, and the subscript i defines each parameter (in a total of 8).  

 

 

The mean values of measured and compensated profiles will overlap, since the compensation normalizes the RNFL for each subject as if each parameter were at the average of that population. 

According to basic statistical considerations, it holds that the compensated RNFL can be defined by Equation 2:  

 

 

 

 

The residuals can be further determined by Equation 3:  

 

 

The sector-wise compensated RNFL thickness measurements were calculated across all subjects. Additionally, we determined the interindividual coefficient of variation (CoV) in 12 clock hour sectors for the measured RNFL and the compensated RNFL, to analyze the reduction of variation obtained by the compensated model. 

Multivariate Model Validation

To validate the multivariate model derived from the training sample, the compensated RNFL according to Equation 1 was applied to the independent validation sample. Variation of the compensated RNFL was compared with the variation when using the conventional approach without compensation. Furthermore, we calculated the CoV in 12 clock hour sectors for the measured and the compensated RNFL across the validation sample. 

Results

Description of Statistical Subsets

Both samples—TS and VS—included a total of 101 subjects each. The demographic data for both samples is presented in Table 1. There were no significant differences between the two samples for any of the variables. 

Table 1

 

Demographic Data of TS and VS

 

Figure 2 illustrates one example of the projection images used both for OD and macula protocols, and the respective result of automatic vessel segmentation and registration. Figure 3 presents a schematic representation of the parameters depending on fovea and OD used for the model development. 

Figure 2

 

Upper: Example of the projection images used for retinal vessel segmentation and registration, centered in macula (left) and OD (right). Lower: Example of the automatic segmentation and registration results.

Figure 3

 

Schematic representation of the parameters used for the multivariate model development. The two solid lines refer to the horizontal axis and the line connecting the OD and fovea centers (marked with X). Dashed lines refer to the major and minor axis of the OD. The angle α refers to the FA; the angle β refers to the ODO.

Moreover, the parameters used for the automated vessel segmentation by applying the Frangi method were s = 1, β = 0.5, and c = 0. 

Multivariate Model Reveals RVD as Having the Strongest Influence in RNFL Profile

The semipartial correlation coefficients of the best fitting multivariate model in the TS in each sector, as selected by AIC, are summarized in Table 2. 

Table 2

 

Training Sample: Statistical Outcome From the Multivariate Regression Analysis (Absolute Values Displayed)

 

Table 2 elucidates how much each parameter influences RNFL profile, by presenting the semipartial correlation coefficient calculated for each parameter, together with the R value for the complete multivariate model. The parameter with most impact on RNFL thickness distribution is RVD, with approximately 90% of the sectors presenting an association with RNFL and an average correlation coefficient of 0.41. 

On average, the R value of the multivariate model was 0.57 (range for the different sectors, 0.37–0.68), which translates into an average R² value of 0.32. This means that, in average, 32% of the intersubject variation might be explained by the reported combination of factors. This percentage might be up to 46% in some sectors (a maximum R value of 0.68 translates into an R² value of 0.46). 

A modified temporal-superior-nasal-inferior-temporal (TSNIT) graph was constructed for each variable based on its individual regression coefficient calculated for each sector (when included on the selected model of that sector). The graphs for each variable are displayed in Figures 4a through 4h. 

Figure 4

 

In the TS, modified TSNIT graphs plotting the regression coefficients (RC) from each of the selected variables of (a) circumpapillary RVD, (b) ODR, (c) ODO, (d) ODA, (e) FD, (f) FA, (g) age, and (h) refractive error.

Figure 4a illustrates the increase in RNFL thickness with higher RVD in a majority of the sectors (237 of 256 sectors), with nonsignificant sectors restricted to a part of the temporal quadrant. Increasing fovea distance (Fig. 4e), as opposed to this, exhibits a main influence in the temporal quadrant, where it increases the RNFL, while in some other parts of the TSNIT curve it decreases the RNFL. The impact of FA (Fig. 4f) extends over two-thirds of the circumference, with the exception of the nasal quadrant. A more positive (less negative) FA—meaning the fovea being located more superior—results in a decrease of RNFL in the superior half and an increase in the inferior half of the temporal quadrant. For the (nasal) superior and (nasal) inferior sectors, we found the exact opposite impact. Regarding the OD parameters, OD area (Fig. 4d) has a substantial positive impact in all quadrants except the temporal, while OD ratio (Fig. 4b) has a meaningful negative role mainly in the nasal and temporal sectors, indicating smaller RNFL thickness with more elongated OD. The impact of OD orientation (Fig. 4c) is limited to the inferior nasal quadrant. Increasing age (Fig. 4g) reduces RNFL mainly in superior and inferior nasal sectors, while refractive error (Fig. 4h) presents a negative slope in the temporal quadrant, combined with a positive slope in most of the remaining sectors. This means that myopic eyes tend to have a thicker RNFL in the temporal quadrant and a thinner RNFL elsewhere compared to hyperopic subjects. 

Compensated RNFL Model Reveals a Clinically Significant Reduction of RNFL Variation

To determine the reduction in variation of RNFL thickness measurements by our compensation model, we applied Equation 1 to all 256 sectors to the measured RNFL of the TS. We calculated the interindividual CoVs of the measured and compensated RNFL thickness values of 12 clock hour sectors, presented in the TSNIT order (Table 3). 

Table 3

 

Training Sample CoVs for Measured and Compensated RNFL Thickness Values

 

On average, the relative reduction in coefficient of variation was 18%. This relative reduction also is graphically visible when comparing the TSNIT profiles of measured and compensated RNFL thickness measurements. Figure 5 illustrates the profile of the TS together with the variation, displayed as the ±2 SD around the mean. The thinner lines refer to measured values of RNFL thickness and thicker lines refer to compensated RNFL thickness. The variation of compensated RNFL clearly is reduced in all regions of the TSNIT profile. 

Figure 5

 

In the TS, TSNIT profiles of the measured and compensated model of RNFL thickness, calculated for the TS according to Equation 1. Dashed gray lines refer to average RNFL measured by FD-OCT ± 2 SD; dashed black lines refer to average values of RNFL compensated ±2 SD; solid line refers to average RNFL measured by FD-OCT.

Successful Validation of the Multivariate Model in a Separate Sample

To validate the described multivariate model and, more specifically, the reduction of CoV by application of the compensation formula, we determined the same presented eight physiological parameters in the independent validation sample of 101 healthy volunteers. We then applied to this sample the same multivariate model and the compensation formula of Equation 1 as derived from the TS. The results obtained for the CoV calculated in 12 clock hour sectors for the measured and the compensated RNFL are presented in Table 4. 

Table 4

 

Validation Sample CoVs for Measured and Compensated RNFL Thickness Values of the VS

 

In the validation sample, the average (18%) and maximum (29%) relative reduction of CoV were similar to those previously described for the TS. The complete average TSNIT profile (±2 SD) for measured and compensated RNFL thickness values is shown in Figure 6. Again, the variation, displayed as the ±2 SD around the mean, clearly is reduced in all regions of the thickness profile of compensated RNFL. 

Figure 6

 

In the VS, TSNIT profiles of the measured and compensated model of RNFL thickness, calculated for the validation sample according to the regression coefficients in Equation 1. Dashed gray lines refer to average RNFL measured by FD-OCT ± 2 SD; dashed black lines refer to average values of RNFL compensated ±2 SD; solid line refers to average RNFL measured by FD-OCT.

Discussion

We present a new multivariate model that facilitates a clinically meaningful reduction of the intersubject variance of RNFL by compensating for variation of some parameters that correlate with RNFL thickness. This model includes physiological parameters automatically extracted from OCT fundus images, such as RVD, OD shape, fovea location, and additionally, age and refractive error. It already has been shown in previous works that circumpapillary RVD has a significant influence on RNFL thickness profile.^23,24 By developing a multivariate model including additional individual characteristics of the healthy population, we obtain a strong interindividual association between the multivariate model and RNFL thickness measurements (mean R = 0.57), as opposed to the more moderate effect we have reported when considering RVD alone (mean R = 0.42).²⁴ Although RVD remains the feature with by far strongest and most widespread impact throughout almost the entire RNFL thickness profile, we have now identified in our multivariate model additional influences of age, refractive error, OD shape, and fovea location. 

Moreover, we have presented a reduction of interindividual variance of RNFL measurements, when calculating a compensated RNFL profile. On average, this compensation reduced the intersubject variability by 18%. Furthermore, this model has been tested in an independent validation sample of 101 healthy subjects and the reduction in CoV (on average 18%) was similar in both samples, which, fundamentally establishes the proof of concept that the interindividual variance of RNFL may be reduced by applying our new method. 

We have discussed previously possible explanations for the correlation of the RVD and the RNFL.^23,24 Briefly, during embryological development, the retinal vessels grow following demands for oxygen and nutrition. In the developing peripapillary retina, both demands are likely largely driven by the developing RNFL.³³ Thus, the locations and thicknesses of the peripapillary retinal blood vessels (reflected by the RVD) should carry an imprint of the RNFL distribution at birth.^34,35 Since the retinal vessel pattern is largely preserved throughout lifetime, this information can be assessed later. 

Concerning refractive error, there is a positive association with the RNFL thickness in three-quarters of the graph, while in the temporal quadrant we have described a negative influence on the RNFL thickness. This means that in the temporal quadrant hyperopic subjects have a thinner RNFL compared to myopic subjects, while in most other locations it is the other way around. The relationship among retinal vessel location, refractive error, and RNFL has been investigated previously and the angle of the major temporal retinal vessel arcades tended to be smaller in myopia, which was associated with increased RNFL thickness in the temporal sector and more temporally shifted peak locations of RNFL.¹⁹ Moreover, the positive regression coefficients in the superior, nasal, and inferior regions found in our model might be partially due to magnification errors. On a hyperopic eye, the image is magnified, meaning that the true measurement circle is, in fact, smaller; thus, corresponding to a thicker RNFL.^11,36 

The impact that fovea distance has on the RNFL has been reported previously¹³ as a probability for having a temporal shift in RNFL profile in subjects where the distance between foveola and optic disc center was increased. This finding fits to our results, where on the temporal quadrant, a larger fovea distance corresponded to thicker RNFL. Regarding the FA, its impact is positive in the temporal inferior and nasal superior quadrants, and negative in the temporal superior and nasal inferior quadrants. This effect is expected, because a more inferiorly located fovea (more negative FA) suggests a rotation of the superior RNFL peak toward the temporal quadrant, and of the inferior RNFL peak away from the temporal quadrant. 

As for the ODA, it has been stated that larger ODs have a thicker RNFL.¹⁵ This effect is largely confirmed by our results. 

One of the limitations of this study concerns the age distribution of included subjects. Both populations consist mostly of young adults, which might explain why age does not show an association with RNFL throughout the majority of the TSNIT curve. This fact does not, however, contradict previous findings correlating RNFL thinning with age,^6,8,37 but only stipulates that this effect is not easily detected in a population of relatively young subjects. This might be a constraint when applying this model to subjects older than 60 years. 

Another limitation is the lack of information on axial length measurements. We did not correct our measurements on the OCT fundus images for magnification effects, due to the lack of axial length measurements. However, it has been described previously that this factor might be used to reduce the intersubject variability of RNFL,¹⁶ as well as it might influence the correlation found between some other parameters and RNFL thickness.^11,19 

Compared to our previous work, we now have developed an automated RVD assessment and a comprehensive multivariate model. To our present knowledge, the combined influence of all individual parameters we have examined has not yet been under careful analysis previously. This work confirms the impact that vessel distribution has on the circumpapillary RNFL profile. Despite the inclusion of seven more parameters in the new model, the association between RVD and RNFL²⁴ was reproduced to a very large extent. 

We have validated our model in an independent VS. In this VS, we have demonstrated that application of the compensation formulas results in a clinically meaningful reduction of interindividual variance of RNFL. This indicates that we are describing really existing physiologic correlations together with a method to exploit these correlations. Thus, we have shown and validated a method that may be used to create normative RNFL data with a narrower range of normal values, which could improve diagnostic separation between early glaucoma and healthy subjects. This, however, remains to be shown in future studies. 

Acknowledgments

Supported by Grant Wiener Wissenschafts-, Forschungs- and Technologiefonds, WWTF – LS11-046. 

Disclosure: I. Pereira, None; H. Resch, None; F. Schwarzhans, None; J. Wu, None; S. Holzer, None; B. Kiss, None; F. Frommlet, None; G. Fischer, None; C. Vass, None