Longitudinal shear waves for elastic characterization of tissues in optical coherence elastography

: In dynamic optical coherence elastography (OCE), surface acoustic waves are the predominant perturbations. They constrain the quantiﬁcation of elastic modulus to the direction of wave propagation only along the surface of tissues, and disregard elasticity gradients along depth. Longitudinal shear waves (LSW), on the other hand, can be generated at the surface of the tissue and propagate through depth with desirable properties for OCE: (1) LSW travel at the shear wave speed and can discriminate elasticity gradients along depth, and (2) the displacement of LSW is longitudinally polarized along the direction of propagation; therefore, it can be measured by a phase-sensitive optical coherence tomography system. In this study, we explore the capabilities of LSW generated by a circular glass plate in contact with a sample using numerical simulations and tissue-mimicking phantom experiments. Results demonstrate the potential of LSW in detecting an elasticity gradient along axial and lateral directions simultaneously. Finally, LSW are used for the elastography of ex vivo mouse brain and demonstrate important implications in in vivo and in situ measurements of local elasticity changes in brain and how they might correlate with the onset and progression of degenerative brain diseases.

approaches [10,11], given the trio of depth penetration of OCT, boundary conditions of tissues, and mechanical excitation wavelengths. SAW propagation has been demonstrated to discriminate elasticity gradients along the propagation direction (usually guided by the tissue surface) with important applications in the study of tissue anisotropy [12], the evaluation of elasticity changes before and after diseases or treatments [13], and the investigation of the relationship between SAW speed and tissue biomarkers such as intraocular pressure in ocular tissue [14]. Nevertheless, the discrimination of elasticity gradient along depth (perpendicular to the tissue surface) is quite limited due to the dispersive nature of SAW [11,15].
Depth-dependent elasticity discrimination of tissues is fundamental for composite and heterogeneous media. For example, human cornea is a highly organized tissue composed of at least five thin layers with differentiated structure, mechanical properties, and physiological functions [16]. Therefore, the elastic characterization of individual corneal layers using non-destructive techniques is of great importance in ophthalmology for better understanding of ocular physiology, diagnosis of ocular pathologies, and monitoring of relevant diseases and treatments [6]. Another application relates how local and global changes in the elastic properties of the brain may correlate with the onset and progression of degenerative brain diseases, such as Alzheimer's disease [17] and Parkinson's disease [18]. Therefore, measuring the heterogeneous distribution of elasticity in the brain will be valuable from a clinical and basic science perspective.
Longitudinal shear waves (LSW) are a special case of shear waves, characterized as polarized slow waves originating in the near-field regime of a perturbation in elastic solids [19]. It has been theoretically studied in [19,20] and utilized in ultrasound elastography [21]. When LSW are generated on the surface of tissues using a force/displacement source vibrating normal to the tissue surface, they can travel along the depth direction (and, importantly, are not guided by the tissue surface as with SAW), which allows them to discriminate depth elasticity gradients via changes in the wave speed. Displacement caused by LSW is longitudinally polarized along the direction of propagation; therefore, it can be measured and detected by Doppler ultrasound or phase-sensitive optical coherence tomography (PhS-OCT). In OCE, preliminary studies have been conducted in phantoms, validating SAW for detecting 1D elasticity gradients along lateral [22] and axial [23] directions. In [23], coaxial excitation of LSW is demonstrated for depth-dependent 1D elasticity characterization in phantoms.
In this study, we explore the capabilities of LSW in detecting 2D elasticity gradients along axial and lateral directions simultaneously, using numerical simulations and tissue-mimicking phantom experiments. LSW are generated by a circular glass plate in contact with the sample. Finally, we demonstrate the use of LSW for the elastography of ex vivo mouse brain tissue with important implications in in vivo and in situ measurements of local elasticity changes in brain for neuromedicine and related clinical applications.

