Open Access
Add to BibSonomy
Abstract
Corneal collagen crosslinking (CXL) is commonly used to prevent or treat keratoconus. Although changes in corneal stiffness induced by CXL surgery can be monitored with non-contact dynamic optical coherence elastography (OCE) by tracking mechanical wave propagation, depth dependent changes are still unclear if the cornea is not crosslinked through the whole depth. Here, phase-decorrelation measurements on optical coherence tomography (OCT) structural images are combined with acoustic micro-tapping (AµT) OCE to explore possible reconstruction of depth-dependent stiffness within crosslinked corneas in an ex vivo human cornea sample. Experimental OCT images are analyzed to define the penetration depth of CXL into the cornea. In a representative ex vivo human cornea sample, crosslinking depth varied from ∼100 µm in the periphery to ∼150 µm in the cornea center and exhibited a sharp in-depth transition between crosslinked and untreated areas. This information was used in an analytical two-layer guided wave propagation model to quantify the stiffness of the treated layer. We also discuss how the elastic moduli of partially CXL-treated cornea layers reflect the effective engineering stiffness of the entire cornea to properly quantify corneal deformation.
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
The cornea is the primary optical element focusing light onto the retina. It contains multiple layers, including epithelium and stroma (Figs. 1(a), (b)). The first acts as a barrier against the external environment and the latter maintains stiffness, transparency and focusing power [1,2]. The microstructure of the stroma is composed of collagen fibrils, arranged in lamellae, lying within a protein rich, hydrated proteoglycan mesh [3,4] (Figs. 1(b), (c)).
Fig. 1. Illustration of the anterior section of the human eye (a), the change in corneal geometry induced by CXL (b), and three-dimensional organization of corneal lamellae (slow-axis NITI symmetry, i.e., symmetry orthogonal to the fibers, or here across lamellae) within the stroma (c). In- and out-of-plane loads are illustrated in (c).
Corneal diseases (such as keratoconus (KC)) and surgical complications from refractive surgeries (such as Laser-Assisted In Situ Keratomileusis (LASIK)) may deform the cornea (ectasia) and alter vision. The prevalence of KC in the general population is estimated to be 1.38 per 1000 [5], and nearly 1 million refractive surgeries are performed each year in the USA. Despite their overall success, however, suboptimal visual outcomes and post-refractive corneal decompensation cannot always be predicted for an individual patient.
Corneal collagen crosslinking is a minimally invasive procedure that can potentially slow the progression of corneal ectasia [6–9]. Ultraviolet (UV) light modifies the microstructure of the cornea soaked in riboflavin and forms additional chemical bonds between collagen fibers in the stroma [10] (Fig. 1(b)). Post-treatment corneas become stiffer and more resistant to enzymatic digestion [11–13]. Although corneal topography (curvature) and thickness maps can be obtained preoperatively, and refractive corrections can be estimated, there is an unmet need to predict corneal decompensation from interventions such as LASIK and CXL therapies. Unfortunately, surgical planning cannot be customized and outcomes (e.g., postoperative corneal ectasia risks) cannot be accurately predicted without quantitatively mapping corneal elasticity. Thus, methods to quantitatively reconstruct corneal elastic moduli are needed.
Ocular response analyzer (ORA - Reichert Technologies) and Dynamic Scheimpflug Analyzer (DSA - Corvis ST – Oculus Opitkgerate GmbH) are the state-of-the-art in clinical measurements of corneal mechanics. They estimate stiffness as the pressure at inward applanation divided by corneal displacement [14–16]. However, measurements induce large corneal deformations that are often clinically unacceptable, require a non-trivial IOP correction in simulations [17] and assume a simple isotropic mechanical model leading to high variability with experimental conditions. Results obtained with the Corvis ST on KC may be contradictory, and some even show no significant change in corneal stiffness pre- and post-CXL surgery [18,19]. In addition, the result is averaged over the entire cornea with no spatial resolution, and the reconstruction is questionable if corneal thickness varies.
Dynamic elastography is a promising tool to probe soft tissue biomechanics. Shear waves can be launched using direct contact excitation [20–23] or radiation force-based techniques [24,25]. By tracking shear wave propagation, either using magnetic resonance imaging (MRI) [26,27] ultrasound [20–22,28] or dynamic phase-sensitive OCT [29–31], one can infer, with an appropriate mechanical description, the linear [32–35], or non-linear [36,37] stiffness moduli of tissue. Optical coherence elastography (OCE) is particularly suited to probe corneal biomechanics non-invasively in a clinical environment [24,30,38–40], as it can be combined with non-contact excitation techniques (for example, using an air-puff or non-contact acoustic micro-tapping (AµT [24])).
Because the cornea is thin and bounded between air and aqueous humor, wave propagation within it is guided, leading to strong geometric dispersion [33,41]. As such, the common approach associating the Rayleigh surface wave group velocity to stiffness [24,30,38–40,42–45] is not appropriate. In addition, the stroma contains collagen lamellae running in-plane across its width. Lamellae make up approximately 90% of tissue thickness and account for most of the cornea’s mechanical structure. They are stacked vertically in approximately 200-500 separate planes [46,47], suggesting an anisotropic mechanical behavior with very different responses to in-plane versus out-of-plane loads (Fig. 1(c)).
To account for this specific architecture we introduced a model of a nearly-incompressible transverse isotropic (NITI) medium [33], in which corneal stiffness is defined by two (in-, ${\mu} $, and out-of-plane, $G$) shear moduli, decoupling tensile/inflation properties from shear responses. Based on this model, we developed an algorithm utilizing guided mechanical waves in a bounded NITI medium to reconstruct both moduli from AµT-OCE. The model was confirmed ex vivo in rabbit [48], porcine [49] and human [32,33] models and in vivo with rabbit models [48].
For ex vivo human corneas, we evidenced that both in- and out-of-plane post-CXL corneal shear moduli experienced an averaged two-fold increase in Young’s modulus and an almost four-fold increase in the out-of-plane shear modulus G [32], assuming the whole thickness was treated. That confirmed that CXL increases inter-corneal lamellae crosslinks but less affects corneal deformational properties, defined by the Young’s modulus that increases less and has more implications for potential refractive changes.
Despite the success in quantifying CXL-induced corneal elastic properties, there is a common situation where the method [32] must be refined. For example, crosslinking is often heterogeneous in depth due to more pronounced riboflavin penetration in the corneal anterior where the solution is applied [50], and because UVA irradiation attenuates and is less effective as it propagates through cornea [51]. Together, this generally produces a clear demarcation line (DL) between treated and untreated regions [52], suggesting a two-layer structure postoperatively.
As noted above, mechanical waves generated in the cornea are typically guided. They occupy the entire thickness of the cornea and, therefore, carry depth-accumulated information. Thus, reconstructing the depth dependence of corneal moduli is difficult with wave based OCE alone without a good estimate of the CXL penetration depth into the cornea.
Blackburn et al. [53] recently introduced a novel metric to track CXL penetration within the cornea using time-resolved OCT. They demonstrated that the phase decorrelation decay rate of the complex OCT signal is reduced in the CXL area, which can be used to distinguish treated from untreated areas postoperatively.
In this paper, we combine the method described in [53] with AµT-OCE measurements to explore possible reconstruction of both in- and out-of-plane corneal elastic moduli over depth in a partially CXL-treated ex vivo human cornea. We developed an analytical model of guided wave propagation accounting for multiple anisotropic layers, each with distinct stiffness moduli and thickness, to properly account for CXL-induced corneal layering. The model and methodology are validated using numerical simulations and phantom studies.
Baseline elastic moduli in an untreated ex vivo cornea can be determined from AµT-OCE prior to CXL so that induced changes in the treated anterior layer can be quantified by fitting the post-CXL wave dispersion dependence in the frequency-wavenumber ($f\textrm{-}k$) domain. We also discuss how the depth-distribution of stiffness affects the effective engineering stiffness of the entire cornea and show that assessing stiffness in both layers is needed to properly predict corneal deformation and quantify surgical outcomes.
2. Method
2.1 Cornea preparation
A corneal-scleral ring stored in Optisol (Chiron Ophtalmics) was obtained from CorneaGen. It came from a 26 year-old donor and was stored for less than 30 days. The corneal-scleral button was mounted on an artificial anterior chamber (Barron, CorzaMedical; see Fig. 2), connected through the inlet port to an elevated bath filled with balanced saline solution (BSS) to apply a controlled pressure mimicking intraocular pressure (IOP) on the anterior segment of the cornea. The outlet port remained closed to allow the IOP to settle at 15 mmHg within the chamber, corresponding to human physiological conditions [54]. CXL followed the Dresden protocol [6]. First, the epithelial membrane was removed. Then, the cornea was soaked in riboflavin for 30 minutes by applying a 50 µL drop of 0.1% riboflavin in 20% dextran solution every two minutes. It was then exposed to 3 mW/cm2 of 370 nm UV light for 30 minutes, while a drop was reapplied every 5 minutes.
Fig. 2. Picture of the experimental set up during UV-CXL. a) Acoustic micro-tapping transducer. b) Artificial anterior chamber with c) inlet port connected to the elevated bath and outlet port closed to control IOP.
2.2 AµT-OCE imaging system
A cylindrically focused air-coupled (AµT) transducer, operating at a 1 MHz frequency, launched mechanical waves in the cornea using acoustic radiation force. A 100 µs linear chirp was sent to the transducer. The chirp, with a bandwidth ranging from 0.95 MHz to 1.05 MHz minimized potential standing wave effects between the transducer and tissue surface. The cylindrical focus generated quasi-planar guided waves within the cornea. More details about AµT can be found in [24]. A spectral domain OCT system with a 46.5 kHz A-scan rate was used to track wave propagation and structural changes [24,30,32,33]. The OCT system operated in M-B mode, where a single push was triggered by the system while 512 consecutive A-scans were taken at a fixed location (M-scan). The M-scan sequence and push excitation were repeated for 256 locations, creating a volume with 256 $x$-samples, 1024 $z$-samples and 512 $t$-samples (see Fig. 3(a)), with an effective imaging range of 6 mm × 1.2 mm × 11 ms.
Fig. 3. Diagram of spectral domain time-resolved OCT and AµT-OCE measurements in a pre-CXL cornea. a) 3D ($x,\; z$ and $t$), wavefield after AµT excitation. The wavefield of the initial time sequence (black dotted region) is used for elastic moduli reconstruction, and data at the end of the sequence are used for phase decorrelation measurements. b) Wavefield ($x\textrm{-}t$ plot) showing the guided wave propagation during the first 5 ms of the acquisition sequence. c) $f\textrm{-}k$ spectrum obtained by 2D-FFT of the $x\textrm{-}t$ plot showing the dispersion dependence of the first anti-symmetric mode ${A_0}$. The red curve indicates the best fit obtained with the NITI model [33]. d) Structural OCT image obtained by averaging the last 7 ms of the raw OCT signal. e) Phase decorrelation function $g(\tau )$ at the location indicated by the red square in d).
The particle velocity along the probe beam direction was obtained from the optical phase difference between two consecutive A-lines at each location [55]. The spatio-temporal ($x\textrm{-}t$) surface signature of the guided wave was computed from an exponentially weighted-average of the particle velocity over the first 180 µm. As shown in Fig. 3(b), the guided wave only propagated during the first 4 ms of the scans. This period was used to determine the stiffness of the material by fitting the computed dispersion curve in the frequency-wavenumber domain ($f\textrm{-}k$) obtained from the 2D Fourier spectrum (Fig. 3(c)). This procedure is detailed in Sections 2.4–2.6. On the other hand, data from the last 7 ms were used to study structural changes with phase decorrelation [53] (see Section 2.7). Raw OCT data before and after CXL are provided in Dataset 1 (Ref. [56]).
2.3 Multi-layer NITI model
Like most biological tissue, the cornea is nearly incompressible. In addition, its structure implies a transverse isotropy [49] and, therefore, its mechanical behavior under small deformation should be described with a NITI model [33]. In Voigt notation, Hook’s law of stress and strain for a NITI material takes the form:
Because the cornea is a nearly incompressible soft tissue, its Young’s modulus does not depend on $\lambda $. Therefore, among the four elastic constants, only G and $\mu$ (respectively the out-of-plane and in-plane shear moduli) are needed to predict corneal deformation under mechanical loading.
The effects of CXL on the cornea depend on depth. Several recent studies showed that postoperative CXL corneas might experience non-uniform crosslinkage with depth. The transition between crosslinked (anterior) and non-crosslinked (posterior) parts tends to be sharp rather than smooth [51,52]. This effect is also observed in our experiments (see Section 2.7). Thus, a two-layer structure is considered an appropriate model to quantify postoperative corneas.
CXL was also shown to change collagen fiber diameter and interfibrillar spacing [50], but nothing suggests a modification of its macroscopic anisotropic organization. Based on this observation, we developed a multi-layer model to predict wave propagation within CXL corneas (see section 1 of Supplement 1) that accounts for any arbitrary number of layers, each with a stress-strain relationship given by Eq. (1) and linked by solid-solid boundary conditions (continuity of normal components of stress and displacement across every interface). Accounting for the external boundary conditions (liquid below and air above the cornea) and the finite thickness of the medium, the dispersion relation for guided waves can be determined directly from stiffness moduli ${G_n}$ and ${{\mu} _n}$ and the thickness ${h_n}$ of each layer. A more detailed description of the multi-layer model is in section 1 of Supplement 1. Although only 2 layers were considered here, the multi-layer model can be used to accurately capture more complicated transitions between crosslinked and non-crosslinked areas.
In an untreated cornea, only the first anti-symmetric mode, referred to as ${A_0}$, typically propagates in the range of frequencies that can be recorded in OCE (usually < 5 kHz). Because a partial-CXL cornea consists, in our approximation, of two horizontally assembled layers, each having a vertically aligned symmetry axis, this symmetry holds for the global material. Thus, only the ${A_0}$-mode is also expected in the partially crosslinked cornea.
2.4 Fitting $f\textrm{-}k$ spectra pre- and post-CXL
Prior to treatment, the cornea was assumed homogeneous, which in our model corresponds to a single NITI layer bounded above by air and below by water. The experimental $f\textrm{-}k$ spectrum (see Fig. 3(c)) was obtained by computing the 2D FFT of the $x\textrm{-}t$ plot. Shear moduli G and $\mu $ in pre-CXL cornea were obtained by fitting the measured $f\textrm{-}k$ spectrum with the analytic dispersion relation for the ${A_0}$-mode [32,33,48,49].
In the CXL-treated cornea, the thickness of both layers can be measured (see Section 3.1), and the posterior layer is assumed to still possess the original (i.e., untreated) elastic properties. Thus, the 2-layer model with known elastic moduli of the bottom (untreated) layer can be used to determine the stiffness of the top layer.
To ensure reliable fitting for all cases, we computed a goodness of fit (GOF) metric $\Phi = \frac{{\mathop \sum \nolimits_f {\chi _{fit}}(f )}}{{\mathop \sum \nolimits_f {\chi _{max}}(f )}}$, where ${\chi _{fit}}(f )$ corresponds to the energy of the 2D spectrum covered by the best analytical solution (one or two layers) at a given frequency f and ${\chi _{max}}(f )$ is the unconstrained maximal energy of the spectrum at frequency f. Based on recent results (see Supplemental Material in [32]), reliable fitting in human ex vivo corneas is associated with values of $\Phi \;>\; 0.9$. An example of a 2D-spectrum and the fitted ${A_0}$-mode obtained with this procedure for the untreated case is shown in Fig. 3(c).
2.5 FEM simulations
We designed finite element (FEM) simulations in OnScale to determine the accuracy of our multi-layer NITI model in reconstructing stiffness along corneal depth. The geometry is shown in Fig. 4(a). Corneal boundary conditions were replicated, with the material bounded above by air and below by water. The speed of sound in all layers (material and water) was fixed to avoid the reflection of compressional waves at the air-liquid and inter-layer boundaries. It also improved the absorption of compressional waves at the absorbing boundaries and, thus, avoided divergent simulations. The simulated transient excitation of broadband elastic waves closely matched that used in AµT experiments. More details about the simulations can be found in [33].
Fig. 4. Finite element simulations to study the effects of a layered structure for post-CXL cornea. a) Geometry of the two-layered material used in simulations, bounded above by air and below by water. b) Top surface spatio-temporal signature ($x\textrm{ - }t$ plot) of the guided wave for a two-layer case with ${G_{ant}} = 296.8\; $kPa, $\mu _{ant} = \; 34.6\; $MPa, ${G_{pos}} = 59.5\; $kPa, ${{\mu} _{pos}} = \; 7.3\; $MPa, ${h_{ant}} = 150\; \mathrm{\mu} \textrm{m}$ and $h_{pos} = 370\; \mathrm{\mu}\textrm{m}$. c) 2D Fourier spectrum of the wave studied in b) showing the main propagating ${A_0}$-mode and, in red, the analytical solution obtained from the 2-layer NITI model with identical parameters.
Based on phase-decorrelation measurements (see Section 2.6), we assumed that post-CXL two layers with distinct thicknesses were formed within the cornea, the top layer being stiffer than the bottom one. Stiffness values assessed from experiments were also used in the simulations. We used the top surface vertical particle velocity of the simulated wave (see Fig. 4(b)), and its associated $f\textrm{-}k$ spectrum (see Fig. 4(c)), to show that the analytical solution obtained from the 2-layer model (see Eq. S39 in Supplement 1, section 1) closely matches the $f\textrm{-}k$ spectrum of numerically simulated wave propagation. This confirms the accuracy of the analytical model in quantifying measurements of corneal elastic moduli and their variation with depth.
2.6 Phantom studies
The appropriateness of the 2-layer model was also validated in phantom experiments. More details about these experiments are given in section 2 of Supplement 1 and raw data are provided in Dataset 1 (Ref. [56]).
Agarose phantoms were prepared using different agar concentrations. Four different phantoms were made. The first two phantoms (Phantom A and Phantom B, Table S1 in Supplement 1, section 2) were homogeneous single layers with different stiffness, in which three (3) and two (2) wt. % of agar were used respectively. Titanium nanoparticles were added to the mixtures (respectively 0.1 and 0.2 wt % $\textrm{Ti}{\textrm{0}_\textrm{2}}$) to create optical scattering and differentiate phantom structures. Both phantoms were isotropic. For all experiments, the phantoms laid on top of a liquid bath, replicating corneal boundary conditions. Their elastic moduli ${ G}$ were measured by fitting the $f\textrm{-}k$ spectra of propagating Lamb waves with a theoretical model of an isotropic layer bounded with air on its top and water at its bottom [33]. The measured moduli were ${{G}_{{{stiff}}}} = 216 \; \textrm{kPa}\pm 2$ kPa and ${{ G}_{{{soft}}}} = {\; }84{\; } \textrm{kPa} \pm 5$ kPa for Phantoms A and B respectively (see Table S1, Supplement 1, section 2). Measured waveforms ($x\textrm{-}t$ plots) in these phantoms, their $f\textrm{-}k$ spectra and fitting curves are shown in Fig. S1.
Two bilayer phantoms were made of the same concentration of agar and $\textrm{Ti}{\textrm{0}_\textrm{2}}$. Phantom C had a stiffer layer on the top of the softer layer. The stiffer solution was brought to 90$^\circ \textrm{C}$ and poured into a petri dish for solidification. Its thickness was about 500 $\mathrm{\mu}$m. Then the softer solution was brought to 90$^\circ \textrm{C}$ and poured on top off the stiff layer. Phantom D had the opposite layer orientation (see Table S1). An example of the OCT structural image of bilayer Phantom C is given in Fig. 5(a). The image-based calculated total thickness was $1000 {\; } \mathrm{\mu} \textrm{m} \pm 9$ $\mathrm{\mu}$m, with approximately 2/3 of the full thickness for the top layer and 1/3 for the bottom layer. Elastic moduli of the layers in the bilayer phantoms were assumed equal to those for single layers (Phantoms A and B).
Fig. 5. Bilayer isotropic phantom experiments. a) structural image of the bilayer phantom with stiffer layer on top. The measured $x\textrm{-}t$ plot is shown in b) and its associated $f\textrm{-}k$ spectrum in c). The best fit solution with the 2-layer model is displayed in c).
The wavefield and associated $f\textrm{-}k$ spectrum in the bilayer Phantom C are shown in Figs. 5(b), (c) respectively. Note that the frequency content in the simulations is a bit broader than in experiments. This is mainly due to a possibly extended push width at the phantom surface provided by the AµT transducer focusing and the influence of material viscosity at frequencies above 3.5 kHz. Clearly, a single ${{ A}_0}$-mode associated with the ‘overall’ phantom is present, with no distinct modes from individual layers. Furthermore, Figs. S3, 4 in section 2 of Supplement 1 clearly show that guided wavefields, group velocity and $f\textrm{-}k$ spectra do not change with depth in bilayer phantoms, i.e., do not depend on the layer. Thus, guided wave propagation completely mixes phantom properties, and data acquired at different depths cannot be used for depth-resolved reconstruction of moduli. A more detailed analysis can also be found in section 2 of Supplement 1. However, this fact simplifies the reconstruction of layer properties with the multi-layer model presented here and described in section 2 of Supplement 1 in detail. Indeed, because the wavefield of the guided ${{ A}_0}$-mode does not change with depth, it can be acquired at any surface or at any depth, or averaged over depth even if the material stiffness varies with depth.
Because the total thickness of the bilayer phantoms is known, and the respective thickness and baseline stiffness of each layer is also known, we used our analytical model and the method introduced above to determine the stiffer layer moduli in Phantoms C and D to mimic the situation with a CXL corneal layer (see Sections below). The value of the softer layer was fixed so that ${{ G}_{{soft}}} = 84\; \textrm{kPa}\; \pm 4$ kPa. Using this method, the stiffer layer was found to be on average ${G}_{{{stiff}}} ^{2{ lay}}={\; }211 \; \textrm{kPa}\; \pm 7$ kPa for the Phantom C shown in Fig. 5, which is in good agreement with the control value ${{ G}_{{{stiff}}}} = 216\; \textrm{kPa} \; \pm 2$ kPa. Measurements with another Phantom D, which has the stiffer layer at the bottom, were also in good agreement with the control value. The results can be found in Table S1 (see Supplement 1, section 2).
Although the phantoms considered here are isotropic, we believe that these studies well justify the approach introduced for anisotropic cornea layers.
2.7 Phase decorrelation OCT (PhD-OCT)
Blackburn et al. [53] have recently introduced a metric to track CXL penetration within the cornea using time-resolved OCT. It was shown that the phase decorrelation decay rate of the complex OCT signal is reduced in CXL areas and can distinguish treated from untreated areas after the procedure.
In our study, the autocorrelation function of the signal g($\tau$) was computed as described in [53]: the correlation of six averaged adjacent pixels within a given A-line is obtained over 15 consecutive samples at 46,500 Hz:
3. Results
3.1 Thickness of the CXL layer
The spatial distribution of the OCT intensity signal (Figs. 6(a), (b), (e)) and phase decorrelation images (Figs. 6(c), (d), (e)) both show clear layering in the treated cornea. The effect of CXL is not homogeneous across the cornea, with a more pronounced effect at the center (∼150 µm) than at the periphery (∼ 100 µm). For the present case, we estimated that about 30% of corneal thickness was treated effectively. The treated cornea is thinner than that prior to CXL (its thickness reduced from 575 µm to 520 µm), as generally observed in the literature [58,59].
Fig. 6. Short time decorrelation before and after CXL. Structural OCT images averaged over the last 7 ms of the OCT scan for a) pre-CXL and b) post-CXL. Maps of decorrelation coefficient Γ for c) pre- and d) post-CXL. e) Profile of OCT intensity and Γ along the red dotted line shown in b) and d) for the CXL cornea.
3.2 Stiffness of CXL-treated corneal layers
AµT-OCE scans, taking approximatively 3 s to acquire and save data, were obtained prior and post-CXL. A space-time x-t plot of the vertically polarized particle velocity (see Figs. 3(a), b) in the untreated cornea was used to compute the $f\textrm{-}k$ spectrum (see Fig. 3(c)), which was then fitted using the procedure detailed in Section 2.4, assuming the cornea as a single homogeneous layer. The fitting routine was also detailed in our recent work [32,48]. Results for the reconstructed in- (${\mu} $) and out-of-plane ($G$) corneal shear moduli pre-CXL are ${\mu} = 7.6\; \textrm{MPa}\; \mp \; ({7,13} )$ MPa and $G = 59.5 \; \textrm{kPa}\; \mp ({5,8} )$ kPa (see Table 1).
Table 1. Corneal thickness ($h$), in- (${\mu} $) and out-of-plane ($G$) elastic moduli and goodness of fit ($\Phi $) pre- and post-CXL.
To reconstruct depth-dependent stiffness moduli after CXL, it is assumed that: i) the thickness of both the anterior and posterior layers can be measured using dynamic OCT from phase-decorrelation and/or intensity variation methods (see Fig. 6, estimated to be ${h_{ant}} = 150\; $µm and ${h_{pos}} = 370$ µm with both methods); ii) the stiffness of the posterior layer remains unchanged after CXL; iii) the effect of CXL is homogeneous in the anterior layer. Fixing known parameters (untreated cornea thicknesses and posterior stiffness moduli), stiffness moduli of the anterior cornea layer can be determined by fitting the wave dispersion curve in the $f\textrm{-}k$ domain with the 2-layer model (Fig. 6(b)). We found an increase in both stiffness moduli ${G_{ant}} = 296.8\; \textrm{kPa}\; \mp\; ({60,\; 92} )$ kPa and ${{\mu} _{ant}} = 34.6 \; \textrm{MPa}\; \mp ({20,23} )$ MPa compared to those for the posterior region ${G_{pos}} = 59.5\; \textrm{kPa}\; \mp ({5,\; 8} )$ kPa and ${\mu} _{pos} = 7.6 \; \textrm{MPa} \; \mp \; ({7,\; 13} )$ MPa (see Figs. 7(c), d). The goodness of fit for the 2-layer model ($\Phi = 0.950$) remained within the range of reliable fitting. The results are summarized in Table 1.
Fig. 7. Reconstruction of moduli in a post-CXL cornea. a) Vertical component of guided wavefield measured in the treated cornea. b) 2D-spectrum computed from a). c) 2D Goodness of fit (GOF) surface when the fit is performed with the 2-layer model to determine the anterior layer elastic moduli. d) Projection of the surface plot for $\Phi = 0.99 \times max(\Phi )$. Error bars are computed as the intersection of the vertical and horizontal directions with the fitted ellipse. The global maximum of $\Phi $ is indicated in c) and d) by the circular marker.
GOF variation for a given sample at a fixed IOP was previously shown to be about 1% [32]. We used this fact to build error bars for the present study, illustrated in Fig. 7(d) for the two-layer fitting procedure. First, we fit the projection of the goodness of fit surface for a value of $\Phi = 0.99 \times max(\Phi )$ (i.e., 1% below the optimal GOF), with an ellipsoidal function that best described the shape of iso-goodness levels. Then, we computed the uncertainty intervals as the intersection of the horizontal and vertical direction with the fitted ellipse. This produces asymmetric uncertainty intervals (particularly for $\mu $), reflecting the asymmetric variation of $\Phi $.
3.3 Mixing rules for effective engineering moduli of layered materials
For clinical diagnosis, the most critical information is how the cornea deforms under physiological loads. In other words, assessing the in-depth distribution of stiffness is only valuable if it can be related to deformation of the layered cornea. This can be achieved, for example, using numerical simulations or by deriving a set of effective moduli that directly relates to the overall stiffness of the composite corneal structure.
A theory for effective moduli of multi-layer materials was developed in the early 1970’s for composite materials. It is now accepted in material science and broadly used in the development of composite structures [60–62]. The derivation of these ‘effective’ material engineering moduli is based on ‘mixing rules’ of stiffness moduli across depth using the following assumptions: i) out-of-plane stresses and in-plane strains are uniform across thickness; ii) in-plane stresses and out-of-plane strains are averaged across thickness based on layer volume fractions. Note that the solution is valid only for low-frequency material deformation, which are, fortunately, usually related to physiologically induced stresses. As such, even if these definitions do not hold for perturbations of any kind, they are appropriate for in vivo corneal response to physiological loads.
Sun et al. [61] have demonstrated that an effective out-of-plane engineering modulus ${G_{eff}}$ can be computed using the inverse mixture rule for out-of-plane material constants of individual material layers:
where $h$ is a total material thickness, ${h_n}$ is a thickness of the ${n^{\textrm{th}}}$ layer and ${G_n}$ is the out-of-plane modulus of this layer. On the other hand, an effective in-plane engineering modulus ${\mu} _{eff}$ can be obtained from the mixture rule: where ${\mu} _n$ is the in-plane modulus of the ${n^{\textrm{th}}}$ layer.Based on the mixture rules described above, the effective low-frequency engineering moduli of a partially treated cornea for our case are computed to be ${G_{eff}} = 77.3 \; \textrm{kPa}\; \mp ({6,10} )$ kPa and ${\mu} _{eff} = 15.4 \;\textrm{MPa} \; \mp \; ({8,11} )$ MPa (see Table 1).
Effective corneal moduli, describing how the entire cornea deforms, may be computed if an overall estimate of surgical outcome is needed, i.e, if the desired effect is reached. For more detailed quantification, and to study longitudinal postoperative changes in the cornea, both pairs of anterior and posterior moduli can be used.
4. Discussion and conclusions
In this study, we combined structural OCT with dynamic AµT-OCE to assess the penetration depth of CXL in the cornea. Analyzing the brightness of structural OCT images and the rate of image decorrelation between consecutive B-scans, we confirm a sharp transition between CXL and untreated cornea layers and measure the thicknesses of treated and untreated layers. This finding suggests a 2-layer medium model for the treated cornea that can be used to reconstruct both in- and out-of-plane elastic moduli in the treated layer.
Numerical simulations and experiments performed in controlled isotropic phantoms confirmed the robustness of the two-layer fitting method (see Supplement 1, section 2). The main feature of ${A_0}$-mode propagation in a heterogeneous in-depth (or layered) medium is that guided waveforms are the same at all depths. Their $f\textrm{-}k$ spectra as well as group velocity are also depth-independent. Thus, attempts to directly translate depth-resolved measurements to mechanical moduli (as computed broadly in the literature) are inaccurate. However, this feature is very helpful for the method considered here because waveforms can be acquired at any surface or depth, or averaged over depth, even if the material stiffness varies with depth.
As we discussed in detail in our previous work [32,48], the ${A_0}$-mode dispersion spectrum is much more sensitive to variations in out-of-plane modulus $G$ than to variations in in-plane modulus ${\mu} $. This results in large error bars in the reconstruction of ${\mu} $ (see Table 1) compared to those for $G$. However, here we used only a single measurement, and reconstruction accuracy can be improved by repeating AµT measurements. Note that the error bar asymmetry comes from the asymmetry of the GOF function (see Figs. 7(c, d)), which was previously described in Refs. [32,48].
Using AµT-OCE, we tracked guided wave propagation in the cornea before and after CXL. Using the NITI model in each layer, we quantified depth-localized corneal stiffening with CXL. Because the ${A_0}$-mode occupies the whole thickness, it yields averaged information about material stiffness (i.e., averaged over its two parts). We determined the effective moduli in the treated corneal layer by fitting the 2D-spectrum with the two-layer model using known thicknesses for each layer and known elastic moduli for the untreated part (obtained from OCE measurements pre-CXL). Effective engineering moduli of the entire cornea could be then calculated using Eqs. (4) and (5), respectively, for out-of- and in-plane moduli.
As explained in Section 3.3, elastic moduli determined for the anterior (CXL-treated) corneal layer can be used to compute effective corneal engineering moduli, where ${{\mu} _{eff}}\; $uses a simple mixture rule whereas ${G_{eff}}$ requires an inverse mixture rule. This implies that even if the treated layer $G$ experiences an almost six-fold increase, its effective increase for the whole tissue is only by a factor of 1.3. On the other hand, ${\mu} $ experiences a local five-fold increase in the anterior layer but the overall in-plane modulus, which is more directly related to deformations in response to physiological loads, increases by a factor of 2. These moduli are particularly important because they can be directly used to estimate the entire corneal deformation under external loads.
Because the mixture rules for engineering effective moduli (Eqs. (4), (5)) assume only low-frequency perturbations, it is interesting to check if they could also describe guided wave behavior in the partially crosslinked cornea when considered as an effective homogeneous single-layer material (see Supplement 1, section 3). We found that the effective engineering mixture rules cannot be applied to quantify ‘effective’ guided wave behavior in the layered medium and, more importantly, the guided wave behavior considered in the ‘effective’ single layer incorrectly describes effective corneal engineering moduli. Fitting the experimentally obtained f-k spectrum (Fig. 7(b)) with the single layer model (Supplement 1, section 3) results in a different (incorrect) set of reconstructed engineering corneal moduli (see Table 2). Thus, measuring the depth of CXL penetration into cornea and implementing the multi-layer guided wave model are both required to accurately assess post-CXL corneal mechanical moduli.
Table 2. Effective moduli measured with use of single or multi-layer models.
Although the relationship between acoustic bulk (longitudinal and shear) wave propagation speeds and mechanical moduli in multi-layer or multi-component media has been reported in several studies [63,64], these relationships have not been explored for guided waves due to their high geometric dispersion. The situation is even more complicated for anisotropic media. The lack of a complete solution to this problem does not affect the goal of this study and is definitely outside its scope. However, we would like to share an interesting observation. We empirically found that both the analytic model and ex vivo experiments in the human cornea sample suggest that a simple mixture rule (Eq. (5)) for both corneal moduli can be applied to approximately compute effective guided wave propagation in a multi-layer NITI medium, as detailed in section 3 in Supplement 1 and in Code 1 (Ref. [65]).
One limitation of this study is the assumption that an untreated cornea is homogeneous in depth. Indeed, several studies indicate that the specific structure of the anterior (compact and interwoven collagen network) and posterior (loosely connected lamellae) stroma makes the anterior cornea stiffer than the posterior. However, no consensus has been reached about the in-depth dependence of corneal elasticity. Reported results varied greatly depending on sample preparation and the methods used to assess corneal properties: indentation [66,67], extensometry [68,69], radiation force [70], Brillouin microscopy [71] and even static OCE [72] methods. Several studies even reported no gradient between anterior and posterior baseline cornea properties (prior CXL) [67,72]. Moreover, all the above-mentioned studies assumed corneal isotropy. Note that the model proposed in this work is not limited to bilayer soft tissues. If a consensus about the in-depth gradient of elasticity in healthy human corneas is reached, this additional information can be taken into account in the proposed model.
Another limitation is that the post-CXL cornea is assumed to contain two homogeneous layers. After corneal CXL, several studies showed changes in stromal keratocytes with cell and collagen fiber shrinkage, chromatin condensation, and apoptotic bodies [73]. Post UVA exposure, riboflavin can produce radicals or singlet oxygen molecules that induce covalent bonds connecting one polymer chain to another, causing in vivo cornea links in collagen fibers [74]. CXL acts unevenly over corneal thickness: UVA intensity is inversely proportional to depth in stromal tissue, and riboflavin concentration decreases linearly with increasing corneal depth in accordance with a diffusion gradient. Furthermore, about 65% of UVA radiation is absorbed in the first 200 µm [51]. The consequence of higher UVA intensity in the anterior is increased stiffening within the first 250 µm. This appears as a DL in slit-lamp examination, defining the transition from cross-linked stroma to untreated tissue [51,75–80]. The DL is typically sharp and relates to the threshold for CXL penetration induced by a combination of multiple different effects of the treatment that are currently not completely understood [81,82]. As such, the assumption of two homogeneous layers in post-CXL corneas is a reasonable starting point. However, if a more gradual transition between CXL corneal layers is evidenced, the multi-layer model introduced here can be used to compute ${A_0}$-mode dispersion for more sophisticated models of CXL (e.g., accounting for a gradient in stiffness or more complex structural changes).
Recent results suggest that reverberant OCE can reconstruct depth-dependent stiffness variations of untreated cornea assuming an isotropic model [23]. It would be interesting to compare it with our method in future studies when reverberant OCE is applied to a NITI material. Note, however, that reverberant OCE is not currently feasible in vivo because it uses contact vibrators. This is why guided wave-based OCE is still the only method capable of in vivo non-contact measurements of corneal anisotropic elasticity, ultimately with sub-mm lateral resolution [83].
Finally, we have shown that phase-sensitive OCT combined with AµT wave excitation can assess both the structure of human cornea and the depth-dependence of moduli due to CXL. These findings are essential for building personalized models of corneal deformation following CXL and, thus, better adapt crosslinking therapy for clinical use and predict its outcome. Further experiments on a larger group of cornea samples are required to generalize the present results.
Funding
National Institutes of Health (R01-AR077560); National Eye Institute (R01-EY026532, R01-EY028753).
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are available in Dataset 1, Ref. [56] and Code 1, Ref. [65].
Supplemental document
See Supplement 1 for supporting content.
References
1. M. Sridhar, “Anatomy of cornea and ocular surface,” Indian J. Ophthalmol. 66(2), 190 (2018). [CrossRef]
2. K. M. Meek and C. Boote, “The organization of collagen in the corneal stroma,” Exp. Eye Res. 78(3), 503–512 (2004). [CrossRef]
3. M. S. Borcherding, L. J. Blacik, R. A. Sittig, J. W. Bizzell, M. Breen, and H. G. Weinstein, “Proteoglycans and collagen fibre organization in human corneoscleral tissue,” Exp. Eye Res. 21(1), 59–70 (1975). [CrossRef]
4. K. M. Meek and C. Knupp, “Corneal structure and transparency,” Prog. Retinal Eye Res. 49, 1–16 (2015). [CrossRef]
5. H. Hashemi, S. Heydarian, E. Hooshmand, M. Saatchi, A. Yekta, M. Aghamirsalim, M. Valadkhan, M. Mortazavi, A. Hashemi, and M. Khabazkhoob, “The prevalence and risk factors for keratoconus: a systematic review and meta-analysis,” Cornea 39(2), 263–270 (2020). [CrossRef]
6. G. Wollensak, E. Spoerl, and T. Seiler, “Riboflavin/ultraviolet-a–induced collagen crosslinking for the treatment of keratoconus,” Am. J. Ophthalmol. 135(5), 620–627 (2003). [CrossRef]
7. P. S. Hersh, R. D. Stulting, D. Muller, D. S. Durrie, R. K. Rajpal, P. S. Binder, E. D. Donnenfeld, D. Durrie, D. Hardten, P. Hersh, F. Price, D. Schanzlin, W. Stark, R. D. Stulting, W. Trattler, and S. Trokel, “United States multicenter clinical trial of corneal collagen crosslinking for keratoconus treatment,” Ophthalmology 124(9), 1259–1270 (2017). [CrossRef]
8. Y. Goldich, A. L. Marcovich, Y. Barkana, Y. Mandel, A. Hirsh, Y. Morad, I. Avni, and D. Zadok, “Clinical and corneal biomechanical changes after collagen cross-linking with riboflavin and UV irradiation in patients with progressive keratoconus: results after 2 years of follow-up,” Cornea 31(6), 609–614 (2012). [CrossRef]
9. G. Wollensak and E. Iomdina, “Long-term biomechanical properties of rabbit cornea after photodynamic collagen crosslinking,” Acta Ophthalmol. 87(1), 48–51 (2009). [CrossRef]
10. P. Maier, T. Reinhard, and M. Kohlhaas, “Corneal collagen cross-linking in the stabilization of keratoconus,” Dtsch Arztebl Int 116(11), 184–190 (2019). [CrossRef]
11. G. Wollensak, E. Spoerl, and T. Seiler, “Stress-strain measurements of human and porcine corneas after riboflavin–ultraviolet-A-induced cross-linking,” J. Cataract Refract. Surg. 29(9), 1780–1785 (2003). [CrossRef]
12. S.-H. Chang, D. Zhou, A. Eliasy, Y.-C. Li, and A. Elsheikh, “Experimental evaluation of stiffening effect induced by UVA/riboflavin corneal cross-linking using intact porcine eye globes,” PLoS One 15(11), e0240724 (2020). [CrossRef]
13. I. M. Beshtawi, C. O’Donnell, and H. Radhakrishnan, “Biomechanical properties of corneal tissue after ultraviolet-A–riboflavin crosslinking,” J. Cataract Refract. Surg. 39(3), 451–462 (2013). [CrossRef]
14. S. Kaushik and S. S. Pandav, “Ocular response analyzer,” J. Curr. Glaucoma Practice 6(1), 17–19 (2012). [CrossRef]
15. W. J. Dupps and C. J. Roberts, “Corneal biomechanics: A decade later,” J. Cataract Refract. Surg. 40(6), 857 (2014). [CrossRef]
16. Y. Zhao, Y. Shen, Z. Yan, M. Tian, J. Zhao, and X. Zhou, “Relationship among corneal stiffness, thickness, and biomechanical parameters measured by Corvis ST, Pentacam and ORA in keratoconus,” Front. Physiol. 10, 740 (2019). [CrossRef]
17. C. J. Roberts, A. M. Mahmoud, J. P. Bons, A. Hossain, A. Elsheikh, R. Vinciguerra, P. Vinciguerra, and R. Ambrósio, “Introduction of two novel stiffness parameters and interpretation of air puff–induced biomechanical deformation parameters with a dynamic Scheimpflug analyzer,” J Refract Surg 33(4), 266–273 (2017). [CrossRef]
18. S. Bak-Nielsen, I. Bach Pedersen, A. Ivarsen, and J. Hjortdal, “The rigidity of corneas before and after corneal cross-linking - as measured by Corvis® ST,” Invest. Ophthalmol. Visual Sci. 54(15), 1613 (2013).
19. M. Sedaghat, M. Naderi, and M. Zarei-Ghanavati, “Biomechanical parameters of the cornea after collagen crosslinking measured by waveform analysis,” J. Cataract Refract. Surg. 36(10), 1728–1731 (2010). [CrossRef]
20. M. W. Urban, C. Pislaru, I. Z. Nenadic, R. R. Kinnick, and J. F. Greenleaf, “Measurement of viscoelastic properties of in vivo swine myocardium using lamb wave dispersion ultrasound vibrometry (LDUV),” IEEE Trans. Med. Imaging 32(2), 247–261 (2013). [CrossRef]
21. Q. Zeng, M. Honarvar, C. Schneider, S. K. Mohammad, J. Lobo, E. H. T. Pang, K. T. Lau, C. Hu, J. Jago, S. R. Erb, R. Rohling, and S. E. Salcudean, “Three-dimensional multi-frequency shear wave absolute vibro-elastography (3D S-WAVE) with a matrix array transducer: implementation and preliminary in vivo study of the liver,” IEEE Trans. Med. Imaging 40(2), 648–660 (2021). [CrossRef]
22. R. Li, X. Qian, C. Gong, J. Zhang, Y. Liu, B. Xu, M. S. Humayun, and Q. Zhou, “Simultaneous assessment of the whole eye biomechanics using ultrasonic elastography,” IEEE Trans. Biomed. Eng. 70(4), 1310–1317 (2023). [CrossRef]
23. F. Zvietcovich, P. Pongchalee, P. Meemon, J. P. Rolland, and K. J. Parker, “Reverberant 3D optical coherence elastography maps the elasticity of individual corneal layers,” Nat. Commun. 10(1), 4895 (2019). [CrossRef]
24. Ł. Ambroziński, S. Song, S. J. Yoon, I. Pelivanov, D. Li, L. Gao, T. T. Shen, R. K. Wang, and M. O’Donnell, “Acoustic micro-tapping for non-contact 4D imaging of tissue elasticity,” Sci. Rep. 6(1), 38967 (2016). [CrossRef]
25. P. Song, H. Zhao, A. Manduca, M. W. Urban, J. F. Greenleaf, and S. Chen, “Comb-push ultrasound shear elastography (CUSE): a novel method for two-dimensional shear elasticity imaging of soft tissues,” IEEE Trans. Med. Imaging 31(9), 1821–1832 (2012). [CrossRef]
26. A. Manduca, T. E. Oliphant, M. A. Dresner, J. L. Mahowald, S. A. Kruse, E. Amromin, J. P. Felmlee, J. F. Greenleaf, and R. L. Ehman, “Magnetic resonance elastography: Non-invasive mapping of tissue elasticity,” Med. Image Anal. 5(4), 237–254 (2001). [CrossRef]
27. O. I. Kwon, C. Park, H. S. Nam, E. J. Woo, J. K. Seo, K. J. Glaser, A. Manduca, and R. L. Ehman, “Shear modulus decomposition algorithm in magnetic resonance elastography,” IEEE Trans. Med. Imaging 28(10), 1526–1533 (2009). [CrossRef]
28. D. C. Mellema, P. Song, R. R. Kinnick, M. W. Urban, J. F. Greenleaf, A. Manduca, and S. Chen, “Probe oscillation shear elastography (PROSE): a high frame-rate method for two-dimensional ultrasound shear wave elastography,” IEEE Trans. Med. Imaging 35(9), 2098–2106 (2016). [CrossRef]
29. H.-C. Liu, B. Gaihre, P. Kijanka, L. Lu, and M. W. Urban, “Acoustic force elastography microscopy,” IEEE Trans. Biomed. Eng. 70(3), 841–852 (2023). [CrossRef]
30. M. A. Kirby, I. Pelivanov, S. Song, Ł. Ambrozinski, S. J. Yoon, L. Gao, D. Li, T. T. Shen, R. K. Wang, and M. O’Donnell, “Optical coherence elastography in ophthalmology,” J. Biomed. Opt. 22(12), 1 (2017). [CrossRef]
31. X. Liang and S. A. Boppart, “Biomechanical properties of in vivo human skin from dynamic optical coherence elastography,” IEEE Trans. Biomed. Eng. 57(4), 953–959 (2010). [CrossRef]
32. M. A. Kirby, I. Pelivanov, G. Regnault, J. J. Pitre, R. T. Wallace, M. O’Donnell, R. K. Wang, and T. T. Shen, “Acoustic micro-tapping optical coherence elastography to quantify corneal collagen cross-linking,” Ophthalmology Science 3(2), 100257 (2023). [CrossRef]
33. J. J. Pitre, M. A. Kirby, D. S. Li, T. T. Shen, R. K. Wang, M. O’Donnell, and I. Pelivanov, “Nearly-incompressible transverse isotropy (NITI) of cornea elasticity: model and experiments with acoustic micro-tapping OCE,” Sci. Rep. 10(1), 12983 (2020). [CrossRef]
34. M. A. Kirby, P. Tang, H.-C. Liou, M. Kuriakose, J. J. Pitre, T. N. Pham, R. E. Ettinger, R. K. Wang, M. O’Donnell, and I. Pelivanov, “Probing elastic anisotropy of human skin in vivo with light using non-contact acoustic micro-tapping OCE and polarization sensitive OCT,” Sci. Rep. 12(1), 3963 (2022). [CrossRef]
35. L. A. Aleman-Castaneda, F. Zvietcovich, and K. J. Parker, “Reverberant elastography for the elastic characterization of anisotropic tissues,” IEEE J. Select. Topics Quantum Electron. 27(4), 1–12 (2021). [CrossRef]
36. H. Latorre-Ossa, J. Gennisson, E. De Brosses, and M. Tanter, “Quantitative imaging of nonlinear shear modulus by combining static elastography and shear wave elastography,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 59(4), 833–839 (2012). [CrossRef]
37. S. Goswami, R. Ahmed, S. Khan, M. M. Doyley, and S. A. McAleavey, “Shear induced non-linear elasticity imaging: elastography for compound deformations,” IEEE Trans. Med. Imaging 39(11), 3559–3570 (2020). [CrossRef]
38. S. Song, Z. Huang, T.-M. Nguyen, E. Y. Wong, B. Arnal, M. O’Donnell, and R. K. Wang, “Shear modulus imaging by direct visualization of propagating shear waves with phase-sensitive optical coherence tomography,” J. Biomed. Opt. 18(12), 1 (2013). [CrossRef]
39. Z. Han, S. R. Aglyamov, J. Li, M. Singh, S. Wang, S. Vantipalli, C. Wu, C. Liu, M. D. Twa, and K. V. Larin, “Quantitative assessment of corneal viscoelasticity using optical coherence elastography and a modified Rayleigh–Lamb equation,” J. Biomed. Opt. 20(2), 020501 (2015). [CrossRef]
40. K. V. Larin and D. D. Sampson, “Optical coherence elastography – OCT at work in tissue biomechanics [Invited],” Biomed. Opt. Express 8(2), 1172 (2017). [CrossRef]
41. I. Pelivanov, L. Gao, J. J. Pitre, M. A. Kirby, S. Song, D. Li, T. T. Shen, R. K. Wang, and M. O’Donnell, “Does group velocity always reflect elastic modulus in shear wave elastography?” J. Biomed. Opt. 24(07), 1 (2019). [CrossRef]
42. Z. Han, J. Li, M. Singh, S. R. Aglyamov, C. Wu, C. Liu, and K. V. Larin, “Analysis of the effects of curvature and thickness on elastic wave velocity in cornea-like structures by finite element modeling and optical coherence elastography,” Appl. Phys. Lett. 106(23), 233702 (2015). [CrossRef]
43. Z. Han, J. Li, M. Singh, C. Wu, C. Liu, R. Raghunathan, S. R. Aglyamov, S. Vantipalli, M. D. Twa, and K. V. Larin, “Optical coherence elastography assessment of corneal viscoelasticity with a modified Rayleigh-Lamb wave model,” J. Mech. Behav. Biomed. Mater. 66, 87–94 (2017). [CrossRef]
44. C. Li, G. Guan, X. Cheng, Z. Huang, and R. K. Wang, “Quantitative elastography provided by surface acoustic waves measured by phase-sensitive optical coherence tomography,” Opt. Lett. 37(4), 722 (2012). [CrossRef]
45. J. Zhu, L. Qi, Y. Miao, T. Ma, C. Dai, Y. Qu, Y. He, Y. Gao, Q. Zhou, and Z. Chen, “3D mapping of elastic modulus using shear wave optical micro-elastography,” Sci Rep 6(1), 35499 (2016). [CrossRef]
46. E. Koudouna, M. Winkler, E. Mikula, T. Juhasz, D. J. Brown, and J. V. Jester, “Evolution of the vertebrate corneal stroma,” Prog. Retinal Eye Res. 64, 65–76 (2018). [CrossRef]
47. M. Winkler, G. Shoa, S. T. Tran, Y. Xie, S. Thomasy, V. K. Raghunathan, C. Murphy, D. J. Brown, and J. V. Jester, “A comparative study of vertebrate corneal structure: the evolution of a refractive lens,” Invest. Ophthalmol. Visual Sci. 56(4), 2764 (2015). [CrossRef]
48. M. A. Kirby, G. Regnault, I. Pelivanov, M. O’Donnell, R. K. Wang, and T. T. Shen, “Noncontact acoustic micro-tapping optical coherence elastography for quantification of corneal anisotropic elasticity: in vivo rabbit study,” Trans. Vis. Sci. Tech. 12(3), 15 (2023). [CrossRef]
49. M. A. Kirby, J. J. Pitre, H.-C. Liou, D. S. Li, R. K. Wang, I. Pelivanov, M. O’Donnell, and T. T. Shen, “Delineating corneal Elastic anisotropy in a porcine model using noncontact OCT elastography and ex vivo mechanical tests,” Ophthalmol. Sci. 1(4), 100058 (2021). [CrossRef]
50. S. Akhtar, A. Smedowski, A. Masmali, A. Alkanaan, A. A. Khan, E. Almutleb, H. K. Mofty, H. I. Al-Debasi, R. Samivel, and T. Almubrad, “Effect of Ultraviolet-A and Riboflavin treatment on the architecture of the center and periphery of normal rat cornea: 7 days post treatment,” Exp. Eye Res. 219, 109064 (2022). [CrossRef]
51. L. Spadea, E. Tonti, and E. Vingolo, “Corneal stromal demarcation line after collagen cross-linking in corneal ectatic diseases: a review of the literature,” OPTH Volume 10, 1803–1810 (2016). [CrossRef]
52. C. Mazzotta, G. Wollensak, F. Raiskup, A. M. Pandolfi, and E. Spoerl, “The meaning of the demarcation line after riboflavin-UVA corneal collagen crosslinking,” Expert Rev. Ophthalmol. 14(2), 115–131 (2019). [CrossRef]
53. B. J. Blackburn, S. Gu, M. R. Ford, V. de Stefano, M. W. Jenkins, W. J. Dupps, and A. M. Rollins, “Noninvasive assessment of corneal crosslinking with phase-decorrelation optical coherence tomography,” Invest. Ophthalmol. Visual Sci. 60(1), 41 (2019). [CrossRef]
54. Y. Shiose, “Intraocular pressure: New perspectives,” Surv. Ophthalmol. 34(6), 413–435 (1990). [CrossRef]
55. R. K. Wang, S. Kirkpatrick, and M. Hinds, “Phase-sensitive optical coherence elastography for mapping tissue microstrains in real time,” Appl. Phys. Lett. 90(16), 164105 (2007). [CrossRef]
56. G. Regnault, M. A. Kirby, R. K. Wang, T. T. Shen, M. O’Donnell, and I. Pelivanov, “Cornea and phantom datasets for: “Possible depth-resolved reconstruction of shear moduli in the cornea following collagen crosslinking (CXL) with optical coherence tomography and elastography”,” figshare (2023), https://doi.org/10.6084/m9.figshare.23929626.
57. B. Chu, “Laser Light Scattering,” Annu Rev Phys Chem. 21(1), 145–174 (1970). [CrossRef]
58. S. A. Greenstein, V. P. Shah, K. L. Fry, and P. S. Hersh, “Corneal thickness changes after corneal collagen crosslinking for keratoconus and corneal ectasia: One-year results,” J. Cataract Refract. Surg. 37(4), 691–700 (2011). [CrossRef]
59. M. M. V. Cordeiro Barbosa, J. B. Barbosa, F. E. Hirai, and A. L. Hofling-Lima, “Effect of cross-linking on corneal thickness in patients with corneal edema,” Cornea 29(6), 613–617 (2010). [CrossRef]
60. P. C. Chou, J. Carleone, and C. M. Hsu, “Elastic constants of layered media,” J. Compos. Mater. 6(1), 80–93 (1972). [CrossRef]
61. C. T. Sun and S. Li, “Three-dimensional effective elastic constants for thick laminates,” J. Composite Mater. 22(7), 629–639 (1988). [CrossRef]
62. A. H. Nayfeh, Wave Propagation in Layered Anisotropic Media: With Applications to Composites, North-Holland Series in Applied Mathematics and Mechanics No. v. 39 (Elsevier, 1995).
63. S. Benjelloun and J.-M. Ghidaglia, “On the sound speed in two-fluid mixtures and the implications for CFD model validation,” European J. Mech. - B/Fluids 90, 152–168 (2021). [CrossRef]
64. A. A. Karabutov, “Laser optoacoustic measurement of paper porosity,” Acoust. Phys. 51(5), 560 (2005). [CrossRef]
65. G. Regnault, M. A. Kirby, R. K. Wang, T. T. Shen, M. O’Donnell, and I. Pelivanov, “Matlab tools for guided waves in composite nearly incompressible materials,” figshare (2023), https://doi.org/10.6084/m9.figshare.23504652.
66. J. M. Dias and N. M. Ziebarth, “Anterior and posterior corneal stroma elasticity assessed using nanoindentation,” Exp. Eye Res. 115, 41–46 (2013). [CrossRef]
67. J. A. Last, S. M. Thomasy, C. R. Croasdale, P. Russell, and C. J. Murphy, “Compliance profile of the human cornea as measured by atomic force microscopy,” Micron 43(12), 1293–1298 (2012). [CrossRef]
68. M. Kohlhaas, E. Spoerl, T. Schilde, G. Unger, C. Wittig, and L. E. Pillunat, “Biomechanical evidence of the distribution of cross-links in corneastreated with riboflavin and ultraviolet A light,” J. Cataract Refract. Surg. 32(2), 279–283 (2006). [CrossRef]
69. J. B. Randleman, D. G. Dawson, H. E. Grossniklaus, B. E. McCarey, and H. F. Edelhauser, “Depth-dependent cohesive tensile strength in human donor corneas: implications for refractive surgery,” J. Refract. Surg. 24(1), 1081597X (2008). [CrossRef]
70. E. R. Mikula, J. V. Jester, and T. Juhasz, “Measurement of an elasticity map in the human cornea,” Invest. Ophthalmol. Visual Sci. 57(7), 3282 (2016). [CrossRef]
71. G. Scarcelli, S. Kling, E. Quijano, R. Pineda, S. Marcos, and S. H. Yun, “Brillouin microscopy of collagen crosslinking: noncontact depth-dependent analysis of corneal elastic modulus,” Invest. Ophthalmol. Visual Sci. 54(2), 1418 (2013). [CrossRef]
72. T. J. Ferguson, S. Singuri, S. Jalaj, M. R. Ford, V. S. De Stefano, I. Seven, and W. J. Dupps, “Depth-resolved corneal biomechanical changes measured via optical coherence elastography following corneal crosslinking,” Trans. Vis. Sci. Tech. 10(5), 7 (2021). [CrossRef]
73. J. S. Dhaliwal and S. C. Kaufman, “Corneal collagen cross-linking: a confocal, electron, and light microscopy study of eye bank corneas,” Cornea 28(1), 62–67 (2009). [CrossRef]
74. P. Kamaev, M. D. Friedman, E. Sherr, and D. Muller, “Photochemical kinetics of corneal cross-linking with riboflavin,” Invest. Ophthalmol. Visual Sci. 53(4), 2360 (2012). [CrossRef]
75. G. Wollensak, E. Spoerl, M. Wilsch, and T. Seiler, “Keratocyte apoptosis after corneal collagen cross-linking using riboflavin/UVA treatment,” Cornea 23(1), 43–49 (2004). [CrossRef]
76. T. Seiler and F. Hafezi, “Corneal cross-linking-induced stromal demarcation line,” Cornea 25(9), 1057–1059 (2006). [CrossRef]
77. G. D. Kymionis, K. I. Tsoulnaras, D. A. Liakopoulos, C. A. Skatharoudi, M. A. Grentzelos, and N. G. Tsakalis, “Corneal stromal demarcation line depth following standard and a modified high intensity corneal cross-linking protocol,” J. Refract. Surg. 32(4), 218–222 (2016). [CrossRef]
78. A. Moramarco, A. Iovieno, A. Sartori, and L. Fontana, “Corneal stromal demarcation line after accelerated crosslinking using continuous and pulsed light,” J. Cataract Refract. Surg. 41(11), 2546–2551 (2015). [CrossRef]
79. Z. Gatzioufas, M. Balidis, and N. Kozeis, “Is the corneal stromal demarcation line depth a true indicator of corneal collagen crosslinking efficacy?” J. Cataract Refract. Surg. 42(5), 804 (2016). [CrossRef]
80. J. B. N. Malta, A. C. Renesto, B. K. Moscovici, H. K. Soong, and M. Campos, “Stromal demarcation line induced by corneal cross-linking in eyes with keratoconus and nonkeratoconic asymmetric topography,” Cornea 34(2), 199–203 (2015). [CrossRef]
81. J. Lin, “Influencing factors relating the demarcation line depth and efficacy of corneal crosslinking,” Invest. Ophthalmol. Visual Sci. 59(12), 5125 (2018). [CrossRef]
82. M. R. Santhiago and J. B. Randleman, “The biology of corneal cross-linking derived from ultraviolet light and riboflavin,” Exp. Eye Res. 202, 108355 (2021). [CrossRef]
83. G. Regnault, M. A. Kirby, M. Kuriakose, T. Shen, R. K. Wang, M. O’Donnell, and I. Pelivanov, “Spatial resolution in optical coherence elastography of bounded media,” Biomed. Opt. Express 13(9), 4851 (2022). [CrossRef]
