Corneal hydration assessment indicator based on terahertz time domain spectroscopy.

Terahertz technology has shown broad prospects for measuring corneal water content, which is an important parameter of ocular health. Based on terahertz time-domain spectroscopy, a new indicator named characteristic ratio (CR) of the sum of low (0.2-0.7 THz) and high (0.7-1.0 THz) frequency spectral intensities, for characterizing corneal hydration is introduced in this work. CR is calculated from the real-time reflection spectra after error elimination of ex vivo human corneal stroma samples which is collected during dehydration under natural conditions (temperature: 22.4 ± 0.3°C; humidity: 20.0 ± 3%). The corresponding relationships between CR and corneal water content are reported. Comparing the linear fitting results with the published similar study, the coefficients of variation of the fitting slope and intercept are 39.4% and 27.6% lower, respectively. This indicates that this approach has the potential to achieve corneal water content in-vivo detection in the future.


Introduction
The stroma constitutes the largest portion (>90%) of cornea which is the outermost structure of the eye [1]. Therefore, many crucial characteristics including physical strength, shape stability and transparency of the cornea are mainly determined by the anatomic and biochemical properties of stroma [2,3]. Corneal tissue water content (CTWC) refers to the percentage of the moisture weight in corneal tissue to the tissue weight after drying. Various ophthalmic diseases, such as corneal edema and Fuch's endothelial dystrophy [4,5] are accompanied by corneal hydration disorders [6][7][8]. In addition, it has been demonstrated that CTWC could significantly affect the tissue ablation rate of most popular keratorefractive surgeries, such as Laser-Assisted in situ Keratomileusis (LASIK), Laser-Assisted Sub-Epithelial Keratectomy (LASEK) and so on [9][10][11][12]. Consequently, CTWC is a vital parameter for clinical ophthalmic diagnoses and cornea-involved surgeries.
Currently, the available clinical diagnosis of CTWC is extrapolated from the central corneal thickness (CCT) measurement by ultrasound or optical coherence tomography (OCT) techniques [13][14][15]. However, this method does not allow for sufficiently early detection of water content disturbances. On the one hand, there are inaccuracies in this indirect mapping relationship derived from the statistical results. On the other hand, it fails to consider the CCT variation among different biological individuals. Moreover, confocal Raman spectroscopy has been reported to measure axial hydration gradient of the rabbit cornea [16]. However, its safety issue still remains to be evaluated.
Recently, terahertz (THz) wave has been proposed as a safe and sensitive method to probe CTWC owing to its low photon energy and high sensitivity to water [17][18][19][20]. Previous reports use THz reflectivity as the key parameter to extract the corneal hydration level. The principle has been verified through the experimental results on ex-vivo porcine corneas [21,22] and in-vivo studies of rabbits [23] and humans [24]. By scanning a focused THz beam on volunteer's eye, Sung S., et al. have obtained the reflectivity maps which represent spatial hydration distribution across human cornea [25][26][27]. But the reflected signal intensity is susceptible to all sorts of factors arising from instruments and samples, which will limit the accuracy and repeatability of measurements and make it challenging to calibrate the relation between CTWC and THz reflectivity.
In this paper, we showed the preliminary results of application of THz time domain spectroscopy for monitoring the long-time dynamics of corneal tissue drying and sensing corneal hydration. By real-time acquisition of the sample reflection spectra during dehydration, the ratios between the sum of the low-frequency band (0.2-0.7 THz) and the high-frequency band (0.7-1.0 THz) spectral intensities could be extracted. The ratio referred to as characteristic ratio (CR), was used as a new indicator for the characterization of CTWC. The essence of the present work lies in that the characteristic absorption peaks of water vapor in atmosphere has been used to eliminate the interferences caused by complex factors during the THz reflection spectra measurement. Finally, the statistical results of the relationship between CR and CTWC with better consistency were obtained. Since this approach was a direct detection of corneal hydration information, it had the potential to be a valuable supplement to CCT measurement.

