Modeling Changes in Corneal Parameters With Age: Implications for Corneal Disease Detection

Journal Pre-proof Modelling changes in corneal parameters with age: implications for corneal disease detection Janelle Tong, Jack Phu, Michael Kalloniatis, Barbara Zangerl PII:

S0002-9394(19)30415-5

DOI:

https://doi.org/10.1016/j.ajo.2019.08.014

Reference:

AJOPHT 11052

To appear in:

American Journal of Ophthalmology

Received Date: 24 February 2019 Revised Date:

19 August 2019

Accepted Date: 19 August 2019

Please cite this article as: Tong J, Phu J, Kalloniatis M, Zangerl B, Modelling changes in corneal parameters with age: implications for corneal disease detection, American Journal of Ophthalmology (2019), doi: https://doi.org/10.1016/j.ajo.2019.08.014. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 The Author(s). Published by Elsevier Inc.

ABSTRACT Purpose: To apply computational methods to model normal age-related changes in corneal parameters and establish their association with demographic factors, thereby providing a framework for improved detection of subclinical corneal ectasia (SCE). Design: Cross-sectional study Methods: 117 healthy participants were enrolled from Centre for Eye Health (Sydney, Australia). Corneal thickness (CT), front surface sagittal curvature (FSSC) and back surface sagittal curvature (BSSC) measurements were extracted from 57 corneal locations using the Pentacam HR from one eye per participant. Cluster analyses were performed to identify locations demonstrating similar variations with age. Age-related changes were modeled using polynomial regression with sliding window methods, and model accuracy was verified with Bland-Altman comparisons. Pearson’s correlations were applied to examine the impacts of demographic factors. Results: Concentric cluster patterns were observed for CT and FSSC, but not BSSC. Sliding window analyses were best fit with quartic and cubic regression models for CT and FSSC/BSSC respectively. CT and FSSC sliding window models had narrower 95% limits of agreement compared to decade-based models (0.015 mm versus 0.017mm and 0.14mm versus 0.27mm respectively), but were wider for BSSC than decade-based models (0.73mm versus 0.54mm). Significant correlations were observed between CT and astigmatism (p = 0.02-0.049) and FSSC and BSSC and gender (p = <0.001-0.049). Conclusions: The developed models robustly described aging variations in CT and FSSC, however other mechanisms appear to contribute to variations in BSSC. These findings and the identified correlations provide a framework that can be applied to future model development and establishment of normal databases to facilitate SCE detection.

Modelling changes in corneal parameters with age: implications for corneal disease detection Janelle Tong1,2, Jack Phu1,2, Michael Kalloniatis1,2, Barbara Zangerl 1,2 1. Centre for Eye Health, University of New South Wales (UNSW), Sydney 2052, NSW Australia 2. School of Optometry and Vision Science, UNSW, Sydney 2052, UNSW Australia

Corresponding author: Dr. Barbara Zangerl Centre for Eye Health, UNSW Sydney 2052, NSW Australia Email: [email protected] Phone: +61 2 8115 0793 Short title: Models of corneal aging changes Supplemental Material available at AJO.com.

INTRODUCTION The cornea forms the primary refractive surface of the eye and plays a pivotal role at maintaining clear focus at the retina. Keratoconus, a progressive condition characterized by corneal protrusion and thinning, can therefore result in significant visual impairment.1, 2 Early detection of keratoconus is essential to identify individuals at risk of developing visual symptoms and those who are unsuitable for refractive surgery.2-5 However, the earliest stage of keratoconus, sometimes termed ‘forme fruste’ or subclinical corneal ectasia (SCE), remains difficult to differentiate from normal corneas, even with the advent of advanced corneal imaging techniques.6, 7 This is likely related to significant inter-individual variation in corneal parameters within the normal population,8-10 with the wide range of normative values contributing to considerable overlap between normal and SCE populations. Normal, age-related changes of corneal thickness (CT) and curvature (CC) are currently not well characterized, with several studies reporting relative corneal thinning and flattening with increasing age10-17 and others having found no such correlations.18, 19 Similarly, consensus on the existence of systematic correlations between corneal parameters and other demographic factors, such as gender, ethnicity and refractive error, is yet to be established.10-13, 16, 18, 20 Measurement noise as a result of normal inter-individual variability is once again likely to contribute to these inconsistencies,8-10 as it may affect the ability of statistical tests to produce replicable and meaningful analyses.21 By clarifying if and how corneal parameters vary with age and other demographic factors, a narrower range of expected normal values could be generated given individual-specific demographic information, which may then facilitate identification of measurements outside of normal limits and improve detection of SCE. Cluster methods have demonstrated success in minimizing inter-individual variability via pooling of individual data points displaying statistically indistinguishable characteristics.21-23 They have been applied in various areas of vision science, including characterization of aging changes in the macular ganglion cell layer (GCL), identification of visual field isocontours, optimization of indices to improve detection of glaucomatous visual fields and characterization of glaucomatous progression patterns.21-26 Given the ability of cluster analyses to meaningfully describe normative biological data by reducing or redistributing measurement noise, applying these principles to corneal measurements may aid in clarifying patterns of age-related change. The present study utilized cluster methods to identify corneal locations showing similar change in CT and CC with age in a normal cohort. From our understanding of corneal variations with increasing eccentricity4 and our previous work on macular GCL aging,21 we hypothesize that normal corneal aging follows similar concentric patterns that may be revealed using a cluster analysis approach. The investigated cluster patterns facilitated the development of regression models to identify patterns of age-related change in corneal parameters, and we sought to validate these models to confirm their robustness. We also examined systematic correlations between CT and CC and demographic factors. This work will improve our understanding of age-related corneal changes and the impacts of demographic factors on variations in corneal measurements, thereby aiding establishment of normal corneal databases and future model development, and has the potential to be applied to algorithms for improved detection of SCE.

1

