Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

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

Open Access Open Access

Abstract

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.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. 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,1216]. 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.

2. Materials and methods

2.1 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.

2.2 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.

 figure: Fig. 1.

Fig. 1. (a) The experimental setup for taking simultaneous weight and THz reflection measurements from the corneal phantom is illustrated. The THz beam was focused at the apex of the stainless-steel sphere. The reference THz-TDS signal from the apex of the metallic sphere is shown in (b) time domain and (c) frequency domain. (d) The water loss of the contact lens during the dehydration experiment is shown as a function of time

Download Full Size | PDF

2.3 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, $\alpha $, 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.

 figure: Fig. 2.

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.

Download Full Size | PDF

2.3.1. Assumptions

We assumed that the hydrogel, methafilicon B, had a density of 1.07 $\textrm{g/c}{\textrm{m}^{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.

2.3.2. 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 $a{x^4} + \; b{x^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.

 figure: Fig. 3.

Fig. 3. (a) The geometry of the lens is shown with a color bar to represent the thickness profile. (b) The cross-section of the contact lens is shown with discretized sections.

Download Full Size | PDF

2.3.3. 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 $\textrm{g/}{\textrm{m}^{3}}$), at each measured time point is determined by,

$$C = \; \frac{{Grams\; of\; Water}}{{Volume\; of\; Hydrogel}} = \; W \cdot {\rho _{Hydrogel}},$$
where ${\rho _{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
$$W = \; \frac{{mass(t )- mas{s_{dry}}}}{{mas{s_{hydrated}}}},$$
where $mas{s_{hydrated}}$ is the mass of a fully hydrated sample and $mas{s_{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 ${\pm} 3\%$ difference between their measured water content and manufacturer value [24,25].

 figure: Fig. 4.

Fig. 4. (a) The diffusion of water (represented by blue arrow) is illustrated outwards from the sample and into the environment. The sample is assumed to have an infinite thickness and $x$ is defined as the axial direction. (b) Eq. (5) was fitted to the linear region of the gravimetric data to extract $\alpha = 3.16 \times {10^{-7}}\,\textrm{m/s}$ and ${D_0} = 1.99 \times {10^{-10}}\,{\textrm{m}^{2}}/\textrm{s}$.

Download Full Size | PDF

Fick’s 2nd law of diffusion governs the transport of water in our model, given by,

$$\frac{{\partial C}}{{\partial t}} = \frac{\partial }{{\partial x}}\left( {{D_0}\frac{{\partial C({x,t} )}}{{\partial x}}} \right),$$
where ${D_0}$ is a constant diffusion coefficient in units of ${\textrm{m}^{2}}\textrm{/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 = \; - \infty $ 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,

$$- {D_0}\frac{{\partial C({x,t} )}}{{\partial x}} = \; \alpha ({{C_{eq}} - {C_s}} )\; \; \textrm{at}\; x = 0\; .$$
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 )$,
$${M_{si}}(t )= \left( {\frac{{{C_{eq}} - {C_0}}}{p}} \right)\left\{ {\exp ({{p^2}{D_0}t} )\textrm{erfc}\left( {p\sqrt {{D_0}t} } \right) - 1 + \frac{2}{{\sqrt \pi }}p\sqrt {{D_0}t} } \right\}\; ,$$
where $p = \alpha /{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.

2.3.4. 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).

 figure: Fig. 5.

Fig. 5. (a) The diffusion of water (represented by blue arrow) is illustrated outwards from the sample and into the environment. The changing thickness is represented by the black arrows. It is important to note that $x = 0$ is now defined at the bottom of the sample rather than the surface, highlighted by a red arrow. (b) Eq. (13) is fitted to the gravimetric data using $\alpha$ from the semi-infinite model. The resulting ${D_{eff}}$ and $k$ values are $3.772 \times {10^{-10}}\,{\textrm{m}^{2}}/\textrm{s}$ and 0.18 respectively.

Download Full Size | PDF

The variable thickness model is also governed by Fick’s 2nd law of diffusion. The diffusion coefficient, $D(C )$, is more realistically represented as a function of concentration,

$$\frac{{\partial C}}{{\partial t}} = \frac{\partial }{{\partial x}}\left( {D(C )\frac{{\partial C}}{{\partial x}}} \right)\; .$$
At the hydrogel-sphere interface, we assumed that there was no diffusion or flux because the evaporative flux only occurs at the surface, therefore,
$$\frac{{\partial C({x,t} )}}{{\partial x}} = 0\; \; \textrm{at}\; x = 0.$$
The surface boundary condition accounts for the thickness of the sample over time, $h(t )$, and is given by,
$$- D(C )\frac{{\partial C({x,t} )}}{{\partial x}} - C\frac{{\partial h(t )}}{{\partial t}} = \alpha ({{C_s} - {C_{eq}}} ),\; $$
where $\alpha $ 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 $\textrm{g/}{\textrm{m}^{3}}$, in the hydrogel, $D(C )$, is represented as an exponential function of water content [21,27],
$$D(C )= {D_{eff}}\exp \left( { - \frac{k}{{C^{\prime}}}} \right),$$
where ${D_{eff}}$ is the effective diffusion coefficient of the unbound water in the hydrogel, $C^{\prime}$ 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,
$$ h(t)=l\left(\frac{1}{S}+\left(\frac{S-1}{S}\right) \frac{\bar{C}(t)}{C_{0}}\right) $$
where l is the starting thickness of the sample and $\bar{C}(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,
$$\frac{{C_n^{t + 1} - C_n^t}}{{{\Delta }t}} = D(C )_n^t\frac{{C_{n + 1}^t - 2C_n^t + \; C_{n - 1}^t}}{{{{({\Delta}x)}^2}}} + \frac{{\partial D(C )}}{{\partial C}}|_n^t{(\frac{{C_{n + 1}^t - C_{n - 1}^t}}{{2{\Delta }x}})^2}\; .$$
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 ${\Delta }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,
$${M_{vt}}(t )= l{C_0} - \; \mathop \smallint \nolimits_0^{h(t )} C({x,t} )dx\; .$$
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_{v{t_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_{v{t_i}}}(t )$ for all sections, given by,
$${M_{total}}(t )= \mathop \sum \limits_{i = 1}^N {M_{v{t_i}}}(t ){A_i},$$
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.

2.4 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.

 figure: Fig. 6.

Fig. 6. The discrete dielectric slabs are shown with the darker slabs representing increasing water content. The complex dielectric constant for each slab can be determined using Bruggeman effective media theory, Eq. (15), where each slab is a binary composite medium of water and hydrogel material. An illustration of the multiple THz reflections is shown by the red arrows and corresponding electric fields, ${E_{r_i}}$. N is the number of layers, here set to 11.

Download Full Size | PDF

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,

$${\tilde{\varepsilon }_w}(\omega )= {\varepsilon _\infty } + \frac{{({\varepsilon _s} - {\varepsilon _2})}}{{1 + i\omega {{\tau }_1}}} + \frac{{({\varepsilon _2} - {\varepsilon _\infty })}}{{1 + i\omega {{\tau }_2}}}\; ,$$
where ${\varepsilon _s}\; $ is the static dielectric constant, ${\varepsilon _2}$ is an intermediate dielectric value, ${\varepsilon _\infty }$ is the high frequency dielectric constant, and ${\tau _1}$ and ${\tau _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,

$${\rho _{w,i}}\left( {\frac{{{{\tilde{\varepsilon }}_i} - {{\tilde{\varepsilon }}_w}}}{{{{\tilde{\varepsilon }}_w} + \; 2{{\tilde{\varepsilon }}_i}}}} \right) + ({1 - {\rho_{w,i}}} )\left( {\frac{{{{\tilde{\varepsilon }}_i} - {{\tilde{\varepsilon }}_h}}}{{{{\tilde{\varepsilon }}_h} + 2{{\tilde{\varepsilon }}_i}}}} \right) = 0\; \forall i \in [{i,N} ],$$
where ${\rho _{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, ${\tilde{\varepsilon }_w}$ is the dielectric constant of water, ${\tilde{\varepsilon }_i}$ is the effective dielectric constant and ${\tilde{\varepsilon }_h}$ is the dielectric constant of the polymer in the hydrogel [5]. The values of ${\rho _{w,1}}$, at the hydrogel-air interface, and ${\rho _{w,N}}$, at the hydrogel-metal interface, are calculated by dividing $C(x,t)$ by the density of water. At $t = 0$, ${\rho _{w,1}} = {\rho _{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 ${\rho _{w,1}}$ decreasing faster than ${\rho _{w,N}}$. The value of ${\tilde{\varepsilon }_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, ${\delta _i}$, by

$${r_i} = \; \frac{{{{\tilde{n}}_{i - 1}} - {{\tilde{n}}_i}}}{{{{\tilde{n}}_{i - 1}} + \tilde{n}{\; _i}}}\forall i \in [{1,N + 1} ]$$
and
$${\delta _i} = \frac{{2\pi }}{\lambda }\; d{\tilde{n}_i}\; \forall i \in [{1,N} ]\; ,\; $$
where d is the uniform layer thickness, $\lambda $ is the wavelength and ${\tilde{n}_i} = \sqrt {{{\tilde{\varepsilon }}_i}} $. The total reflected signal for each layer can be calculated by solving the multiple internal reflections in a geometric series [5,18,3537]. The aggregate reflection coefficient, ${{\Gamma }_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,
$${\Gamma_i} = \frac{{{r_i} + {{\Gamma}_{i + 1}}{e^{ - 2j{\delta _i}}}}}{{1 + {r_i}{{\Gamma}_{i + 1}}{e^{ - 2j{\delta _i}}}}},\; \; {{\Gamma}_{N + 1}} = \; \; {r_{N + 1}}\; \forall i \in [{1,N + 1} ]\; .\; \; $$
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, $\tilde{n}$. 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,
$${E_{model}}(t )= {{{\cal F}}^{ - 1}}\{{{{\Gamma}_1} \cdot {E_{ref}}(f )} \}.$$

3. Results

3.1 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.

 figure: Fig. 7.

Fig. 7. The THz time-domain traces are shown over the course of the experiment. Each trace represents a single time point and corresponds to a water content shown in the color bar.

Download Full Size | PDF

3.2 Modeling results

As shown in Fig. 4(b), the first step in the diffusion model was to use the semi-infinite assumption and determine $\alpha $ 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 $\alpha $ 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.

 figure: Fig. 8.

Fig. 8. The hydration profile is shown as a function of axial location in the hydrogel for the (a) thinnest and (b) thickest sections within our 2.7 mm diameter beam spot at the focus. Time is represented by the color bar.

Download Full Size | PDF

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,

$$correction\; factor = \frac{{2{d_{model\; }}{n_{air}}}}{c}\; ,\; $$
where c is the speed of light in a vacuum and ${n_{air}}$ is the refractive index of ambient air.

 figure: Fig. 9.

Fig. 9. The simulated THz time-domain traces (black) and the corresponding measured THz signal (red) is shown at multiple time points.

Download Full Size | PDF

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.

4. 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.

5. 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.

Acknowledgments

The authors would like to thank the Department of Ophthalmology at Stony Brook University for helping us acquire the necessary samples.

Disclosures

The authors declare that there are no conflicts of interest related to this work

References

1. J. Fischbarg and D. M. Maurice, “An update on corneal hydration control,” Exp. Eye Res. 78(3), 537–541 (2004). [CrossRef]  

2. Y. S. Rabinowitz, “Keratoconus,” Surv. Ophthalmol. 42(4), 297–319 (1998). [CrossRef]  

3. A. P. Adamis, V. Filatov, B. J. Tripathi, and R. C. Tripathi, “Fuchs’ endothelial dystrophy of the cornea,” Surv. Ophthalmol. 38(2), 149–168 (1993). [CrossRef]  

4. J. Ytteborg and C. H. Dohlman, “Corneal edema and intraocular pressure. II. Clinical results,” Arch Ophthalmol. 74(4), 477–484 (1965). [CrossRef]  

5. Z. D. Taylor, J. Garritano, S. Sung, N. Bajwa, D. B. Bennett, B. Nowroozi, P. Tewari, J. Sayre, J.-P. Hubschman, S. Deng, E. R. Brown, and W. S. Grundfest, “THz and mm-wave sensing of corneal tissue water content: electromagnetic modeling and analysis,” IEEE Trans. Terahertz Sci. Technol. 5(2), 170–183 (2015). [CrossRef]  

6. Z. D. Taylor, R. S. Singh, D. B. Bennett, P. Tewari, C. P. Kealey, N. Bajwa, M. O. Culjat, A. Stojadinovic, H. Lee, J.-P. Hubschman, E. R. Brown, and W. S. Grundfest, “THz medical imaging: in vivo hydration sensing,” IEEE Trans. Terahertz Sci. Technol. 1(1), 201–219 (2011). [CrossRef]  

7. Z. Taylor, R. Singh, M. Culjat, J. Suen, W. Grundfest, and E. Brown, “THz imaging based on water-concentration contrast - art. no. 69490D,” Proc. SPIE 6949, 69490D (2008). [CrossRef]  

8. Z. D. Taylor, J. Garritano, S. Sung, N. Bajwa, D. B. Bennett, B. Nowroozi, P. Tewari, J. W. Sayre, J. Hubschman, S. X. Deng, E. R. Brown, and W. S. Grundfest, “THz and mm-wave sensing of corneal tissue water content: in vivo sensing and imaging results,” IEEE Trans. Terahertz Sci. Technol. 5(2), 184–196 (2015). [CrossRef]  

9. L. Guo, X. Wang, P. Han, W. Sun, S. Feng, J. Ye, and Y. Zhang, “Observation of dehydration dynamics in biological tissues with terahertz digital holography [Invited],” Appl. Opt. 56(13), F173–F178 (2017). [CrossRef]  

10. Q. Sun, R. I. Stantchev, J. Wang, E. P. Parrott, A. Cottenden, T. W. Chiu, A. T. Ahuja, and E. Pickwell-MacPherson, “In vivo estimation of water diffusivity in occluded human skin using terahertz reflection spectroscopy,” J. Biophotonics 12(2), e201800145 (2019). [CrossRef]  

11. Q. Sun, E. P. Parrott, Y. He, and E. Pickwell-MacPherson, “In vivo THz imaging of human skin: Accounting for occlusion effects,” J. Biophotonics 11(2), e201700111 (2018). [CrossRef]  

12. D. Bennett, Z. Taylor, P. Tewari, S. Sung, A. Maccabi, R. Singh, M. Culjat, W. Grundfest, J. P. Hubschman, and E. Brown, “Assessment of corneal hydration sensing in the terahertz band: in vivo results at 100 GHz,” J. Biomed. Opt. 17(9), 0970081 (2012). [CrossRef]  

13. S. Sung, S. Chantra, N. Bajwa, R. McCurdy, G. Kerezyte, J. Garritano, J.-P. Hubschman, W. Grundfest, S. X. Deng, and Z. Taylor, “Direct measurement of corneal tissue water content by reflection imaging at Terahertz Frequencies,” Invest. Ophthalmol. Visual Sci. 56, 1644 (2015).

14. D. B. Bennett, Z. D. Taylor, P. Tewari, R. S. Singh, M. O. Culjat, W. S. Grundfest, D. J. Sassoon, R. D. Johnson, J. P. Hubschman, and E. R. Brown, “Terahertz sensing in corneal tissues,” J. Biomed. Opt. 16(5), 057003 (2011). [CrossRef]  

15. R. S. Singh, P. Tewari, J. L. Bourges, J. P. Hubschman, D. B. Bennett, Z. D. Taylor, H. Lee, E. R. Brown, W. S. Grundfest, and M. O. Culjat, “Terahertz sensing of corneal hydration,” Conf. Proc. IEEE Eng. Med. Biol. Soc. 2010, 3021–3024 (2010).

16. S. Sung, S. Selvin, N. Bajwa, S. Chantra, B. Nowroozi, J. Garritano, J. Goell, A. Li, S. X. Deng, E. Brown, W. S. Grundfest, and Z. D. Taylor, “THz imaging system for in vivo human cornea,” IEEE Trans. Terahertz Sci. Technol. 8(1), 27–37 (2018). [CrossRef]  

17. S. Sung, S. Dabironezare, N. Llombart, S. Selvin, N. Bajwa, S. Chantra, B. Nowroozi, J. Garritano, J. Goell, A. Li, S. X. Deng, E. Brown, W. S. Grundfest, and Z. D. Taylor, “Optical system design for noncontact, normal incidence, THz imaging of in vivo human cornea,” IEEE Trans. Terahertz Sci. Technol. 8(1), 1–12 (2018). [CrossRef]  

18. D. B. Bennett, W. Li, Z. D. Taylor, W. S. Grundfest, and E. R. Brown, “Stratified media model for terahertz reflectometry of the skin,” IEEE Sens. J. 11(5), 1253–1262 (2011). [CrossRef]  

19. Z. Gu and P. Alexandridis, “Drying of poloxamer hydrogel films,” J. Pharm. Sci. 93(6), 1454–1470 (2004). [CrossRef]  

20. N. J. Bauer, J. P. Wicksted, F. H. Jongsma, W. F. March, F. Hendrikse, and M. Motamedi, “Noninvasive assessment of the hydration gradient across the cornea using confocal Raman spectroscopy,” Invest. Ophthalmol. Vis. Sci. 39, 831–835 (1998).

21. H. P. Blandin, J. C. David, J. M. Vergnaud, J. P. Illien, and M. Malizewicz, “Modelling of drying of coatings: Effect of the thickness, temperature and concentration of solvent,” Prof. Org. Coat. 15(2), 163–172 (1987). [CrossRef]  

22. Contamac, “Contaflex FDA contact lens”, retrieved Feb. 07, 2019.

23. S. J. Vincent, D. Alonso-Caneiro, H. Kricancic, and M. J. Collins, “Scleral contact lens thickness profiles: The relationship between average and centre lens thickness,” Cont. Lens Anterior Eye 42(1), 55–62 (2019). [CrossRef]  

24. N. Efron, P. B. Morgan, I. D. Cameron, N. A. Brennan, and M. Goodwin, “Oxygen permeability and water content of silicone hydrogel contact lens materials,” Optom. Vis. Sci. 84(4), E328–E337 (2007). [CrossRef]  

25. P. McConville and J. M. Pope, “A comparison of water binding and mobility in contact lens hydrogels from NMR measurements of the water self-diffusion coefficient,” Polymer 41(26), 9081–9088 (2000). [CrossRef]  

26. J. Crank, The mathematics of diffusion (Oxford university press, 1979).

27. L. Ion and J. M. Vergnaud, “Process of drying a polymeric paint by diffusion-evaporation and shrinkage. Determination of the concentration-dependent diffusivity,” Polym. Test. 14(5), 479–487 (1995). [CrossRef]  

28. W. A. Douthwaite, “Contact Lens Optics and Lens Design,” Contact Lens Optics and Lens Design (2006).

29. T. Ahmed, I. V. Belova, and G. E. Murch, “Finite difference solution of the diffusion equation and calculation of the interdiffusion coefficient using the Sauer-Freise and Hall Methods in binary systems,” Procedia Eng. 105, 570–575 (2015). [CrossRef]  

30. C. J. Budd and G. J. Collins, “An invariant moving mesh scheme for the nonlinear diffusion equation,” Appl. Numer. Math. 26(1-2), 23–39 (1998). [CrossRef]  

31. M. H. Arbab, T. C. Dickey, D. P. Winebrenner, A. Chen, M. B. Klein, and P. D. Mourad, “Terahertz reflectometry of burn wounds in a rat model,” Biomed. Opt. Express 2(8), 2339–2347 (2011). [CrossRef]  

32. J. T. Kindt and C. A. Schmuttenmaer, “Far-Infrared dielectric properties of polar liquids probed by femtosecond terahertz pulse spectroscopy,” J. Phys. Chem. 100(24), 10373–10379 (1996). [CrossRef]  

33. E. Pickwell-MacPherson, B. Cole, A. Fitzgerald, V. Wallace, and M. Pepper, “Simulation of terahertz pulse propagation in biological systems,” Appl. Phys. Lett. 84(12), 2190–2192 (2004). [CrossRef]  

34. Y. S. Jin, G. J. Kim, and S. G. Jeon, “Terahertz dielectric properties of polymers,” J. Korean. Phys. Soc. 49, 513–517 (2006).

35. G. P. Kniffin and L. M. Zurk, “Model-based material parameter estimation for terahertz reflection spectroscopy,” IEEE Trans. Terahertz Sci. Technol. 2(2), 231–241 (2012). [CrossRef]  

36. P. Lecaruyer, E. Maillart, M. Canva, and J. Rolland, “Generalization of the Rouard method to an absorbing thin-film stack and application to surface plasmon resonance,” Appl. Opt. 45(33), 8419–8423 (2006). [CrossRef]  

37. E. W. Max Born, Principles of Optics (Sixth Edition) (Pergamon, 1980).

38. X. Cheng, S. J. Petsche, and P. M. Pinsky, “A structural model for the in vivo human cornea including collagen-swelling interaction,” J. R. Soc. Interface 12(109), 20150241 (2015). [CrossRef]  

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (9)

Fig. 1.
Fig. 1. (a) The experimental setup for taking simultaneous weight and THz reflection measurements from the corneal phantom is illustrated. The THz beam was focused at the apex of the stainless-steel sphere. The reference THz-TDS signal from the apex of the metallic sphere is shown in (b) time domain and (c) frequency domain. (d) The water loss of the contact lens during the dehydration experiment is shown as a function of time
Fig. 2.
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.
Fig. 3.
Fig. 3. (a) The geometry of the lens is shown with a color bar to represent the thickness profile. (b) The cross-section of the contact lens is shown with discretized sections.
Fig. 4.
Fig. 4. (a) The diffusion of water (represented by blue arrow) is illustrated outwards from the sample and into the environment. The sample is assumed to have an infinite thickness and $x$ is defined as the axial direction. (b) Eq. (5) was fitted to the linear region of the gravimetric data to extract $\alpha = 3.16 \times {10^{-7}}\,\textrm{m/s}$ and ${D_0} = 1.99 \times {10^{-10}}\,{\textrm{m}^{2}}/\textrm{s}$.
Fig. 5.
Fig. 5. (a) The diffusion of water (represented by blue arrow) is illustrated outwards from the sample and into the environment. The changing thickness is represented by the black arrows. It is important to note that $x = 0$ is now defined at the bottom of the sample rather than the surface, highlighted by a red arrow. (b) Eq. (13) is fitted to the gravimetric data using $\alpha$ from the semi-infinite model. The resulting ${D_{eff}}$ and $k$ values are $3.772 \times {10^{-10}}\,{\textrm{m}^{2}}/\textrm{s}$ and 0.18 respectively.
Fig. 6.
Fig. 6. The discrete dielectric slabs are shown with the darker slabs representing increasing water content. The complex dielectric constant for each slab can be determined using Bruggeman effective media theory, Eq. (15), where each slab is a binary composite medium of water and hydrogel material. An illustration of the multiple THz reflections is shown by the red arrows and corresponding electric fields, ${E_{r_i}}$. N is the number of layers, here set to 11.
Fig. 7.
Fig. 7. The THz time-domain traces are shown over the course of the experiment. Each trace represents a single time point and corresponds to a water content shown in the color bar.
Fig. 8.
Fig. 8. The hydration profile is shown as a function of axial location in the hydrogel for the (a) thinnest and (b) thickest sections within our 2.7 mm diameter beam spot at the focus. Time is represented by the color bar.
Fig. 9.
Fig. 9. The simulated THz time-domain traces (black) and the corresponding measured THz signal (red) is shown at multiple time points.

Equations (20)

Equations on this page are rendered with MathJax. Learn more.

C = G r a m s o f W a t e r V o l u m e o f H y d r o g e l = W ρ H y d r o g e l ,
W = m a s s ( t ) m a s s d r y m a s s h y d r a t e d ,
C t = x ( D 0 C ( x , t ) x ) ,
D 0 C ( x , t ) x = α ( C e q C s ) at x = 0 .
M s i ( t ) = ( C e q C 0 p ) { exp ( p 2 D 0 t ) erfc ( p D 0 t ) 1 + 2 π p D 0 t } ,
C t = x ( D ( C ) C x ) .
C ( x , t ) x = 0 at x = 0.
D ( C ) C ( x , t ) x C h ( t ) t = α ( C s C e q ) ,
D ( C ) = D e f f exp ( k C ) ,
h ( t ) = l ( 1 S + ( S 1 S ) C ¯ ( t ) C 0 )
C n t + 1 C n t Δ t = D ( C ) n t C n + 1 t 2 C n t + C n 1 t ( Δ x ) 2 + D ( C ) C | n t ( C n + 1 t C n 1 t 2 Δ x ) 2 .
M v t ( t ) = l C 0 0 h ( t ) C ( x , t ) d x .
M t o t a l ( t ) = i = 1 N M v t i ( t ) A i ,
ε ~ w ( ω ) = ε + ( ε s ε 2 ) 1 + i ω τ 1 + ( ε 2 ε ) 1 + i ω τ 2 ,
ρ w , i ( ε ~ i ε ~ w ε ~ w + 2 ε ~ i ) + ( 1 ρ w , i ) ( ε ~ i ε ~ h ε ~ h + 2 ε ~ i ) = 0 i [ i , N ] ,
r i = n ~ i 1 n ~ i n ~ i 1 + n ~ i i [ 1 , N + 1 ]
δ i = 2 π λ d n ~ i i [ 1 , N ] ,
Γ i = r i + Γ i + 1 e 2 j δ i 1 + r i Γ i + 1 e 2 j δ i , Γ N + 1 = r N + 1 i [ 1 , N + 1 ] .
E m o d e l ( t ) = F 1 { Γ 1 E r e f ( f ) } .
c o r r e c t i o n f a c t o r = 2 d m o d e l n a i r c ,
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.