Noninvasive Assessment of Corneal Crosslinking With Phase-Decorrelation Optical Coherence Tomography

Abstract

Purpose: There is strong evidence that abnormalities in corneal biomechanical play a causal role in corneal ectasias, such as keratoconus. Additionally, corneal crosslinking (CXL) treatment, which halts progression of keratoconus, directly appeals to corneal biomechanics. However, existing methods of corneal biomechanical assessment have various drawbacks: dependence on IOP, long acquisition times, or limited resolution. Here, we present a method that may avoid these limitations by using optical coherence tomography (OCT) to detect the endogenous random motion within the cornea, which can be associated with stromal crosslinking.

Methods: Phase-decorrelation OCT (PhD-OCT), based in the theory of dynamic light scattering, is a method to spatially resolve endogenous random motion by calculating the decorrelation rate, Γ, of the temporally evolving complex-valued OCT signal. PhD-OCT images of ex vivo porcine globes were recorded during CXL and control protocols. In addition, human patients were imaged with PhD-OCT using a clinical OCT system.

Results: In both the porcine cornea and the human cornea, crosslinking results in a reduction of Γ (P < 0.0001), indicating more crosslinks. This effect was repeatable in ex vivo porcine corneas (change in average Γ = −41.55 ± 9.64%, n = 5) and not seen after sham treatments (change in average Γ = 2.83 ± 12.56%, n = 5). No dependence of PhD-OCT on IOP was found, and correctable effects were caused by variations in signal-to-noise ratio, hydration, and motion.

Conclusions: PhD-OCT may be a useful and readily translatable tool for investigating biomechanical properties of the cornea and for enhancing the diagnosis and treatment of patients.

Corneal ectasia causes impaired quality of life¹ and is a leading indication for corneal transplantation.² Keratoconus, the most common form of corneal ectasia, affects approximately 1 in 375 individuals³ and has a predilection for progression in the second and third decades of life.⁴ Another type of ectasia, postsurgical ectasia, is a rare (0.04% to 0.6%) but serious complication of common refractive surgeries.^5–8 Although the etiology of ectasia is still incompletely understood, localized biomechanical abnormalities are believed to play a central role in the development of ectatic diseases like keratoconus and postrefractive surgery ectasia.^9,10 Computational modeling simulations^11,12 and experimental tissue studies that directly manipulate corneal material properties with collegenase¹³ have shown that regional weakening to levels similar to keratoconic explants^14,15 can reproduce the topographic characteristics of corneal ectasia without corneal thinning or other morphologic changes as a precondition. Given strong evidence of the causal role of corneal biomechanical abnormalities in these conditions, there is an equally strong rationale for developing clinically feasible, spatially sensitive biomechanical measurement methods. 

Further, the importance of understanding corneal biomechanics in ectasia is underscored by current treatment. The best practice for treating keratoconus is to stiffen the cornea through corneal crosslinking (CXL). The physical rationale behind CXL, that increasing molecular bonds halts the disease, aligns with the theory that weakened mechanical properties in the cornea cause ectasias. Proxy variables such as corneal curvature change are currently used to approximate treatment effect. The ability to readily and directly measure mechanical properties would improve our knowledge of disease progression and our ability to assess treatments that could reduce residual refractive error. 

Although there is growing interest in measuring cornea biomechanics to better understand ectasias and related treatments, even the most promising clinical techniques for corneal biomechanical measurement have notable limitations. For instance, optical coherence elastography (OCE) assesses mechanical properties by mechanically perturbing the cornea with a translating lens,¹⁶ acoustic vibration source,^17,18 or air-puff¹⁹ while imaging the resulting tissue response with optical coherence tomography (OCT). OCE provides spatially resolved two-dimensional (2D) or 3D maps of corneal displacement. However, all these methods require some degree of patient contact or perturbation, are influenced by IOP,²⁰ and the measured quantities are displacements that ultimately depend on analytical simplifications or complex inverse finite element methods to be related to mechanical property estimates. 

Another investigational methodology is Brillouin microscopy, which detects frequency shifts of light due to inelastic scattering in the tissue that correlate to local stiffness.²¹ This method can produce spatially resolved, 3D Brillouin-shift images of the cornea without the need for corneal contact or macro-deformation.²² However, Brillouin systems are large and expensive, do not leverage current clinical imaging technology, and involve lengthy acquisition times that make clinical translation and large-scale adoption more difficult. 

Finally, the ocular response analyzer (ORA; Reichert, Inc., Buffalo, NY, USA) and the Corvis STL (Oculus, Wetzlar, Germany) are commercially available devices that monitor the corneal deformation and recovery responses to a metered air puff with either infrared surface reflections or a single cross-sectional Scheimpflug image, respectively, to infer corneal properties.^23,24 These methods are attractive because they are readily available, inexpensive, and offer very short acquisition times (20 to 30 ms). However, these approaches are limited when measuring corneal properties due to confounders such as IOP²⁵ and contributions from extracorneal structures like the sclera,²⁶ as well as because the relationship between various response metrics and traditional biomechanical properties is not yet worked out. Most importantly, these methods assess bulk responses and do not provide spatially resolved information across or through the cornea. For this reason, ORA and CorvisST may underperform in ectasia risk screening applications where higher sensitivity is required to detect early, localized corneal biomechanical abnormalities.^27,28 

It appears that there is still an unmet clinical need for a high-resolution, high-speed, and noncontact method of determining mechanical properties of the whole cornea. In this study, we present a new method, which we call phase-decorrelation OCT (PhD-OCT), to determine corneal properties associated with crosslinking. OCT is well established as a useful imaging tool in the field of ophthalmology. PhD-OCT requires no modification to the hardware of existing clinical ophthalmic anterior-segment OCT systems. This paper provides introductory characterization of the PhD-OCT measurement and its sensitivity to suspected confounding factors that arise from both the imaging system and the tissue sample. Additionally, this work studied the feasibility of PhD-OCT for differentiating the crosslinking state of the cornea, as well as feasibility for in vivo clinical use. 