METHODS Participant Recruitment This was a cross-sectional study where data were retrospectively collected from patients attending the Centre for Eye Health (CFEH, Sydney, Australia) for corneal assessment or anterior eye imaging from January 2017 to December 2017. CFEH is an intermediate-tier, referral-based clinic primarily servicing patients referred by community optometrists in the metropolitan Sydney region. The clinical protocol for corneal assessment at CFEH includes thorough ocular and medical history, visual acuity measurement with habitual spectacle correction, slitlamp examination, Pentacam HR Scheimpflug imaging (Oculus, Germany), Medmont E300 corneal topography (Medmont, Australia), Irx3 wavefront aberrometry (Imagine Eyes, France) and Cirrus or Spectralis anterior optical coherence tomography (Carl Zeiss Meditec, Inc., USA and Heidelberg Engineering, Heidelberg, Germany). Refractive error information was obtained from referral data, while self-reported demographic information and medical history were extracted from intake forms completed prior to examination. Signed consent to use clinical and demographic data for research purposes was obtained for each participant as per ethics protocols approved by the University of New South Wales Australia Human Research Ethics Advisory (UNSW Australia HREA) panel, and this study adhered to the tenets of the Declaration of Helsinki throughout its duration. The following criteria determined eligibility for inclusion in this study: absence of corneal ectasia in both eyes confirmed by two experienced clinicians based on clinical examination and normal imaging results, including corneal topography and tomography scans of adequate quality that did not show patterns suggestive of ectasia, visual acuity better than 20/25 (logMAR > 0.1), spherical equivalent refractive error between +3.00 diopters and -6.00 diopters and refractive astigmatism <2.00 diopters, for participants up to 60 years of age based on the critical age range for the detection of ectatic changes.27, 28 As suitable threshold values for summary parameters, such as the Belin/Ambrósio enhanced ectasia D value, to distinguish normal and SCE eyes with high sensitivity and specificity are yet to be established in the literature,7 we did not rely on such parameters in isolation to confirm absence of keratoconus. Additionally, as corneal astigmatism is an instrument-specific metric and corneal topography and tomography scans were manually inspected in a holistic manner for signs suggestive of ectasia, corneal astigmatism was not identified as an individual inclusion/exclusion criterion in addition to refractive astigmatism. Exclusion criteria were defined as any history of contact lens wear or ocular surgery, a history of diabetes or other systemic neuropathies, current pregnancy or hormone therapy, including contraceptive pill and hormone replacement therapy, and an unavailable complete medical history. Pentacam Measurements Corneal measurements were obtained from the Pentacam HR. For eligible participants, scan quality was manually assessed by an experienced clinician (JT); those that did not show an ‘OK’ quality score according to the Examination Quality Specifications provided by the instrument review software (Pentacam Version 6.08r27, Oculus, Germany) or with data missing as a result of imaging artefacts, most commonly shadowing from eyelids, were excluded from analyses. One eye per participant was included at random where both eyes met inclusion criteria, otherwise the eye meeting inclusion criteria was included in this study. CT and CC maps were directly exported from the instrument review software. CT, front surface sagittal curvature (FSSC) and back surface sagittal curvature (BSSC) measurements were extracted from the default 57-point locations on the ‘large color map’ setting as per the instrument review software, covering a total central corneal diameter of 8.00mm (Figure 1). While Montalbán et al.11 reported significant correlations between mean front and back corneal radius in normal eyes, similar relationships have not been established for pointwise measurements, and given the importance of the posterior corneal surface in early detection of corneal ectasia3, 29 we chose to consider BSSC measurements obtained in addition to FSSC measurements. FSSC and BSSC were chosen over corneal elevation measurements as they do not rely on reference sphere calculations that may vary between 2

examinations.30, 31 All left eye scans were transposed to right eye format by horizontally transposing data across the vertical midline, to allow for direct location-specific comparisons between scans. Cluster Methods Measurements were grouped by participant age in decades (≤10 years, 11-20 years etc. up to 51-60 years, Table 1), and average CT, FSSC and BSSC were calculated for each decade bracket at each of the 57-point locations. Cluster analyses were performed using SPSS Statistics Version 25.0 (IBM Corporation, New York, NY, USA). Hierarchical clustering utilizes an agglomerative technique to group statistically similar measurements into successively larger clusters, which were defined through within-groups linkage and squared Euclidean distance to compute distances between clusters.22 The separability of the initial clusters obtained using hierarchical clustering was calculated using the d’ statistic22: ݀ᇱ =

|‫ݔ‬ଵ − ‫ݔ‬ଶ | ඥ0.5 × ሺߪଵ ଶ + ߪଶ ଶ ሻ

Clusters were systematically merged, starting with the cluster pair with the lowest d’ value, until all comparisons reached our threshold for statistical separability set at d’ >1 at a final number of clusters k. Unlike hierarchical clustering, k-means clustering involves the assignment of points into a pre-determined number of clusters, which is chosen prior to the algorithm being performed. To enable direct comparisons with hierarchical cluster patterns, the desired number of clusters was set at k and separability of clusters was calculated in an identical fashion. Regression Analysis Average CT, FSSC and BSSC were calculated per decade bracket and cluster and fitted with linear and second order polynomial (quadratic) regression curves, and sum-of-squares F tests were performed using Graphpad Prism Version 7.04 (La Jolla, CA, USA) to compare goodness of fit of these models to the dataset. Sliding window analyses were also performed using decade windows, with measurements averaged over each window (i.e. ≤10 years, 1120 years...51-60 years, then ≤11 years, 12-21 years…52-60 years and so forth), to minimize the potential bias of using semi-arbitrary decade groupings while maintaining the advantages of pooling similar measurements in reducing measurement noise. The results of the sliding window analyses were fitted with the most appropriate polynomial regression curves using Graphpad Prism Version 7.04. To determine whether regression models varied between clusters, one-way ANOVA with Tukey’s multiple comparison tests of the regression coefficients was performed for each cluster pattern and corneal parameter. The points of inflection for each polynomial regression model were calculated using Matlab Version R2017b (Mathworks, Natick, MA, USA) to identify time points where significant points of change in the investigated corneal parameters occurred. Regression Model Validation To verify the accuracy of the derived polynomial models, participants in the 31-40 year bracket were selected as a 35-year-old age-normative cohort (Table 1), approximating the mean age of the entire cohort (mean 32.9 years, Table 1). The regression curves derived from the sliding window models and models derived from single averaged measurements per decade bracket (hereby referred to as decade-based models) were utilized to convert the 57-point measurements for CT, FSSC and BSSC measurements from all other participants to 35-year-old equivalent data. Bland-Altman comparisons were then performed between the age-corrected and age-normative cohorts. Statistical Analyses Statistical analyses were performed using SPSS Statistics Version 25.0 and GraphPad Prism Version 7.04. Normality of distribution of CT, FSSC and BSSC at each location was determined using the D’Agostino and Pearson normality test. Correlations between point 3

