Age-Related Changes in Corneal Deformation Dynamics Utilizing Scheimpflug Imaging

Purpose To study age-related changes in corneal deformation response to air-puff applanation tonometry. Methods Fifty healthy subjects were recruited for a prospective study and divided into two equal age groups (≤ 28 and ≥ 50 years old). Up to three measurements by a corneal deformation analyser based on the Scheimpflug principle were performed on the left eye of each subject. Raw Scheimpflug images were used to extract changes in anterior and posterior corneal profiles, which were further modelled by an orthogonal series of Chebyshev polynomial functions. Time series of the polynomial coefficients of even order exhibited a dynamic behavior in which three distinct stages were recognized. A bilinear function was used to model the first and the third stage of corneal dynamics. Slope parameters of the bilinear fit were then tested between the two age groups using Wilcoxon rank sum test and two-way non-parametric ANOVA (Friedman) test. Results Statistically significant changes (Wilcoxon test, P<0.05) between the age groups were observed in the phase of the second applanation dynamics for the posterior corneal profile. In a two-way comparison, in which the corneal profile was used as a dependent variable, statistically significant changes (ANOVA/Friedman test, P = 0.017) between the groups were also observed for that phase. Conclusion Corneal biomechanics depend on age. The changes in corneal deformation dynamics, which correspond to mostly free return of the cornea to its original shape after the air pulse, indicate that the age related differences in corneal biomechanics are subtle but observable with high speed imaging.

Introduction systemic disease and contact lens wear. For each subject, up to three measurements of corneal deformation response (in terms of high speed Scheimpflug camera recordings) were carried out on left eyes only. Each acquired series of images was saved for further analysis. To assess uniformity of subjects in age groups, central corneal thickness (CCT) and the intraocular pressure (IOP) registered by Scheimpflug analyser was noted.
All measurements were performed in succession, allowing about one minute break between each acquisition. To take into account diurnal variations in subjects IOP [26] all measurements were conducted in the mornings between 10 am and 12 pm.
The study adhered to the tenets of the Declaration of Helsinki and was approved by the Ethics Committee of the Medical University of Wroclaw (KB 503/2011 agreement). Written informed consent was obtained from all participants included in this study. Patient records were anonymized and de-identified prior to analysis.