Methods

OCT Devices

To acquire the needed OCT data for phase decorrelation processing, an M-B scan pattern is used. Our M-B scans consist of M-scans captured in one spatial location for approximately 10 ms before moving to the next location in the cross-sectional B-scan. Two different OCT systems were used in this study. The first is a custom-made, spectral-domain, benchtop system with a central wavelength of 1310 nm, an A-line camera rate of 47 kHz (Sensors Unlimited, Princeton, NJ, USA), and an axial and lateral resolution of 15 and 20 μm, respectively. Further details on the linear-in-k OCT spectrometer are presented in a prior publication.²⁹ A second scanner, (Bioptigen Envisu C-class; Leica, Wetzlar, Germany), FDA approved for patient use, was used to image patients at the Cole Eye Institute at the Cleveland Clinic. This is also a spectral-domain OCT system, but with a central wavelength of 840 nm and an A-line rate of 35 kHz. 

PhD-OCT

The principle of PhD-OCT is based on the theory of dynamic light scattering (DLS) by particles undergoing Brownian motion.³⁰ When laser light illuminates randomly moving particles suspended in fluid, the amplitude and phase of the scattered light field evolve in a stable manner over time. A time-dependent complex-valued autocorrelation of the signal, called g⁽¹⁾, follows an exponential decay (Equation 1).³¹ The decay constant Γ is directly related to the Brownian diffusion coefficient, D (Equation 2).³¹ Γ is the metric reported for tissues by PhD-OCT. If the size and shape of the scattering particles are known, Γ may be directly linked to the viscosity of the medium by the Stokes-Einstein equation (Equation 3).³²  

\(\def\upalpha{\unicode[Times]{x3B1}}\)\(\def\upbeta{\unicode[Times]{x3B2}}\)\(\def\upgamma{\unicode[Times]{x3B3}}\)\(\def\updelta{\unicode[Times]{x3B4}}\)\(\def\upvarepsilon{\unicode[Times]{x3B5}}\)\(\def\upzeta{\unicode[Times]{x3B6}}\)\(\def\upeta{\unicode[Times]{x3B7}}\)\(\def\uptheta{\unicode[Times]{x3B8}}\)\(\def\upiota{\unicode[Times]{x3B9}}\)\(\def\upkappa{\unicode[Times]{x3BA}}\)\(\def\uplambda{\unicode[Times]{x3BB}}\)\(\def\upmu{\unicode[Times]{x3BC}}\)\(\def\upnu{\unicode[Times]{x3BD}}\)\(\def\upxi{\unicode[Times]{x3BE}}\)\(\def\upomicron{\unicode[Times]{x3BF}}\)\(\def\uppi{\unicode[Times]{x3C0}}\)\(\def\uprho{\unicode[Times]{x3C1}}\)\(\def\upsigma{\unicode[Times]{x3C3}}\)\(\def\uptau{\unicode[Times]{x3C4}}\)\(\def\upupsilon{\unicode[Times]{x3C5}}\)\(\def\upphi{\unicode[Times]{x3C6}}\)\(\def\upchi{\unicode[Times]{x3C7}}\)\(\def\uppsy{\unicode[Times]{x3C8}}\)\(\def\upomega{\unicode[Times]{x3C9}}\)\(\def\bialpha{\boldsymbol{\alpha}}\)\(\def\bibeta{\boldsymbol{\beta}}\)\(\def\bigamma{\boldsymbol{\gamma}}\)\(\def\bidelta{\boldsymbol{\delta}}\)\(\def\bivarepsilon{\boldsymbol{\varepsilon}}\)\(\def\bizeta{\boldsymbol{\zeta}}\)\(\def\bieta{\boldsymbol{\eta}}\)\(\def\bitheta{\boldsymbol{\theta}}\)\(\def\biiota{\boldsymbol{\iota}}\)\(\def\bikappa{\boldsymbol{\kappa}}\)\(\def\bilambda{\boldsymbol{\lambda}}\)\(\def\bimu{\boldsymbol{\mu}}\)\(\def\binu{\boldsymbol{\nu}}\)\(\def\bixi{\boldsymbol{\xi}}\)\(\def\biomicron{\boldsymbol{\micron}}\)\(\def\bipi{\boldsymbol{\pi}}\)\(\def\birho{\boldsymbol{\rho}}\)\(\def\bisigma{\boldsymbol{\sigma}}\)\(\def\bitau{\boldsymbol{\tau}}\)\(\def\biupsilon{\boldsymbol{\upsilon}}\)\(\def\biphi{\boldsymbol{\phi}}\)\(\def\bichi{\boldsymbol{\chi}}\)\(\def\bipsy{\boldsymbol{\psy}}\)\(\def\biomega{\boldsymbol{\omega}}\)\(\def\bupalpha{\unicode[Times]{x1D6C2}}\)\(\def\bupbeta{\unicode[Times]{x1D6C3}}\)\(\def\bupgamma{\unicode[Times]{x1D6C4}}\)\(\def\bupdelta{\unicode[Times]{x1D6C5}}\)\(\def\bupepsilon{\unicode[Times]{x1D6C6}}\)\(\def\bupvarepsilon{\unicode[Times]{x1D6DC}}\)\(\def\bupzeta{\unicode[Times]{x1D6C7}}\)\(\def\bupeta{\unicode[Times]{x1D6C8}}\)\(\def\buptheta{\unicode[Times]{x1D6C9}}\)\(\def\bupiota{\unicode[Times]{x1D6CA}}\)\(\def\bupkappa{\unicode[Times]{x1D6CB}}\)\(\def\buplambda{\unicode[Times]{x1D6CC}}\)\(\def\bupmu{\unicode[Times]{x1D6CD}}\)\(\def\bupnu{\unicode[Times]{x1D6CE}}\)\(\def\bupxi{\unicode[Times]{x1D6CF}}\)\(\def\bupomicron{\unicode[Times]{x1D6D0}}\)\(\def\buppi{\unicode[Times]{x1D6D1}}\)\(\def\buprho{\unicode[Times]{x1D6D2}}\)\(\def\bupsigma{\unicode[Times]{x1D6D4}}\)\(\def\buptau{\unicode[Times]{x1D6D5}}\)\(\def\bupupsilon{\unicode[Times]{x1D6D6}}\)\(\def\bupphi{\unicode[Times]{x1D6D7}}\)\(\def\bupchi{\unicode[Times]{x1D6D8}}\)\(\def\buppsy{\unicode[Times]{x1D6D9}}\)\(\def\bupomega{\unicode[Times]{x1D6DA}}\)\(\def\bupvartheta{\unicode[Times]{x1D6DD}}\)\(\def\bGamma{\bf{\Gamma}}\)\(\def\bDelta{\bf{\Delta}}\)\(\def\bTheta{\bf{\Theta}}\)\(\def\bLambda{\bf{\Lambda}}\)\(\def\bXi{\bf{\Xi}}\)\(\def\bPi{\bf{\Pi}}\)\(\def\bSigma{\bf{\Sigma}}\)\(\def\bUpsilon{\bf{\Upsilon}}\)\(\def\bPhi{\bf{\Phi}}\)\(\def\bPsi{\bf{\Psi}}\)\(\def\bOmega{\bf{\Omega}}\)\(\def\iGamma{\unicode[Times]{x1D6E4}}\)\(\def\iDelta{\unicode[Times]{x1D6E5}}\)\(\def\iTheta{\unicode[Times]{x1D6E9}}\)\(\def\iLambda{\unicode[Times]{x1D6EC}}\)\(\def\iXi{\unicode[Times]{x1D6EF}}\)\(\def\iPi{\unicode[Times]{x1D6F1}}\)\(\def\iSigma{\unicode[Times]{x1D6F4}}\)\(\def\iUpsilon{\unicode[Times]{x1D6F6}}\)\(\def\iPhi{\unicode[Times]{x1D6F7}}\)\(\def\iPsi{\unicode[Times]{x1D6F9}}\)\(\def\iOmega{\unicode[Times]{x1D6FA}}\)\(\def\biGamma{\unicode[Times]{x1D71E}}\)\(\def\biDelta{\unicode[Times]{x1D71F}}\)\(\def\biTheta{\unicode[Times]{x1D723}}\)\(\def\biLambda{\unicode[Times]{x1D726}}\)\(\def\biXi{\unicode[Times]{x1D729}}\)\(\def\biPi{\unicode[Times]{x1D72B}}\)\(\def\biSigma{\unicode[Times]{x1D72E}}\)\(\def\biUpsilon{\unicode[Times]{x1D730}}\)\(\def\biPhi{\unicode[Times]{x1D731}}\)\(\def\biPsi{\unicode[Times]{x1D733}}\)\(\def\biOmega{\unicode[Times]{x1D734}}\)\begin{equation}\tag{1}{g^{\left( 1 \right)}}\left( \tau \right) = \left\langle {{{{{\left\langle {E\left( t \right){E^*}\left( {t + \tau } \right)} \right\rangle }_{pixels}}} \over {\sqrt {{{\left\langle {E\left( t \right){E^*}\left( t \right)} \right\rangle }_{pixels}}} \sqrt {{{\left\langle {E\left( {t + \tau } \right){E^*}\left( {t + \tau } \right)} \right\rangle }_{pixels}}} }}} \right\rangle = {e^{ - \Gamma \cdot \tau }} \approx 1 - \Gamma \cdot \tau \end{equation}

 

