The biomechanical properties of the cornea can be significantly altered by several ocular diseases, such as keratoconus.¹ Hence, assessing corneal biomechanical properties are required for detecting the onset and monitoring progression of the ocular diseases as well as assessing the effects of therapeutic interventions. Previous studies have shown the feasibility of commercially available instruments such as the ocular response analyzer and the CorVis,² and emerging techniques such as Brillouin microscopy³ for estimating corneal biomechanical properties. Although these methods can provide important measurements that can reflect some corneal biomechanical features, quantitative assessment of the corneal viscoelasticity is still a challenge.

Optical coherence elastography (OCE) is an emerging technique with a great potential for noninvasive measurements of the local biomechanical properties of tissues with high spatial and temporal resolutions.⁴ Similar to other elastographic techniques, such as magnetic resonance elastography⁵ and ultrasound elastography,⁶ OCE combines the corresponding imaging modality with mechanical loading. We previously used OCE in different applications, such as soft tissue tumor detection,⁷ cornea,⁸ and cardiac muscle⁹ elasticity estimation. However, quantification of tissue viscoelastic properties from OCE measurements requires the selection of a proper model that can accurately map the parameters of measured elastic waves to, e.g., Young’s modulus.

In this study, we demonstrate, for the first time to the best of our knowledge, application of a modified Rayleigh–Lamb frequency equation (RLFE) to quantify the viscoelastic properties of porcine corneas from OCE measurements. Based on the temporal displacement profiles of a focused air-pulse-induced elastic wave, phase velocities over a range of angular frequencies were calculated and fitted using the RLFE to extract the corneal elasticity and viscosity. A validation study was performed on agar phantoms before the RLFE method was applied to the porcine corneas.

The home-built OCE system is composed of a focused air-pulse delivery system¹⁰ and a phase-stabilized swept source optical coherence tomography (PhS-SSOCT) system.¹¹ The system utilized a broadband swept laser source (HSL2000, Santec, Inc., California) with a central wavelength of $\sim 1310\mathrm{nm}$, bandwidth of $\sim 150\mathrm{nm}$, scan rate of 30 kHz, and experimentally measured phase stability of $\sim 40\mathrm{nm}$ during cornea experiments. The focused air-pulse delivery system was comprised of a controller with a signal input for synchronization, a solenoid-controlled air gate, and an air-pulse port with a flat edge and diameter of $\sim 150\mu \mathrm{m}$. The system is capable of delivering a short duration focused air-pulse ($\le 1\mathrm{ms}$) with a Gaussian profile to induce a small amplitude deformation (order of $\mu \mathrm{m}$) within the cornea tissue.

The air-pulse was positioned $\sim 2\mathrm{mm}$ away from the apex of the cornea to ensure that the elastic wave propagated across the apex and $\sim 250\mu \mathrm{m}$ away from the surface of the cornea. Imaging of the elastic wave propagation was performed by synchronizing 501 M-mode images at successive positions across $\sim 6.1\mathrm{mm}$ tissue surface.¹² The phase velocity of the elastic wave was calculated by performing a fast Fourier transform on the elastic wave temporal displacement profiles at each measurement position and imaged depth, producing a phase shift, $\mathrm{\Delta }\theta$, at a specific angular frequency, $\omega$, and depth. For each depth, the distance, $\mathrm{\Delta }r$, of the elastic wave propagation was calculated, taking the curvature of the cornea into account.¹³ The phase velocity of the elastic wave at the angular frequency, $\omega$, is expressed as: $c_p=\omega \mathrm{\Delta }r/\mathrm{\Delta }\theta$, which was obtained by least-squares linear fitting. The resulting phase velocities at each frequency were then averaged over the imaged depths in the sample.

Validation experiments were first performed on 2% agar phantoms ($w/w$, $n=3$) with a cylindrical shape: diameter $D=50\mathrm{mm}$ and height $H=11\mathrm{mm}$. The phantoms were assumed to be homogeneous and isotropic with a negligible viscosity.¹⁴ The agar phantoms were placed on a ring support such that the central part of the bottom surface was free of the contact. The excitation position was at the center of the top surface of the phantom and all OCE measurements were taken in the central region. Assuming free boundary conditions, the Rayleigh-Lamb characteristic equation for an anti-symmetric wave has the form¹⁵

In Eq. (1), $h$ is half the thickness of the phantom ($h=(1/2)H=5.5\mathrm{mm}$), and $k=\omega /c_p$ is the wave number at angular frequency $\omega$. For Eq. (1), the following relationships can be expressed

By solving Eq. (1), we obtained the numerical relationship between $\omega$ and $c_p$, which was then compared with the OCE experimental data to quantify the Young’s modulus. Elasticity values obtained from OCE measurements were then validated on agar phantoms (same concentration 2% $w/w$, $n=3$) using uniaxial mechanical compressional tests (Model 5943, Instron Corp., Massachusetts).

