Investigation of water diffusion dynamics in corneal phantoms using terahertz time-domain spectroscopy.

Perturbation of normal corneal water content is a common manifestation of many eye diseases. Terahertz (THz) imaging has the potential to serve as a clinical tool for screening and diagnosing such corneal diseases. In this study, we first investigate the diffusive properties of a corneal phantom using simultaneous THz time-domain spectroscopy (THz-TDS) and gravimetric measurements. We will then utilize a variable-thickness diffusion model combined with a stratified composite-media model to simulate changes in thickness, hydration profile, and the THz-TDS signal as a function of time. The simulated THz-TDS signals show very good agreement with the reflection measurements. Results show that the THz-TDS technique can be used to understand water diffusion dynamics in corneal phantoms as a step towards future in vivo quantitative hydration sensing.


Introduction
The human cornea is a transparent and avascular tissue that transmits light to the retina. Transparency of the cornea is maintained by the critical spacing and organization of collagen fibers and relative state of dehydration. The endothelium is the innermost layer of the cornea and employs a tightly controlled aqueous transport system to supply nutrients and metabolites to the cornea and maintain corneal clarity [1]. Intraocular surgery and corneal injury can cause dysregulation of this transport system and would result in a perturbation of normal corneal water content. This change in water content is also a common manifestation of corneal diseases such as keratoconus [2], Fuch's dystrophy [3] and dry eye syndrome. In current clinical practice, corneal water content can be determined by first using pachymetry to measure corneal thickness and then calculating water content using a linear relationship determined by Ytteborg and Dohlman in 1965 [4]. Recent in vivo studies implemented Terahertz (THz) and mm-wave methods to determine corneal hydration compared to the calculation method based on pachymetry [5]. These results suggested that a linear function may not be adequate in explaining the relationship between corneal hydration and corneal thickness.
THz imaging has been shown to be a suitable method for assessing hydration of biological tissue due to its sensitivity to small changes in water content. Taylor et al. showed that THz reflectivity has a linear relationship to water content using polypropylene filter paper with varying saturations of water [6]. Other materials, such as polyester/cellulose clean room wipes, were also used to observe dehydration process using THz reflectivity and a precision scale with 0.01 g resolution [7]. The sensitivity of THz reflectivity to fluctuations in a sample's water content between 75% and 85% water was characterized using Noise Equivalent Reflectivity Difference (NERD) for a value of 0.1204% [8]. Guo et al. showed how a coherent THz holography system can be used to understand the dehydration dynamics of pork, mutton and cattle tissue [9]. Furthermore, THz spectroscopy has been used in conjunction with water diffusion models to estimate the reflectivity of occluded human skin [10,11]. Therefore, THz Time-Domain Spectroscopy (THz-TDS) and imaging are promising tools for probing the water content of in vivo and ex vivo biological tissue.
The high homogeneity and relatively low physiological variability of corneal tissue make THz radiation a fitting tool for hydration measurements [5,6,[12][13][14][15][16]. In addition to a sensitivity analysis, the study by Taylor et al. showed the feasibility of THz imaging for cornea analysis by performing in vivo experiments on dehydrating healthy rabbit corneas through a mylar window [8]. More recently, efforts have been made to omit the need for an imaging window and develop a non-contact and normal-incidence THz imaging system designed specifically for in vivo and clinical assessment of corneal hydration [16,17].
In this study, we aim to develop THz-TDS experimental methods and theoretical diffusion models to investigate the water diffusion dynamics in corneal phantoms. We utilized an analytical balance (1 mg accuracy) and a non-contact THz-TDS system to observe the dehydration dynamics and THz reflectivity of a contact lens sample. Previous THz dehydration studies relied on either overall hydration measurements or literature values for in vivo hydration profiles with constant sample thickness to simulate the THz response of the sample [5,7,9,18]. We created a diffusion model, based on Gu et al. [19], to predict a thickness and hydration profile at each measured time point in the drying process. The predicted hydration profiles were then input in a stratified composite-media model to simulate THz-TDS measurements. The modeled THz-TDS signal showed very good agreement with the measured data. Coupling a dynamic diffusion model with THz measurements can be extended to the analysis of water content in other diffusive samples and biological tissues.

Gravimetric measurements
Corneal phantoms with tightly controlled material parameters are challenging to fabricate because of their spherical geometry and sub-millimeter thickness. Contact lenses, however, possess a spherical shape with a radius of curvature near the size of the human cornea. Moreover, the water content and dimensions (which depend on the diopter rating) of contact lenses are precise and standardized because of a well-regulated manufacturing process. Most contact lenses are created using hydrogel materials with lower water contents than the cornea (which is approximately 75% to 85%) [8,20]. Insofar as the primary experimental concerns are a relative study of water content in a spherical-cap shape with variable thickness, contact lenses will work well as corneal phantoms.
The contact lens used in this study was a Hydrasoft model with −20 diopter power, 8.6 mm base curvature, 14.2 mm diameter, and a reported initial water content of 55% (CooperVision Inc, Lake Forest, CA, USA). The edge thickness of the contact lens was determined by a cross-section of the center of a hydrated lens taken using a stereoscope.

Terahertz system
THz measurements were taken with an ASynchronous OPtical Sampling (ASOPS) THz System (Menlo Systems Inc. Newton, NJ, USA) which uses two 1560 nm femtosecond fiber lasers to simultaneously pump an InGaAs/InAlAs THz emitter and probe a LT InGaAs/InAlAs receiver at slightly offset repetition rates. The ASOPS parameters were set to a sampling rate of 100 MHz and a difference frequency of 100 Hz. An ASOPS measurement system was ideal for this study because it allowed for rapid time-domain scans with high signal-to-noise ratios. Each measurement consisted of an average of 100 time-domain scans and were taken at five second intervals. Figure 1(a) shows the contact lens mounted on a stainless-steel sphere (J.W. Winco, Inc. New Berlin, Wisconsin, USA) with a diameter of 16.002 mm. The diameter of the sphere was chosen to be close to the base curvature of the contact lens so that the mounted sample would be secure and stationary during the experiment. The mounted sample was secured to the measuring pan of a 1 mg accuracy laboratory analytical balance (Ohaus Inc, Parsippany, NJ, USA). The collimated THz beam was focused onto the center apex of the stainless-steel sphere using a high-density polyethylene (HDPE) f-θ lens. The choice of an f-θ lens, which is part of a collocated raster scanning reflection setup developed in our laboratory, is not essential for the experiments conducted in this study. As shown in Fig. 1(b-c), a reference THz measurement was obtained from the stainless-steel sphere before mounting the contact lens. The THz-TDS sample measurements were obtained from the apex of the sphere during the dehydration experiment. The signal-to-noise ratio of the reference was about 60 dB, and the data acquisition time for the averaged 100 spectra was 1 second, given the described ASOPS settings.

Water-loss model
Changes in a hydrogel's water content are typically accompanied by a change in its physical dimensions, i.e. its thickness [21]. Hence, we implemented a hydrogel diffusion model that calculates the thickness and water content profiles using initial parameters from the manufacturer, initial geometry and gravimetric measurements to determine water loss [19]. The flowchart in Fig. 2 shows the division of this model into two parts. The first part, describing the semi-infinite thickness case, is a simplified approach with an analytical solution that provides an estimate for the mass transfer coefficient, α, and the constant diffusion coefficient, D 0 , which will be defined in Section 2.3.3. These parameters are then entered in to the second part, which represents the hydrogel as a slab with a finite and time-dependent thickness. This variable-thickness model is a more realistic scenario without an analytical solution but can be solved numerically to extract the thickness and the water content profile. These variables are then used in the stratified media equations to simulate the dielectric response of the sample.

THz-TDS Reference Signal
Sample Geometry  Fig. 2. The flowchart for the two models used in this study is shown along with the parameters extracted from each step. This flowchart provides an overview of the models before describing them in detail in later sections.

Assumptions
We assumed that the hydrogel, methafilicon B, had a density of 1.07 g/cm 3 and was stable over the course of the drying process [22]. Moreover, environmental conditions, such as temperature and humidity, were assumed to be stable over the course of the experiment. The center thickness of the sample was reported by the manufacturer. Finally, the surface area of the contact lens, as probed by the THz beam, was assumed to be constant during the experiment.

Geometry
Optometrists are provided with the base curvature from contact lens manufacturers to properly fit a patient's eye. The front curvature of the contact lens, however, is typically proprietary information. The curvature of the front surface is needed for determining the thickness profile and the surface area of the contact lens. We approximated the geometry of the lens using the center thickness, density, and base curvature, which were provided by the manufacturer, and our measurement of the initial weight. Given the symmetry of the sample, we chose a compound quadratic polynomial of the form ax 4 + bx 2 + c to represent the upper surface curve. The constants in this equation were optimized to give the known values of center and edge thickness, as well as the total volume of the sample. The resulting geometry, as shown in Fig. 3, agrees with general thickness profiles of the negative power contact lenses taken with optical coherence tomography (OCT) [23], as well as an evaluation of the sample's cross-section from a stereoscope. The thickness at the center of the contact lens was about 60 µm, at the thickest region was 635 µm, and at the edges tapered to 200 µm.

Semi-infinite model
As illustrated in Fig. 4(a), the dehydration of the sample is initially modeled by a one-dimensional diffusion process through a semi-infinite medium. The water concentration, C, (in units of g/m 3 ), at each measured time point is determined by, where ρ hydrogel is the density of the hydrogel and W is the water content (% mass). The initial concentration, C 0 , was determined using Eq. (1). The water content was defined as the mass fraction of water in the hydrogel where mass hydrated is the mass of a fully hydrated sample and mass dry is the mass of a dried sample. At the final time point, the sample is equilibrated with the atmosphere. We found the initial hydration to be 57.8%, which was close to the manufacturer reported water content (55%). Other studies also reported a ±3% difference between their measured water content and manufacturer value [24,25]. Fick's 2 nd law of diffusion governs the transport of water in our model, given by, where D 0 is a constant diffusion coefficient in units of m 2 /s, and the concentration, C, is a function of both axial position, x, and time, t [26]. In this scenario, it is assumed that the hydrogel is an infinitely thick slab, where the water content, W, at x = −∞ is always equal to the initial water content described in Fig. 4(a). At the surface of the sample, x = 0, the boundary condition for Eq. (3) is given by, The rate of evaporative flux is governed by D 0 and the mass transfer coefficient, α. C eq is the water concentration of the hydrogel at equilibrium, which is set equal to the water concentration of the surrounding environment (∼0.8% by mass). C s is the concentration at the surface boundary of the film. The analytical solution for the semi-infinite model is found by solving Eq. (3), given the initial and boundary conditions, for water loss per unit area, M si (t), where p = α/D 0 [26]. As shown in Fig. 4(b), M si (t) was fitted to the linear region of the measured water loss per unit area. The surface area was determined by the rotation integral of the polynomial expression from Section 2.3.2. The analytical solution begins to diverge from the measurements once the water concentration at the bottom of the film deviates from the initial water content value, which invalidates the semi-infinite thickness assumption. Thus, we fit the values of α and D 0 to the linear region of the water loss data in Fig. 4(b), where our assumption is valid [19]. These variables serve as inputs for the second part of the model.

Variable-thickness model
The variable-thickness model considers a more realistic scenario where the water content of the sample is finite and influences the thickness. As water leaves the sample, its thickness decreases, as shown in Fig. 5(a).  The variable thickness model is also governed by Fick's 2 nd law of diffusion. The diffusion coefficient, D(C), is more realistically represented as a function of concentration, At the hydrogel-sphere interface, we assumed that there was no diffusion or flux because the evaporative flux only occurs at the surface, therefore, The surface boundary condition accounts for the thickness of the sample over time, h(t), and is given by, where α was obtained from the semi-infinite model. C s is the concentration at the surface or at x = h(t). The diffusion coefficient of water, in units of g/m 3 , in the hydrogel, D(C), is represented as an exponential function of water content [21,27], where D eff is the effective diffusion coefficient of the unbound water in the hydrogel, C is the dimensionless concentration, or C/ C 0 , and k is a scaling constant [25]. The optimized D 0 value from the semi-infinite model was used as a starting point for D eff in the fitting. Unlike D 0 , which approximated the diffusion process as a constant value, D(C) describes this process over a range of hydrations. The thickness of the contact, h(t), is given by, where l is the starting thickness of the sample andC(t) is the average concentration. The expansion factor, S, describes the degree of swelling for similar contact lenses [28]. Because there is no analytical solution to the variable thickness model, it was numerically solved using the Finite-Difference Time-Domain (FDTD) method. The FDTD version of the Fick's equation [29] is given by, The subscript n represents the mesh points along the x axis and the superscript t denotes steps in time. To account for the decreasing thickness as the hydrogel dries, we implemented a invariant moving mesh scheme [30]. This scheme decreases the spatial intervals by ∆x to match the thickness changes over time. The total quantity of water lost per unit area for the variable thickness model, M vt (t), is found by using, These equations assume that the diffusing medium has a uniform thickness in the lateral direction. As shown in Fig. 3(a), the thickness of contact lenses can vary by hundreds of micrometers over small distances, especially at high power. Thus, the thickness profile of the contact lens was discretized into 60 sections, illustrated by the red lines in Fig. 3(b), where each discrete section has a uniform thickness. For each section, the water loss per unit area was modeled using Eqs. (5)- (12) and was given by M vt i (t), where i describes each section. The final M total (t) from Fig. 5(b) was determined by taking the weighted sum of all M vt i (t) for all sections, given by, where A i is the percentage of total surface area for each section and N is the total number of sections. The D eff and k values from Eq. (9), were obtained by minimizing the error between M total (t) and the measured water loss per unit area. D eff and k values were kept constant between each section. The resulting fit, shown in Fig. 5(b), had an excellent agreement with the measured water loss data.

Electromagnetic model
The theoretical model for the dielectric response of the corneal phantom in this study is based on the stratified media model from Taylor et al. and Bennet et al. [8,18]. Modeling dehydrating tissue as a stack of composite dielectric media is a more realistic representation of our measurements because the water loss at the surface of the sample is dominated by evaporation and the water loss at the deeper sections is dominated by diffusion toward the surface. Consequently, the water content at each axial position will vary, so a model with discrete dielectric slabs of varying water content, represented in Fig. 6, is an adequate way to model THz propagation through our sample. Each discrete dielectric slab is modeled as a binary composite medium of water and the hydrogel polymer. The dielectric constant of water is given by the double-Debye equation, where ε s is the static dielectric constant, ε 2 is an intermediate dielectric value, ε ∞ is the high frequency dielectric constant, and τ 1 and τ 2 are relaxation time constants that represent hydrogen bond processes and structural rearrangements [31]. The double-Debye equation has been empirically shown to match the real dielectric constant of water [32,33].
The effective dielectric constant was determined using Bruggeman effective media theory for each layer, i, given by, where ρ w,i is the volume fraction of water of each layer, N is the number of layers, which is set equal to 11 in each section in this study,ε w is the dielectric constant of water,ε i is the effective dielectric constant andε h is the dielectric constant of the polymer in the hydrogel [5]. The values of ρ w,1 , at the hydrogel-air interface, and ρ w,N , at the hydrogel-metal interface, are calculated by dividing C(x, t) by the density of water. At t = 0, ρ w,1 = ρ w,N = 0.61, which is in agreement with the manufacturer's water concentration specification. Our variable thickness model then predicts the water loss with ρ w,1 decreasing faster than ρ w,N . The value ofε h was found by fitting a simulated THz-TDS signal to the measurement with the lowest water content. The resulting index of refraction was similar to that of PMMA, a polymer with a similar structure, but with slightly higher complex index of refraction [34]. The larger complex value can be attributed to the hydroxy groups in the hydrogel, thereby making it more polar than PMMA. Material parameters of our samples are challenging to compute because of the nonuniform thickness profile of the contact lens within the beam spot size.
Once the dielectric properties of each layer were determined, the propagation of the THz wave with normal incidence at each layer was calculated using the Fresnel reflection coefficient, r i , and the effective electrical length, δ i , by and where d is the uniform layer thickness, λ is the wavelength andñ i = √ε i . The total reflected signal for each layer can be calculated by solving the multiple internal reflections in a geometric series [5,18,[35][36][37]. The aggregate reflection coefficient, Γ i , at the air-hydrogel surface was then calculated using a recursive function that accounts for all reflections from each layer. This method, commonly referred to as Rouard Method [36], starts at the deepest layer, N + 1, and provides a relationship for the aggregate reflection coefficient, given by, To account for the nonuniform thickness of the contact lens, we used a weighted average of the thickness and hydration profiles of 25 discretized sections within the 2.7 mm beam focus spot, where the weights were based on the Gaussian shape of the THz beam. The recursive function described in Eq. (18) considers any changes in the phase due to both propagation through each layer and the complex index of refraction,ñ. Therefore, the final time-domain result, E model (t), can be acquired by convolving the aggregate reflection coefficient by the reference pulse from the stainless-steel sphere, E ref (f ), and taking the inverse Fourier transform given by,

Water-loss measurements
Assuming the hydrogel material is conserved, and that weight-loss only reflects evaporation, Fig. 1(b) shows the overall water loss as a function of time. The extinction coefficient of water is relatively high in the THz regime [32], so we expect reduced attenuation in the reflected THz signal as the sample dries. The real part of the index of refraction of water at most THz frequencies is higher than that of the hydrogel, so we expect a larger reflection coefficient amplitude when our sample has a higher water content. This behavior is shown by the time domain traces in Fig. 7, where with high water content (blue) the amplitude of the first peak is larger than the same value at lower water content (yellow). However, as the sample dries, the reflection from the hydrogel-metal interface becomes larger. This observation is attributed to the reduced attenuation of the THz beam propagating through the sample and the increased Fresnel transmission coefficient at the air-hydrogel interface.

Modeling results
As shown in Fig. 4(b), the first step in the diffusion model was to use the semi-infinite assumption and determine α and D 0 by fitting the results of Eq. (5) and the linear region of the weight measurements. Unlike the variable-thickness model, the concentration at the sphere-hydrogel interface is kept constant at C 0 . As a result, the water loss in the semi-infinite model continues past the measurements as shown in Fig. 4(b). The extracted values for α and D 0 were then used as inputs into the variable-thickness diffusion model, shown in Fig. 5(a). Extracted values for the constants, D eff and k, were determined by minimizing the error between M total (t), given by Eq. (13), and the measured water loss per unit area. After implementing both parts of the model, we obtained excellent agreement (R 2 = 0.998) with the measured total hydration loss. Examples of the hydration profiles for the thinnest and the thickest sections within our beam spot size are shown in Fig. 8(a-b), respectively. Each discretized radial section had a different thickness, and therefore resulted in a distinct hydration profile. Most of the water loss occurred at the early time points, which can be seen by the sparse amount of lines that become denser as time progresses.
The THz-TDS traces at various time points, along with their corresponding simulated results, are shown in Fig. 9. A time shift correction factor was applied to the simulated results to account for the change in time-of-arrival after placing the contact lens on the stainless-steel sphere. This correction factor was based on the modeled thickness, d model , and calculated by, where c is the speed of light in a vacuum and n air is the refractive index of ambient air. The reflection from the hydrogel-sphere interface is almost invisible at the initial timepoints because most of the THz signal has been attenuated by the higher water content. At later time  points, when the water content decreases, the reflection from the hydrogel-ball interface begins to appear. Simulated THz signals show very good agreement with the measured data, particularly at later time points.

Discussion
In this study, we aimed to build THz experimental methods and diffusion models to gain a detailed understanding of water loss in a corneal phantom. Our model provides a simple one-dimensional representation of the water diffusion in the axial direction. In more complicated models, the diffusion of water could be considered in both the axial and the perpendicular directions. Future work will include expanding our model to include diffusion in the perpendicular directions and tracking it using spatial mapping of a spherical target, based on the experimental techniques in [16,17]. Furthermore, these results can be generalized to other tissue types or non-biological samples. Translating our findings into real-world clinical applications would require an extension of our model that incorporates the transport phenomena in the corneal endothelium, similar to the model described in [38]. Additionally, the reflected signal from the backside of the in vivo corneal samples will be weaker. The dielectric properties of the corneal stroma may also deviate from effective media models due to the active and passive transport through the endothelial layer. However, if these challenges are met, a model of a living human cornea could provide significant insights into the mechanisms of corneal endothelial transport.

Conclusion
In this study, the diffusion dynamics of a corneal phantom were investigated by simultaneously obtaining gravimetric and THz-TDS measurements. The THz reflection results of the dehydrating sample showed very good agreement with the combination of a two-step hydrogel diffusion model and a stratified composite-media model. The diffusion model in this study utilized the water loss measurements obtained from an analytical balance to determine an axial hydration profile at each time point. The axial hydration profile was then used in the stratified composite-media model to determine an aggregate THz reflection coefficient. The reflection coefficient was convolved with the reference measurement to obtain the simulated THz reflection signal ultimately in the time-domain, which showed very good agreement with the experimental results. Future work would involve the expansion of our one-dimensional diffusion model to three dimensions, and inclusion of the active and passive transport processes in the corneal endothelial layer to achieve a more realistic model for in vivo studies.

Funding
National Institute of General Medical Sciences (R01GM112693); Stony Brook University.