measurements of CT, FSSC and BSSC and ethnicity, spherical equivalent refractive error and refractive astigmatism were analyzed using Pearson’s correlations, while correlations with gender were analyzed using unpaired Welch’s t-tests, as previous studies have reported correlations between corneal parameters and these factors.10-13, 16, 18, 20 RESULTS Study Cohort Of 117 participants included in this study, 51% were female (Table 1). The most prevalent self-reported ethnicities were Caucasian (53%) and East Asian (23%); other less commonly reported ethnicities included Indian (5%), Middle Eastern (4%), Central or South American (3%), Pacific Islander (3%), Indigenous Australian (2%) and mixed descent (5%). Additionally, only 8% of the entire cohort exhibited against-the-rule corneal astigmatism; although this sample was too small to investigate prevalence with increasing age, 56% of this sample with against-the-rule astigmatism were within the oldest age bracket (51-60 years), suggesting that the previously reported association of increasing against-the-rule astigmatism with age may also be present in the studied cohort.10, 32 Meanwhile, 82% of the cohort with corneal astigmatism equal to or greater than 0.75 diopters demonstrated withthe-rule corneal astigmatism (Supplementary Table 1). Cluster Analysis and Regression Models Hierarchical clustering identified 7, 8 and 7 statistically separable clusters for CT, FSSC and BSSC measurements respectively. Each cluster represents groups of test locations demonstrating statistically similar age-related change, which are conveying using the colored cluster patterns in Figure 2. Concentric patterns of age-related change were observed for CT and FSSC measurements for both hierarchical and k-means cluster analyses, however concentricity was not maintained in cluster patterns describing BSSC measurements. Regression analyses were applied to average measurements pooled per decade bracket and per cluster according to the obtained cluster patterns, to characterize age-related variation (Figure 3A). Change per decade for clustered CT measurements was best fitted using linear regression, whereby the central clusters showed a gradual increase in CT and the peripheral clusters showed a gradual decrease with age, however the goodness of fit of these models was superior for central clusters, as indicated by higher coefficients of determination (R2 = 0.002-0.62, Supplementary Tables 2 and 3). For FSSC and BSSC measurements, linear and quadratic regression models alternately provided superior fits for different clusters (R2 = 0.05-0.98 and R2 = 0.006-0.85 respectively, Figures 3B and 3C, Supplementary Tables 2 and 3). Accounting for the possibility that decade-based models artificially reflect aging changes due to the semi-arbitrary decade brackets, sliding window analyses using decade windows were performed. In such analyses, fourth degree polynomial (quartic) regression models provided the most appropriate fit for CT measurements (R2 = 0.53-0.75, Figure 4A and Supplementary Tables 4 and 5), and FSSC data were best fitted with third degree polynomial (cubic) regression models (R2 = 0.42-0.74, Figure 4B and Supplementary Tables 4 and 5) across all clusters and cluster patterns. One-way ANOVA of the regression coefficients for CT and FSSC models showed no statistically significant difference between clusters (p = 0.115>0.999), indicating that these parameters changed in a similar fashion with age regardless of location. The 3 stationary points derived from CT quartic models were consistent across all clusters (mean 17.2 ± 0.3 years, 32.6 ± 0.9 years and 51.7 ± 1.2 years, Table 2), as were the 2 stationary points derived from FSSC cubic models (mean 21.4 ± 2.5 years and 51.0 ± 4.0 years, Table 2). For BSSC data, although a cubic regression model best described age-related changes, regression patterns varied more markedly between clusters (Figure 4C and Supplementary Tables 4 and 5). This was reflected in the statistically significant difference in regression coefficients between cluster 5 and all other clusters in regression models derived from kmeans clustering (p = 0.0006-0.044), although no statistically significant difference between 4

clusters was noted for models based on hierarchical clustering (p = 0.116->0.999). Differences between clusters were also evident in the widely varying stationary points; while most clusters displayed 2 stationary points (mean 22.2 ± 8.2 years and 42.4 ± 5.2 years, Table 2), several displayed only 1 stationary point (mean 45.8 ± 6.7 years) and 1 peripheral cluster showed no stationary points, indicating a continuous increase in BSSC in this location with age. Model Validation using Age-Normative Comparisons Bland-Altman comparisons were performed to verify the accuracy of the sliding windowbased and decade-based models in describing age-related variation. The widths of the 95% limits of agreement (LoA) intervals were expressed as percentages of the range of CT, FSSC and BSSC values to facilitate comparisons between the different parameters. 95% LoA interval widths were similar between CT and FSSC, but larger in BSSC measurements (Table 3 and Figure 5). Comparisons between sliding window and decade-based models suggested that the sliding window model based on hierarchical clusters most accurately described variations in CT, and both hierarchical and k-means clustering-based sliding window models better predicted changes in FSSC, as indicated by narrower 95% LoA and bias closer to zero. In contrast, decade-based models better described changes in BSSC with age. Correlations with Demographic Factors Small but statistically significant correlations were observed between inferocentral CT measurements and magnitude of refractive astigmatism (r = 0.18 to 0.20, p = 0.03-0.048, Figure 6), which was subsequently confirmed with the observed significant correlations between CT measurements and corneal astigmatism (r = 0.18 to 0.21, p = 0.02-0.049). Significant higher CC measurements in females compared to males were also observed in FSSC measurements across the entire cornea (mean difference -0.12 to -0.22mm, p = <0.0001-0.048, Figure 6) and BSSC measurements at the mid-peripheral and peripheral temporal cornea (mean difference -0.11 to -0.16mm at flagged locations, p = 0.0003-0.049). Additionally, scattered peripheral BSSC measurements at the nasal cornea showed significant correlations with spherical equivalent refractive error (r = 0.19 to 0.23, p = 0.010.046) and ethnicity (r = -0.18 to -0.21, p = 0.03-0.047). All other correlations did not return statistically significant results. Gender-Specific Regression Analysis Due to the widespread differences in FSSC between genders, sliding window analyses of the FSSC data were repeated with the cohort separated by gender. While FSSC followed a cubic pattern of change in the female cohort (R2 = 0.23-0.54, Figure 7), similar to the overall cohort, age-related change in the male cohort defaulted to a linear regression pattern, particularly in the more peripheral clusters, and resultant regression models provided generally poorer fits (R2 = 0.009-0.57, Figure 7 and Supplementary Tables 6 and 7).T-tests with Welch’s correction comparing coefficients of regression at corresponding clusters indicated significantly different rates of change with age between genders for all clusters (p = <0.0001-0.015), excluding 1 paracentral cluster in both hierarchical and k-means cluster maps (p = 0.515-0.617, Cluster 2 in Figure 2B). DISCUSSION Practical Implications Diagnosis of SCE is a persistent clinical problem, in spite of advances in corneal imaging technologies developed to detect ectasia with greater accuracy, with several global parameters showing only moderate sensitivity and specificity at 70.6% to 82.3% and 61.8% to 73.2% respectively.33, 34 The expected large inter-individual variability within normative cohorts is likely to contribute significantly to the evident overlap between normal and SCE measurements. By characterizing the influences of confounding factors resulting in such variability in greater detail, a narrower range of values considered within normal range in the 5

