Reverberant 3D optical coherence elastography maps the elasticity of individual corneal layers

The elasticity mapping of individual layers in the cornea using non-destructive elastography techniques advances diagnosis and monitoring of ocular diseases and treatments in ophthalmology. However, transient Lamb waves, currently used in most dynamic optical coherence and ultrasound elastography techniques, diminish the translation of wave speed into shear/Young’s modulus. Here, we present reverberant 3D optical coherence elastography (Rev3D-OCE), a novel approach leveraging the physical properties of diffuse fields in detecting elasticity gradients not only in the lateral direction, but also along the depth axis of the cornea. A Monte Carlo analysis, finite element simulations, and experiments in layered phantoms are conducted to validate the technique and to characterize the axial elastography resolution. Experiments in ex vivo porcine cornea at different intraocular pressures reveal that Rev3D-OCE enables the elastic characterization of single layers that matches the anatomical description of corneal layers with unprecedented contrast in the dynamic OCE field.

Processing pipeline for the Rev3D-OCE method. Structural and motion volumes provided by the PhS-OCT system are used for detecting topographic surfaces ( , for = 0, 1, 2, … , 299) in cornea and projections of particle velocity OCT ( , ) on them. Then, temporal and spatial filtering is applied to each projection frame � OCT ( , ). A temporal Fourier transform is applied to each spatial location for the calculation of magnitude ( ) and phase ( ). Subsequently the complex matrix ( ) = �π/2 ( )e i ( ) is formed for the local calculation of auto-correlation in a 2 × 2 mm 2 window. The average auto-correlation profile is fitted to Equation 2 of the main manuscript for the local estimation of wave number * . Shear wave speed maps s ( ) = 0 / ( ) * are calculated for each surface in cornea by moving the window through the entire projection frame, and by fixing the a priori known excitation frequency 0 = 2π 0 , where 0 = 2 kHz used during experiments. Finally, the compounding of s ( ) calculated at each surface is represented in a speed color-coded 3D volume.