Spectral measurements
The experiments were carried out by our home-made terahertz time-domain spectroscopy (THz-TDS) system [28], which uses two commercial photo conductive antenna (TERA15-FC, Menlosystems, Germany) as the emitter and receiver, respectively. Corneal stromata excised from the healthy human body during femtosecond laser small incision corneal refractive surgery (i.e., small incision lenticule extraction, SMILE) were used as specimens. A total of four corneal samples were all stored in the Ringer's solution to keep their optical characteristics intact [29]. Each sample is approximately 6.6 mm in diameter and 104-138µm in thickness when excising from the body. During the measurement, the corneal samples were all exposed to constant natural environment (temperature: 22.4 ± 0.3°C; humidity: 20.0 ± 3%). The hemispherical design of the sample holder shown in Fig. 1 mainly took account of that the radius of curvature is between 7.5 and 8.0 mm at the 3 mm central optical zone of human cornea and the surface is almost spherical [3]. Meanwhile, the Ringer's solution filled the central groove, so that sample's back surface could come in direct contact with it. The purpose of this was to exclude the effect of the Fabry-Perot interference caused by large differences of refractive index between interfaces. Therefore, it was possible to avoid the introduction of new interference factor to spectral intensity. It is worth mentioning that phase ambiguity introduced by multiple reflections, i.e., F-P effect, as discussed in Ref. [30] could also be safely ruled out. In addition, in order to minimize the misalignment when changing samples, the position of the sample holder was carefully referenced by using a 3-D translation stage with spatial precision of 10 µm.
The terahertz radiation generated by the transmitter (Tx) was collimated into a parallel beam via an off-axis profile mirror. After being deflected by the mirror, a TPX lens with focal length of 5 cm focused terahertz wave on the surface of sample. Then the sample reflection signal passed through another set of identical TPX lens and mirror, and was finally collected by the receiver (Rx). A series of spectra collected during dehydration needed to be processed at first to eliminate the effects of interferences. After this, CR corresponding to each spectrum at every dehydration time could be calculated.

Elimination of interference factor effects
The original spectra shown in Fig. 2 might be subjected to the intensity modulations of various factors, such as changes of corneal thickness and surface curvature, fluctuations of femtosecond laser power and so on. Those factors that affected the spectral intensities would interfere with the subsequent calculation and extraction of corneal hydration information contained in the spectrum, so they might be considered as sources of error in the study. Considering the impact analysis of each factors, on the one hand, limitations on changes in corneal thickness during dehydration made it very difficult to use the reference mirror as calibration. On the other hand, it requires extensive investigation to confirm whether the final influence of various factors could be regarded as a simple superposition. The complicated interference sources would also greatly increase the complexity and workload of the measurement. Hence, the spectral treatment proposed here did not tediously analyze every possible influencing factor that might exist in the experiment. Rather, in view of that when the ambient humidity was constant during measurement, the value of absorbing peak of water vapor in atmosphere contained in the spectrum should be constant.
In order to facilitate the analysis, the sample reflectance spectrum could be equivalently regarded as the modulated result of the transmission spectrum of air (S air (f )). Modulation factors included collagen (the main constituent of the corneal stroma except water, S collagen (f )), moisture in cornea (S H 2 O (f )), and combined effects of interferences (β(f )), while f all denoted THz frequency. Therefore, the measured reflection spectrum (S(f )) could be understood in terms of the following Eq. (1): Then, it was assumed that the sample was in a dry air environment (S dryair (f )) which was completely free of water vapor. Therefore, the actual measured signal could be regarded as a result of the transmission of a gas chamber filled with water vapor and the reflection of the sample surface. The virtual water vapor chamber set here, as shown by the black dotted box in Fig. 1, represented equivalently the effect of experimental constant temperature and humidity environment on signals (S vapor (f )). And Eq. (1) could be rewritten as: Based on the above analysis, the intensity of the characteristic absorption peak of water vapor in spectrum should remain unchanged theoretically in the constant temperature and humidity environment. Then the variation of intensity existed in the actual situation reflected the modulation effect of β(f ). That is, it could be considered as a reference to guide the elimination of complicated interferences.
Since interferences caused by certain factors were sensitive to frequency, multiple water vapor characteristic absorption peaks were required to correct the intensity interferences of different frequency segments. However, it was restricted to the number of water vapor absorption peaks in the frequency range of 0.2-1.0 THz, which band had a good signal-to-noise ratio in the sample reflection spectra. The fluctuations of two water-vapor peak intensities at 0.56 THz and 0.76 THz were finally chosen to express the average interference levels of the low frequency band (0.2-0.7 THz) and the high frequency band (0.7-1.0 THz), expressed as β(f l ) and β(f h ), respectively. For the purpose of extracting β(f l ) and β(f h ), two "reference values" of the water-vapor peak intensities, denoted by β 0 (0.56THz) and β 0 (0.76THz), need to be set firstly.
Equation (3) describes that the mean of intensity values of the water vapor absorption peaks at 0.56 THz and 0.76 THz in all spectra was computed separately as the "reference value", where N is the total number of spectra obtained from the experiments of four samples, and S i (f ) represents one spectrum corresponding to each sample at a certain moment. Thereafter, β 0 (0.56T) and β 0 (0.76T) could be used to obtain the intensity deviation multiples of low frequency band and high frequency band in each spectrum. These two multiples were not only the intensity fluctuation levels of the two water-vapor peaks at 0.56 THz and 0.76 THz, but also represented the average interference impact ratio of the two frequency intervals, that is, β(f l ) and β(f h ): Finally according to Eq. (2) by utilizing the β(f l ) and β(f h ) to correct the intensities of the two bands, the spectrum after interference elimination (S (f )) were expressed as Eq. (5) and displayed in Fig. 3. Note that the change of spectral intensity (at 0.56 THz and 0.76 THz) caused by moisture variation in the sample was negligible when compared to the characteristic absorption of water vapor [21].