context of specific demographic information could theoretically be calculated, which has potential to be applied to methods facilitating more accurate detection of SCE. Using hierarchical and k-means cluster algorithms, we confirmed overall concentric patterns of aging for CT and FSSC measurements. Previous studies have analyzed CT and CC data pooled in concentric rings,11, 35, 36 based on the assumption that concentric corneal locations demonstrate similar properties. The cluster patterns derived from statistical clustering algorithms identify locations that are more appropriate to pool, thereby providing an improved framework to perform further analyses in the context of reduced measurement noise. In a similar nature, sliding window methods can successfully filter out pertinent information from noise. The higher order polynomial regression models obtained from sliding window analyses (Figure 4) were not replicated in decade-based models (Figure 3), and CT and FSSC sliding window models showed greater consistency across clusters and consequently less discrepancy with the age-normative cohort. This suggests that models utilizing single averaged measurements to represent data within semi-arbitrarily defined age-brackets, as is commonly performed in cohort studies,16, 17, 32, 37 may not describe corneal aging changes with sufficient precision. Interestingly, BSSC measurements did not exhibit concentric cluster patterns, but several mid-peripheral and peripheral points mirrored patterns of change identified in central locations, and models based on these cluster patterns showed poor predictive accuracy. Collectively, these locations did not show significant correlations with the investigated demographic factors (Figure 6). As investigation of other factors that could impact on BSSC, including location-specific variations in anterior chamber depth and corneal biomechanical properties,37 were outside the scope of the present study, future work exploring such factors would be worthwhile to aid development of more accurate models describing BSSC variations. Changes in Corneal Thickness with Astigmatism As the majority of the cohort with corneal astigmatism greater than 0.75 diopters demonstrated with-the-rule astigmatism, we postulate that the small but statistically significant correlations between higher magnitudes of astigmatism and thinner inferior paracentral CT were due to irregular tear film pooling secondary to relative steepening along the vertical meridian. While the tear film is purportedly excluded from Pentacam-derived CT measurements,38 the air-tear film interface forms the anterior reflecting surface captured by the Pentacam, and therefore inconsistencies in the tear film may interfere with accurate detection of the corneal epithelial surface.39, 40 Alternately, there may be increased variability in inferior CT measurements in individuals with moderate with-the-rule astigmatism; while Huang et al.41 did not report increased variability in inferior CT, given that their cohort was not stratified by magnitude of astigmatism it is difficult to ascertain its possible role on measurement variability. The similar location of minimum CT in SCE and established keratoconus42, 43 suggests that in the presence of moderate with-the-rule astigmatism, reduced paracentral CT should not be utilized as a diagnostic criterion for SCE in isolation. Rather, as proposed by Ambrósio et al.,4 rates of change in CT measurements from central and peripheral locations may be more pertinent in distinguishing SCE, and the identified cluster patterns have the potential to be applied in future analyses characterizing rates of CT change between normals and SCE eyes. Factors Potentially Contributing to Variations in Corneal Parameters The widespread and significant differences in CC measurements between genders were consistent with previous studies reporting significantly flatter keratometry readings in males,10, 44, 45 and significant differences between gender-specific FSSC regression models were also observed (Figure 7). While there were no significant correlations between CT and gender, the location of the two peak stationary points calculated from the CT quartic models 6

