Optimizing Structure–Function Relationship by Maximizing Correspondence Between Glaucomatous Visual Fields and Mathematical Retinal Nerve Fiber Models

Abstract

Purpose.: To introduce a method to optimize structural retinal nerve fiber layer (RNFL) models based on glaucomatous visual field data and to show how such an optimized model can be used to reduce noise in visual fields while probably preserving clinically important features.

Methods.: Correlation coefficients between age-adjusted deviation values of pairs of visual field test locations were calculated from 103 visual fields of eyes with moderate glaucomatous damage. Distances between those test locations were defined for various parameters of a mathematical RNFL model. Then, the correspondence between the structural and functional data was defined by the spread, or variance, of the correlation coefficients for all distances. The model parameters that minimized this spread constituted the optimized model. To reduce noise in visual fields, the optimized model was used to smooth visual field data according to the RNFL's structure. The resulting fields were compared with visual fields that were smoothed based on the regular testing grid.

Results.: The optimal parameters for the RNFL model reduced the variance of the correlation coefficients by 78% and were well within the range of parameters previously determined from fundus photographs. Smoothing the visual fields based on the optimized RNFL model strongly reduced noise while keeping important features.

Conclusions.: Mathematic RNFL models can be optimized based on visual field data, resulting in a strong structure–function relationship. Taking the RNFL's shape, as defined by such an optimized model, into account when smoothing visual fields results in better noise reduction while preserving important details.

Introduction

Standard automated perimetry (SAP) is an essential tool in the diagnosis of glaucoma and in monitoring glaucomatous progression. ¹ By determining local sensitivity values at various test locations (commonly called retinal sensitivity), it produces an estimate of the visual field. However, the reproducibility of such visual field sensitivity measurements is rather poor, due to high levels of noise caused by, for example, eye movement, blinking, or imprecisions of the measuring device, and most importantly the variability of the patient's response. ^2–5 A number of different factors have been shown to affect visual field measurements. ⁶ While noise may be reduced by using summary parameters, such as the mean defect (MD) or visual field index (VFI), these parameters hide important details of the visual field, such as its spatial structure. Noise may also be reduced by exploiting spatial and/or temporal relationships in the measurements. ² These methods aim to average or smooth points that are highly correlated, while points that are uncorrelated are excluded in this average. 

The spatial correlation in visual fields does not need to depend primarily on the distance between test locations as measured on the grid (i.e., Euclidean distance). Instead, this correlation is primarily determined by the structural position of the ganglion cells and their axons. ² The presence of wedge-shaped retinal nerve fiber layer (RNFL) defects in many glaucomatous eyes ⁷ underlines this. In images of these wedge-shaped defects, there is a strong correlation between the points within the wedge. However, points within the defect are largely uncorrelated with points outside the defect, regardless of their Euclidean distance. The same applies to locations within the visual field: in the case of an arcuate scotoma, locations within the scotoma are highly correlated, while neighboring points at the edge of the scotoma are largely uncorrelated. The Euclidean distance, which is for instance used in Gaussian filters, is not well suited to capture this characteristic of the correlation between test locations. Therefore, we would like to find a different measure that defines distance in a way that is more coherent with the correlations found in visual fields and what is known about the RNFL. 

It has been recognized previously by other authors that using conventional uniform or Gaussian filters may blur defects in visual fields ⁸ and is therefore of limited use. ⁹ Several authors developed improved filters that incorporate the RNFL structure. Crabb et al. ⁹ built a filter based on point-wise multiple regression of each test location on all other test locations, where the regression coefficients function as weights for the estimation of the sensitivity value at that first location. Gardiner et al. ² pursued a similar approach by applying regression with some constraints on the covariances to derive weights. Betz-Stablein et al. ¹⁰ used weights that were based on adjacency in the regular grid of test locations as well as physical adjacency to incorporate the spatial structure of glaucomatous damage into their model. 

Mathematic models of the RNFL morphology can be used to define distances (other than the Euclidean distance) between points that may then be used (as a covariate) to model the correlation between these points (e.g., the correlation between points is inversely related to their modeled distance). However, these models often contain free parameters (shape parameters) that allow them to optimize their shape, either to individual eyes or to a population of eyes, as in the present study. Previously, the parameter optimization was done by first tracing nerve fibers in fundus photographs and then selecting the shape parameters that resulted in the best match between observed and modeled nerve fibers. ^11,12 This is a very time-consuming approach and its accuracy may depend on the operator. ¹³  