Theoretical derivations on reverberant shear wave fields
As stated in the manuscript, a reverberant field can be understood as the superposition of plane shear waves traveling in a random direction 1 . Shear waves are transversal body-type waves; therefore, the particle velocity produced by these perturbations is perpendicular to the direction of propagation. For a given 3D point in the Cartesian system, three orthogonal vectors in the spherical coordinate system are found: �, � , and � ( Figure 1a in the main manuscript). Then, the corresponding particle velocity vector field ( , ), at position and time , in a reverberant chamber produced by plane waves propagating with a wave number and radial frequency 0 is modeled as where the index represents a realization of the random unit vector � describing the direction of wave propagation, and the index represents a realization of the random unit vector � describing a direction of particle velocity parallel to the disk formed by the basis vectors � and � defined within a realization of . Then, � • � = 0. Finally, is an independent, identically distributed random variable describing the magnitude of particle velocity within a realization of . The summation over is understood to be taken over the 4π solid angle, and the summation over is taken over a 2π angle within the disk.
Typically, optical coherence tomography (OCT) systems measure particle velocity or displacement along one single axis that we call the sensor axis. For simplicity in the mathematical derivations, we choose the x-axis � to be the sensor axis. Then, ( , ) = ( , ) • � describes the particle velocity scalar field along � , and where = � • � is a scalar random variable. Then, using the spherical coordinate system (see definition of angles in Figure 1a in the main manuscript), we can represent � = cos( ) � + sin( ) � , where � = cos( ) cos( )� + cos( ) sin( )� − sin( ) � , and � = −sin( )� + cos( )� + 0 � . Here, is defined as the angle between � and �, defining a realization of � within the disk formed by the basis vectors � and �. Therefore, we find that = � • � = sin( ) cos( ) cos( ) − cos( ) sin( ). ( By taking the auto-correlation function of Equation 2, , in space and time, we obtain: where {•} represents an ensemble average and the asterisk represents conjugation. The product of the two series will include vanishing cross terms of the form: Since and are independent, and realizations of are uncorrelated, Equation 4b becomes: where = 〈 2 〉 is the expected value of 2 over both q and l realizations extracted from the curly braces since and { , � } are independent random variables. In an ideal diffuse field, the ensemble or spatial averaging will assign equal weighting to all directions � of incident shear waves in the 4π solid angle, and, given a realization of � , equal weithting in the direction of particle velocity � in the 2π angle within the disk formed by the basis vectors � and � ( Figure  1a in the main manuscript). Then, the average of the summation over discrete directions of each incident wave becomes the average over all directions in the polar coordinate system 2 . By substituting Equation 3 into Equation 6b and averaging in the polar coordinate system we have: where = sin( ) . For solving Equation 7, two cases are considered: the auto-correlation direction (1) perpendicular, and (2) parallel to the sensor axis. For case 1, we first consider Δ = Δ along the z-axis. Then, the exponential term � • Δ in Equation 7 becomes: where Δ cos ( ) is the projection of � in Δ using spherical coordinates. Next, solving the triple integral in Equation 7 gives: where 0 and 1 are spherical Bessel functions of the first kind of zero and first order, respectively. It is observed that the spatial and temporal component of (∆ , ∆ ) are separable and can be redefined as (∆ ) (∆ ). An identical result (when replacing the subscripts 'z' in Equation 9 with 'y') can be found when the y-axis is selected as the auto-correlation axis, Δ = Δ , since it is perpendicular to the sensor axis � . Finally, in order to be consistent with the main manuscript, we can redefine the sensor axis to be � and the correlation axis to be � or � and change the subscripts accordingly since the ideal diffused field is isotropic. Then, accounting for the axis redefinition and neglecting the temporal component of Equation 9 we obtain Equation 2 of the main manuscript.
For case 2, we redefine the sensor axis to be � and consider the auto-correlation direction Δ = Δ along the z-axis (parallel to the sensor axis) for simplicity in the mathematical derivations. Then, where = � • � is the projection of � in the new sensor axis � . Following the same derivation flow as in Equation 3, we find that = − sin( ) sin( ). Then, taking the autocorrelation of Equation 10 along the z-axis and using � • Δ = Δ cos ( ) in the exponential term, we have Solving the triple integral in Equation 11 gives: which is identical to Equation 3 of the main manuscript after neglecting the temporal component.

Spectral-domain phase-sensitive optical coherence tomography (PhS-OCT) system
The spectral-domain PhS-OCT used in the experiments was a custom-built system as shown in Supplementary Figure 3. The system was spectrometer-based with a light source coming from a superluminescent light emitting diode )EXS210045-01, EXALOS, Schlieren, Switzerland( with a central wavelength of 1307 nm and a full-width half-maximum )FWHM( spectral width of 100 nm. The average output power was 12 mW. The theoretical axial point-spread function )i.e. depth resolution limit( was computed to be approximately 8 µm. The interferometer (see Supplementary Figure 3a) was designed based on a free-space Michelson interferometer configuration. The output from the light source was delivered to the interferometer through a single mode optical fiber. At the interferometer, the input light was collimated to a beam size of about 7 mm. The collimated light beam was split by a beam-splitter cube )BS015, Thorlabs, Newton, NJ, USA(, having 50:50 power split ratio. Half of the power of the light beam was delivered to a sample arm, consisting of dual-axes galvanometric mirrors )GVS012, Thorlabs, USA(, an objective lens )LSM05, Thorlabs, USA(, and a sample holder. The objective lens has a measured lateral resolution of about 15 µm, a working distance of about 90 mm, and a maximum scanning field of view of about 25 mm × 25 mm. The objective lens focused a light beam onto a sample and collected back-scattered light from the sample. The collected back-scattered light was, then, delivered back to the beam splitter to be combined with the reference light beam. An additional half-power of the light beam was propagated along the reference beam path which includes a dispersion compensator (LSM05DC, Thorlabs, USA) and a retro-reflector prism (PS975-C, Thorlabs, USA) mounted on a length-adjustable lens tube (SM1NR1, Thorlabs, USA). The reference light beam was reflected back and combined with the sample beam at the beam splitter. The combined light beam was coupled to another single mode fiber and delivered to a custom built spectrometer 3 .
The spectrometer (see Supplementary Figure 3a) consisted of a reflective collimator, a reflective grating, a focusing mirror, and a line-array sensor. Light output at the other end of the fiber was placed at the focal point of a 90° off-axis parabolic (OAP) mirror with a reflective focal length of 101.6 mm (MPD249-M01, Thorlabs, USA) to produce a collimated beam of about 20 mm beam size. The collimated beam was directed to a reflective grating of 600 lines mm -1 (GR50-0613, Thorlabs, USA). A dispersed light beam was then focused by a 45° OAP mirror with a reflective focal length of 202.3 mm (MPD284-M01, Thorlabs, USA) to the line-array sensor. The line-array sensor was an InGaAs high speed line-scan camera, consisting of 2048 pixels of 10 µm pitch (2048R, Sensor Unlimited, Princeton, NJ, USA) 3 . The spectral interference signal captured by the spectrometer was transferred to a computer for signal processing to obtain the depth profile signal. The signal processing consisted of interpolation of the captured spectrum to a linear wave number space, fast Fourier transformation, and logarithm mapping and display. The imaging depth, measured from a depth position of 10 dB amplitude drop, was approximately 5 mm. An axial resolution of the system was measured from the full width at FWHM of the depth profiles to be between 15-20 µm in air. The capturing of spectral interference signals was synchronized with the scanning of the light beam on the sample to produce 2D, 3D, or 4D OCT datasets. The Doppler phase shift detection scheme was implemented to monitor the propagation pattern of vibrational waves 4 . The data acquisition and OCT signal processing were implemented in the Labview program (Labview, National Instruments, Austin, TX, USA).