are consistent with those derived from FSSC cubic models. Overall, this calls into question the potential role of biological or biochemical differences, such as hormonal variations, in modulating CT and FSSC. The presence of estrogen, progesterone and androgen receptors in the human cornea certainly renders this possible,46-48 with suggested mechanisms including estrogen-induced fluid retention and estrogen- and progesterone-induced negative regulation of collagenolysis.49, 50 Studies have reported significant corneal steepening in the second and third trimesters of pregnancy,51, 52 however other studies did not note significant variations in CC at different stages of pregnancy,53-55 and while Kiely et al.56 described fluctuations in keratometric readings at different stages of the menstrual cycle, these results were subsequently not replicable.53, 55, 57, 58 Similarly, the existence of CT fluctuations at different stages of pregnancy and the menstrual cycle are controversial.52, 54, 55, 57, 59-61 Further work investigating the contribution of hormonal variations on corneal variations would provide additional insight to aid earlier detection of changes outside expected normal ranges. Anatomical variations including axial length, corneal diameter and height have been proposed as alternate explanations for gender-specific differences in CC, as in the absence of significant refractive error, reduced corneal power is required to maintain focus at the retina of a larger, longer eye.62-65 While these measurements were not obtained in the present study, high refractive error formed part of the exclusion criteria and no significant correlations between spherical equivalent refractive error and CC were observed, hence large variations in these parameters that could have significantly influenced FSSC findings would not be expected in the present cohort. Furthermore, although axial length may contribute to absolute differences in CC between genders62, 63 and corneal diameter and height have been reported to show significant correlations with CC,63, 65-67 these associations are unlikely to explain the gender-specific differences in FSSC regression, as studies have not suggest different rates of change between genders.57, 64, 65, 67-69 Similarly, while we did not observe significant differences in CT between genders, the quartic aging patterns observed in CT have not been observed in variations in axial length or corneal diameter with age.62, 64, 67, 68 It is therefore unlikely that these anatomical factors are solely responsible for the observed patterns of aging in CC and CT. Limitations The current study relied on self-reported medical history from intake forms, which poses the issue of participants not disclosing their full medical history. Future studies incorporating access to comprehensive medical records would be valuable to verify the findings of the current study. Moreover, although a range of ethnicities were included in the present study, due to the demographics of patients attending CFEH only 22% of participants reported ethnicities other than Caucasian or East Asian, with ≤5% of participants in each additional category. This has probably contributed to the insignificant correlations observed between ethnicity and the investigated corneal parameters, contrary to a previous study,20 and studies incorporating a more diverse participant base may reveal additional systematic correlations. Lastly, while a longitudinal design may be more appropriate to characterize age-related regression, following individual participants across the desired age range would be impractical to perform. CONCLUSION The present study investigated the utility of cluster algorithms and sliding window analyses to minimize noise and aid identification of age-related regression patterns in corneal parameters. Cluster analysis has demonstrated that CT and FSSC follow concentric patterns of aging, and the robust quartic and cubic regression models describing aging variations in CT and FSSC respectively suggest that these parameters are influenced by hormonal variations. In contrast, the non-concentric patterns of aging and variable regression models obtained for BSSC measurements indicate that mechanisms other than or in addition to those investigated contribute to BSSC variations. In conjunction with the systematic correlations observed with astigmatism and gender, these models provide a framework that 7

should be considered in establishment of normal corneal databases and development of future models describing normal variations in corneal parameters, which may be applied in future directions to improve detection of SCE.

8

ACKNOWLEDGEMENTS AND DISCLOSURES The authors thank Dr. Maria Markoulli for her advice and manuscript feedback. Some of the preliminary results in the current manuscript were presented to the American Academy of Optometry annual meeting, November 9 2018, San Antonio, and further results as described within the manuscript were presented at the Association for Research in Vision and Ophthalmology (ARVO) annual meeting, April-May 2019, Vancouver. JT, JP, MK and BZ are named as inventors on a patent application relating to the presented body of work. This work was supported by Guide Dogs NSW/ACT through their support to clinical service delivery at CFEH. The funding body had no role in the conceptualization or writing of the paper. Financial Disclosures J. Tong, P; J. Phu, P; M Kalloniatis, P; B. Zangerl, P.

9

REFERENCES 1. Vazirani J, Basu S. Keratoconus: current perspectives. Clin Ophthalmol 2013;7:201930. 2. Mas Tur V, MacGregor C, Jayaswal R, O'Brart D, Maycock N. A review of keratoconus: Diagnosis, pathophysiology, and genetics. Surv Ophthalmol 2017;62(6):77083. 3. Gomes JAP, Tan D, Rapuano CJ, et al. Global consensus on keratoconus and ectatic diseases. Cornea 2005;34(4):359-69. 4. Ambrósio R, Jr., Alonso RS, Luz A, Coca Velarde LG. Corneal-thickness spatial profile and corneal-volume distribution: tomographic indices to detect keratoconus. J Cataract Refract Surg 2006;32(11):1851-9. 5. Fernández Pérez J, Valero Marcos A, Martínez Peña FJ. Early diagnosis of keratoconus: what difference is it making? Br J Ophthalmol 2014;98(11):1465-6. 6. Buhren J, Kook D, Yoon G, Kohnen T. Detection of subclinical keratoconus by using corneal anterior and posterior surface aberrations and thickness spatial profiles. Invest Ophthalmol Vis Sci 2010;51(7):3424-32. 7. Shetty R, Rao H, Khamar P, et al. Keratoconus screening indices and their diagnostic ability to distinguish normal from ectatic corneas. Am J Ophthalmol 2017;181:1408. 8. Ambrósio R, Jr, Klyce SD, Wilson SE. Corneal topographic and pachymetric screening of keratorefractive patients. J Refract Surg 2003;19(1):24-9. 9. Flanagan GW, Binder PS. Precision of flap measurements for laser in situ keratomileusis in 4428 eyes. J Refract Surg 2003;19(2):113-23. 10. Orucoglu F, Akman M, Onal S. Analysis of age, refractive error and gender related changes of the cornea and the anterior segment of the eye with Scheimpflug imaging. Cont Lens Anterior Eye 2015;38(5):345-50. 11. Montalbán R, Piñero DP, Javaloy J, Alió JL. Scheimpflug photography-based clinical characterization of the correlation of the corneal shape between the anterior and posterior corneal surfaces in the normal human eye. J Cataract Refract Surg 2012;38(11):1925-33. 12. Goto T, Klyce SD, Zheng X, Maeda N, Kuroda T, Ide C. Gender- and age-related differences in corneal topography. Cornea 2001;20(3):270-6. 13. Cho P, Lam C. Factors affecting the central corneal thickness of Hong KongChinese. Curr Eye Res 2009;18(5):368-74. 14. Ueno Y, Hiraoka T, Miyazaki M, Ito M, Oshika T. Corneal thickness profile and posterior corneal astigmatism in normal corneas. Ophthalmology 2015;122(6):1072-8. 15. Sanchis-Gimeno JA, Lleó-Pérez A, Alonso L, Rahhal MS. Caucasian emmetropic aged subjects have reduced corneal thickness values. Emmetropia, CCT and age. Int Ophthalmol 2004;25(4):243-6. 16. Galgauskas S, Juodkaite G, Tutkuviene J. Age-related changes in central corneal thickness in normal eyes among the adult Lithuanian population. Dove Press 2014;9:114551. 17. Ueno Y, Hiraoka T, Beheregaray S, Miyazaki M, Ito M, Oshika T. Age-related changes in anterior, posterior, and total corneal astigmatism. J Refract Surg 2014;30(3):1927. 18. Prasad A, Fry K, Hersh PS. Relationship of age and refraction to central corneal thickness. Cornea 2011;30(5):553-5. 19. Rüfer F, Schröder A, Bader C, Erb C. Age-related changes in central and peripheral corneal thickness. Determination of normal values with the Orbscan II topography system. Cornea 2007;26(1):1-5. 20. Aghaian E, Choe JE, Lin S, Stamper RL. Central corneal thickness of Caucasians, Chinese, Hispanics, Filipinos, African Americans, and Japanese in a glaucoma clinic. Ophthalmology 2004;111(12):2211-9. 21. Yoshioka N, Zangerl B, Nivison-Smith L, et al. Pattern recognition analysis of agerelated retinal ganglion cell signatures in the human eye. Invest Ophthalmol Vis Sci 2017;58(7):3086-99. 10