In this paper, we introduce a method to optimize mathematic RNFL models based on visual field data of moderately glaucomatous eyes, thereby establishing a relation between structural and functional measurements. The proposed method does not require manual tracing of nerve fiber bundles and does not depend on user-selected parameters. Furthermore, we show how an optimized model can be used to improve the estimation of threshold deviation values, compared with methods that do not take RNFL morphology into account. 

Methods

In this section, the RNFL model is described, followed by a description of the data which was used for the optimization and evaluation of our methods. Then, a novel method to optimize the relationship between the structural RNFL model and the functional visual field data is introduced. This method results in an optimized RNFL model that is finally used to reduce noise in visual fields based on RNFL morphology. 

Model and Data

A mathematic model of the RNFL provides a description of the course of retinal nerve fiber bundles in a mathematic formula or set of formulas. Several of these mathematical RNFL models exist. ^11,12,14 However, the complexity of these models may be very different. For example, an axial model is very simple but quickly becomes inaccurate with increasing distance from the optic nerve head (ONH), while a model with many parameters may be able to largely capture normal morphologic variation between different eyes. With such an RNFL model, the distance between two test locations may be defined in a different way than the conventional Euclidean distance. For example, the distance may be measured along or across fibers. 

The model used in this paper describes the path of the nerve fiber bundles by a single equation that contains two parameters ¹¹ :  

 

 

Here, a test location in the visual field grid, given by its x and y coordinates, is defined by α and β, representing the angle of that point measured at the ONH and the fovea, respectively. A and B are two shape parameters that may be adjusted to adapt the model to the actual nerve fiber bundles. A more detailed description and interpretation of α, β, and the parameters A and B can be found elsewhere.¹¹ In Equation 1, the value of Ψ selects a specific nerve fiber bundle. The angle α^* at which this nerve fiber bundle meets the ONH can be derived from Equation 1 by setting y to 0 and β to Display Formula , yielding  

 

 

Airaksinen et al. ¹¹ optimized the shape parameters A and B on observed nerve fiber bundle patterns in fundus photographs. 

In this study, a generic procedure is presented to find optimal parameters for a population of eyes based on visual field data. The model described above is used as an example. 

To be able to match an RNFL model to functional data, visual field data needs to contain sufficient spatial structure. This is not the case for the visual fields of both healthy eyes and advanced glaucomatous eyes; they are either completely normal or largely defective. Visual fields of moderately glaucomatous eyes (with a mean defect between −6 dB and −12 dB) typically contain at least one visible focal defect. These focal defects provide the structure required to find the correspondence between the RNFL model and the functional data. Similar data selection strategies were employed in earlier work that defined structure–function maps. ¹⁵  

For our analysis, we used data from an ongoing study at the Rotterdam Eye Hospital. In that study, SAP was performed on a Humphrey Field Analyzer II (HFA; Carl Zeiss Meditec AG, Jena, Germany), with the white-on-white 24-2, Full-Threshold or SITA-Standard program. The study adhered to the tenets of the Declaration of Helsinki and informed consent was obtained from all participants. More details on inclusion and exclusion criteria may be found elsewhere. ¹⁶ For the present study, all visual fields measured in the years 1998 through 2011 with an MD between −6 dB and −12 dB were selected. If a patient had more than one field with an MD in that range, the field from the mean time point was included. In all, we used 103 visual fields from 103 patients, which is a subset of the data used by Bryan et al. ¹⁷ All data that was used in our analyses are freely available at http://orgids.com. Since our goal was to relate defects in the visual field (and not sensitivity values) to the RNFL model and the deviation from age adjusted normal threshold values represent the pattern of loss better than the raw threshold sensitivities, we used those threshold deviations for our analysis. The raw data were exported from the HFA and analyzed with the open source software R, version 3.0.0. ¹⁸ The observations directly above and below the ONH were not used in our analysis, leaving 52 test locations per eye. 