\begin{equation}\tag{2}\Gamma = {\left( {{{4 \pi n} \over \lambda }} \right)^2}D\end{equation}

 

\begin{equation}\tag{3}D = {{{k_B}T} \over {6 \pi \eta R}}\end{equation}

 

In these equations, E is the electric field, Γ is the decorrelation coefficient, n is index of refraction of the solution (assumed to be 1.333 for water unless noted otherwise), λ is the wavelength of the light, D is the diffusion coefficient, k_B is Boltzman's constant, T is temperature, η is the viscosity of the solution, and R is the radius of the spherical scattering particles. 

Successive publications in the field of DLS have made progress toward studying biological tissues and using OCT to take advantage of its high contrast and spatial resolution. In 1998, Boas et al. used low-coherence interferometry to spatially resolve Brownian motion.³³ In 2010, Kalkman et al. furthered this objective by mapping the diffusion coefficient in two dimensions with spectral-domain OCT.³⁴ In 2012, Lee et al. demonstrated the full integration of DLS theory and OCT by deriving a field autocorrelation function and producing high-resolution, 3D maps of flow and diffusion within blood vessels in vivo.³⁵ 

Other researchers have used similar approaches of using the decorrelation of real- or complex-valued OCT data to measure tissue properties. Chhetri et al. and Blackmon et al. published methods of using scattering from exogenous particles to obtain absolute measurements of tissue properties based on diffusion.^36,37 Other work has used endogenous scattering by tissue, such as holographic tissue dynamics spectroscopy,³⁸ field-based dynamic light scattering microscopy,³⁹ and complex differential variance for tissue thermal ablation.⁴⁰ All of this work supports the feasibility of OCT with motility contrast as a tool to measure tissue properties. However, clinical mobility-contrast imaging of the cornea, a relatively stable tissue, imposes additional challenges in requiring a short acquisition time and robustness in regard to patient motion. No methods in the literature have addressed these challenges. 