22. Phu J, Khuu SK, Nivison-Smith L, et al. Pattern recognition analysis reveals unique contrast sensitivity isocontours using static perimetry thresholds across the visual field. Invest Ophthalmol Vis Sci 2017;58(11):4863-76. 23. Phu J, Khuu SK, Bui BV, Kalloniatis M. Application of pattern recognition analysis to optimize hemifield asymmetry patterns for early detection of glaucoma. Transl Vis Sci Technol 2018;7(5):3. 24. Bae HW, Rho S, Lee HS, et al. Hierarchical cluster analysis of progression patterns in open-angle glaucoma patients with medical treatment. Invest Ophthalmol Vis Sci 2014;55(5):3231-6. 25. Bae HW, Ji Y, Lee HS, et al. A hierarchical cluster analysis of normal-tension glaucoma using spectral-domain optical coherence tomography parameters. J Glaucoma 2015;24(4):328-33. 26. Asaoka R. Mapping glaucoma patients' 30-2 and 10-2 visual fields reveals clusters of test points damaged in the 10-2 grid that are not sampled in the sparse 30-2 grid. PLoS One 2014;9(6):e98525. 27. Ertan A, Muftuoglu O. Keratoconus clinical findings according to different age and gender groups. Cornea 2008;27(10):1109-13. 28. Romero-Jiménez M, Santodomingo-Rubido J, Wolffsohn JS. Keratoconus: a review. Cont Lens Anterior Eye 2010;33(4):157-66. 29. Belin MW, Ambrósio R, Jr. Scheimpflug imaging for keratoconus and ectatic disease. Indian J Ophthalmol 2013;61(8):401-6. 30. Belin MW, Khachikian SS. An introduction to understanding elevation-based topography: how elevation data are displayed - a review. Clin Exp Ophthalmol 2009;37(1):14-29. 31. Labiris G, Giarmoukakis A, Rampotas A, et al. Variability of different reference bodies in normal, keratoconus, and collagen crosslinked corneas. Eur J Ophthalmol 2015;25(6):468-73. 32. Hayashi K, Hayashi H, Hayashi F. Topographic analysis of the changes in corneal shape due to aging. Cornea 1995;14(5):527-32. 33. Hashemi H, Beiranvand A, Yekta A, Maleki A, Yazdani N, Khabazkhoob M. Pentacam top indices for diagnosing subclinical and definite keratoconus. J Curr Ophthalmol 2016;28(1):21-6. 34. Steinberg J, Casagrande MK, Frings A, et al. Screening for subclinical keratoconus using swept-source Fourier domain anterior segment optical coherence tomography. Cornea 2015;34(11):1413-9. 35. Hashemi H, Asgari S, Emamian MH, Mehravaran S, Fotouhi A. Five year changes in central and peripheral corneal thickness: The Shahroud Eye Cohort Study. Cont Lens Anterior Eye 2016;39(5):331-5. 36. Hashemi H, Asgari S, Emamian MH, Mehravaran S, Fotouhi A. Age-related changes in corneal curvature and shape: The Shahroud Eye Cohort Study. Cornea 2015;34(11):1456-8. 37. Sharifipour F, Panahi-Bazaz M, Bidar R, Idani A, Cheraghian B. Age-related variations in corneal biomechanical properties. J Curr Ophthalmol 2016;28(3):117-22. 38. Zhuang H, Zhou X, Xu J. A novel method for pachymetry mapping of human precorneal tear film using Pentacam with fluorescein. Invest Ophthalmol Vis Sci 2010;51(1):156-9. 39. Barkana Y, Gerber Y, Elbaz U, et al. Central corneal thickness measurement with the Pentacam Scheimpflug system, optical low-coherence reflectometry pachymeter, and ultrasound pachymetry. J Cataract Refract Surg 2005;31(9):1729-35. 40. Lee J-H, Kim JH, Kim SW. Repeatability of central corneal thickness measurement using rotating Scheimpflug camera in dry and normal eyes. Eye Contact Lens 2018;44(6):S29-S32. 41. Huang J, Ding X, Savini G, et al. A Comparison between Scheimpflug imaging and optical coherence tomography in measuring corneal thickness. Ophthalmology 2013;120(10):1951-8. 11