Sample preparation and characterization
Tissue-mimicking phantoms were made of gelatin powder of different concentrations for simulating softer (3% gelatin) and stiffer (5% gelatin) tissues. In all cases, the same concentration of 0.5% intralipid powder was used to provide optical scattering to the phantoms (measured refractive index of phantoms ∼1. 35). Five different elasticity configurations were tested (where L = left, R = right, T = top, and B = bottom): uniform 3%, uniform 5%, [L3% − R5%] horizontally distributed halves, and [T3% − B5%] or [T5% − B3%] vertically distributed layers. The average thickness of the top layer was measured to be ∼0.8 mm and the bottom layer is considered as a semi-infinite media. The approximate size of each phantom was 80 × 80 mm in the lateral extent, and 30 mm along depth.
The Young's modulus of each phantom concentration was measured by taking three cylindrical samples (n = 3) of ∼40 mm in diameter, and conducting a stress-relaxation compression test in a MTS Q-test/5 universal testing machine (MTS, Eden Prairie, Minnesota, USA) at a strain of 5% for a total time of 600 s. The stress-time plots obtained were fitted to a Kelvin-Voigt fractional derivative rheological model for the frequency-dependent Young' s modulus (E) estimation [24]. For the excitation frequencies used in the experiments (1 kHz and 2 kHz), Young's moduli were found to be similar (assuming negligible non-linear elastic behavior) for each concentration: E 3% = 3.95 ± 0.23 kPa and E 5% = 10.57 ± 0.68 kPa, for 3% and 5% concentrations, respectively. Assuming that the phantom is an isotropic and homogeneous material of density of ∼1 g/cm3, and a Poisson's ratio of ∼0.5, the shear wave speed can be calculated [1] and reported as: c 3% s = 1.15 ± 0.03 m/s, and c 5% s = 1.88 ± 0.06 m/s, for 3% and 5% concentrations, respectively. Brains were procured from euthanized (no more than one hour) 8-week old male C57B1/6 mice (Charles River Laboratories, Wilmington, MA, USA) and placed in a 2% concentration agarose gel solution of artificial cerebrospinal fluid, which slows the desiccation of brain tissue during the OCT imaging process. The temperature of the agarose solution was controlled not to exceed 37 • C before solidification to ensure all sides of the brain were covered with a separation between brain and gel surfaces not exceeding 0.1 mm for optimal OCT imaging. No optical scatterers were added in the coupling gel to ensure deeper OCT penetration in brain tissue. During OCE experiments, a harmonic loading glass coverslip is placed on the top of the agarose gel volume containing the mouse brain. The gel in between the coverslip and brain tissue (thickness ∼0.1 mm) serves as an elastic coupling medium. The shear wave speed in the 2% concentration agarose gel was measured to be 2.67 ± 0.05 m/s, which is in the range of previously reported shear wave speed in mouse brain [25] and avoids the creation of strong wave reflections.

Experimental setup
The schematic of the experimental setup is shown in Fig. 1. Here, a PhS-OCT system is used for the motion detection of LSW generated in the sample via a synchronized mechanical excitation using a piezoelectric actuator attached to a circular glass coverslip. The PhS-OCT system is implemented with a swept-source laser (HSL-2100-WR, Santec, Aichi, Japan) with a center wavelength of 1318 nm and a full-width half-maximum (FWHM) bandwidth of 125 nm. The imaging depth was measured to be 5 mm (-10 dB sensitive fall-off). The optical lateral resolution is approximately 20 µm, and the FWHM of the axial point spread function after dispersion compensation is 10 µm. The A-line rate is 20 kHz, which provides a sampling time T s = 50µs. A galvo-scanner was used to change the lateral position of the beam in the sample. The synchronized control of the galvo-scanner, mechanical excitation, and the OCT data acquisition were conducted using a LabVIEW (Version 14.0f1, National Instruments Corporation, Austin, TX, USA) platform connected to a work station. The mechanical excitation system begins with a function generator (AFG320, Tektronix, Beaverton, OR, USA) output signal connected to an ultra-low noise power amplifier (PDu150, PiezoDrive, Callaghan, NSW, Australia) feeding a small piezoelectric actuator (SA030305, PiezoDrive, Callaghan, NSW, Australia) of 3 × 3 × 5 mm. The actuator is fixed to the top surface of a circular glass coverslip (Thermo Fisher Scientific, Waltham, MA, USA). The bottom surface of the glass coverslip is touching the sample. The OCT laser beam passes through the glass while vibrating along the axial direction for the generation and measurement of LSW as described in Fig. 2. Distortion in the phase information of OCT A-lines generated by the movement of the glass plate is modeled and compensated for (see Section 2.3). The signal set in the function generator for experiments in phantoms and ex vivo brain tissue is a 1 kHz (5 cycles) or 2 kHz (10 cycles) harmonic signal. The region of interest (ROI) was defined to be 10 mm × 2.5 mm along the lateral and depth axis of the cross-sectional B-mode plane, respectively. Fig. 1. Experimental setup using a PhS-OCT imaging system synchronized with mechanical excitation using a piezoelectric actuator attached to a circular glass coverslip. Two types of samples are shown: gelatin tissue-mimicking phantoms, and embedded mouse brain tissue. Fig. 2. Configuration of the sample, glass coverslip, and the OCT scanning probe for analyzing waves generated in the sample and particle velocities measured by OCT.

