Possible depth-resolved reconstruction of shear moduli in the cornea following collagen crosslinking (CXL) with optical coherence tomography and elastography

: Collagen crosslinking of the cornea (CXL) is commonly employed to prevent or treat keratoconus. Although the global change of stiffness before and after treatment can be monitored with OCE, the depth dependence of this change is still unclear. Here we propose to combine phase-decorrelation measurement applied to OCT structural images and acoustic micro-tapping optical coherence elastography (AµT-OCE) to explore possible reconstruction of stiffness with depth within crosslinked corneas in an ex vivo human cornea sample. The analysis of experimental OCT images can help define the penetration depth of CXL into the cornea, which varies from ∼ 100μm in the periphery to ∼ 150μm in the central area and demonstrates a sharp transition between areas. This information was used in numerical simulations of the mechanical guided wave propagation in the two-layer medium with different elastic moduli corresponding to those for CXL and normal corneas. The simulation supports our hypothesis that ‘averaged’ elastic moduli can be reconstructed from guided wave propagation in the cornea with partial penetration of CXL, which may then be used to decouple the moduli in each layer. However, additional studies are required to define and characterize the reconstruction algorithm for the more general case of arbitrary CXL penetration depth and its possible variation over the cornea.


Introduction
At the interface between air and the inner eye, the cornea provides both protection and the primary optical element focusing light onto the retina. It contains multiple layers, including the epithelium and the stroma. The first acts as a barrier against the external environment and the latter maintains its 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].
Corneal diseases (such as keratoconus (KC)) and surgical complications from refractive surgeries (such as LASIK) may deform the cornea (ectasia) and reduce 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 (CXL) is a minimally invasive procedure that can potentially slow the progression of corneal ectasia [6][7][8][9]. UV light modifies the microstructure of cornea soaked in riboflavin and forms additional chemical bonds between collagen fibers in the stroma [10]. Post-treatment corneas become stiffer and more resistant to enzymatic digestion [11][12][13]. Although corneal topography (curvature map) and thickness map can be obtained prior to surgery, and needed 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. post-surgery corneal ectasia risks) cannot be accurately predicted without quantitatively mapping corneal elasticity. Thus, methods that can quantitatively reconstruct corneal elastic moduli are in high demand.
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][15][16]. Over 100 papers are published annually on the topic, and the Journal of Cataract and Refractive Surgery devoted an entire issue to the topic [15]. 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.
There is no consensus in the literature on the elastic model for cornea and the stiffness range even for healthy subjects. The most common model assumes an incompressible, isotropic, and linear elastic solid, where a single parameter, the Young's modulus (or equivalently the shear modulus = /3), defines elasticity. Destructive mechanical tests can determine ex vivo, with reported values for human cornea (at low-strain) of 800 kPa to 4.7 MPa for tensile loading [20][21][22][23][24], and 100 kPa to 3 MPa for inflation loading [25][26]. Note that the destructive nature of mechanical tests precludes their clinical translation.
Optical coherence elastography (OCE) is a promising tool to probe corneal biomechanics. Dynamic OCE can potentially deliver a non-contact and non-invasive measurement in a clinical environment [27][28][29][30][31] as it excites mechanical waves in the cornea (for example, using an airpuff or acoustic micro-tapping (AT [27]) and tracks them using phase-sensitive optical coherence tomography (OCT). Analyzing the shear wave propagation speed or dispersion leads to an estimate of the shear modulus. OCE studies reported corneal shear moduli in the range of 1.8 -52.3 kPa [32,33], in close agreement with values obtained from torsional and rheometry testing of ex vivo cornea (2.5 -47.3 kPa) [34][35][36]. However, all shear-based methods produce moduli differing by 1-2 orders of magnitude from those reported by tensile and inflation tests.
For years the Rayleigh wave group velocity was used to reconstruct corneal elasticity [31, 38, 27, 28-32, 39, 40], resulting in a 2 orders of magnitude mismatch in Young's modulus compared to tensile and inflation tests. We clearly showed that this approach is inappropriate [41] since the cornea is a layer bounded by air at the top and aqueous humor on the bottom. Thus, guided waves propagate with strong geometric dispersion [41,42].
Recently, we hypothesized that corneal anisotropy is the primary cause of these discrepancies [42]. Corneal microstructure supports this hypothesis. The stroma contains collagen lamellae running in-plane across its width. Lamellae make up approximately 90% of tissue thickness and account for the majority of the cornea's mechanical structure. They are stacked vertically in approximately 200-500 separate planes [43,44], suggesting an anisotropic mechanical behavior with very different responses to in-plane versus out-of-plane loads. We introduced a model of a nearly-incompressible transverse isotropic (NITI) medium [42] to explain the discrepancy. It is defined by two (in-, , and out-of-plane, ) shear moduli, decoupling tensile/inflation properties from shear responses. Based on this model, an algorithm utilizing guided mechanical waves in a bounded NITI medium to reconstruct both moduli from AT-OCE measurements has been developed.
In [45] we showed that both in-and out-of-plane post-CXL corneal shear moduli increase compared to their pre-surgery values. It is interesting that at the end of the CXL procedure corneal Young's modulus changes on average by about 2 times whereas the out-of-plane shear modulus changes by almost 4 times. Indeed, the CXL procedure is intended to better crosslink corneal lamellae, which it does. However, deformational properties are defined by the Young's modulus, which does not change nearly as much and has implications for potential refractive changes.
Note, however, that the effects of crosslinking are always homogeneous in depth as the penetration of riboflavin is more pronounced in the anterior region of the cornea, where the solution is applied during soaking [46]. Riboflavin penetration is purposely limited to minimize damage to the endothelium and ensure the biological integrity of the cornea [47].
As noted above, mechanical waves generated in the cornea are usually not bulk shear or surface Rayleigh waves; they are guided waves. Reconstructing corneal moduli using local estimates of the propagation speed at different depths leads to large errors. In fact, guided waves occupy the entire depth of the cornea and, therefore, average information over depth. Thus, reconstructing the depth dependence of corneal moduli is practically impossible without a good estimate of CXL penetration within the cornea.
Blackburn et al. [48] have 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 CXL areas and can be used to distinguish treated and untreated areas after a CXL procedure.
In this letter, we combine the methods described in [48] with AT-OCE measurements to explore possible reconstruction of both in-and out-of-plane corneal elastic moduli over depth. In addition, we use numerical simulations of guided wave propagation in a two-layer medium consisting of a CXL cornea layer on the top of the untreated cornea to demonstrate how guidedwave characteristics can be used to reconstruct corneal elasticity in both parts of the medium. Although our results show promise in monitoring the effects of CXL using a two-layer model, a more detailed study on model-based reconstruction of the in-depth distribution of elastic properties in CXL cornea is required to generalize our findings.

Cornea preparation
One corneal-scleral ring, stored in Optisol (Chiron Ophtalmics) was obtained from CorneaGen. This sample from a 26 years-old donor was stored for less than 30 days. The corneal-scleral button is mounted on an artificial anterior chamber (Barron, CorzaMedical; see figure 1), connected through the inlet port to an elevated bath filled with balanced saline solution (BSS) to apply a controlled pressure on the anterior segment of the cornea. The outlet port remained closed to allow the pressure to settle at 15 mmHg within the chamber, corresponding to human physiological conditions [49]. Crosslinking was performed following the Dresden protocol [6]. First, the epithelial membrane was removed from the sample. Then, the cornea was soaked in riboflavin for 30 min, by applying a 50µL drop of 0.1% riboflavin in 20% dextran solution every two minutes. The cornea was then exposed to a 3mW/cm 2 of 370 nm ultra-violet (UV) light for 30 minutes, while a drop was re-applied every 5 minutes.

AµT-OCE imaging system
A spectral domain OCT system with a 46.5 kHz frame rate was used to track wave propagation and structural changes within the cornea. As described in previous studies [27,31,42,45], a cylindrically focused air-coupled transducer, operating at 1MHz, generates a spatio-temporal sharp displacement at the surface of the cornea, using reflection based radiation force. Because of the cylindrical geometry of the transducer, the push is line-shaped and generates planar guided waves within the bounded tissue. The OCT system operated in M-B mode. A single push is triggered by the system while 512 consecutive A-scans are taken at a fixed location (Mscan). The M-scan sequence and push excitation are repeated for 256 locations, creating a threedimensional volume with 256 -samples, 1024 -samples and 512 -samples (see figure 2.a)), with an effective imaging range of 6mm × 1.2mm × 11ms. The particle axial vibration is obtained from the optical phase difference between two consecutive A-line scans at each location [51]. The spatio-temporal ( -) surface signature of the guided wave is computed from an exponential weighted-average of the vertical velocity over the first 180µm of the anterior part of the cornea. As shown in Figure 2.b), the guided wave only propagates during the first 4ms of the scans. This first part will be used to determine the stiffness of the material by fitting the measured dispersion curve in the frequency-wavenumber domain ( -), obtained from the 2D Fourier spectrum (see Figure 2.c)). This procedure is detailed in section 2.3. On the other hand, data from the last 7ms are used to study the structural changes with phase decorrelation (see section 2.4).