42. Romero-Jiménez M, Santodomingo-Rubido J, González-Méijome J-M. The thinnest, steepest, and maximum elevation corneal locations in noncontact and contact lens wearers in keratoconus. Cornea 2013;32(3):332-7. 43. Temstet C, Sandali O, Bouheraoua N, et al. Corneal epithelial thickness mapping using Fourier-domain optical coherence tomography for detection of form fruste keratoconus. J Cataract Refract Surg 2015;41(4):812-20. 44. Scholz K, Messner A, Eppig T, Bruenner H, Langenbucher A. Topography-based assessment of anterior corneal curvature and asphericity as a function of age, sex, and refractive status. J Cataract Refract Surg 2009;35(6):1046-54. 45. Eysteinsson T, Jonasson F, Sasaki H, et al. Central corneal thickness, radius of the corneal curvature and intraocular pressure in normal subjects using non-contact techniques: Reykjavik Eye Study. Acta Ophthalmol Scand 2002;80(1):11-5. 46. Gupta PD, Johar K, Sr., Nagpal K, Vasavada AR. Sex hormone receptors in the human eye. Surv Ophthalmol 2005;50(3):274-84. 47. Wickham LA, Gao J, Toda I, Rocha EM, Ono M, Sullivan DA. Identification of androgen, estrogen and progesterone receptor mRNAs in the eye. Acta Ophthalmol Scand 2000;78(2):146-53. 48. Suzuki T, Kinoshita Y, Tachibana M, et al. Expression of sex steroid hormone receptors in human cornea. Curr Eye Res 2009;22(1):28-33. 49. Sato T, Ito A, Mori Y, Yamashita K, Hayakawa T, Nagase H. Hormonal regulation of collagenolysis in uterine cervical fibroblasts. Biochem J 1991;275(3):645-50. 50. Zhou H, Kimura K, Orita T, Nishida T, Sonoda K-H. Inhibition by female sex hormones of collagen degradation by corneal fibroblasts. Mol Vis 2011;17:3415-22. 51. Park SB, Lindahl KJ, Temnycky GO, Aquavella JV. The effect of pregnancy on corneal curvature. CLAO J 1982;18(4):256-9. 52. Ataş M, Duru N, Ulusoy DM, et al. Evaluation of anterior segment parameters during and after pregnancy. Cont Lens Anterior Eye 2014;37(6):447-50. 53. Bahadir Kilavuzoglu AE, Cosar CB, Bildirici I, Cetin O, Ozbasli E. Estrogen- and progesterone-induced variation in corneal parameters according to hormonal status. Eye Contact Lens 2018;44(5):S179-84. 54. Efe YK, Ugurbas SC, Alpay A, Ugurbas SH. The course of corneal and intraocular pressure changes during pregnancy. Can J Ophthalmol 2012;47(2):150-4. 55. Manges TD, Banaitis DA, Roth N, Yolton RL. Changes in optometric findings during pregnancy. Am J Optom Physiol Opt 1987;64(3):159-66. 56. Kiely PM, Carney LG, Smith G. Menstrual cycle variations of corneal topography and thickness. Am J Optom Physiol Opt 1983;60(10):822-9. 57. Hashemi H, Mehravaran S, Rezvan F. Changes in corneal thickness, curvature, and anterior chamber depth during the menstrual cycle. Can J Ophthalmol 2010;45(1):67-70. 58. Oliver KM, Walsh G, Tomlinson A, McFadyen A, Hemenger RP. Effect of the menstrual cycle on corneal curvature. Ophthal Physiol Opt 1996;16(6):467-73. 59. Giuffrè G, Di Rosa L, Fiorino F, Bubella DM, Lodato G. Variations in central corneal thickness during the menstrual cycle in women. Cornea 2007;26(2):144-6. 60. Goldich Y, Cooper M, Barkana Y, et al. Ocular anterior segment changes in pregnancy. J Cataract Refract Surg 2014;40(11):1868-71. 61. Ziai N, Ory SJ, Khan AR, Brubaker RF. β-human chorionic gonadotropin, progesterone, and aqueous dynamics during pregnancy. Arch Ophthalmol 1994;112(6):8016. 62. Hashemi H, Khabazkhoob M, Miraftab M, et al. The distribution of axial length, anterior chamber depth, lens thickness, and vitreous chamber depth in an adult population of Shahroud, Iran. BMC Ophthalmol 2012;12:1-8. 63. Hoffmann PC, Hütz WW. Analysis of biometry and prevalence data for corneal astigmatism in 23,239 eyes. J Cataract Refract Surg 2010;36(9):1479-85. 64. Kim J-H, Kim M, Lee S-J, et al. Age-related differences in ocular biometry in adult Korean population. BMC Ophthalmol 2016;16(1):146.

12

65. Wong TY, Foster PJ, Johnson GJ, Klein BEK, Seah SKL. The relationship between ocular dimensions and refraction with adult stature: The Tanjong Pagar Survey. Ophthalmol Vis Sci 2001;42(6):1237-42. 66. Jonuscheit S, Doughty MJ, Martín R, Rio-Cristobal A. Relationship between corneal thickness and radius to body height. Optom Vis Sci 2016;94(3):380-6. 67. Hashemi H, Khabazkhoob M, Emamian MH, Shariati M, Yekta A, Fotouhi A. Whiteto-white corneal diameter distribution in an adult population. J Curr Ophthalmol 2015;27(12):21-4. 68. Grosvenor T. Reduction in axial length with age: an emmetropizing mechanism for the adult eye? Am J Optom Physiol Opt 1987;64(9):657-63. 69. Ninn-Penderson K. Relationships between preoperative astigmatism and corneal optical power, axial length, intraocular pressure, gender, and patient age. J Refract Surg 1996;12(4):472-82.

13

FIGURE CAPTIONS Figure 1. The 57 locations from which CT, FSSC and BSSC measurements were obtained from Pentacam HR scans. Figure 2. Cluster patterns derived from hierarchical and k-means clustering for CT, FSSC and BSSC measurements. Each cluster, denoted by different colors, indicates locations demonstrating statistically similar change with age. Locations outlined in black indicate those that showed variable cluster assignment between hierarchical and k-means patterns. ‘N’ denotes nasal locations. Figure 3. Preferred regression models describing age-related change in CT, FSSC and BSSC, with data pooled as per the derived cluster patterns and per decade bracket, with each color representing a different cluster. Error bars indicate 1 standard deviation per age group and cluster, and stars indicate clusters where a quadratic regression model was preferred over linear regression using sum-of-squares F test at p<0.05. As models were similar between those derived from hierarchical and k-means clustering for each parameter, only those derived from hierarchical clustering are shown for clarity. ‘N’ denotes nasal locations. Figure 4. Polynomial regression models describing age-related change in CT, FSSC and BSSC, obtained using sliding window analyses and in hierarchical cluster patterns. Only models derived from hierarchical clustering with data pooled as per clusters 1 and 7, denoted by orange and green locations respectively, are shown for clarity. ‘N’ denotes nasal locations. Figure 5. Bland-Altman plots showing differences in CT, FSSC and BSSC between a 35 year old normative cohort and the remainder of the study population converted to a 35 year old equivalent (Table 1), with age-correction performed using the polynomial regression models obtained using sliding window analyses (Column 1) and decade-based models (Column 2). To aid comparisons between parameters, for CT measurements the difference between the age-normative and age-equivalent cohorts is expressed in millimeters. The dashed lines indicate the overall bias of the age-corrected measurements, and the dotted lines show the 95% limits of agreement. For clarity, only Bland-Altman plots for models derived from hierarchical clustering are shown. ‘N’ denotes nasal locations. Figure 6. Locations showing statistically significant Pearson’s correlations between CT measurements and refractive astigmatism (blue) and corneal astigmatism (orange), and BSSC measurements spherical equivalent refractive error (purple) and ethnicity (green). Locations showing statistically significant differences with Welch’s t-tests in FSSC and BSSC between genders are also shown (fuchsia). Circles with more than 1 color indicate locations where multiple significant correlations were observed. ‘N’ denotes nasal locations. Figure 7. Cubic models obtained using sliding window analyses for front surface sagittal curvature measurements, with data pooled as per the derived cluster maps and the study cohort separated by gender. For clarity, only models derived from hierarchical clustering with data pooled as per clusters 1 and 7, denoted by orange and green locations respectively, are shown for clarity. ‘N’ denotes nasal locations.

