Evaluation of Structure-Function Relationships in Longitudinal Changes of Glaucoma using the Spectralis OCT Follow-Up Mode

The detection of glaucoma progression is an essential part of glaucoma management. Subjectivity of standard automated perimetry (SAP) prevents the accurate evaluation of progression, thus the detection of structural changes by optical coherence tomography (OCT) is attracting attention. Despite its objectivity, there is controversy about the appropriateness of the use of OCT, because many previous studies have indicated OCT results may not reflect the deterioration of visual field. A reason for this dissociation may be the test-retest variability of OCT, a major cause of which is misplacement of the measurement location. Recent advantages of spectral-domain OCT (SD-OCT), especially Spectralis OCT with an eye-tracking system (follow-up mode) enable measurement at approximately the same location as previous examinations. In addition to utilizing Spectralis follow-up mode, we introduced structure-function relationship map and nonlinear relationship between SAP and OCT results in considering structure-function relationship in longitudinal changes. The introduction of these two ideas in our study population improved the correlation between the SAP and OCT (R = 0.589 at most). The results of this study support the practical use of OCT in glaucoma progression but also stress the importance of focus on the corresponding focal changes and the consideration of disease severity.

However, their report included a limited number of patients and did not compare the correlation with visual field (VF) progression. Another reason is the complex relationship between structure and function. A number of cross-sectional studies have revealed a detailed correlation between sensitivity in the sectional visual field and the thickness of the corresponding sectorial cpRNFLT 23 . The most famous idea is the Garway-Heath map ( Fig. 2A) 24 , which associates OCT measurements and SAP results. We have reported another structure-function relationship map (Nakanishi map) based on clinical data collected at our glaucoma clinic 25 .
Another attempt to clarify the structure-function relationship in glaucoma is discussed in a review by Hood et al. 26 . Their argument consists of two issues. The first is how to calculate the sensitivity of sectorial visual field. They insist that the mean of the deviations should be calculated after each value is converted to the anti-log (Fig. 2B). The second concerns the nonlinear relationship between SAP and OCT values (Fig. 2C). In their cross-sectional clinical study, the exponential conversion of SAP values (or the logarithmic conversion of the sectorial cpRNFLT values) resulted in a linear relationship between SAP testing and OCT measurement 27 .
In this study, we investigated whether the application of the ideas and hypotheses (Fig. 2) discussed above could improve the consistency between longitudinal changes in SAP testing and OCT measurement. We also evaluated the ideas and hypotheses that improved the structure-function relationship in longitudinal changes in glaucoma.  ; the rate of global cpRNFLT change was −0.95 ± 0.07 µm/year (95% CI, −1.08 to −0.81 µm/year). The linear mixed models detected significant negative rates of change in 6 sectors of the Spectralis cpRNFLT. We also analyzed the linear mixed models on the sectors of sensitivity thresholds and total deviations in SAP and of cpRNFLT in Spectralis OCT defined by the structure-function relationship maps. These results are shown in Supplemental Tables 1-3. In all sectors, the significant longitudinal changes were detected.
To evaluate the correlation of longitudinal changes between VF and cpRNFL thickness, the baseline and rate of change of each sector of SAP or cpRNFLT in each eye were estimated by BLUPs, followed by visualization of the correlation by scatter plots and calculation of the Pearson's correlation efficient using the BLUPs of the VF and cpRNFLT sectors. The scatter plots of the BLUPs of the SAP MD and the global cpRNFLT are shown in Supplemantal Fig. 1. The correlation coefficients of the estimated baseline and the rate of change were 0.548 (moderate relationship) and 0.358 (weak relationship), respectively. To increase the correlation between VF and cpRNFLT, we introduced three methods (Fig. 2): the structure-function relationship map (Garway-Heath map and Nakanishi map), the arithmetic mean of the sector VF, and a nonlinear model of VF and cpRNFLT. The results of scatter plot analysis after introducing these methods are shown in Fig. 3 and Supplemental Figs. 2-5. In terms of structure-function relationship map, analysis of the rates of change ( Fig. 3 and Supplemental Fig. 3B) revealed that the correlation coefficients between corresponding sectors in the Garway-heath map were significantly higher than those between non-corresponding sectors (P < 0.01, Wilcoxon test).
The baselines of the scatter plots (Supplemental Figs 2 and 3A) revealed that the introduction of a nonlinear model to either SAP or cpRNFLT straightened the logarithmic correlation curve. Although in the analysis of rates of change ( Fig. 3 and Supplemental Fig. 3B) it did not increase the correlation coefficients significantly (P = 0.077, Wilcoxon test), nonlinear models succeeded in improving the structure-functional correlation of longitudinal changes especially in the nasal inferior sector of the SAP TD with a nonlinear model of the Garway-Heath map (temporal superior sector of cpRNFLT in the linear model); the correlation coefficient of longitudinal evaluations in this sector of these models was 0.589. In the superior sector of the SAP (inferior sector of the cpRNFLT), the highest correlation coefficient was 0.538 between superior SAP sensitivity threshold (nonlinear model; Nakanishi map) and temporal inferior cpRNFLT (linear model). The scatter plots with the highest correlation coefficient in the superior or inferior sector in each structure-function relationship map are shown in Fig. 4. In terms of calculation of sectorial sensitivity of VF, arithmetic mean did not significantly improve the correlation compared with geometrical (conventional) mean in Fig. 3 and Supplemental Fig. 3B (P = 0.37, Wilcoxon test).