Optimizing the Structure–Function Relationship

Optimizing the structure–function relationship requires an objective function that expresses the strength of this relationship, based on some of its properties. Our approach is based on two assumed properties: first, the correlation coefficient between two test locations is inversely related to the distance between those two locations. Close test locations should, thus, have a high correlation, and test locations far apart should show low correlations. Second, test locations with similar distances should have similar correlations. 

The first property was implemented by modeling the relationship between any kind of distance and correlation as a nonparametric function that is forced to be monotonically decreasing. This function was defined by a penalized B-spline model with monotonicity constraints. ¹⁹ In a B-spline approach a weighted sum of a set of piecewise defined polynomials is used to approximate and/or interpolate the data. A more detailed description as well as a graphic illustration may be found elsewhere. ¹⁹ The second property was then evaluated by calculating the root mean square error (RMSE) of the fit of this function to the correlation data. Consequently, the objective function to optimize (i.e., minimize) is given by this RMSE. 

Visual field defects in glaucoma are commonly nerve fiber bundle defects, which implies that test locations within a bundle show similarly reduced sensitivity values, and hence are highly correlated. Correspondingly, test locations within the same nerve fiber bundle should therefore have a small distance, in agreement with the first property defined above. For the shape model used, the angle of the nerve fiber bundle at the ONH fulfills this criterion. The distance between two test locations was thus defined by the angular difference between the two corresponding nerve fiber bundles at the ONH. 

In the conventional representation of a visual field in x and y coordinates, the distance between test locations is given by the Euclidean distance. Figure 1 shows a scatter plot of the Pearson correlation coefficients as a function of the Euclidean distance, for all pairs of visual field test locations, similar to what has been shown before. ²⁰ The distances are quantized by the spacing of the test grid. The red line indicates the monotonically decreasing penalized B-spline fit, which had a corresponding RMSE of 0.3. For the whole range of measured Euclidean distances, the variability of the correlation coefficients was rather large, as illustrated by the large spread of the data around the red line. 

Figure 1

 

Scatter plot of correlation coefficients versus Euclidean distances.

When using the RNFL model to express the distance between test locations by the angular difference of the nerve fiber bundles at the ONH, a scatter plot like the one in Figure 2 results. For this example, A = 0.8 and B = 0.015 were selected; the course of the nerve fiber bundles for all 52 test locations are shown in the inset graph of Figure 2. The red line again indicates the result of the fit of the correlation curve. Compared with the Euclidean distance, the angular distance of this model resulted in a reduced variation of the data points around this line (RMSE = 0.23). 

Figure 2

 

Example scatter plot of correlation coefficients versus ONH angle differences. Inlay shows layout of the RNFL model with parameters A = 0.8 and B = 0.015.

Because the correlation coefficients are based solely on the test locations, the values of the correlation coefficients are independent of the type of distance that was used. Instead, changing the distance measurement results in different distance values for those correlation coefficients. Graphically, this results in moving points along the horizontal axis. 

The negative correlation coefficients in our data are due to the selection of moderately glaucomatous eyes. The visual fields of these eyes typically show considerable defects in one hemifield, but they have not progressed to such an extent that they show very large defects in both hemifields. As a result, if one hemifield shows a scotoma, the other hemifield is likely to be less affected. This inverse relationship is expressed by negative correlation coefficients. 

Our proposed method now optimizes the values for A and B (or other shape parameters in case of a different RNFL model) such that the RMSE is minimized. This optimized model shows the least variation of correlation coefficients around the fitted curve. 

Application: Model-Based Smoothing for Noise Reduction

Estimated visual fields are noisy representations of the true sensitivity at given retinal locations. To reduce noise, correlated measurements may be combined. This is often done by smoothing, which may be done temporally or spatially. The correlation between the measurements implies that they are not independent. Smoothing then reduces the noise, which is independently distributed, and retains the correlated signal. Correlation of measurements implies that their noise-free values can be represented by a smaller number of values. A noise-free representation of the 52 visual field test locations thus requires fewer than 52 parameters and the same holds for a smoothed visual field. A smoothed field is therefore considered to be a simpler representation than the observed visual field. It approximates the 52 observed thresholds deviation values but does not exactly reproduce them. 