Derivation of corneal tissue water content
Due to the limitation of the system space, the electronic balance could not be used to measure the sample weight in real time in this experiment. Since the change of optical path could be estimated according to the time delay of the reflected signal and the geometric path of THz beam shown in Fig. 4(a), it was possible to infer the geometry thickness variation of the sample during dehydration. Before the THz time-domain spectrum measurement, an electronic balance (accuracy: 0.0001g) was used to record the wet weight m 0 of corneal sample at the initial time t 0 .The sample was then placed into the THz optical path, and a set of THz time-domain signals with corresponding spectral data were collected every 40-60 seconds. As the cornea dehydrated, the variation in thickness caused the translation of the THz time domain signal in the time window. The time delay ∆t is difference value of each time-domain signal's position compared with the signal of t 0 , as shown in Fig. 4(b). And then the THz optical path difference could be calculated through multiplying ∆t by the speed of light c.The half of the optical path difference was the right-angle side length in the triangle ABC, so the height BD of the triangle could be calculated, that is, the amount of change in corneal thickness ∆h. Further, based on the empirical relationship between corneal thickness and its water content, a derivation from the amount of change in signal's delay time to the reduction of corneal water content could be achieved. The above derivation process was expressed by the following Eq. (6): where ∆CTWC is the variation of corneal water content, and 40%/mm is the approximate slope of the linear fit of the relationship between CCT and CTWC within the physiologically relevant range [31]. Therefore, it was only necessary to measure the dry weight m dry of the sample after dehydrating in a 90°C drying oven for 3 hours. Thus, the water content of sample under different hydration conditions could be derived according to signal time domain information by Eq. (7): In order to prove the feasibility of the above processing method, experiments on the rehydrated samples to obtain the relationships between variation of corneal water content and dehydration time was supplemented by means of weighing method. Under the premise of the same temperature and humidity environment, the rehydrated corneal stroma sample was taken out from the Ringer's solution and placed on the electronic balance weighing platform. The weight reading was recorded every 20 seconds, and the relationship between the dehydration time t and the CTWC (%) was finally obtained through the Eq. (8): where m(t) represents the net weight of the sample at time t. In Fig. 5, data obtained by weighing method were compared with the results of time delay estimation method, and the two results had good consistency. Moreover, The differences between weighing method and time delay estimation for Cornea1 ( Fig. 5(a)) and Cornea4 (Fig. 5(d)) may be attributed to two aspects. One is the inaccuracy of the derivation relationship between CCT and CTWC in time delay estimation. The other is that multiple hydration may inevitably affect the structure and function of corneal stroma to varying degrees, thereby the dehydration rate as well. The thickness variations has been estimated according to the broadening of the signal pulse waveform with respect to the reference pulse waveform reflected by a flat gold mirror. It indicates an averaged thickness variation of less than 12 µm. It could lead to approximately a CTWC measurement uncertainty of 0.5% according to the results shown in Fig. 5. Hence, the effect of the thickness variations in the corneal samples within the THz spot would not play a significant role in our work.

Computation of characteristic ratio
Both results of published research [21,23] and our experiment of this work (Fig. 3) reflected the fact that THz sensing to corneal hydration was more sensitive in the low frequency band. So, the low-frequency stage of spectrum could be regarded as the characteristic frequency band mainly carrying the sample hydration information, while the high-frequency stage could be seen as the reference frequency band mainly containing information about components other than water. After the spectral interference elimination, there is: Then for each spectrum, two intensity sums were obtained according to the frequency segment, and finally the characteristic ratio, CR, was received by dividing S i (f l ) by S i (f h ).
According to Eq. (10), for the same system performance (mainly referred to the energy collection efficiency of high frequency and low frequency) and constant experimental environment, the impact of S vapor (f ) and S dryair (f ) could be represented by a constant term k: Thus, it could be clearly fould out from the simplified Eq. (12) that CR was an applicable parameter which was defined to express the state of corneal hydration:

Results
The availability of the spectral processing proposed above was demonstrated and judged firstly through the experimental results by considering two types of interference factors respectively.
(1) Geometric deformation of the sample. Variation of corneal stroma thickness resulted from dehydration performed as the change in the height of signal reflection in measurement. A mirror was placed on a platform which could achieve one-dimensional translation in the vertical direction, and then the mirror was raised in steps of 0.04 mm (total translation of 0.68 mm). Eighteen sets of reflection spectra were collected, and some of them were plotted in Fig. 6(a). It was not difficult to find from original spectra that this factor ought to belong to the frequency-sensitive type.
(2) Power fluctuation of the femtosecond pulsed laser. Because the dehydration process of cornea exposed to the natural condition would continue for a period of time, the comparability between different spectra might be decreased when femtosecond laser power was instable. Tuning the femtosecond optical power by means of adjusting fiber coupling efficiency, a series of air transmission spectra were acquired shown in Fig. 7(a). The results calculated according to the above spectral processing method were plotted in Fig. 7(b).
The results in Fig. 7 confirm one of the advantages of our method: it is independent on the input laser power fluctuation (also the power of THz wave). However, it is so crucial in conventional method that a stable reference signal is required.
So overall, the data processing method of this study had a good elimination effect for different types of factors. Based on analyses above, the spectral calculation results of four corneas were plotted in Fig. 8, with water content derived from the method of time delay estimation. Then the four sets of discrete data were linearly fitted respectively. The calculation of coefficients of variation (COV) was performed on the slope and intercept of the fitting results in order to quantify the improvement of our results. Comparison with the published relationships of terahertz reflectivity and corneal water content was shown in Table 1. As displayed in Table 1, our results had a 39.4% reduction in the coefficients of variation of fitting slope and a 27.6% decrease in intercept (both compared with the better one of the previously published results), while it also ensured a high goodness of fit. The calculation of COV reflected the better comparability and consistency of the results for different samples in this study. By using the fitting results of Fig. 8, the corresponding estimated CTWC values fluctuated from 74.27% to 88.95% in the case of unprocessed data and from 73.19% to 79.89% in the case of processed data indicated in Fig. 6(b), when the geometric displacement was considered as the main source of the measurement error. Similarly, the corresponding fluctuation of the deduced CTWC values changed from 72.70% to 74.40% in the case of unprocessed data and varied from 72.62% to 75.83% in the case of processed data shown in Fig. 7(b).
By referring to the measurement results of Ref. [32], the effect of vapor concentration fluctuation on our measurement has been carefully evaluated that it may play a negligible role  Table 1, although significant improvement has been achieved by our work, there still appears to be residual inter-sample variation in Fig. 8. The cause might be the water content estimation errors among individual samples as shown in Fig. 5, which requires more statistics study in future.
Besides, in view of the ultimate goal of achieving detecting corneal water content in vivo, the relational expression for CTWC inversion calculation must be based on a larger statistical sample size and a more detailed classification of population features in the future. As such, its reliability and applicability would be guaranteed.

Summary
A new indicator for corneal water content assessment based on terahertz time domain spectroscopy was proposed in this study. The reflectance spectrum of the stroma excised from human cornea during dehydration in natural conditions (temperature: 22.4 ± 0.3°C; humidity: 20.0 ± 3%) was systematically acquired in real time. For the data processing, we took advantage of the characteristic absorption peaks of water vapor in atmosphere to represent the effect of interferences generated by multifarious sources, and finally obtained a good elimination effect. By extracting CR which contained the corneal hydration information from the spectrum, the relationship between CR and corneal water content was finally obtained and then linearly fitted. The coefficients of variation of slope and intercept were reduced respectively by 39.4% and 27.6% when compared to the previously published results. This meant that the relationships between different samples were more comparable and consistent. Consistency was a requirement that must be met of the established relationship for CTWC inversion calculation in in-vivo study. In summary, based on the results of this research, we had the potential to continue to develop its broad application prospects of non-invasive CTWC detection in the future with its advantages of non-contact, higher comparability and consistency between samples with no complicated calibration.