Acquisition scheme
The MB-mode acquisition approach, as described in Zvietcovich et al., was used for data acquisition [26]. It consists of triggering the excitation (via a burst of the 1 kHz or 2 kHz continuous sinusoidal waveform) and acquiring 100 A-lines during a 5 ms period at a single position within the lateral axis of the ROI. Subsequently, the galvanometer controlling the x-axis changes position to scan an adjacent lateral location. The excitation and acquisition process are then repeated for different positions. The process is stopped when the x-axis has covered the 10 mm lateral length of the ROI. The number of samples taken in the x-axis is 200. The total acquisition time is 1 second per experiment. The 3D data (2D spatial, 1D temporal) is then reorganized into 100 2D spatial frames of 10 × 2.5 mm, where each frame corresponds to a virtual instantaneous snapshot of sample motion, and the apparent frame rate is 20 kHz.
Methodologically, in a medium with refractive index n, the depth-dependent phase difference of a given A-line at two consecutive instants t 0 and t 1 (t 0 < t 1 ) for a given lateral x 0 position, is related to particle velocity v z (z) in the axial direction by where φ is the phase of the A-line signal acquired with the PhS-OCT system, n is the refractive index of the medium, λ 0 is the center wavelength of the laser, T s and is the time sampling resolution. In practice, given the experimental setup, three media with different refractive indexes are concatenated coaxially as shown in Fig. 2: coupling air, glass, and the sample. Then, ∆φ (z) measurements (and, therefore, v z (z)) can be affected by the motion of surfaces within the transition of materials. An appropriate model to compensate for those effects and estimate the true phase difference ∆φ T (z) produced by any scatterer within the sample is presented as where ∆φ (z), ∆φ (s 1 ), and ∆φ (s 2 ) are phase differences measured at the sample, the boundary between the sample and the bottom surface of the glass coverslip, and the boundary between the top surface of glass coverslip and air, respectively. Also, n 1 , n 2 , and n 3 are the refractive indexes of air, the glass coverslip, and the sample, respectively. The derivation of Eq. (2) is presented in the Appendix Section. Assuming the glass coverslip is rigid, then ∆φ (s 2 ) = ∆φ (s 1 ), and Equation 2 reduces to which is the same correction factor found in [27] for a two-medium configuration (air and sample). We select Loupas' algorithm [28] for the accurate estimation of ∆φ (z), which increased the signal-to-noise ratio compared to other approaches. For the corrected phase, we use n 1 ≈ 1 (air), and n 3 ≈ 1.35 (gelatin phantom), and the ∆φ (s 1 ) measurement at the top surface of the glass coverslip. Finally, spatio-temporal particle velocity fields v z (x, z; t) are obtained within the ROI for each experiment.

Numerical simulation
Numerical simulations of LSW propagation produced by an axisymmetric, uniformly distributed transient displacement excitation applied on the top of a 3D solid volume were conducted using finite elements in Abaqus/CAE version 6.14-1 (Dassault Systems, Velizy-Villacoublay, France). The 3D solid elastic media of 16 × 16 × 10 mm size is subjected to a spatially-uniform (disk shape) and temporally-transient (one cycle of a 1 kHz signal) displacement field on the top surface, and zero displacement and rotation on the bottom and lateral surfaces of the solid (see Fig. 3(a)). The spatial distribution of the displacement field in the disk was softened at the borders in order to avoid drastic changes that may produce numerical errors to propagate (ringing effect). Similarly, the transient 1 kHz signal was softened at the initial transition from zero displacement to the first cycle. The solid was meshed with an approximate grid size of 33 µm and using linear and hexahedral dominant elements (C3D8R). The type of simulation was selected to be dynamic-implicit with a minimum temporal discretization of 1.65 µs and analyzed during a time range of 3 ms.  Table 1.
Two linear elastic materials whose properties mimicked the mechanical measurements of the 3% and 5% gelatin concentration phantoms described in Section 2.1 were chosen to characterize four different configurations: uniform 3%, uniform 5%, [L3% − R5%] horizontally-distributed halves, and [T3% − B5%] vertically-distributed layers as shown in Fig. 3(b). In all cases, the diameter of the displacement disk in the top surface of the solid is 9 mm. Table 1 shows the selected density, Poisson's ratio, Young's modulus, and the equivalent shear wave speed for each designated material. Finally, for the uniform 5% material case, LSW is simulated for different disk diameters of displacement distribution: {4, 5, 6, 7, 8, 9} mm in order to evaluate its effect on the generation of LSW wavefronts. Simulated particle velocity v F E M z (x, y, z; t) estimated along the z-axis is evaluated at the xz-plane crossing through the center of the volume (justified by the symmetric assumptions of the model) and organized in 2D spatial frames along time, generating v F E M z (x, z; t) for the calculation of LSW speed. Table 1. Elastic material parameters of 3% and 5% materials defined in Abaqus/CAE version 6.14-1.