PhD-OCT differs from prior techniques in two ways: (1) PhD-OCT interrogates a much shorter maximum time lag, just 0.3 ms, than previous techniques, and (2) the decorrelation is fit to the first-order Taylor expansion of the exponential decay. To start, OCT data are acquired in the form of an M-B scan. This means that the imaging beam is allowed to remain in place for 10 ms (meanwhile acquiring an M-scan of 500 A-Lines) before moving to the next lateral position. These M-scans are acquired in as many lateral positions as required to construct the desired B-scan cross section. The raw spectral data for each A-Line has the nonvarying spectral background signal subtracted and is then Fourier transformed. The Fourier transform of a real-valued signal yields a complex-valued resulting signal. To calculate each PhD-OCT data point in the cross section, the corresponding M-scan is analyzed. Starting at the first A-line in the M-scan, the complex decorrelation curve (Equation 1) of six axially adjacent pixels is calculated over 15 A-lines. The number of adjacent A-lines analyzed in the correlation is referred to as “maximum time lag.” This process is repeated starting at the 2nd, 3rd,…, 485th A-line in the M-scan, resulting in 485 14-point curves that may be used to calculate the slope, Γ. The average slope, Γ, of all these decorrelation curves is calculated. In the case of fluid samples with known scatterer diameter, Equation 2 is used to calculate the diffusion coefficient, and Equation 3 is used to calculate the viscosity of the solution. In Figure 1A, the exponentially decaying correlation measured in a fluid sample undergoing Brownian motion is shown. In Figure 1B, the computation of the decay coefficient by a linear fit, in the form of y = −Γx + B, to the initial part of the decay curve is shown. 

Figure 1

 

Phase decorrelation analysis. A and B show the g¹(t) decorrelation (Equation 1) taken from a 3-mL sample of 25% glycerol (w/v) solution with 0.025% v/v 0.5-μm-diameter polystyrene beads (Magsphere, Inc., Pasadena, CA, USA) The data were acquired in an M-scan of 5000 A-lines. To calculate a single g¹(t) curve, the correlation of six adjacent pixels within one A-line was computed as a function of time across successive A-lines. The nth point in the curve is the median correlation value between all pairs of A-lines separated in time by n/(A-line rate) seconds. To produce A and B, 20 g¹(t) curves from the same sample were calculated, and the mean and SD of all 20 curves was calculated. (A) g¹(t) is shown from a relative time lag of 0 (1 A-line) to 53 ms (2499 A-lines). The data are fit to an exponential function. (B) The same data are shown from a relative time lag of 0 (1 A-line) to 0.3 ms (15 A-lines). The data are fit to a linear function in the form y = −Γx + B.

In the case of the cornea, we hypothesize that Γ will be inversely related to the degree of collagen confinement. In a cornea treated with crosslinking therapy, the collagen becomes more confined,⁴¹ and a corresponding decrease in random motion is expected. Thus, we hypothesize that increased crosslinking of corneal tissue should be related to a decrease in Γ. This is investigated in this study. Furthermore, in a healthy cornea, the collagen is highly confined by the normal microstructural arrangement,⁴² whereas in a keratoconic cornea, the collagen is less confined.^43,44 Thus, we hypothesize that PhD-OCT may be able to detect early changes in collagen microstructure associated with keratoconus. This will be investigated in future studies. 

Glycerol and Gelatin Phantoms

The first experiments were undertaken to verify the basic premise of the technique: that is, as a derivation of DLS principles, PhD-OCT could accurately assess Brownian dynamics of samples. Engineered phantoms were created to characterize the range of sensitivity of PhD-OCT to random small motions in fluid and gel. To assess fluid samples, 1 mL glycerol phantoms was created with between 0% and 70% glycerol in water (v/v), with 0.025% (w/v) polystyrene beads (0.5 μm diameter, PSF-500; Magsphere, Inc., Pasadena, CA, USA) and maintained at room temperature. Each sample was imaged with 10 M-scans, each of 5000 A-lines, and spaced 200 μm apart. Every M-scan was considered as an independent measurement. With the assumptions that the particles are spherical and uniform, the viscosity of a given sample could be calculated using Equations 1 to 3, and then the decorrelation values were averaged within each sample. 

The gelatin phantoms were manufactured by heating mass/volume proportions of water and gelatin from porcine skin (G-2500 Type A; Sigma-Aldrich Corp., St. Louis, MO, USA) until dissolved and then adding 0.025% w/v of the same polystyrene beads and casting 1-cm-diameter cylinders. Samples were maintained at 4°C until imaging. Each sample was imaged with 10 M-scans (1 M-scan is 5000 A-lines), each spaced 100 μm apart. Every M-scan was considered an independent measurement. Γ was calculated using Equations 1 and 2, and then the values were averaged in depth within each sample. 

A second group of experiments aimed to characterize the dependence of PhD-OCT on the density of scatterers. This is needed because the development of the speckle field is necessary to DLS. This was accomplished by varying the concentration of polystyrene beads over two orders of magnitude from a w/v bead concentration from 0.0038% to 0.22%. From the images, corresponding scattering coefficients of 0.13 and 0.18 mm⁻¹ were calculated. This range was selected to approximate variability that may be seen in corneas, as bovine corneas have been reported to have a nominal scattering coefficient of 0.191 mm⁻¹ at 1310 nm.⁴⁵ The beads were suspended in a 25% glycerol solution. Each sample was measured in 10 different locations to reduce noise. 

Porcine Cornea Preparation

Whole porcine globes were obtained from a local abattoir, delivered via courier in a cooler, and used within 6 hours postmortem. Upon arrival, they were transferred to a 4°C refrigerator. Before use, the cornea was de-epithelialized with the back of a scalpel, and Optisol GS (Chiron Ophthalmics, Irvine, CA, USA) was applied to the anterior surface continuously for 30 minutes at room temperature to normalize hydration. Then, the globe was examined under the OCT imaging system and observed for signs of preexisting abnormalities. For the experiments described below, any cornea with abnormalities apparent under visual or OCT inspection was rejected. Finally, a 21-gauge needle was inserted into the posterior segment. The needle was attached to a gravity-fed pressure system with an adjustable-height 500-mL saline reservoir and digital in-line pressure transducer attached to a continuous pressure monitor (Siemens SC 7000; Siemens, Berlin, Germany). Before each experiment, all tubes were checked for air bubbles, the pressure monitor was zeroed to atmospheric pressure, and the IOP was increased to 10 mm Hg. 

