Assessment of corneal properties based on statistical modeling of OCT speckle

: A new approach to assess the properties of the corneal micro-structure in vivo based on the statistical modeling of speckle obtained from Optical Coherence Tomography (OCT) is presented. A number of statistical models were proposed to ﬁt the corneal speckle data obtained from OCT raw image. Short-term changes in corneal properties were studied by inducing corneal swelling whereas age-related changes were observed analyzing data of sixty-ﬁve subjects aged between twenty-four and seventy-three years. Generalized Gamma distribution has shown to be the best model, in terms of the Akaike’s Information Criterion, to ﬁt the OCT corneal speckle. Its parameters have shown statistically signiﬁcant diﬀerences (Kruskal-Wallis, p < 0.001) for short and age-related corneal changes. In addition, it was observed that age-related changes inﬂuence the corneal biomechanical behaviour when corneal swelling is induced. This study shows that Generalized Gamma distribution can be utilized to model corneal speckle in OCT in vivo providing complementary quantiﬁed information where micro-structure of corneal tissue is of essence.


Introduction
The cornea is the major part of the eye responsible for about two thirds of its total refractive power [1]. Corneal refraction is dependent on its properties and a change in its shape or structure contributes to a variation in optical properties of the eye. From the optical point of view, some of these changes may be corrected by the crystalline lens in early ages but become gradually significant when the crystalline lens starts to lose its power of accommodation [2]. Alterations of the corneal shape and structure may affect not only its optical properties but also contribute to changes in its biomechanical properties [3]. Changes to corneal tissue may have pathological origin as in the case of keratoconus [4] or be a result of corneal surgery [5,6], which is nowadays commonly used to compensate refractive errors. Ageing is another natural consequence which contributes to corneal changes during lifetime. The cornea is found to increase its stiffness as the characteristics of viscoelastic behaviour decrease with age [7,8]. Besides these irreversible structural changes, corneal tissue also experiences daily alterations when subjected to hypoxia, i.e., reduced amount of oxygen provided to corneal epithelium. Hypoxia induces corneal swelling and consequently, corneal edema, which may be noticed not only after few hours of sleep but also, in some cases, due to long contact lens wear [9,10]. Corneal edema is usually associated with an increase in central corneal thickness (CCT) which reflects the macro-structural point of view. However, from the micro-structural perspective, it is hard to point a reliable parameter describing how the biomechanical corneal properties are influenced.
Understanding the dynamics of corneal biomechanics is of high interest in vision science to evaluate the results of refractive surgery or the progression of keratoconus and is a factor to consider in glaucoma management. Glaucoma is a concern for irreversible blindness, as the world's population ages [11,12] while corneal biomechanics may significantly interfere with the measurement of Intraocular Pressure (IOP) leading to inaccurate tonometric readings. Therefore, further studies are needed to characterize corneal properties in order to gain better understanding of its mechanical behaviour in short and long periods of time.
Optical Coherence Tomography (OCT) has become a popular technology to visualize and analyze the anterior segment of the eye. It is particularly attractive for its capability to provide, in real time, cross sectional images of tissue structure in vivo [13]. Interferometry, the measurement technique on which OCT is based, gives rise to speckle, which has been considered a form of noise that degrades the quality of the OCT image. Speckle noise diminishes the image contrast which makes the boundaries between structures difficult to resolve. Therefore, in order to have smoother and better-outlined borders, speckle is usually suppressed for different OCT applications such as, intravascular imaging [14], skin imaging [15] or retinal layer segmentation [16][17][18]. However, as Schmitt et al. showed [19], OCT speckle can also be signal-carrying and be used to infer about the micro-structure of the tissue. In fact, the potential of the speckle has already been consolidated in ultrasound for the purpose of characterizing and identifying abnormalities in different kind of tissues such as crystalline lens [20], myocardium [21] or in the blood [22]. Generally speaking, the statistics of speckle signals are investigated by fitting probability distributions to the signal envelope intensities. The spatial distribution of signal-carrying and signal-degrading of the speckle have demonstrated to be sensitive to morphological changes in tissue.
With regards to OCT there are just few studies in the literature. Farhat et al. explored OCT backscatter signals to differentiate between modes of cell death in vitro [23]. Grzywacz et al. [24] proposed an automatic statistical procedure based on OCT speckle for a retina with diabetic retinopathy. More recently, Amini et al.
[25] suggested a normal-Laplace distribution for statistical modeling of retinal OCT images allowing precise identification of each layer based on the probability density function (pdf) of each OCT intra-retinal layer. As far as we know there are no attempts to statistically model the corneal speckle obtained from OCT, besides our early works [26]. Therefore, in this study, different models to predict the corneal quality are presented and discussed. Our aim is to demonstrate that age-related and short-term changes in corneal micro-structural properties can be assessed based on the analysis of OCT speckle.
The remainder of the paper is organized as follows. In Section II the proposed models, image acquisition and age-related and daily corneal changes are described. The best model to fit the corneal speckle and its applicability is presented in Section III. In Section IV, the obtained results are discussed and lastly, conclusions are given in Section V.