NITI model and fitting
The microstructure of the cornea implies a transverse isotropy in the direction normal to the corneal interface [52]. Like most biological tissue, the cornea is nearly-incompressible. Consequently, an exact description of its behavior under external loading requires a NITI (for Nearly Incompressible Transversely Isotropic) model [42]. In Voigt notation, the Hook's law of stress and strain for a TI material takes the form: where denotes engineering stress, denotes engineering strain, denotes shear stress, = 2 denotes shear strains, the subscript , and refer to the Cartesian axes and , , and are four independent elastic constants. In previous studies [52,53], we have demonstrated that , which accounts for tissue tensile anisotropy, cannot be determined from guided wave propagation but that the in-plane Young's modulus can be approximated as = ≅ 3 assuming tensile isotropy ( = 0). Thus, among the four elastic constants, only and , respectively the out-of-plane and in-plane shear moduli, are needed to predict deformations of a nearly incompressible tissue under mechanical loading. Accounting for the appropriate boundary conditions (water below and air above) and the finite thickness of the medium, the dispersion relationship of guided waves can be determined directly from and . Note that in cornea, only the first anti-symmetric mode, referred to as 0 , propagates in the range of frequencies that can be recorded in elastography (typically < 5kHz).
The experimentalspectrum (see Figure 2.c)) is obtained by computing the 2D-FFT of the --plot. Because the boundary conditions of the cornea on the anterior chamber highly differ with physiological ones (the cornea is clamped here), a structural resonance of the cornea is observed in the low-frequency regime (below 1kHz). To remove this artifact, we applied a temporal super-Gaussian filter (∼ 0.5ms FWHM) that follows the maximum of vibration velocity prior to computing the Fourier transform. The shear moduli and are obtained from fitting the measured --spectrum with the analytical dispersion relationship of the 0 mode.
To ensure reliable fitting, we computed a goodness of fit metric = ∑ ( ) corresponds to the energy of the 2D spectrum covered by the best analytical solution at a given frequency and ( ) is the unconstrained maximal energy of the spectrum at frequency . Based on recent results (see supplemental material in [12]), reliable fitting in human ex vivo corneas are associated with values of > 0. 9 . More details about the complete fitting procedure can be found in [45]. An example of a 2D-spectrum and the fitted 0 mode obtained with this procedure are shown in Figure 2.c).