In this paper, smoothing is performed by two-dimensional (2D) penalized B-splines. ^19,21 In brief, this method represents the field by a relatively large number of B-splines, which are then smoothed up according to some smoothing parameters. The choice of these smoothing parameters determines the complexity of the smoothed visual field. One way to express complexity is by the effective dimension (ED), which can be calculated for one-dimensional smoothing by ²¹  

 

 

Here, B is the B-spline matrix (also called design matrix), λ is the smoothing parameter, and D is the difference matrix used to penalize differences in adjacent coefficients. Since B and D are predefined for all fields, the ED only depends on the smoothing parameter used. Related concepts are effective or equivalent degrees of freedom. ²² More details, as well as the extension to EDs for 2D penalized B-splines, may be found elsewhere. ²¹  

The resulting difference between the original, nonsmoothed visual field and the smoothed visual field was expressed by the RMSE of all threshold deviation values. (Note that this is a different RMSE that was used to quantify the structure–function relationship when optimizing the shape parameters of the RNFL model.) Different smoothing methods (e.g., with or without taking the RNFL structure into account) with the same complexity resulted in different errors between the observed and smoothed fields. The method with the smallest RMSE for a certain ED may be considered the best, providing the most accurate representation of the field for all models with the same complexity. Alternatively, the method with the smallest ED for a certain RMSE may be considered optimal, because it is the least complex representation of all models with the same error. 

For smoothing according to the RNFL model, both the angle (for smoothing across bundles) and the length (for smoothing along a bundle) are used to determine which of the splines that represent the visual field are close to each other and should therefore have similar coefficients. The angle was given explicitly by Equation 2 while the length was numerically approximated from the model equations. ¹¹  

The reproducibility of estimated visual field sensitivities is limited, which sets a lower bound on the RMSE. Representing the field with a higher accuracy than this RMSE only amounts to modeling noise. Based on earlier described estimates of the sensitivity-dependent variability, ³ the average noise level for visual fields in our data set was determined to be approximately 4 dB. Note that this should be considered only as a rough indication of what accuracy may be considered reasonable; the variability of single visual field points can be considerably worse. 

In this paper, we consider three models for smoothing fields: isotropic x-y smoothing, anisotropic x-y smoothing, and smoothing across and along nerve fibers as given by our optimized RNFL model. Here, isotropic means that the same amount of smoothing was done along the x-axis (λ _x) and y-axis (λ _y). For anisotropic smoothing, this may be different. Isotropic smoothing for the RNFL model was not considered meaningful; smoothing across (λ _across) or along (λ _along) fibers operates on very different scales (ONH angle versus length). 

For the three models we evaluated the RMSE and ED for a range of smoothing parameters. The sets of optimal combinations of RMSE and ED were compared between the three models. 

Results

Optimizing the Structure–Function Relationship

To find optimal parameters for the RNFL model, a grid search was done on the shape parameters of the model: A was evaluated from 0.1 to 0.8 with a step size of 0.1, and B was evaluated from 0 to 0.04 with a step size of 0.005. For each combination of shape parameters, the monotonic penalized B-spline model was fitted and the corresponding RMSE was computed. Figure 3 shows a contour plot of the RMSE, with the minimum error at A = 0.4 and B = 0.02. Figure 4 shows the correlation-distance plot and nerve fiber layout of this optimized model. Again, the red line indicates the monotone fitted line. The RMSE was calculated to be 0.14 for the optimized model. 

Figure 3

 

Estimated contour plot of the RMSEs for the shape parameters A and B. The determined optimum was A = 0.4 and B = 0.02, as indicated by the cross.

Figure 4

 

Scatter plot of correlation coefficients versus ONH angle differences for all combinations of two test locations in the visual field for the optimized RNFL model.

Application: Model-Based Smoothing for Noise Reduction