Methods
In order to statistically model the corneal OCT speckle and to analyse long-term (years) and short-term (hours) changes in corneal properties, two separated studies were performed. The first study, focused on corneal long-term changes (age-related) whereas the second study, explored short-term (daily) changes in corneal properties which may occur due corneal swelling.

Age-related changes study
For the purpose of studying long-term changes of corneal properties related to age, OCT data from sixty-five subjects with healthy corneas, aged between 24 and 73 years were divided into three age groups containing approximately the same number of subjects. The first group corresponds to 25 images of subjects up to 25 (24.4 ± 0.5, range between 24 and 25) years old, second group includes 20 images of 20 subjects between 26 and 50 (31.3 ± 4.6, range between 27 and 45) years old and finally, the third group corresponds to 20 images of 20 subjects older than 52 (61.2 ± 8.4, range between 52 and 73) years old. The data of each subject corresponds to one single eye chosen randomly. Unlike the daily changes of corneal tissue, age-related changes are not expected to change significantly during a period of a day or even months. Therefore one single image per subject was considered sufficient for the purpose of studying age-related changes.

Short-term changes study
Four subjects (two male and two female volunteers of 26, 53, 29 and 45 years old) were recruited for this part of study. Baseline measurements of both eyes were obtained by a OCT and a biometer (IOLMaster 700, Carl Zeiss Meditec AG, Germany) before inserting a contact lens in the right eye. Contact lens fitting was checked to ensure adequate corneal coverage and movement. The right eye was patched for three hours with an ophthalmic eye patch and a smooth tissue between the patch and the eye lid was added to ensure complete eye closure. Normal blinking was allowed for the control (left) eye. Measurements by OCT and IOLMaster were taken immediately after the removal of the eye patch and the contact lens, after 60 and 120 minutes later. Firstly, the measurements were taken in the right eye and then in the left eye and assessment by OCT was followed by IOLMaster at all time points. For each time point, 10 measurements by OCT and 6 measurements by IOLMaster were taken. Similarly to other studies where corneal edema was induced [27, 28], Hydroxyethyl-methacrylate soft contact lenses were used. The lenses had a water content of 60%, a diameter of 14.2 mm with an average oxygen permeability of 27 × 10 −11 (cm 2 /sec)(mlO 2 /ml · mmHg) at 35 • C and the base curve was 8.7 mm.