Corneal Imaging

Unless otherwise specified, all corneal cross sections are the result of a 10-mm M-B scan. Each scan consists of 100 evenly spaced M-scans of 1000 A-lines. Regions of specular reflection and low signal strength were avoided or omitted when calculating the average Γ value within each corneal cross section. All B-scans and M-B-scans of corneas were taken parallel to the vertical meridian, offset laterally from the corneal apex by approximately 1 mm to avoid specular reflection. 

Signal-to-Noise Ratio Experiment

In OCT systems, sources of noise include thermal, shot, and quantization noise from at the camera, as well as noise from the laser light source and galvanometer jitter and other vibrations. For our initial experiments in PhD-OCT, we were most concerned with the effect of noise from the camera. This is because in imaged regions with low signal intensity, this noise may overwhelm the phase decorrelation analysis. To examine this, three porcine globes were prepared as described. To investigate the effect of signal-to-noise ratio (SNR) on Γ, neutral density (ND) filters of varying optical densities were placed in the sample arm of the benchtop OCT system. One M-B scan was acquired for each globe for each ND filter. To measure the signal intensity, the percent change in intensity value from the segmented sample region of each cross-sectional image was taken and compared with the case when no ND filter was in place. Average Γ from each cross-sectional image was measured, and the percent change from the no-filter case was calculated. 

Hydration Experiment

Although corneal hydration is well controlled in healthy patients, it is of interest to observe the effects of hydration on the PhD-OCT measurement, because in ex vivo experiments, the hydration must be externally controlled. To vary the hydration, six porcine globes were prepared as described previously. Photrexra Viscous (20% dextran and 1.46 mg/mL riboflavin 5′-phosphate; Avedro, Waltham, MA, USA), a solution used clinically to titrate hydration of corneas during the CXL procedure, was applied dropwise to cover the corneal surface and replenished at 5-minute intervals for 30 minutes. Approximately 2 minutes after every application of dextran, OCT and PhD-OCT images were acquired of the cornea. From the OCT reflectance images, the central corneal thicknesses for each cross-sectional image were calculated with a custom automatic segmentation script (MATLAB 2016a; The MathWorks, Natick, MA, USA). The central corneal thickness measurement is taken as a surrogate metric for hydration state. From the PhD-OCT images, average Γ from each cross-sectional image was measured, and the percent change between the t = 0 image and each successive image was calculated. 

IOP Experiment

To examine the effect of elevated IOP on the corneal PhD-OCT images, six porcine globes were imaged. Each globe was prepared as described above, except the IOP was progressively increased from 5 to 25 mm Hg in 5-mm Hg steps. After each increase in IOP, the globe was allowed to equilibrate for 2 to 5 minutes as pressure was maintained before being imaged. Average Γ from each cross-sectional image was measured. 

Motion

To understand the artifacts that may be created by a patient's movement during clinical imaging, controlled motion experiments were conducted. 

Lateral

Three porcine globes were sequentially placed on a motorized, calibrated translation stage (T-LSR300A; Zaber Technologies, Inc., Vancouver, Canada) under the OCT sample arm. Each cornea was applanated with a weighted and stabilized glass coverslip. IOP was then raised to 10 mm Hg. As the sample was imaged, the stage was translated laterally at a variety of speeds between 0 and 1 mm/s. Five M-scans of 5000 A-lines were acquired during the displacement. 

Axial

A porcine globe was placed in the imaging area and the IOP was then raised to 10 mm Hg. The central portion of the cornea was imaged (M-scan, 5000 A-lines) at 47 kHz, whereas the reference arm was moved manually back and forth. This manually applied motion captured a range of displacement speeds sufficient to analyze motion dependence from 0 to 1 mm/s, as measured by Doppler OCT. 

CXL

An experimental adaptation of the Dresden crosslinking protocol⁴⁶ was followed to induce crosslinking in the ex vivo porcine cornea samples. After preparation as described above, a 1.46 mg/mL riboflavin 5′-phosphate ophthalmic solution (Photrexra; Avedro, Waltham, MA, USA) was applied continuously to the corneal surface for 30 minutes in a dark environment. To conserve the clinical grade riboflavin, a 1-cm retaining ring was used to prevent excess riboflavin from immediately running off the cornea. The isotonic riboflavin formulation was used in favor of the hypertonic 20% dextran solution to avoid dehydrating the cornea. For the CXL treatment, the corneas were then exposed to 2.9-mW/cm² 360-nm UV light from a UV-A diode (PriaVision, Inc., Menlo Park, CA, USA) for 30 minutes while the riboflavin application was continued. For the sham treatment, the corneas were kept in the dark environment without UV exposure while riboflavin exposure was continued for 30 minutes. The experiments were conducted pairwise such that one sham-treated eye and one CXL-treated eye were completed concurrently. The same CXL protocol was used when monitoring the progress of crosslinking during UV exposure, with the exception that the UV source was briefly removed to allow the cornea to be imaged once every 10 minutes during UV exposure. For all experiments, the radiance was calibrated using a UV-A light meter (handheld UV-A meter, HHUVA1; Omega Engineering, Inc., Stamford, CT, USA). In calculating and displaying the data, only information from the anterior third of the cornea was used, as this is where most of the crosslinking takes place.⁴⁷ 

Clinical Imaging