14

Table 1. Demographic characteristics of the study population. The 35-year-old normative cohort utilized in Bland-Altman comparisons was the 31-40 age bracket, and the remaining participants were included in the age-corrected cohort. Cohort

n

Age ± SD, y

All participants ≤10 11-20

117 32.90 ± 16.00 9 9.21 ± 0.67 25 16.32 ± 2.69

21-30 31-40 41-50

24 19 19

51-60

21

Gender, M:F

SE ± SD, Diopters

57:60

Eye tested, OD:OS 65:52

5:4 12:13

4:5 12:13

26.03 ± 3.32 34.54 ± 3.29 46.27 ± 2.45

16:8 11:8 7:12

14:10 10:9 11:8

57.05 ± 3.17

6:15

14:7

-2.15 ± 1.24 -0.99 ± 1.52 -1.60 ± 1.98 -0.55 ± 0.82 +0.02 ± 1.31 +0.58 ± 1.51 -0.72 ± 1.81

-0.69 ± 1.69

Astigmatism ± SD, Diopters -0.66 ± 0.64 -0.58 ± 0.71 -0.53 ± 0.61 -0.97 ± 0.75 -0.67 ± 0.66 -0.57 ± 0.56 -0.58 ± 0.54

Age-corrected 98 32.58 ± 46:52 55:43 -0.66 ± 0.65 cohort 17.42 n, number of participants; y, years; SD, standard deviation; M, male; F, female; OD, right eye; OS, left eye; SE, spherical equivalent

Table 2. Coefficients of determination (R2) for the polynomial regression models derived from sliding window analyses for the investigated corneal parameters and stationary points calculated from each regression model. N/A indicates clusters where regression models for that particular cluster did not have an additional stationary point or no stationary points were derived Hierarchical K-means Point Point Point Point Point Point 2 2 R R 1 2 3 1 2 3 Corneal thickness Cluster 1 0.71 17.2 31.4 53.3 0.53 17.2 31.9 52.9 Cluster 2 0.70 17.4 31.8 52.9 0.71 17.3 31.5 53.2 Cluster 3 0.75 17.6 32.3 51.5 0.71 17.3 32 52.3 Cluster 4 0.70 17.4 32.6 51.4 0.75 17.6 32.3 51.5 Cluster 5 0.71 16.4 33.1 50.8 0.69 17.3 32.2 52.1 Cluster 6 0.69 17.1 33.0 51.9 0.69 16.8 33.4 50.8 Cluster 7 0.74 17.4 34.3 49.6 0.74 17.4 34.3 49.6 Front surface sagittal curvature Cluster 1 0.43 23.6 47.3 0.43 23.6 47.4 Cluster 2 0.49 24.7 44.8 0.49 24.6 44.9 Cluster 3 0.67 32.2 50.8 0.46 23.1 51.1 Cluster 4 0.57 22.6 51.9 0.62 22.7 52.2 Cluster 5 0.63 21.1 53.9 0.60 21.8 53.6 Cluster 6 0.70 23.1 53.8 0.74 22.8 54.7 Cluster 7 0.42 25.5 47.4 0.45 25.2 48.0 Cluster 8 0.62 24.7 57.4 0.62 24.7 57.4 Back surface sagittal curvature Cluster 1 0.67 32.0 36.0 0.63 30.3 34.4 Cluster 2 0.64 39.2 N/A 0.52 9.5 40.6 Cluster 3 0.30 18.3 41.3 0.30 17.9 42.2 Cluster 4 0.67 29.6 49.7 0.67 29.6 49.7 Cluster 5 0.47 17.4 43.1 0.66 15.3 44.2 Cluster 6 0.72 45.7 N/A 0.64 52.5 N/A Cluster 7 0.69 N/A N/A 0.72 N/A N/A 2 R , coefficient of determination

Table 3. Parameters of Bland-Altman comparisons between corneal measurements for a 35 year old normative cohort and the cohort age-corrected to 35 year equivalent (Table 1) using the polynomial regression models derived from sliding window and decade-based models. The width of the 95% limits of agreement interval was expressed in both absolute units (mm) and percentage of the range of corneal thickness, front surface sagittal curvature and back surface sagittal curvature measurements. Hierarchical Bias ± SD, mm 95% LoA interval width, mm (%)

K-Means Bias ± SD, mm 95% LoA interval width, mm (%)

Sliding window models Corneal thickness -0.002 ± 0.004 0.015 (3.27) -0.003 ± 0.004 Front surface sagittal -0.04 ± 0.04 0.14 (4.59) -0.04 ± 0.04 curvature Back surface sagittal -0.11 ± 0.19 0.73 (17.85) -0.12 ± 0.18 curvature Decade-based models Corneal thickness 0.005 ± 0.004 0.017 (3.79) 0.005 ± 0.004 Front surface sagittal -0.09 ± 0.07 0.27 (9.22) -0.04 ± 0.09 curvature Back surface sagittal -0.03 ± 0.14 0.54 (11.76) -0.03 ± 0.14 curvature SD, standard deviation; LoA, limits of agreement; mm, millimeters

0.018 (3.88) 0.14 (4.59) 0.70 (17.11)

0.017 (3.79) 0.34 (11.60) 0.54 (11.76)