Discussion
In this study, we investigated the correlation of the longitudinal changes between the SAP indices and the cpRN-FLT measured using Spectralis Follow-up mode. Comparing our previous report 19 , the correlation coefficiency of the rate of change between SAP MD and cpRNFLT became higher by introducing Follow-up mode (R = 0.358, with Follow-up mode; R = 0.16 without Follow-up mode). Moreover, the introduction of the structure-function relationship map and nonlinear relationship between the VF and cpRNFLT improved the correlation between the SAP indices and cpRNFLT in corresponding sectors (R = 0.589 or 0.538; compare Fig. 4 with Supplemental Fig. 1). These results justify not only the application of those two ideas (structure-function relationship map and nonlinear relationship in Fig. 2) in considering the structure-function relationship in longitudinal studies but also practical use of OCT in evaluating the glaucoma progression: progressive thinning of cpRNFL indicate deterioration of visual field defects. The usefulness of introducing structure-function relationship map was expected for the following reasons; the focal structural or functional change in glaucoma is not necessarily followed by changes in other regions, as glaucoma is characterized by focal damage of the visual field and retinal nerve fiber layer (nasal step, Bjerrum scotoma, focal rim thinning or nerve fiber layer defect [NFLD]), not by their global deterioration 1,28,29 . Previous studies of the structure-function relationship in glaucoma progression 11,12 that compared the global indices revealed a relationship; the mean deviation in Humphrey 24-2 only corresponded with the temporal side of the cpRNFL. Thus, recent studies have also focused on progressive RNFL thinning as determined by event analysis (Guided Progression Analysis [GPA]), using the RNFL thickness map acquired by SD-OCT 30,31 . Although this strategy is informative for the discrete evaluation of progression, the detailed relationship between longitudinal changes of VF and cpRNFL remains to be determined (We will discuss the problem of event analysis later).
We used two types of structure-function relationship map in this study: the Garway-Heath map and the Nakanishi map. While the nasal inferior sector of the Garway-Heath map was correlated most strongly with the temporal superior sector of the cpRNFLT, the strongest correlation was in the superior sector of visual field came from combination with the superior sector of VF and the temporal inferior (not inferior) sector of cpRNFLT in Nakanishi map. This result replicated the correlation in the original cross-sectional report from Nakanishi et al. 25 . Hood et al. depicted the schematic model of the glaucomatous damage of the macula, indicating that retinal ganglion cells within "vulnerable region" in the inferior visual field including macula project to the inferior quadrant of the disc (corresponding to the temporal inferior sector of cpRNFLT in Nakanishi map) 32 . Our results may support Hood's schematic model even in the longitudinal study. While the correlation differed according to the structure-function relationship map, we could not determine the superiority of either map because this study included cases that were also included when creating the Nakanishi map. As described in Table 1, the participants in this study were myopic like those in the studies from Asia and may have different characteristics from those who contributed to induce Garway-Heath map.
The application of the nonlinear relationship between the VF and cpRNFLT also improved the correlation of the evaluations of the longitudinal changes in some sectors, indicating that the disagreement between the VF and RNFLT in previous longitudinal studies may have been due to the hypothesis that the longitudinal relationship between the VF and the RNFLT is linear. The nonlinear relationship implies that the simple application of event analysis to both VF and RNFLT may mislead judgment. In the detection of glaucoma progression, the baseline of the VF or RNFLT should be considered. The thinning of thick RNFL rarely affects the visual field, but slight changes in thin RNFL may lead to remarkable VF deterioration. The famous continuum proposed by Weinreb insisted that structural changes precede functional changes 33 . Abe et al. also showed in a longitudinal study that the detection rate of progression by SD-OCT (but not by SAP) was higher in eyes with less severe disease. An increase in disease severity increased the chance of detection by SAP, but not SD-OCT 34 . The nonlinear model between SAP and cpRNFLT supported their hypothesis and reports. We should consider the disease severity regardless of the use of SAP or SD-OCT to detect glaucoma progression.
We also evaluated a way to calculate the average sensitivity or deviation of the SAP sector, but the difference was not significant. According to Weber-Fechner's law 35 , SAP indices should be calculated on a logarithmic scale, as visual acuity is. This means Weber-Fechner's law also supports the conventional method of the use of the mean, total deviation, or MD slope. Our results could not determine the correctness of Weber-Fechner's law or  37 , because in this study changes in treatment were not prohibited. The introduction of a nonlinear model in this study was merely to improve the structure-function relationship in terms of longitudinal changes. This study had several limitations. The first was the small sample size of both cases (patients) and experiments. To avoid the incorrect evaluation of longitudinal changes, we applied linear mixed effect models and BLUPs 38 , but a larger cohort may reveal the precise and detailed relationship between VF and RNFLT. The second was the exclusion of cases with ERM or retinoschisis from the study. As the strategy of this study relied on the hypothesis that decrease of RNFLT directly results in VF progression, we had to exclude other factors that could influence RNFLT without changing the VF. As stated in the first section of the Results, many glaucoma cases are accompanied by ERM or retinoschisis. Determination of the optimal evaluation method of the structural longitudinal changes in those cases is goal for future research. The third was the modeling used in this study. As previous studies have suggested 26 , cpRNFLT is never less than 30 µm even in advanced cases, as components other than the nerve fiber layer (vessels and glial tissues) remain 39 . Correlation of this floor effect to the model was not possible, as in some cases the sectorial thickness of the cpRNFL measured by the Spectralis OCT decreased to <10 µm. Another cross-sectional study could not detect the floor effect 40 , so we must verify its validity in future studies.
In conclusion, we evaluated the correlation between structural and functional changes of glaucoma in a longitudinal study. The results of this study support the practical use of OCT in glaucoma progression but also stress the importance of focus on the changes in the corresponding sectors and the consideration of disease severity. Future studies will reveal that other modalities to measure the structural changes are also correlated with functional changes or with those detected by cpRNFLT.