Phase Decorrelation OCT (PhD-OCT)
Blackburn et al. [48] have 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 CXL areas and can be used to distinguish treated and untreated areas of the cornea after a procedure. The complex-valued autocorrelation of the signal ( ) is computed over 15 consecutive sample at 46,500 Hz, for six consecutive pixels within a given A-line: which is expected to follow an exponential decay [54]: where Γ is the decorrelation coefficient that is inversely proportional to the Brownian diffusion coefficient [54], meaning that the more coherent the material, the smaller the decorrelation coefficient. This procedure is performed starting at the , + 1, + 2, … A-line, where is the first time-sample used for the phase-decorrelation ( ( ) = 4ms). The decorrelation coefficient Γ is then computed using the averaged ( ) over the number of starting points by fitting with a first order polynomial (see Figure 2.e)): < ( ) > = − Γ. , where <> denotes the average over the number of starting points. In crosslinked parts of the cornea (anterior), the tissue stiffens, which implies that Γ should be smaller than in the untreated region (posterior). For post-processing, we rejected all fits were b < 0.95, corresponding in general to peripheral regions where the SNR (Signal to Noise Ratio) is too low.

Corneal stiffness before and after CXL
A single scan of AµT-OCE, taking approximatively 3s to acquire and save, was taken before and after CXL. The signature of the vertically polarized velocity was used to compute the -spectrum, which was fitted using the procedure detailed in section 2.3. The results for the reconstructed stiffness before and after CXL are shown in Table 1. As generally observed in the literature [55,56], the thickness is slightly reduced following crosslinking. The measured stiffnesses increased for from 61.8kPa to 132.2kPa and for from 5.3MPa to 8.1MPa. Such increase in both stiffness moduli following CXL is in good agreement with previously reported results for ex vivo human corneas [45]. The goodness of fit metric is = 0.96 before CXL and = 0.95 after CXL, which confirms the robustness of the fit and the validity of the results.