LSW speed estimation
Local speeds of LSW are calculated using two methods: (1) time-of-flight approach for transient excitations and (2) correlation approach for harmonic excitations. The time-of-flight approach is based on space-time waveform tracking along the propagation direction. For a given x-axis position, the local group velocity of the wavefront c LSW is approximated by the slope of space-time curves within a window (W) along depth by where z w is the location of the window in the z-axis, and ∆z is the distance traveled by the wavefront during time ∆t within the window along depth. Equation (4) can be adapted to a least square regression approach as described in [26] for reducing noise and increasing accuracy in the estimation of wave speed. The sequential estimation c LSW (z w ) of along depth for each lateral position until the entire ROI is covered generates a 2D shear wave speed map (SWSM), also called an elastogram.
The correlation approach (for harmonic excitation), also referred to as Hoyt's method in [26,29], leverages the complex-valued spatial particle velocity field v C plx z (x, z) obtained by: (1) Calculating the Fourier transform of the time signal at a given (x 0 , z 0 ) position: (2) evaluating the transformed signal at the angular frequency Defining the z-axis as the LSW propagation direction, for a given lateral position x 0 , the first lag of the auto-correlation function R of v C plx z using a kernel window of size L along depth is: where z n = n∆ z , for n = 0, 1, 2, ..., N − L − 1; where N is the number of samples along the z-axis, and ∆ z is the depth sampling resolution. Then, the local spatial wavenumber k can be estimated as: where Re {·} and Im {·} refer to the real and imaginary parts of a complex number. Subsequently, the LSW speed can be calculated as c LSW (x 0 , z n ) = ω 0 /k (x 0 , z n ). The sequential estimation of c LSW (x 0 , z n ) along depth for each lateral position until the entire ROI is covered generates a 2D SWSM. Motion information (particle velocity) coming from pixels with a low structural OCT signal were masked out and, therefore, not taken into account for the LSW speed calculation.