All 103 fields were smoothed with different amounts of smoothing for isotropic and anisotropic x-y smoothing (λ _x and λ _y) and for smoothing along and across the nerve fibers as given by the optimized RNFL model (λ _along and λ _across). Different smoothing parameters resulted in a different complexity (expressed by ED) and representation error (expressed by RMSE). For x-y smoothing, the results are presented in Figure 5A. The solid and dashed lines mark the sets of optimal smoothing parameters for anisotropic and isotropic smoothing, respectively. Points on those lines are optimal in the sense that one criterion (ED or RMSE) cannot be improved without adversely affecting the other one. The results for smoothing according to the optimized RNFL model are shown in Figure 5B. Again, the solid line marks the optimal smoothing parameters. 

Figure 5

 

Calculated RMSE and ED for each combination of smoothing parameters λ _x and λ _y for x-y smoothing (A) and λ _along and λ _across for RNFL-based smoothing (B).

In Figure 6, the optimal sets of smoothing parameters for all three methods are displayed together with the estimated noise level of 4 dB. For x-y smoothing, anisotropic smoothing produces significantly better results (Wilcoxon signed-rank test, P < 0.0008 for each value of the ED ≤ 20) in terms of RMSE and ED. The best method, however, was smoothing based on the optimized RNFL model. Over the relevant range of ED (cf. estimated noise level) its optimal set lay below the two other optimal (P < 0.025). That means that RNFL-based smoothing can produce a less complex representation of a visual field with the same error or a representation with a smaller error for the same complexity. It also implies that the RNFL-based smoothing method provided the least complex representation for reproducing a visual field up to the level of its reproducibility: ED is 14.2 for the RNFL-based model compared with ED equaling 16.0 for anisotropic x-y smoothing, and ED equaling 18.0 for isotropic x-y smoothing. 

Figure 6

 

Comparison of the optimal set of smoothing parameters for all three smoothing methods. The horizontal line indicates the RMSE that corresponds to the approximated noise level of the measurements.

An example of an unprocessed visual field and the same field after smoothing based on the three different methods is shown in Figure 7. The smoothing parameters for all methods were chosen such that the same complexity (∼15 ED) resulted. All smoothed visual fields showed reduced spatial variability compared with the original field, but retained different levels of detail. Notice that the paracentral superior scotoma in the original field was largely retained by the optimized RNFL method (Fig. 7D), but it was smoothed away when applying smoothing based on the Euclidean distance. Clinically, smoothing away potentially important scotomas obviously is undesirable. 

Figure 7

 

Example of an unprocessed visual field of a moderately glaucomatous eye (A) and the same field after isotropic x-y smoothing (B), anisotropic x-y smoothing (C), and smoothing based on the optimized RNFL model (D).

Discussion

We have introduced a general method to optimize structural RNFL models based on functional data. This method can be applied to any model that associates test locations in the visual field with angles at the ONH. The same method can also optimize models that relate test locations to sectors at the ONH, but these discretized models may lead to oversimplified representations of the visual field, with insufficient spatial resolution. 

Introducing the difference between ONH angles as defined by the optimized RNFL model instead of the Euclidean distance as defined by the 24-2 test grid reduced the unexplained variance in the correlation coefficients by approximately 78%. This large decrease shows the strong structure–function relationship that exists between the morphology of the RNFL and the visual field threshold deviation values. In addition, the optimal shape parameters found with our method (A = 0.4 and B = 0.02) are well within the range of the shape parameters found by fitting the RNFL model to fundus photos. ¹¹  

It has been previously shown that a filter to smooth visual fields should probably include the anatomical structure of the RNFL, because filters that are based on the grid structure of the measurements can blur scotomas. ² This concurs with our results: when smoothing based on the x-y coordinates is applied to the visual field, either by isotropic smoothing (Fig. 7B) or anisotropic smoothing (Fig. 7C), the relatively small defect in the central visual field is blurred. Optimized RNFL models provide a great opportunity to incorporate the structure of the RNFL in the smoothing procedure. This is clearly illustrated in Figure 7D, which shows that this defect is preserved when smoothing based on the RNFL model. Due to the defect's location close to the fovea, such a defect should be considered important and it should therefore not be degraded by processing of the visual field. 

In our analysis, penalized B-splines were applied to perform smoothing, because of their flexibility and ease of implementation for different assumptions on the structure of the smoothing (e.g., grid-like x-y coordinates or coordinates based on the RNFL model). Smoothing of visual fields could be further improved by including weights in the fitting process. Such weights may, for example, be based on the variability of measurements for different sensitivity levels ^3,23 or on the density of ganglion cells over the measured area of the retina. ^24,25 It is important to note that our approach does not take into account that the strength of relation between points with the same distance is probably not uniform over the whole field, depending on the density of the ganglion cells. 