Depth dependance of CXL
The results for phase decorrelation are shown in Figure 3. For the untreated case (Figures 3.a,  c), both the OCT intensity and Γ are constant in directions normal to the interface. In the center, the higher scattering of laser light slightly influences the results, suggesting that Γ depends on the signal to noise ratio of the system. The CXL area can be seen in both structural OCT and Γ maps after CXL (Figures 3.b, d)). As shown in Figure 3.e), the drop of OCT intensity is associated with an increase of Γ, which demonstrates that both metrics can be used to identify the CXL area. The thickness of the CXL layer varies from ∼100μm in the periphery to ∼150μm in the central area of the cornea, which agrees with recent observation of lateral changes in the effect of crosslinking in corneas [46].

Numerical simulations
We designed a finite element (FEM) simulation to study the layering effect of CXL (section 3.2) on the reconstructed stiffness using the NITI model. The geometry is shown in Figure 4.a). Boundary conditions of the cornea were replicated so that the material is bounded above by air and below by water. The speed of sound in all layers (material and water) was fixed to avoid reflection of compressional waves at the different boundaries. It also improved the absorption of waves at the absorbing boundaries and, thus, avoided divergence of the simulations. We used transient excitation mimicking AµT experiments to generate broadband elastic waves within the material. More details about the simulations can be found in [42,50]. Based on the phasedecorrelation measurements (see Section 3.2), we assumed that after CXL, two layers with distinct stiffnesses but with identical thicknesses (each 0.25mm) were formed within the cornea, the top layer being stiffer than the bottom one.
We fixed 1 = 200kPa, 1 = 10MPa, 2 = 50kPa and 2 = 4MPa, and a total thickness ℎ = 0.5mm. Similar to OCE experiments, we used the top surface vertical signature of the guided wave (Figure 4.b) to compute its 2D -k spectrum (Figure 4.d). We also performed simulations in a single-layer medium with moduli corresponding to that averaged over the two-layer structure, i.e. 1 = 125kPa, 1 = 7MPa, 2 = 125kPa and 2 = 7MPa. The computedplot (Figure 4.c) for the homogeneous layer with averaged parameters is nearly identical to that for the two-layer medium. Furthermore, an analytically calculated dispersion curve for the case of averaged parameters fits well the 2D spectrum computed for the two-layer case with = 0.961 (red line in Figure 4.d). Thus, the simulation supports our hypothesis that 'averaged' moduli can be reconstructed from guided wave propagation in the cornea with partial penetration of CXL. These preliminary results show that defining the CXL penetration depth from structural OCT images may help decouple elastic moduli in each cornea layer. However, additional studies are required to define and characterize the reconstruction algorithm for the more general case of arbitrary CXL penetration depth and its possible variation over the cornea.

Discussion and Conclusions
In this study we combined structural OCT imaging with dynamic OCE to assess the penetration depth of CXL treatment in the cornea. Analyzing brightness of structural OCT images and image decorrelation between consecutive B-scans, we conclude that there is a sharp transition between CXL and untreated cornea. This finding allowed us to suggest a model of a two-part medium for the treated cornea, where cornea elastic properties remain constant in each part. This finding can be confirmed with destructive measurements of both layers, but it is outside the scope of the current study.
Both experimental results and numerical simulations show that guided waves deliver 'averaged' or 'effective' elastic properties of the cornea. When the thickness of the CXL layer can be measured, the two-layer model can be used to reconstruct mechanical moduli in both parts. Note, however, that such reconstructions may not be simple and require additional simulations, experimental studies, and depth-resolved mechanical tests.