LSW propagation in numerical simulations
The effect of plate diameter on LSW propagation was evaluated using a uniform 5% material property in simulations. Figure 4 shows a sequence of motion frames at instants {1.2, 1.8, 2.4} ms when plate diameters of 9 mm, 6 mm, and 4 mm are excited with one cycle of a 1 kHz signal. The LSW wavefronts are propagating along the depth direction from the surface of the medium. In all cases, the edge of the plate in the surface of the medium produces shear waves propagating outwards and inwards, and towards depth as expected in the theoretical description of Figure 2. While the outward propagating Rayleigh and shear waves fade with distance from the source, the inward propagating shear wave can interact with the LSW and perturb the estimation of speed and amplitude. In the propagation sequence of Figure 2, a dashed border box shows that the inward propagating shear wave is approaching the center of the LSW faster when the diameter of the disk is smaller. This will have important implications for selecting the appropriate plate diameter for excitation of LSW.
The displacement amplitude and speed of the LSW are analyzed from space-time representations extracted at the center of the plate along depth (position x 1 = 8 mm in Fig. 4)   the disk diameter influences the spatio-temporal effect of the inward propagating shear wave (dashed border box) in distorting the amplitude and speed of the LSW wavefront. Displacement of the LSW main peak before the interference of the inward propagating shear wave (∼1.2 ms in Fig. 5(d)) tends to decrease when the plate diameter is increased, accounting for the geometrical spreading of the initial LSW wavefront.
For the simulation, we ensured that all disk diameters were excited with the same initial displacement in the sample. Later in time, the inward propagating shear wave interferes with the LSW wavefront producing an initial displacement increment as shown in Fig. 5(d). This effect can be more dramatically observed when the speed of the LSW is analyzed in Fig. 5(e). The perturbation caused by the shear wave deviates the LSW speed from the nominal value of 1.943 m/s. This near field effect occurs closer to the surface as the plate diameter is decreased and it fades during propagation (as in the far field regime). When the plate diameter is increased, the perturbation occurs farther in depth and then disappears as the LSW propagates deeper in the material. Since the measured OCT depth of field typically ranges from 1 mm to 3mm, a larger diameter plate is desired to avoid the near field effect.
In Fig. 6, the potential of utilizing LSW for the [L3% − R5%] characterization of the horizontally-distributed materials is explored. Figure 6(a) shows a sequence of motion frames at instants {1.2, 1.8, 2.4} ms when a plate of 9 mm diameter is used. The main LSW wavefront is propagating from the surface of the media toward depth at a different speed in each half. In the boundary of both halves, the refraction of waves is stronger as it propagates, limiting the capabilities of detecting sharp changes in LSW speed between both halves. Figures 6(b,c) show space-time representations of wave propagation at lateral positions x 1 = 6 mm (3% material) and x 2 = 10 mm (5% material), respectively. The first perturbation detected corresponds to the pressure waves, which are also called P-waves or longitudinal fast waves because their speed is at least two orders of magnitude faster in nearly incompressible solids such as biological tissue. The second perturbation corresponds to the LSW propagation and its slope corresponds to the wave speed. The time-of-flight approach for speed estimation in both (3% and 5%) material sides is reported in Fig. 6(d) when the entire propagation extent is considered. Speeds of 1.18 m/s (3% material) and 1.94 m/s (5% material) are consistent with the material's ground truth definitions reported in Table 1. Small discrepancies between estimations and ground truth speed values are attributed to the refraction effects in the boundary.  Fig. 7, the potential of utilizing LSW for the characterization of the layered [T5% − B3%] vertically distributed media is explored. Figure 7(a) shows a sequence of motion frames at instants {1.2, 1.8, 2.4} ms when a plate of 9 mm diameter is used. The main LSW wavefront is propagating along the depth direction from the surface of the media at different speeds in each layer. In the boundary of both layers, the refraction of waves occurs, limiting the capabilities of detecting sharp changes in LSW speed between both layers. Figure 7(b) shows the space-time representation of wave propagation at lateral position x 1 = 8 mm. The first perturbation detected corresponds to the pressure wave, as in previous cases. The second perturbation describes the LSW propagation and how its slope changes in each layer. The time-of-flight approach for speed estimation in both the 5% (top layer) and 3% (bottom layer) materials is reported in Fig.  7(c), when the entire propagation extent is considered. Speeds of 1.08 m/s (3% material) and 1.75 m/s (5% material) are consistent with material ground truth definitions reported in Table  1. Speed underestimation in both layers is attributed to the interference of waves coming from the plate edges that rapidly collapse the boundary between the 5% and 3% layers. In Fig. 7(d), depth-dependent speed estimations in both layers are reported when using the time-of-flight approach in a small moving window of 0.1 mm. LSW speed variability can be explained by the repetitive refraction of waves produced at both layer boundaries (or a classical resonance effect) as well as the propagation of numerical error during simulations. The transition from T5% to B3% (z-axis) was fitted to a sigmoid function for OCE resolution characterization described as where c 3% and c 5% are the average LSW speeds in the 3% and 5% media, respectively; z 0 and τ represent the location and width of the transition, respectively. Then, we define the OCE resolution along the z-axis as the distance from 20% to 80% of the transition, also called R 2080 [26]. We found R 2080 = 0.109 mm as reported in Fig. 7(d).