Figure 1 shows the agar phantom elasticity assessment by RLFE and validation by mechanical testing. Figure 1(a) presents a typical temporal displacement profile from the agar phantoms as measured by OCE. The duration of displacement [$\sim 2\mathrm{ms}$ as marked in Fig. 1(a)] was assumed as a half wavelength for determining the lower limit of frequencies to be used in the RLFE. For the agar phantoms, angular frequencies lower than $1500\mathrm{rad}/\mathrm{s}$ were due to low frequency noise and were not used for calculations. Figure 1(b) depicts the phase velocities calculated by the RLFE and as measured by OCE. From the RLFE and phase velocity data, the Young’s modulus of the 2% agar phantoms was estimated at 160 kPa. Deviations of the OCE data from the RLFE may be due to the fact that the Lamb wave model is strictly based on a thin layer geometry, to which the phantoms do not strictly adhere. Figure 1(c) shows the Young’s moduli of the agar phantoms as estimated by the RLFE and measured by uniaxial mechanical compression testing. The RLFE results agreed well with the mechanical testing measurement, which demonstrates the feasibility of this method for quantitative elasticity assessment.

Fig. 1

(a) A typical temporal displacement profile of a single point on a 2% agar phantom. (b) Optical coherence elastography (OCE) measurements of the elastic wave dispersion for 2% agar phantom fitted with Rayleigh–Lamb frequency Eq. (1). (c) Quantitative results of the Young’s modulus of 2% agar phantom assessed by Rayleigh–Lamb frequency equation (RLFE) and mechanical testing.

Next, experiments were performed on porcine corneas ($n=3$) in whole eyeball configuration with artificially controlled intraocular pressure, $\mathrm{IOP}=20\mathrm{mmHg}$. The average thickness of the corneas was $H=0.9\mathrm{mm}$ (de-epithelialized). The RLFE used for agar phantoms was modified because the cornea has two distinct features. First, the viscosity of the cornea must be taken into account. We assumed a Kelvin–Voigt viscoelastic model, which led to a complex dynamic shear modulus of the form: ${\mu }_D=\mu +i\eta \omega$, where $\eta$ is the shear viscosity and $i$ is the imaginary unit. Second, the fluid effect caused by the aqueous humor must be considered. Although the cornea has a convex shape, it was treated as a flat disk for simplification (the effect of the corneal curvature on the estimated values of $E$ and $\eta$ is currently the topic of our investigation). In addition, the boundary conditions were assumed as zero stress on the corneal anterior surface and equal stresses and vertical displacements between the corneal tissue and aqueous humor on the posterior surface. For this case, the Lamb wave cannot be considered as symmetric or anti-symmetric. Based on these assumptions, the RLFE for the cornea was then modified as

Here, $d=(1/2)H=0.45\mathrm{mm}$ is half of the porcine corneal thickness, ${\rho }_F=1000\mathrm{kg}/{\mathrm{m}}^3$ is the fluid density, $\rho =1062\mathrm{kg}/{\mathrm{m}}^3$ is the cornea density,¹⁶ and $c_F=1500\mathrm{m}/\mathrm{s}$ is the speed of sound in water.

Figure 2 presents the corneal viscoelasticity assessment by the modified RLFE from OCE measurements. A typical temporal displacement profile as measured by OCE from a porcine cornea is shown in Fig. 2(a). Similar to the agar phantoms, the displacement duration of $\sim 2.5\mathrm{ms}$ marked in Fig. 2(a) was used to determine the lower frequency limit of $\sim 1250\mathrm{rad}/\mathrm{s}$ to be adopted in the modified RLFE. Figure 2(b) shows the estimated Young’s modulus ($E=60\mathrm{kPa}$) and shear viscosity ($\eta =0.33\mathrm{Pa}·\mathrm{s}$) of three porcine corneas at an artificially controlled IOP of 20 mmHg, which were estimated by fitting the phase velocities from the OCE measurements with the modified RLFE numerical results. The estimated Young’s modulus is of the same order as measured by atomic force microscopy,¹⁷ which suggests that the proposed method could provide accurate quantitative viscoelastic characterization of the cornea under the proper assumptions.

Fig. 2

(a) A typical temporal displacement profile from a porcine cornea. (b) Rayleigh–Lamb wave dispersion for porcine corneas ($n=3$) fitted with modified Rayleigh–Lamb frequency Eq. (4), with estimated Young’s modulus ($E=60\mathrm{kPa}$) and shear viscosity ($\eta =0.33\mathrm{Pa}·\mathrm{s}$).

To summarize, in this study, the viscoelasticity of porcine corneas was quantitatively assessed by fitting the phase velocities of an elastic wave measured by PhS-SSOCE to a modified Rayleigh-Lamb wave model, which accounted for the fluid effect on the corneal posterior surface. RLFE validation was conducted on 2% agar phantoms by uniaxial mechanical compression testing, which verified the accuracy of this method. The Young’s modulus of porcine corneas at 20 mmHg IOP was estimated to be $\sim 60\mathrm{kPa}$ and the shear viscosity as $\sim 0.33\mathrm{Pa}·\mathrm{s}$. The results indicate that the combination of PhS-SSOCE and the modified Rayleigh-Lamb characteristic equation may be potentially useful for assessing corneal viscoelasticity in vivo.