Methods
This prospective longitudinal study adhered to the tenets of the Declaration of Helsinki and was approved by the Institutional Review Board and Ethics Committee of the Kyoto University Graduate School of Medicine. Study subjects were prospectively enrolled at Kyoto University Hospital between April 2008 and October 2016, and informed consent was obtained from all patients. All patients included in this study had a normal angle on gonioscopy and a Snellen equivalent best-corrected visual acuity (BCVA) of 20/20 or better to ensure high imaging quality. Subjects also had a typical glaucomatous VF defect (as seen on SAP using the Swedish interactive threshold algorithm [SITA] 24-2) with corresponding optic disc changes (narrowing of the neuroretinal rim) and/or thinning of the retinal nerve fiber layer (RNFL), based on the Ocular Hypertension Treatment Study (OHTS) criteria 41 . Eyes were excluded from analysis if they had a cataract that could affect visual function, vitreoretinal disease, pathologic myopia with patchy chorioretinal atrophy, a lacquer crack lesion or choroidal neovascularization, uveitis, or prior ocular surgery including laser therapy. The use and changes of topical IOP-lowering agents was permitted during this study. Patients were also excluded from participation if they had a neurological disease, diabetes mellitus, or any other systemic disease that might affect the eye or the VF (such as a cerebrovascular event, uncontrolled hypertension, or blood disorders). Data from patients who developed an epiretinal membrane or retinoschisis on the macula or around the optic disc during follow-up were excluded because they could influence retinal thickness measurements. To avoid the introduction of selection bias, data from patients who underwent ocular surgery during the follow-up period were included until the date of surgery. Clinical examinations. All subjects underwent a comprehensive baseline ophthalmic examination, including measurement of IOP with a Goldman applanation tonometer and uncorrected and best-corrected visual acuity with a Landolt chart at 5 m. Subjects also underwent slit-lamp examination, gonioscopy, stereoscopic optic disc photography (3-Dx simultaneous stereo disc camera; Nidek, Gamagori, Japan), SAP with a Humphrey Visual Field Analyzer (Carl Zeiss-Meditec, Dublin, CA, USA) using the 24-2 SITA standard testing protocol (HFA + 24-2 SITA), and OCT examinations (Spectralis HRA + OCT system; Heidelberg Engineering, Heidelberg, Germany). All subjects had at least four VF and OCT examinations during the follow-up period. The first and last VF and OCT examinations were performed within 3 months.
Visual field examinations. Only reliable VF tests were used for analysis. The reliabilities were defined as the fixation loss (<20%), false-positive error rate, and false-negative rate (<33%). A glaucomatous VF defect was defined as (1) glaucoma hemifield test results outside normal limits; (2) more than three significant (P < 0.05) and one highly significant (P < 0.01) non-edge contiguous points on the same side of the horizontal meridian as in the pattern deviation plot; or (3) a pattern standard deviation (PSD) < 5% in an otherwise normal VF. All VF findings were confirmed in at least two consecutive testing sessions. To compensate for the subject's learning process, the first VF test was not included in the analysis if it was the first VF test the subject had taken.
Spectral-domain optical coherence tomography imaging. The Spectralis HRA + OCT system was used to evaluate the cpRNFLT. The eye tracking system of this instrument allows accurate averaging of up to 100 B-scans (7 μm axial resolution) acquired at an identical location, which can efficiently reduce speckle noise, and enables measurement at approximately the same location during every examination. For cpRNFL imaging, we performed a 3.46 mm diameter circular scan centered on the optic disc and consisting of 1536 A-scans. A total of 16 scans were acquired and averaged to obtain the final scan used in analysis. For repeated examinations of cpRNFLT, only the measurements taken in follow-up mode were used, as they measured the same location as the first examination.