LSW propagation in tissue-mimicking phantoms
In Fig. 8, LSW propagation is explored in tissue-mimicking phantoms of uniform 3% and 5% gel concentrations as a 9 mm glass coverslip is continuously excited at 1 kHz. Motion frames taken at 1 ms instants show differing wavelengths between the 3% and 5% materials ( Fig. 8(a)). The correlation approach was used to generate the 2D LSW speed elastograms as shown in Fig. 8 Fig. 9 as a glass coverslip of 4 mm or 9 mm diameter is excited at 1 kHz and 2 kHz. In Fig. 9(a), motion frames captured at 1 ms instants demonstrate a Rayleigh wave propagating outwards from the glass coverslip borders along the phantom surface, and LSW propagating along depth from the phantom surface. Differentiated wavelength patterns in all cases between the 3% and 5% halves are evident.  Smaller wavelength patterns are found when the excitation frequency is increased, allowing for the elastic characterization of smaller components (which corresponds to improvement in OCE resolution). Refraction and reflection effects are found in the boundary between both halves as predicted in numerical simulations, creating strong artifacts in the LSW speed calculation shown in Fig. 9(b) using the correlation approach. In the 9 mm plate the size of artifacts is 1 mm approximately, corrupting the average estimation of speed and lateral OCE resolution in the transition of both material halves. Average LSW speeds of 0.76 ± 0.11 m/s and 1.68 ± 0.19 m/s are reported for the L3% and R5% material cases, respectively, considering all excitation frequencies and disk diameters. The average transition from L3% to R5% (x-axis) in the regions of strongest presence of artifacts was fitted to a sigmoid function (Eq. (7) along the x-axis) for OCE resolution characterization. We found R 2080 = 0.78 mm for the 9 mm plate diameter, and R 2080 = 0.46 mm for the 4 mm plate diameter, averaged over all frequencies. Discrepancies between R 2080 estimations for the 1 kHz and 2 kHz cases (for a given plate diameter) are weak since the presence of artifacts dominates the transition compared to the LSW wavelength.  Fig. 10 when a 4 mm and 9 mm diameter glass coverslip is exited at 2 kHz. In Fig. 10(a), motion snapshots obtained at 1 ms instants show the changes in LSW wavelength pattern in each layer when traveling along depth from the top surface. Propagation of Rayleigh waves in the surface and reflections of LSW produced in the boundary between layers are present as predicted in numerical simulations. Figure 10(b) shows 2D LSW elastograms calculated using the correlation approach from Figure 10a. Average LSW speeds for all cases including uniform, and horizontally-and vertically-distributed material layers are reported in Fig. 11. Finally, more defined LSW speed transitions between layers are found in vertically distributed media when compared to horizontally distributed materials due to the directionality of reflections produced in the boundary. The average transitions from T3% to B5% and T5% to B3%, choosing regions of strongest presence of artifacts, and averaging for both disk diameter cases, are reported to be R 2080 = 0.101 mm and R 2080 = 0.128 mm, respectively, by applying Eq. (7) along the z-axis. Discrepancies in LSW speed values in phantoms when compared to mechanical measurements are attributed to reflections produced in the boundaries and the interference of Rayleigh and shear waves generated by the harmonic excitation.

LSW propagation in ex vivo mouse brain
Results from the ex vivo mouse brain experiment are demonstrated in Fig. 12. The glass coverslip was placed at the top surface of the agar phantom where the brain was suspended, covering sections of cerebral cortex, midbrain (superior and inferior colliculus), and cerebellum (see Fig. 12(a)). The 9 mm diameter coverslip was harmonically excited at 2 kHz. An en face structural image (showing a transverse cut) taken at a depth of 1.2 mm is shown in Fig. 12(b), which demonstrates finer details within the cerebral cortex. In Fig. 12(c), a motion snapshot taken at the same depth at time 1.8 ms shows Rayleigh waves propagating outwards and parallel to the coverslip orientation, and a large LSW aberrated by the brain structure. In order to characterize LSW propagation along depth, a cross-section along the xz-plane (y 1 = 6.2 mm) is analyzed through time in Fig. 12(d). Here, the LSW wavefront can be tracked at time instants {1.2, 1.8, 2.4} ms.
The correlation approach is used for LSW speed estimation in the selected cross-section ( Figure 12(e)-top) and a 2D speed elastogram is reported in Fig. Overall, these results demonstrate differentiated elasticity values between the cerebral cortex (stiffer) and cerebellum/midbrain (softer) regions. Frontal cuts in Figs. 12(e) and 12(f) intersect with the sagittal cut in Fig. 12(g) in two different regions as shown in Fig. 12(a). LSW speed values are reported within ±0.1 mm from each intersection: 1.09 ± 0.29 m/s at x 1 in Fig.12(e)bottom, 2.99 ± 0.80 m/s at x 1 in Fig. 12(f)-bottom, 0.81 ± 0.45 m/s at y 1 Fig. 12(g)-bottom, and 2.91 ± 0.22 m/s at y 2 in Fig. 12(g)-bottom. We found that LSW speeds at the intersections of each frontal and sagittal cut agree reasonably well within the error bounds.