Image acquisition
The SOCT Copernicus HR (Optopol, Poland) was used to acquire the image of the anterior eye, specifically the corneal layers and part of the anterior chamber. The signal source was a super luminescent diode with a center wavelength of 850 nanometres. The axial and lateral resolutions of the system were 3 and 15 micrometres in tissue, respectively. The images were sampled through 1024 A-scans (columns) with 848 pixels per A-scan. With the help of the graphic user interface of SOCT Copernicus HR, all corneal images were carefully acquired approximately at the same axial position within the bands of instrument's depth of focus. Hence, taking into account the fact that the imaged cornea has a width of about a quarter of the depth of focus listed in instrument's technical specifications, we assume that images were acquired in the focal plane.
A ROI (250 rows by 450 columns) was selected to study long-term changes, as shown in Fig. 1. Such criteria was adopted because age-related changes affect the entire corneal stroma [29]. For the purpose of studying short-term changes, a region of interest (ROI) consisting of 150 rows by 450 columns was selected. This ROI was limited to middle and posterior corneal stroma because the anterior part is highly resistant to hydration even in extreme cases [30] and consequently is less susceptible to corneal swelling. Regarding the lateral dimensions, ROI was empirically limited to half of the available field of view because the OCT signal from the corneal stroma depends on the incident angle of the probe beam. Using a limited area, the effects of corneal curvature and its inter-subject variability are reduced.
An automatic routine was used to ensure that all ROIs remained at the same position. First, the original image was cropped resulting with the image presented in Fig. 1. Then, corneal edges were sharpened using the unsharp masking method of Matlab (Mathworks, Natick, MA). A threshold array based on the mean value of each A-scan was created. The threshold was used to detect the epithelium (upper layer) and the endothelium (lower layer). Both were fitted using a polynomial of third order in a least-squares sense. The Bowman's layer was detected finding the highest peak between the epithelium and endothelium and consequently approximated by a polynomial of third order. Finally, a ROI between the Bowman's capsule and the endothelium was selected accordingly the established dimensions, considering the null differentiation at epithelium curvature as a reference to the central horizontal position. For the purpose of maximizing the amount of signal-carrying data, it was important to ensure that no direct central light reflection was present in the images. However, it is also necessary to ensure that all images are acquired approximately at the same corneal position. Therefore, in the first instance, corneal apex was detected finding simultaneously the specular reflection in the 0 and 90 degree scans. Then, the instrument was moved with respect to the eye one single step of a motor to the inferior/nasal direction and a horizontal scan was obtained. That single step was estimated at 150 micrometres using test images. In addition, IOLMaster 700 was used to measure the corneal thickness at the corneal apex. To date, IOLMaster 500 has been the gold standard in optical biometry. However, the agreement between the gold standard and the IOLMaster 700 has shown to be excellent [31]. The B-scan images of OCT, in principle, could also be used to measure the corneal thickness in a similar way to that of IOLMaster [32]. However, the IOLMaster 700 is equipped with an integrated swept source OCT technology that helps ensuring repeatable axial measurements of the eye including that of CCT.