Structure-function relationship map.
To more accurately evaluate the sectorial relationship between SAP and cpRNFL thickness in OCT, we introduced the "structure-function relationship map." In this study, we used the Garway-Heath map 24 (most famous and widespread in glaucoma society) and the Nakanishi map (produced from a clinical study at our hospital; Fig. 2A) 25 . To calculate the sectorial cpRNFL thickness, we exported 768 values of the cpRNFLT along the 360-degree OCT scan using RNFL export software (Heidelberg Engineering). The arithmetic means of the values corresponding to each sector were used as the sectorial thickness of cpRNFL. For the Nakanishi map, in which the fovea-optic nerve head (ONH) center angle should be considered, the averaged values were chosen after adjustment for the influence of the fovea-ONH angle measured using the built-in software. In SAP, the sectorial threshold sensitivities or total deviations (TDs) were also calculated from the means of the values corresponding to the defined sectors in the map. In this study, the means of the VF values were calculated on an apostilb scale both arithmetically (anti-log TD) and geometrically (arithmetically on a dB scale; conventional TD), as it remains controversial whether the total deviation values should be anti-logged before averaging 26 . Both types of means were analyzed.
Nonlinear relationship between structure and function. When comparing the longitudinal changes of the VF and OCT results, the nonlinear relationship between structure and function in the cross section must be considered. To correct the nonlinear relationship to a linear one, the sectorial means of the VF values (threshold sensitivities and total deviations) were converted to exponential values (1/Lambert [1/L] scale) using the following formula: dB 10 log (1/L) log (1/L) (1) 10 1 259 Similarly, the sectorial means of cpRNFL thickness were converted to logarithmic values (defined as LogRNFLT in this study): LogRNFLT 10 log (cpRNFLT) log (cpRNFLT) These exponential VF values (1/L) and logarithmic OCT values (LogRNFL) were also used to analyze the longitudinal change in each case as nonlinear models. The unconverted values were used in linear models in the next analysis (described in the next subsection).

Statistical analysis.
The longitudinal time trends of the outcome measures were evaluated using linear mixed models fitted with random intercepts and coefficients at both the subject and eye levels. Linear mixed-effects modeling allowed the determination of correlations among repeated measurements and the use of two eyes from the same subject 38 . Linear mixed-effects models can also manage data sets with multiple missing data points or with high variation in examination times. The following equation describes the corrections applied to our data. where Y ijt is the individual measurement at visit t; β 0 and β 1 are the fixed-effects coefficients; ζ 0j and ζ 1j are the random patient effects associated with the intercept and time slope; and ζ 0i|j and ζ 1i|j are the random effects associated with including both eyes of a single subject. As we evaluated the simple relationship of the longitudinal changes between the VF values and cpRNFL thickness, as we usually do in clinical situations, we did not enter other fixed effects such as age, axial length, or imaging quality.
To evaluate the individual heterogeneity of baselines and time trends in each examined eye, estimates of best linear unbiased predictors (BLUPs) were calculated and analyzed in this linear mixed model. The relationships between VF and OCT cpRNFL in longitudinal changes were evaluated by correlation coefficients calculated from the BLUPs of the VF values (arithmetic or geometrical means × linear or exponential values = 4 patterns) and OCT values (linear or logarithmic values). To visualize the strength of the correlation coefficients, the color of the dots in the scatter plot varied according to the strength of the correlation coefficient ( Fig. 3; Supplemental Figs 2-5).
All P-values presented are two-sided values. The statistical significance was defined as P < 0.05. The differences of correlation coefficients in each model (structure-function relationship map, nonlinear model and calculation of mean deviation. Also see Fig. 2) were compared using the Wilcoxon signed rank test. All analyses were performed using R ver. 3.3.2 (R Foundation for Statistical Computing, Vienna, Austria) and SAS ver. 9.4 (SAS Institute, Inc, Cary, NC) statistical software.

Data Availability Statement
We uploaded the raw data used in this work as Supplement Dataset with the manuscript.