Data analyses and statistical methods
The raw Scheimpflug images were numerically processed using a custom-written program in Matlab (Math Works, Inc., Natick, MA, USA). Segmentation algorithms including Otsu's thresholding method [27] for edge detection of the anterior and posterior corneal contours were applied for every frame. Because of increasing noise at the cornea periphery, two values of data trimming were considered. This corresponded to cutting out either 5% or 10% of the extracted profile data on both sides. Final analysis was performed for a 10% cut that showed to be more robust in terms of the subsequent parametric modelling of the corneal edge data. In the preliminary analysis on a smaller group of subjects [28], standard polynomial approximation with the optimally set order of six (in the Akaike Information Criterion sense [29]) was performed for detected edges. The results of this analysis were not entirely conclusive when the extended set of data was used. Subsequently, this was followed with Chebyshev Changes in corneal shape for a few sample frames. The first frame refers to the normal cornea state, the next frames are associated with the first applanation state (n = 31), the concavity state (n = 74), and the second applanation state (n = 95) while the last frame corresponds to return of the cornea to its initial state.
doi:10.1371/journal.pone.0140093.g001 polynomial approach as they form a set of orthogonal and complete series of basis functions in which estimated coefficient is independent of the other in the series [30]. More importantly, the boundary behavior of Chebyshev polynomials at the edges of the extracted edges is more stable. Time varying polynomial Chebyshev coefficients a 0 (t),a 1 (t),. . .,a n (t), were estimated for a given time instant t using a least squares procedure from a model: p c ðx; tÞ ¼ a n ðtÞT n ðxÞ þ a nÀ1 ðtÞT nÀ1 ðxÞ þ . . . þ a 1 ðtÞT 1 ðxÞ þ a 0 ðtÞT 0 ðxÞ þ εðtÞ; Where T k (x), k = 0,1,. . .,n, calculated recursively, is the k-th Chebyshev polynomial with x representing the horizontal image axis in pixels, and ε (t) denotes the measurement and modeling error. Similarly as in the case of standard polynomials, Akaike Information Criterion was used to determine the optimal model order, whose median for all considered measurements was equal to six. Fig 3 shows an example of the timevarying Chebyshev polynomial coefficients, a 6 (t), describing changes in anterior and posterior corneal surfaces for two selected subjects from the young and older group. Similar time variation of corneal deformation was observed for all even order Chebyshev polynomial coefficients, i.e., for a 2 (t) and a 4 (t).
For each modelled measurement consisting of about 32 ms time series of the sixth order Chebyshev polynomial coefficient a three-stage process was observed. The first stage comprises the time from the first recorded image to the maximum concavity state (from about 7 to 11 ms from the start of the recording). The second phase describes the corneal oscillation period [24] (from about 11 to 16 ms) while the third stage corresponds to outgoing concavity and return of the cornea to its initial state (corresponds to the rest of the 32 ms recording). In this study, the first and the third stage of that process is considered.
To separate the slow and fast changes in corneal deformation, each of the considered two stages has been quantitatively described by a constrained bilinear fit. This results in further division of each of the stages into two separate phases. Unlike in a piece-wise linear regression, in the constrained bilinear fit the estimate of the second line is conditioned on the estimate of the first line resulting in a model without discontinuities [31]. An example of the constrained bilinear fitting is shown in Fig 4. The performance of bilinear estimation procedure has been investigated using nonparametric bootstrap in which residuals between the considered phase time series and the bi-linear model were resampled. The bootstrap estimated standard deviation of the intercept was small and amounted to less than a sample point which is equivalent to less than 0.24 ms. This indicates that the proposed procedure of bilinear fitting is justified.
As a result, we obtain eight slope coefficients m A1 , m A2 , m A3 , m A4 , m P1 , m P2 , m P3 , m P4 , where subscripts A and P denote the anterior and posterior surface, respectively and the numeral indicates the line order. The first line (m A1 or m P1 ) describes the initial dynamic state of the cornea surface (anterior or posterior, respectively) before the first applanation. The second and third lines (m A2 or m P2 and m A3 or m P3 ) indicate the ingoing and outgoing corneal dynamics, respectively. The last line (m A4 or m P4 ) describes the dynamics of the return to the original state of cornea. These parameters were considered for further statistical analysis. The same procedure was applied to both groups of subjects.
Since normality of the data was rejected (Jarque-Bera test, P<0.05), all slope coefficients were tested in relation to subject age with Wilcoxon rank sum test for equal medians. Additionally, of interest were relationships between the four sets of parameters (i.e., m A1 vs. m P1 , m A2 vs. m P2 , m A3 vs. m P3 , and m A4 vs. m P4 ). These data was tested with Jarque-Bera normality test and the Bartlett's test for equal variances. Since in majority of cases normality and equal variance criteria were not fulfilled, non-parametric two-way ANOVA (Friedman) test was applied to examine the data for age dependence. The significance level for all tests was set at α = 0.05. All statistical analyses were performed in MATLAB.

Results
Baseline characterization of both groups is shown in Table 1 (age, gender, CCT, and IOP). The group average IOP and CCT were (mean ± S.D.) 15.2 ± 1.9 mmHg and 559.4 ± 35.1 μm, respectively for the young group, and 14.8 ± 2.3 mmHg and 554.9 ± 36.8 μm, respectively for the older group. No statistically significant differences between age groups in CCT and IOP were observed (p = 0.39 and p = 0.66, respectively).
Considering age-related changes in the slope coefficients (see Table 2), statistically significant differences were observed only for the m P3 parameter (Wilcoxon rank sum test, p = 0.009, p = 0.007, and p = 0.012 for a 6 (t), a 4 (t), a 2 (t), respectively). Table 3 presents a summary of statistical test results for two-way comparison, in which the corneal profile was used as a dependent variable. Statistically significant differences (two-way ANOVA (Friedman)) were observed for m A3 vs. m P3 (p = 0.016). Table 3 includes also the results of differential analysis where statistically significant changes were found for dm 3 = m A3 −m P3 (p = 0.017 and p = 0.005 for a 6 (t) and a 4 (t), respectively).