Other approaches to include physiological information on the RNFL in smoothing filters have been published previously. ^2,9,10,20 As mentioned above, and in contrast to other authors, ^14,26 we do not aim to produce an individualized RNFL map, but a general model that can provide a sufficient approximation to the RNFL structure, as shown in Figure 4. Our optimized RNFL model captures the variability in the individual RNFL structures in the RMSE, which contains the interindividual variability as well as the error made by approximating the real RNFL structure with the model. Our RNFL model differs from other models for instance in simplifying assumptions like symmetric hemispheres and that the ONH and fovea lie on a horizontal line. It has been shown previously that several ocular parameters may have significant influence on the layout of the RNFL. ²⁷ Although many of those parameters are currently not available for our study, it would be possible to include information on an individual eye's distance and angle between fovea and ONH into the model we used. Another factor to consider when comparing structural and functional measurements near the fovea is the displacement of ganglion cells. ²⁸ We estimated that the resulting additional error (∼1.7° for the four points nearest to the fovea) was insignificant in our analysis, but correction may be required when analyzing correlation of structure and function for an individual or when looking at denser grids. ²⁹  

Several other studies use the map by Garway-Heath et al. ¹⁵ that relates visual field test locations to ONH angles based on fundus photos. Recently, it has been shown that tracing nerve fiber bundles in such photos is associated with poor reproducibility, suggesting that the resulting nerve fiber maps have a high variability (Denniss J, et al. IOVS 2013;54:ARVO E-Abstract 1883). An optimized mathematical RNFL model, as proposed in the present work, does establish the same kind of relationship, but for continuous rather than discretized ONH angles, and can, therefore, also be used in these and other filters. Furthermore, ONH angles derived from an optimized RNFL model, as done in this study, relate any location in the visual field to an angle at the ONH, without being restricted to a certain pattern of test locations. 

Several previous studies explored structure–function relationships by calculating correlations or regression models between visual field sensitivities and measurements of the RNFL thickness, as determined by optical coherence tomography ^30–32 or by a Heidelberg Retina Tomograph. ³³ Our method is able to establish and optimize such a relationship without additional measurements, thereby showing a clear and strong structure–function relationship between the RNFL morphology and the visual field threshold deviation values. 

In the present work, we introduced a method to adapt a mathematic model of the RNFL morphology to visual field data. This new method clearly establishes a strong and significant relationship between structure and function, without requiring laborious tracing of RNFL bundles in fundus photographs. We also showed how this optimized model could be used to reduce noise in visual fields. Despite its large amount of noise, SAP plays an important part in the clinical management of glaucoma. Our results showed that including information on the morphology of the RNFL may reduce noise without affecting important defects in the visual field and may therefore be an important new aid in interpreting and processing visual field data in clinical glaucoma care. Ultimately, processing visual fields in any way, including smoothing, can only be considered beneficial if it increases the contrast-to-noise ratio, which is defined as the ratio of the useful range of the measurement (contrast) and its noise. ³⁴ Temporal averaging of repeated visual fields is the most straightforward method to achieve this, but it requires too many resources for clinical use. However, it may be used in a study setup to produce a ground truth. Future studies are needed to quantify the effect of filtering visual fields based on the presented approach and its contribution to the assessment of glaucoma and glaucomatous progression. 

Acknowledgments

Supported by grants from Stichting Wetenschappelijk Onderzoek Het Oogziekenhuis (Prof. Dr. H.J. Flieringa) and Stichting Combined Ophthalmic Research Rotterdam. 

Disclosure: N.S. Erler, Carl Zeiss Meditec (F), Heidelberg Engineering (F); S.R. Bryan, None; P.H.C. Eilers, None; E.M.E.H. Lesaffre, None; H.G. Lemij, Carl Zeiss Meditec (F), Heidelberg Engineering (F); K.A. Vermeer, Carl Zeiss Meditec (F), Heidelberg Engineering (F)