Patients undergoing either CXL or LASIK were recruited and imaged either pre- or postoperatively as an adjunct to the laboratory experiments to pilot the use of the approach in a clinical population of interest. Patients were imaged with a commercially available anterior segment scanner (Bioptigen Envisu C-class; Leica, Wetzlar, Germany). The left and right eyes of each patient were scanned while the imaging probe was firmly clamped to a slit-lamp stand in a normal examination room using the 12-mm telecentric objective lens for anterior imaging (Fig. 2). Patients were asked to rest their heads in a chin-rest and head-strap to reduce motion. The research was approved by the Cleveland Clinic institutional review board (IRB #13-213) and followed the tenets of the Declaration of Helsinki. Informed consent was obtained from the subjects after explanation of the nature and possible consequences of the study. Each scan consisted of 100 M-scans of 300 A-lines each, spaced over 10 mm, for a total acquisition time of just under 3 seconds. 

Figure 2

 

Photograph of the Bioptigen Envisu (C-class; Leica, Wetzlar, Germany) OCT scanner using the 12-mm telecentric objective lens for anterior imaging configured to acquire clinical data for the PhD-OCT analysis.

Image Processing and Analysis

After acquisition, data were processed with a custom-designed script (MATLAB 2016a; The MathWorks). Portions of this script related to phase-decorrelation processing are provided in the Supplementary Materials S1. This script reads in the data and calculates the complex decorrelation coefficient, g⁽¹⁾ (Equation 1), which is then used to calculate Γ for each row of the M-scans constituting the M-B scans. For corneal samples, the script also includes automatic segmentation algorithms to select the corneal stroma from the image, as well as motion- and intensity-correcting algorithms, as described below. These algorithms were defined based on the results of the experiments described above. 

Correcting for Motion

In the case of axial motion, the axial velocity is readily found by computing the Doppler OCT image from acquired M-scans. However, OCT is not sensitive to motion in the lateral direction, yet lateral motion also impacts the PhD-OCT measurement. To address this, the correction algorithm assumes that, over a few milliseconds, a patient's involuntary motion may be regarded as a linear trajectory in 3D space. In this case, the ratio between axial and lateral motion is constant and the axial motion may be used as a proxy measure of the total velocity. Therefore, the clinically obtained images are dynamically self-corrected for motion artifacts by looking for linear correlations between the Doppler OCT signal and Γ over successive milliseconds-long blocks of time, and linearly corrected as needed. Additionally, to address large amount of bulk motion, the algorithm corrects for discontinuities in the M-B scan cornea profile by applying multipixel shifts to each column of the image to align the corneal surface in each column to fall on the quadratic trend line of the uncorrected corneal surface profile. 

Correcting for SNR

To account for the effect of SNR, the experimental results for SNR dependence in corneal images are fit with a second-order polynomial, which is then used to correct for the anticipated decrease of Γ in low-SNR regions of corneal images. 

Results

In Figure 3A, measured viscosity of various glycerol-water mixtures (N = 5) with 0.025% polystyrene beads are in excellent agreement with expected absolute viscosity values.⁴⁸ Each measurement consists of 5000 A-lines. Viscosity was calculated from the first 0.2 ms of each decorrelation curve. In Figure 3B, hydrogel phantoms were prepared with varying concentration of porcine gelatin (Sigma-Aldrich Corp.). For contrast, 0.025% w/v beads, as described previously, were added. Each sample was imaged 10 times. The expected linear trend is evident (r² = 0.97). These data support the idea that PhD-OCT is sensitive to small changes in sample material properties. 

Figure 3

 

Glycerol phantoms. Controlled phantoms were used to preliminarily characterize the phase decorrelation computation. (A) Measured absolute viscosity of various glycerol-water mixtures with 0.025% polystyrene plotted with expected values⁴⁸ (dashed line). Values are plotted as mean ± SD (N = 5). (B) Gelatin phantoms. Varying concentrations of gelatin phantoms were imaged. The mean and SD Γ values are plotted with a linear fit (r² = 0.971).

The second group of experiments aimed to characterize important possible confounding factors, such as IOP and hydration, specific to the cornea. As needed, each of these factors could be corrected for or ignored. In Figure 4, the concentration of beads in a 25% glycerol solution was varied to mimic tissues of varying scattering coefficient and speckle development. On the left axis, D is shown to trend toward an asymptote as the bead concentration increases (dotted line is exponential trend line, r² = 0.82). On the right axis, the scattering coefficient is shown to proportionally increase with logscale bead concentration (dotted line is linear trend line, r² = 0.99). In corneas, the scattering density is high and similar between samples, so small changes should not confound the measurement. 

Figure 4

 

Apparent diffusion coefficient as a function of scatterer concentration. The concentration of beads in a 25% glycerol solution was varied. The observed diffusion coefficient is shown to trend toward an asymptote as the bead concentration increases. (Dotted line is exponential trend line, r² = 0.82.) On the right axis (red), the scattering coefficient is shown to proportionally increase with logscale bead concentration (dotted line is linear trend line, r² = 0.99). Representative B-scan sections of varying bead concentrations are inset to show speckle development.

SNR was found to be a confounding factor in the decorrelation analysis of cornea samples (Fig. 5A). However, the trend (r² = 0.996) between SNR and Γ may be corrected for because the SNR at each point of the image is known. The relationship between SNR and Γ, in which Γ decreases with decreasing SNR, may be explained by the fact that a perfectly random signal has an autocorrelation of 0 at all times. Therefore, increasing noise in a signal will “flatten” the autocorrelation curve and reduce the calculated Γ. 

Figure 5

 

Characterization of confounding factors in cornea samples. (A) As the incident light intensity on the sample is reduced, decreases due to the decrease in SNR. A quadratic trend line was fit to the data (n = 3; r² = 0.94). (B) As hydration is varied, with change in central corneal thickness taken as a surrogate measure, Γ tends to decrease (r² = 0.512). (C) As IOP of the porcine globe is varying across a physiologic range, the effect (r² = 0.004) on Γ is observed to be small compared with other variations in the measurement.

It has been previously reported that the hydration of the cornea has an effect on its biomechanical properties. Some studies^49–53 have detected a decrease in stiffness with increasing hydration, whereas other studies^54,55 have the opposite effect. These discrepancies could be caused by a combination of (1) not fully disambiguating the effects of corneal thickness from the corneal stiffness measurement and (2) the effect of probing a viscoelastic tissue in different spatio-temporal regimes. PhD-OCT shows a decrease (r² = 0.47) in Γ as the cornea dehydrates and the central corneal thickness decreases (Fig. 5B). Typically, a decrease in Γ is associated with a stiffer tissue. 

Many other methods of measuring corneal biomechanics have a dependence on the IOP because IOP opposes deforming forces alongside corneal stiffness. However, when IOP was varied over a physiologic range, the correlation in Γ was very small (r² = 0.49) compared with other variations in the measurements and natural variance between eyes (Fig. 5C). 

The effect of lateral and axial motion was investigated in Figure 6 because detecting or correcting for motion is a critical step to imaging human subjects, where bulk motion may appear as artifact in the PhD-OCT measurement of Γ. As the lateral displacement speed of the cornea sample increased (Fig. 6A), Γ increased as well. Additionally, the effect of this increase is amplified when longer maximum lag times are considered. Three samples were translated at 0.1 mm/s, 0.5 mm/s, and 1 mm/s while 5000 A-lines were acquired. The mean and SD of the average Γ within the cornea are shown. In Figure 6B, purely axial motion was obtained by continuously translating the reference mirror at a variable speed as 5000 A-lines were acquired. Γ was calculated with 1.5, 2.3, and 3.2 ms of maximum lag. A second-order polynomial was fit to Γ versus axial speed for each maximum lag. 

Figure 6

 

Effect of motion on Γ. (A) Increased lateral motion from a translation stage increases the induced motion artifact observed in porcine cornea (N = 3), as does increased maximum lag time. The mean and SD of the average Γ within the cornea are shown. Points are staggered for visualization. (B) Purely axial motion was obtained by continuously translating the reference mirror while imaging a porcine cornea (N = 1). This motion, along with lag time, corresponds to increased measured D. The mean and SD of the average Γ within the cornea are shown. Points are staggered for visualization.

The crosslinking protocol as described above was applied to ex vivo porcine eyes to demonstrate feasibility of monitoring the procedure in human patients. These experiments yielded clear and expected results. In Figures 7A and 7B, examples are shown of a cornea before and after the CXL protocol, as well as a cornea before and after the sham treatment, respectively. White arrows indicate the location of the demarcation line, which is visible in OCT reflectance images and is associated with the edge of the acellular region created by CXL.⁵⁶ In the PhD-OCT images, the same region corresponds to a band of decreased Γ. The sham treatment exhibits no demarcation line or band of decreased Γ. The results from the whole experiment, with five eyes in each group, are summarized in Figure 7C, where there is a significant difference (P < 0.0001) between the CXL-treated corneas and the sham-treated corneas. In Figure 7D, there is a clear relationship between the increasing UV irradiance of the cornea and the decreasing average Γ, suggesting that PhD-OCT has sufficient sensitivity to monitor CXL treatment in real time, potentially allowing feedback for clinicians to titrate UV irradiance to effect. 

Figure 7

 

Corneal crosslinking of ex vivo porcine cornea. Examples of OCT reflectance images (upper images) and PhD-OCT images (lower images) of ex vivo porcine corneas before and after (A) a crosslinking procedure and (B) a sham procedure. The black dotted line on the PhD-OCT images indicates the segmentation of the anterior third of the cornea. Field of view is 4.5 mm lateral by 1.2 mm axial. (C) The average Γ within the anterior third of the corneas was computed for samples before and after CXL (n = 5) and sham (n = 5) procedures. The relative postprocedure change in Γ was calculated for each sample and plotted. The two groups are significantly different (P < 0.001). (D) For three additional ex vivo porcine globes, corneal imaging was conducted during the crosslinking procedure, at the start of UV irradiation, and every 10 minutes thereafter. The average Γ of the anterior third is plotted. Note that in the post-CXL example, in the OCT reflectance image, the corneal demarcation line is highlighted with white arrows. This region corresponds to the region of decreased mobility in the PhD-OCT image.

Data gathered from the benchtop experiments described above was used to design correction algorithms to reduce artifacts in images acquired from human subjects. In Figure 8A, the top image shows the result of phase-decorrelation processing on raw data from a clinical system. Several artifacts are noticeable: (1) geometrical distortion from motion during imaging, (2) distortion of Γ due to motion, and (3) SNR artifact in the periphery of the cornea. All of these issues are mitigated by applying the corrections described above to the image, as seen in the bottom image. Figure 8B shows two examples of normal corneas and two examples of corneas 3 months after the patients underwent a crosslinking procedure. In these two corneas, a band of decreased D in the anterior third is clearly visible and very similar to those observed in the ex vivo porcine corneas (Fig. 6). 

Figure 8

 

Images acquired with a clinical OCT system. Field of view is 3.7 mm lateral by 1 mm axial. (A) To generate interpretable images, corrections were applied to account for SNR and patient motion. An example of a clinically obtained image before and after corrections have been applied is shown. (B) Examples of four unrelated scans are shown. Two pre-LASIK patients were imaged and found to have relatively homogenous corneas. However, two patients at 3 months after CXL exhibit patterns of decreased Γ in the anterior third, very similar to the ex vivo porcine eyes post-CXL.

Discussion

These results characterize the PhD-OCT method for corneal imaging, including imaging-dependent factors, such as SNR, and physiologic factors, such as hydration. The presented images of treated corneas demonstrate clinical feasibility and utility, where PhD-OCT has a number of advantages over existing technologies. First, spatial resolution (both axial and lateral) of material properties of soft tissue has been demonstrated. This is important because biomechanical abnormalities in corneal degenerations such as keratoconus are not evenly distributed in space.^{9,10,12,22,57} This can limit the sensitivity of current measurement methods for early diagnosis. Mapping of cornea biomechanical properties also represents an opportunity to tailor a patient's crosslinking treatment to the most affected areas of the cornea. It may also be useful for observing localized surgical effects and wound healing in the cornea. Second, no perturbation of the cornea is required. With this come the benefits of decreased instrument and measurement complexity, increased patient comfort, and absence of dependence on intraocular pressure. Third, the acquisition time is short, allowing the scan to fit into clinical workflow. Fourth, because this method can be implemented using existing clinical OCT systems, without additional hardware, it may be relatively straightforward to translate to clinical use compared with novel imaging systems. 

A key to the PhD-OCT method is the analysis of correlation decay over very short time frames. There are several benefits to this approach. First, the total acquisition time decreases, which is particularly desirable for clinical imaging. Second, it allows the simplification of the parameter fit to a line by Taylor expansion, as opposed to the ill-conditioned problem of fitting an exponential decay with an offset. Numerical simulations (data not shown) demonstrated that the linear fit results in a narrower parameter distribution and improved accuracy than exponential decay given the same acquired data. Finally, the shorter interrogation time yields a result less affected by motion (Fig. 7). 

Figures 3A and 3B demonstrate the large dynamic range of PhD-OCT measurements across fluid and hydrogel samples. In the case of fluid samples, the results corresponded very closely with theoretical values, with no calibration performed. Figure 4 suggests that it is valid to compare measurements from the same tissue type, even if scatterer properties vary, as long as the speckle field is sufficiently developed, as is the case in corneas. 

The preliminary investigation of confounding factors presented in Figures 5 and 7 conform to expectations. In Figure 5A, SNR is shown to decrease the apparent Γ. This result may seem counterintuitive compared with the effect of noise in similar OCT techniques, such as complex differential variance OCT,⁵⁸ but is expected because the addition of random noise will tend to flatten the slope of the g⁽¹⁾ curve. In Figure 5B, the trend of dehydration making the cornea stiffer corresponds with other recent results in cornea biomechanics literature.⁵⁹ Additionally, the apparent lack of dependence on IOP (Fig. 5C) may be expected because the simple models of Brownian motion and dynamic light scattering, as discussed above, do not include pressure as a factor in determining random displacements. 

Last, Figure 6 is a characterization of motion artifacts. As is evident in Figure 6A, increasing lateral motion causes a linear increase in the apparent Γ, which is expected because any motion should result in a decrease in the time-dependent autocorrelation curve and thus an increase in Γ. The same is true of axial motion, shown in Figure 7B. This characterization is important for the development of motion-correction algorithms whereby detected motion during the imaging time may be used to correct for motion-induced artifacts in the PhD-OCT measurement. Additionally, processing the data with different maximum lag times highlights the importance of lag time in reducing motion artifacts. In both cases, a longer maximum lag time means that motion has a larger influence on Γ. 

The results demonstrating sensitivity to corneal crosslinking (Fig. 7) indicate that the method is capable of distinguishing crosslinked eyes from noncrosslinked eyes. Figure 7D further suggests that this sensitivity directly corresponds to degree of crosslinking. Additionally, Figure 7A raises questions about the correspondence of crosslinking depth to demarcation line. To the authors' knowledge, how the demarcation line correlates with actual crosslinking, or the etiology of the demarcation line itself, is not fully demonstrated in the literature. The relationship between the demarcation line and PhD-OCT could provide useful information about the microstructural properties of the stroma during crosslinking. 

The clinical images shown in Figure 8 are encouraging. Additional clinical data are needed to determine whether meaningful conclusions may be drawn from individual corneal images. We expect that PhD-OCT may have several points of impact in the clinic. First, PhD-OCT, as demonstrated here, may be able to give clinicians real-time feedback on CXL procedures to facilitate tailoring of procedures to achieve a desired effect. Second, because it is hypothesized^9,10 that biomechanical perturbations occur early on in the progression of corneal ectasias and before irreversible geometrical distortion is caused, it is plausible that PhD-OCT could be used to screen for early-stage keratoconus if sufficient sensitivity and specificity are demonstrated in future clinical studies. Groups that are particularly at risk could be screened and appropriate actions could be taken to correct or stabilize an ectasia before the patient's vision is impaired. Third, it could be more broadly used to characterize the corneal integrity of prospective candidates for LASIK or other refractive surgeries. Currently, spatially resolved corneal biomechanical property assessment is not available in the clinical workflow of determining whether a patient is suitable for refractive surgery.⁶⁰ The absence of such screening technology could contribute to the number of patients experiencing postoperative structural instability.⁶⁰ 

Acknowledgments

The authors thank the J.H. Routh company for providing samples and the staff of the Cleveland Clinic Cole Eye Institute for assistance with the clinical images. 

Research reported in this manuscript was supported by the National Eye Institute, the National Heart, Lung, and Blood Institute, and the National Institute of Biomedical Imaging and Bioengineering, of the National Institutes of Health under award numbers R01EY028667, R01HL083048, R01HL126747, T32EB007509, T32EY007157, C06RR12463, R01EY023381, IP30EY025585, an Unrestricted Grant Award from Research to Prevent Blindness to the Department of Ophthalmology at Cole Eye Institute (RPB1508DM), Foundation Fighting Blindness Center Grant to the Cole Eye Institute (CCMM08120584CCF). The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. 

Disclosure: B.J. Blackburn, P; S. Gu, P; M.R. Ford, None; V. de Stefano, None; M.W. Jenkins, P; W.J. Dupps Jr, Avedro Inc. (C), Alcon (C), P; A.M. Rollins, P