Speckle modeling
The statistical literature contains a large number of probabilistic models that could be applied to fitting the corneal OCT speckle. So far, our study has focused on the most common models which have been used to fit light or ultrasound speckle such as, Generalized Gamma ( The Generalized Gamma is a three parameter distribution which includes many well known distributions as special cases [37]. Its probability density function may be written as where a is the scale parameter, v and p are the shape parameters, x is the pixel intensity and Γ represents the conventional Gamma function. Generalized Gamma distribution reduces to Gamma distribution for f (x; a, v, 1) and Weibull distribution considering v = p. Also with two-parameters, Nakagami distribution is related to the Gamma distribution since given a random variable from the Nakagami distribution X ∼ f (x; m, ω), by setting a = m, v = ω/m and taking the square root of a Gamma variable Y , X = √ Y . Finally, it is possible to reduce GG distribution to a one-parameter distribution defining, f (x, d √ 2, 1, 2), d > 0, which is known as Rayleigh distribution.
The K-distribution is a re-parametrization of the usual form of the family of gamma distributions. It combines two gamma distributions and can be written in the following way, where α = ϕ − L, β = L + ϕ − 1, ξ = Lϕ/µ, and K α is a modified Bessel function of the second kind of order α. The lognormal distribution is a two parameter distribution whose logarithm has a normal distribution. Its probability density function is defined as, where µ is the logarithmic mean and σ is the logarithmic standard deviation. Finally, the Rician distribution has the density function, with non-centrality parameter s ≥ 0 and scale parameter ϑ ≥ 0, for x > 0. I 0 is the zero-order modified Bessel function of the first kind. If x has a Rician distribution with parameters s and ϑ, then (x/ϑ) 2 has a non-central chi-square distribution with two degrees of freedom and non-centrality parameter (s/ϑ) 2 . The parameters (θ) of the presented models were estimated using the method of maximum likelihood [38]. One issue that may need to be considered in modeling the speckle statistics is the presence of the background noise. In general, the pdf of the acquired ROI data corresponds to the convolution of the OCT speckle pdf and that of the background noise, when the speckle and noise are assumed to be independent. Hence, only in the case where the pdf of the background noise approximates a delta function one could assume equality of the acquired ROI data pdf to the OCT speckle pdf. For this, a requirement of high signal to noise ratio is necessary.

Model selection
In model selection, the competing models are fitted, evaluated and ranked. To test whether the kernel density function (Epanechnikov) of the raw data statistically differs from the candidate model, a two-sample Kolmogorov-Smirnov test at level of significance of 5% was used. After selecting which of the competing models does not present a statistically significant difference, goodness of fit and model complexity were balanced.
Goodness of fit describes the discrepancy between observed values and the values expected. In general, a better goodness of fit is obtained by more complex models as they are described by a larger number of parameters. However, increasing complexity also increases the risk of over-fitting the data. Therefore, following the principle of parsimony, models with a smaller number of parameters are preferred. Different models with different combinations of goodness of fit and complexity are found in literature. The comparison of two nested models using the ratio of the likelihood (likelihood ratio test) is one of commonly used criteria for model selection [39]. It compares how much more likely it is that the data comes from a simpler model compared to a more complex one without considering any information about the number of parameters in each model. Akaike's information criterion (AIC) is an extension of the likelihood ratio test, which includes goodness of fit and the model complexity [40]. The AIC value is calculated for each model k by using the following expression, where L k and p k denote the likelihood and the number of parameters associated with model k.
The model with the lowest AIC value is selected as the model that best fits the raw data.

Statement and statistical analysis
In both studies, all subjects were treated in accordance with the principles of the Declaration of Helsinki and informed consent was obtained from all subjects after the goals of the research and consequences of participation had been discussed. All the statistics reported in this paper came from the respective ROI images of each study group. To test the statistical differences between the three age groups, the non-parametric ANOVA (Kruskal-Wallis) test with 0.05 level of significance was used.

Results
The results presented in this study are divided in two parts. The first one, explores which distribution is the most appropriated to fit the corneal speckle from OCT and the second, focuses on its application to analyse age-related and short-term changes that occur in corneal structure. Figure 2 shows the cumulative distribution function (cdf) of the raw data approximated by the Epanechnikov kernel function and the corresponding cdf of each competing model for a subject selected at random. The cdf by itself does not allow to conclude about which models are more appropriate to fit the raw data. Therefore Kolmogorov-Smirnov test was used to test the statistical difference between the fitted data and the corneal raw data. Average and standard deviation of the supremum function of two-sample Kolmogorov-Smirnov statistics is shown in Fig. 3. The plotted data corresponds to ten measurements of four subjects aged between 26 and 53 years in normal conditions (without corneal swelling). The critical value for a significance level of 0.05, is represented by the dashed horizontal line. Values above the critical value mean there is a statistically significant difference between the compared data samples. Raw data modeled by Generalized Gamma, Weibull, Nakagami and Rayleigh distribution did not present statistically significant differences when compared with raw data approximated by the    Epanechnikov kernel function. Therefore, the best model considering the trade-off between goodness of fit and the model complexity was chosen according to the Akaike's Information Criterion. Figure 4 shows the AIC value for each of these models. Despite the highest penalty applied to Generalized Gamma, as a three-parameter distribution, it presented the lowest AIC value which means it is the best in that sense proposed model to fit the corneal speckle. For this reason, the Generalized Gamma distribution was chosen to study whether corneal OCT speckle is influenced by long and short-term changes. Figure 5 shows example of corneal OCT images for subjects with different ages and the respective central corneal thickness. It is well known that CCT changes with age [41] however it is also an anatomically dependent feature. Bigger eye globes may indicate thicker corneas whereas thinner corneas may correspond to smaller eyes. Hereupon, subjects of Fig. 5a and 5d have the same CCT but the latter is twice older. Therefore, we would expect different structural properties due ageing. Table 2 shows the variation of the Generalized Gamma parameters for each age-group. The GG a and GG p decrease with age whereas the GG v increases with age.  The probability density function corresponding to the mean values of table 2 for each age-group is shown in Fig. 6. Figure 7 shows examples of corneal OCT images before inducing swelling (a), immediately after the removal of the eye patch (b), and after 60 (c) and 120 (d) minutes. Table 1 shows the central corneal thickness and Generalized Gamma parameters before patching the right eye (OD) and at the moment of removing the eye patch (0 minutes), 60 minutes and 120 minutes later.

Age-related changes and swelling analysis
The values of the control eye (OS) are also presented for comparison. The last column indicates the p-value between groups based on the non-parametric ANOVA (Kruskal-Wallis) test. A statistically significant difference was observed in the right eye for all parameters (p < 0.001). The control eye (left) did not present statistically significant differences for CCT and GG v , p > 0.05. However, GG a and GG p (shape parameters) showed a slight but significant change over four hours, with p = 0.005 and p = 0.016, respectively. The reasons for such difference are discussed in the next section. Figure 8 shows the GG pdf modeled with the mean values from table 1 for the patched eye. The different lines express how GG pdf changes in the corneal swelling. Despite achieving statistically significant changes of the corneal speckle for each of the subjects, it was not clear how those short-term changes are affected by long-term changes related to age. The shape parameter GG p has showed the highest variation of 25.5% from the baseline after inducing corneal swelling and therefore it was considered the most sensitive parameter to analyze the corneal structure. In Fig. 9 the mean value of GG p and the respective linear regression for each time-point after removing the eye patch are shown, considering four subjects with 25, 29, 45 and 53 years old. The line slope shows how fast each subject recovered from corneal swelling from the micro-structural point of view. According to the obtained results, age seems to be an important factor since younger subjects recovered faster than the older subjects. A similar analyzis was done for macro-structural point of view based on the CCT, however there was no relation observed between the slope and the age.

Discussion
The K-S test did not reject (p > 0.05) the hypothesis of equal distribution for the Generalized Gamma, Weibull, Nakagami and Rayleigh distribution. These distributions belong to the gamma family and as showed in the Methods section, the last three distributions (i.e., W, N, R) can be derived from the GG distribution. A more complex model is expected to fit the data better than a model with less parameters, however it is also more susceptible to over-fitting. Nonetheless, Akaike's Information Criterion has showed that despite the highest penalty being applied to GG distribution, it is the best distribution to fit the corneal OCT speckle data. Similar results were obtained by Ali et al. [42] and Seevaratnam et al. [43] who chose the GG distribution as the best distribution to fit the speckle in OCT skin images and to quantify temperature changes in tissue-mimicking fluid phantoms, respectively.
Tunis et al. [44] showed that GG parameters are sensitive to structural changes within cells. Although their study was based on ultrasound analysis, a similar behaviour in the speckle development may be expected in OCT [23]. It has been suggested that the scale parameter, GG a , measures the average backscattered power and can be related to the average scatterer cross-section. The average backscattered power is dependent on both the geometrical cross-section and scatterer organization and, therefore, the GG a should reflect both. On the other hand, the shape parameters, GG v and GG p , were related to the scatterer density and the ratio of the two shape parameters, can be used as an estimate of effective scatterer number density [44]. The speckle resulting from the interaction of light with the sample scatters changes with the properties of the light source, the propagating beam, and the aperture of the detector [19]. Therefore, for different OCT systems with different light source properties, we may expect a different profile of the speckle intensity histogram. A change in lateral and axial resolution of the beam will probably change the interpretation of the scatter density and the respective back cross-section mentioned by Tunis et al. [44]. Consequently, the Generalized Gamma distribution may not be the most appropriate probability density function to model the corneal speckle of all OCT instruments. In this study, short and age-related corneal changes were studied by modeling and analyzing speckle statistics using Generalized Gamma distribution. Short-term corneal changes were obtained by patching the right eye of each subject for three hours, which may correspond to the swelling that occurs during sleeping time. Corneal swelling occurs when the amount of oxygen falls below an adequate level and a series of acute corneal responses occur, including an increase in stromal lactate, a reduction in intercellular pH, and an increase in corneal hydration. These acute responses are reversible when normal oxygen is restored and corneal tissue returns to its natural state. Corneal swelling has been studied from the macro-structural point of view, where CCT is taken as main reference. However, it is still unknown how the micro-structure is affected in vivo. The studies described in literature regarding to corneal structural properties are limited to ex vivo analysis which may not entirely correspond to in vivo corneal properties. Some authors may argue that corneal hysteresis (CH) obtained by the Ocular Response Analyser (ORA) [ Our results show that corneal speckle based on OCT imaging implied micro-structural changes that occur in corneal tissue. A statistically significant difference for all GG parameters among the different time-groups was observed. The scale parameter decreased significantly after corneal swelling whereas both shape parameters increased. Corneal hydration of the posterior stroma may contribute to an increase of the inter-fibrillar space reducing the backscatter cross section of the speckle and consequently increasing the scatter density. In this specific case, what may occur is not an increase of the number of scatters but a higher differentiation of the fibrils due the increase of the inter-fibrillar space. The second shape parameter has shown to be the most sensitive to these structural changes (25.5%), approximately three times more than the variation in CCT (7.8%). Variation of GG p with age was observed. This may be of importance to calibrate IOP values, mainly when subjects have the same CCT but different ages. Certainly, we would expect an older subject to have stiffer cornea than a young subject even despite having the same CCT. However, so far, there was no way to measure such differences in vivo. In addition it was seen that unlike the CCT, GG p did not reach the baseline values after two hours of de-swelling. These results suggest that macro and micro-structure are not directly correlated and do not recover at the same rate. A statistically significant difference was also observed in the control eye for both shape parameters. However, their mean value did not present any apparent trend. Such statistical difference may have an origin in two different sources. Firstly, daily changes may occur and consequently affect the corneal micro-structure and secondly, the fibril structural organisation is slightly different from the corneal apex to the limbal region and the acquired images may not exactly correspond to the same cross-section.
A larger ROI including all stroma was used to analyze the influence of age-related changes on the signal envelope statistics. The information in that ROI came mainly from the collagen fibrils embedded in an extracellular matrix rich in proteoglycans, glycoproteins and keratocytes. The micro-structure of the stroma suffers from age-related changes including an increase in the cross-sectional areas of both collagen fibrils and fibrillar molecules which may be due an increase in non-enzymatic cross-linking [7,8]. In addition, there is also a decrease in the inter-fibrillar spacing, related to changes in the proteoglycan composition of the inter-fibrillar matrix. These structural changes in collagen fibrils play an important role in corneal biomechanics contributing to tissue stiffening with age. This may be an important issue for glaucoma management since a few studies have reported an overestimation of intraocular pressure as a consequence of the ageing role in corneal biomechanics [49]. In our study, age-related changes are reflected in corneal speckle statistics. It was seen that age-related changes narrow the distribution around its centre and cause a shift along the intensity axis in the direction to a lower intensity. This is understood as an increase of the backscatter cross-section and a reduction of the scatter density which may originate, essentially, from the reduction of inter-fibrillar spacing.
In addition to the biomechanical and elastic properties of the tissue, speckle is influenced by the system characteristics, such as size and temporal coherence of the light source, multiple scattering and phase aberrations of the propagating beam, and the aperture of the detector [50]. All of these variables certainly contribute to a different speckle behaviour in optical coherence tomography of living tissue. In any case, it is straightforward to extend the proposed method to other instruments where access to raw images is provided.

Conclusion
In this study, we explored OCT speckle from the signal-carrying perspective to analyze short-term and age-related changes in corneal structure. The presented experimental findings show that Generalized Gamma distribution is a good model to fit the corneal speckle. The parameters of Generalized Gamma distribution can be utilized to quantitatively evaluate the structure and elastic properties of the cornea in-vivo. This is of high interest since there is no standard method to infer about the corneal micro-structural properties. Moreover, the proposed technique of statistically modeling OCT speckle provides complementary information to better understand alterations of the collagen framework in-vivo, which can be further applied in refractive surgery, ocular tonometry, corneal transplantation and other areas of ophthalmology where quantified information on micro-structure of corneal tissue is of essence.

Funding
This work was supported by the Marie Curie Innovative Training Networks, Ageing Eye, AGEYE, 608049.