Discussion
The propagation of LSW for the 2D elastic characterization of heterogeneous media is explored in this paper as an alternative to improve the discrimination between different elastic structures and layers, beyond the limited results possible using surface acoustic waves. Numerical simulation of LSW propagation using a vibrating plate for characterizing horizontally-and vertically-gradient materials established a starting point for: (1) understanding the detectability of these waves using axial-displacement sensitive systems (such PhS-OCT), (2) learning the properties of these waves by analyzing speed and amplitude changes based on the medium's elasticity distribution, and (3) recognizing the inconveniences caused by the near field effects, the interference with other types of waves (i.e. shear, pressure, or Rayleigh waves), and undesired reflections produced from the tissue boundaries. The reason for choosing a plate instead of any other actuator shape for the generation of LSW stems from the desire for a quasi-planar extended wavefront to avoid interference with other wave types. Another group [23] investigated the formation of LSW using a ring actuator instead of a plate for 1D elasticity measurements.
The propagation of LSW in tissue-mimicking phantoms has been demonstrated for 2D elastic characterization of elastic horizontally-and vertically-gradient materials. Aberrations of the wavefronts caused by reflections in the material boundaries are more dramatic in the [L3%− R5%] compared to the [T3% − B5%] or [T5% − B3%] cases since reflections created in the former case tend to propagate perpendicular to the boundary plane. Such reflections in the [L3% − R5%] case produce strong artifacts in the 2D LSW speed maps and, therefore, discrepancies between The color bar represents LSW speed in m/s. Finally, (f) and (g) show B-mode (top) and 2D LSW elastogram (bottom) taken at the cross-sectional xz-plane (or frontal cut) (y 2 = 2.7 mm) and yz-plane (or sagittal cut) (x 1 = 3.1 mm), respectively. The color bar represents LSW speed in m/s. measured LSW speeds and mechanical measurement result. Similarly, the boundary between dissimilar regions is strongly dominated by the presence of artifacts. Reflections can be minimized by the use of directional filters at the cost of compromising resolution [30].
The resolution of LSW-based elastography along the axial direction is limited by the LSW wavelength, and any reflections produced in the boundaries. When the boundary between different layers in materials is oriented perpendicular to the axial axis, the direction of the reflected and refracted waves tends to propagate axially. This case is very similar to Rayleigh waves traveling laterally in heterogeneous media and it has been intensively studied in [26]. Then, stiffer samples (higher Young's modulus) may produce higher LSW speeds and, therefore, longer wavelengths which could diminish the elastography resolution. This effect can be compensated for by increasing the mechanical excitation frequency in order to reduce the LSW wavelength. This alternative is constrained by the attenuation of waves in the sample (which in tissues tends to increase monotonically with the excitation frequency) and the spatial and temporal resolution capabilities of the OCT imaging system.
We conclude that LSW-based OCE has an axial elastography resolution approximately 5× greater than the lateral resolution, with an R 2080 on the order of 0.1 mm at 2 kHz harmonic excitation. Our OCE resolution measurements using a sigmoid function to estimate transitions between two differentiated materials accounts for the joint effect of OCT imaging resolution, motion and wave speed estimation kernels, and excitation wavelength. A more formal resolution study [31], separating the individual effect of each component on the total resolution is left for future work. LSW-based OCE has a unique advantage over other dynamic OCE techniques such as (1) SAW elastography, and (2) reverberant wave fields. On the one hand, we demonstrated that LSW provides of a good discrimination of elasticity gradient along depth which very limited, if not impossible, using SAW due to its dispersive nature. On the other hand, while LSW accounts for a very specific type of wave whose speed is the same as shear wave speed (desired since such speed can be converted into shear modulus in a straight forward fashion [1]), reverberant fields leverage in the harmonic excitation of tissues that could include, besides shear waves, other SAW branches such as Rayleigh or Lamb waves (with speeds that are not easily converted into shear modulus). When dealing with tissues with highly irregular boundary conditions and heterogeneous elastic properties, LSW-based OCE is limited by the wave mode conversion and strong reflections propagating in different directions. Then, reverberant wave fields may be established and, combined with LSW excitation methods, localized shear wave speed can still be recovered by using the appropriate estimator [32]. Therefore, the implementation of reverberant-LSW OCE could potentially improve the lateral resolution of LSW-based OCE. This will be material for a future work.
An experiment involving mouse brain elastography using LSW was performed and demonstrated that the cerebellum and midbrain were found to be softer than the cerebral cortex. This is consistent with various magnetic resonance elastography studies which have shown similar findings in both mice and humans [25,[33][34][35]. In these studies, a general trend for shear storage modulus is from 4 kPa to 16 kPa for the cerebral cortex, and from 1 kPa to 4 kPa for the cerebellum/midbrain. Our results show an average shear modulus of 9.36 kPa for the cerebral cortex, and 1.10 kPa for the cerebellum/midbrain, which are consistent with the aforementioned trend. In addition, the 2D LSW speed map of the cerebral cortex in Fig. 12(f)-bottom shows variations that may be attributed to the heterogeneity of the tissue type, which has also been previously reported using magnetic resonance elastography [35]. Artifacts created by strong wave reflections produced by the brain's irregular boundary conditions may affect the quality of speed estimation in Fig. 12(f)-bottom.
Finally, the results of LSW for the elastography of ex vivo mouse brain tissue indicates a promising method that can be implemented for in vivo studies. Pathological changes in neurodegenerative diseases are hypothesized to be correlated with local and global elastic changes in the brain. Preliminary in situ OCE studies have been performed to characterize elasticity changes in mouse brain [36]. The use of a transparent glass coverslip as an observation and protection window after craniotomy is a common practice in optical brain studies. It is ideal for studying mouse models since this procedure keeps the brain isolated and at normal intracranial pressure [37]. Even though OCT signals from brain scans can be obtained at an appropriate signal-to-noise ratio when imaging through a cranial window and glass coverslip, the skull/brain architecture prevents the brain from receiving localized mechanical excitation, which is fundamental for OCE applications. In this paper, we demonstrated the potential of producing LSW using the same glass coverslip for mouse brain studies, enabling the possibility of in situ and in vivo studies using OCE. Future research will focus on in vivo OCE of mouse brain using LSW, taking into account the boundary conditions involving the skull, cerebrospinal fluid, and the glass coverslip.

Conclusions
In this paper, we investigated the capabilities of LSW in characterizing 2D elastic material gradients simultaneously along axial (vertically-distributed layers) and lateral directions (horizontallydistributed halves) using a PhS-OCT imaging system for OCE applications. Numerical simulations were used to explore the detectability, properties, and inconveniences of LSW. Subsequently, experimental results in tissue-mimicking phantoms were corroborated with numerical simulations and confirmed the capabilities of LSW to characterize heterogeneous materials in OCE. Waves were generated with a circular glass plate in contact with the sample and vibrated at a continuous sinusoidal excitation while OCT signals were acquired in the sample passing through the glass. A motion compensation model is proposed to correct for the phase errors/artifacts due to glass motion. LSW speeds estimated in phantoms agreed with mechanical measurements. Finally, we demonstrated the use of LSW for the elastography of ex vivo mouse brain tissue, in particular observing speeds in the superficial cerebral cortex, midbrain, and cerebellum. Results enable the implementation of LSW-OCE for in vivo mouse brain studies that could be important for neuromedicine studies and clinical applications.

Appendix
Given the experimental setup shown in Fig. 2, and extending the derivation from [27], the total equivalent optical path length (OPL) from the point scatterer of interest (located at depth z in Fig. 2) and the fiber coupler, where the reference beam and the test beam interfere in Fig. 1, is OPL 0 (z) = L 1 n 1 + L 2 n 2 + L 3 n 3 , for time t = 0 when there is no motion. The roundtrip path of the light was already considered in Equation 1. After a time t = T s , motion at s 1 , s 2 , and z ( Figure 2) produces local displacements T s v s 1 , T s v s 2 , and T s v z , respectively (motion towards the z-axis is defined to be positive). Then, the equivalent OPL becomes: Then, the OPL difference, OPL 1 (z) − OPL 0 (z), between t = 0 and t = T s is given by: In the limiting cases, when the scatterer in z is located at any of the surfaces, then lim ∆OPL (z) = (n 1 − n 2 ) T s v s 1 + (n 2 ) T s v s 2 , ∆OPL (s 1 ) = lim Then, using Eq. (1) in Eqs. (11) and (12) to relate particle velocities v s 1 and v s 2 with phase differences ∆φ (s 1 ) and ∆φ (s 2 ), respectively, and using ∆OPL (z) = ∆φ (z) λ 0 /(4π) in Equation 10, we found that ∆φ (z) = 4πn 3 λ 0 T s v z + 1 − n 3 n 2 ∆φ (s 2 ) + n 3 n 2 − n 3 n 1 ∆φ (s 1 ) .
∆φ (z) is the apparent phase difference produced by the scatterer in z and measured by the OCT system; however, the true displacement T s v z generated by the scatterer produces a phase difference ∆φ T (z) = (4πn 3 /λ 0 ) T s v z (using Equation 1), which we call true phase difference. Finally, rearranging Equation 13, and isolating ∆φ T (z) we found Equation 2.