Discussion
High-speed Scheimpflug images of the left eyes of 50 healthy volunteers in two age groups were used in the study. All subjects exhibited normal IOP values [32] and normal central corneal thickness values [33], which were not statistically different between the groups. The dynamics of the corneal deformation response was assessed in relation to subject's age. We proposed the time-varying Chebyshev polynomial based model, in which higher even order terms exhibited a distinct three-stage dynamic behavior. Those even terms mostly capture corneal asphericity changes undergoing during cornea applanation. Likewise, the odd terms of the polynomial fit mostly capture the rotational movement of the eye globe. Focus was made on the first and the third stage of corneal deformation dynamics which include the times of the first and second corneal applanation. Those two stages were further optimally sequenced (by fitting a constrained bilinear model) into phases that were linearly modelled. The first phase represents the pre-applanation cornea state in which the applied air pressure is sufficiently small to cause only linear deformation of cornea. The second phase includes the first applanation and the peak air pressure value. The third phase corresponds to the second applanation of  Table 3. Results (p-values) of two-way ANOVA (Friedman) for testing age-related differences in slope parameters of the bilinear models. Subscripts A and P correspond to the anterior and posterior corneal profiles, respectively. the cornea while the fourth phase is related to the return of corneal surface to its original state where again linear approximation to cornea deformation can be used. Considering each surface separately, parameters of bilinear model showed age-related differences, particularly in the third phase of the posterior corneal profile where those differences were statistically significant. This confirms the preliminary results reported in [28]. On the other hand, when the anterior and posterior surfaces were considered as dependent factors in the analysis, statistically significant age-related changes were evident in the second and third phases for the 5% trimming of the profile data but not for one corresponding to the 10% trim. This suggests that the 5% trim may not be sufficient or that changes in corneal asphericity may play some role in this result. On the other hand, differential analysis, which carries more sensitive information, indicates statistically significant differences in the third phase for the 10% trim but not for the 5% trim thus confirming the results of previous analysis.
It is worth noting that in the first and the fourth phase of the corneal deformation dynamics the applied outer pressure is relatively low to cause any substantial deformation of the corneal profile and to show any age-related differences. In the second phase the outer pressure applied to the cornea is substantially higher than that of the IOP and it includes the first applanation.
Our results indicate that the corneal deformation dynamics in that phase do not change with age. We conclude that the pressure applied to the cornea is sufficiently high for the corneal biomechanical properties to play a substantial role in that dynamics. Finally, age-related changes in biomechanical parameters of the cornea are evident in the third phase of corneal deformation dynamics where the outer pressure earlier applied to the cornea is being released.
It is now well established that the material properties of cornea depend on age [8][9][10][11][12][13][14]. Elsheikh et al. [4] showed statistically significant differences in corneal stiffness in relation to age. Tonnu et al. [14] demonstrated that subject age has a differential effect on the IOP measurements made by the Goldmann applanation tonometer (GAT) and ocular blood flow tonograph (OBF) compared to the handheld tonometer. Klein et al. [34] suggested that age could be associated with an overestimation of IOP value. Kotecha et al. [8] found the effect of age on IOP measurement suggesting an age-related corneal biomechanical change that may induce measurement error additional to that of CCT. The observed in our study age-related changes in corneal deformation confirm those earlier ex-and in-vivo studies and provide a new insight into particular phases of corneal deformation dynamics.
Various age-related changes in the biomechanics of ocular components have been reported in the literature [35][36][37][38][39]. Full knowledge of the aging processes in the cornea could bring an important insight into diagnosis and treatment of eye diseases. Numerical analyses of the corneal deformation dynamics incorporating subject's age can be used to build a biomechanical model of the cornea that subsequently could help those endeavours.
Supporting Information S1 Raw Data. An Excel file containing the bi-linear fitting slope parameters values for all even order time-varying Chebyshev polynomial coefficients analysis (a 6 (t), a 4 (t) and a 2 (t)). (XLSX)