1. Introduction

Mechanical properties of biological tissues are unavoidably influenced by the progression of pathophysiological processes. Therefore, non-invasive measurements of these properties (i.e. elasticity and viscosity), also known as elastography, are useful for disease diagnostics, surveillance, and treatment [1]. Examples of applications in different tissue types include cornea, skin, liver, breast, and brain [1, 2]. Since the 1980’s, elastography has been implemented in various imaging modalities such as ultrasound, magnetic resonance, and optical coherence tomography (OCT) [1, 3]. In particular, OCT-based elastography, also called optical coherence elastography (OCE), has been intensely studied in the last two decades due to its favorable axial and lateral resolution (∼15 μm) in a variety of different biological tissues [4]. Much of the OCE work has been studied and categorized in review papers from years 2015 [3], 2016 [2], and 2017 [4–6]. Here, a taxonomy of OCE methods is conducted according to the following categories: spatial and temporal extent of loading, excitation approaches, and tissue-specific applications.

Dynamic OCE methods, an important subset of OCE techniques, leverage the propagation and measurement of mechanical waves for tissue viscoelastic characterization, which is usually estimated from the displacement and speed of such waves [7–9]. Surface acoustic waves (SAW), also referred to as surface guided waves, are the predominant perturbation in most dynamic OCE 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.

2. Materials and methods

2.1. 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%, $\left[L3\%-R5\%\right]$ horizontally distributed halves, and $\left[T3\%-B5\%\right]$ or $\left[T5\%-B3\%\right]$ 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\pm 0.23$ kPa and $E_{5\%}=10.57\pm 0.68$ kPa, for 33% = 1.15 $\pm 0.03$ m/s, and $c_s^{5\%}=1.88\pm 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\pm 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.

2.2. 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\mu$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.

 

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 mousebrain tissue.

Download Full Size | PDF

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\times 3\times 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. 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.

Download Full Size | PDF

2.3. 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₀ and t₁ ($t_0<t_1$) for a given lateral x₀ position, is related to particle velocity $v_z\left(z\right)$ in the axial direction by

is the refractive index of the medium, λ₀ 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, $\mathrm{\Delta }\varphi \left(z\right)$ measurements (and, therefore, $v_z\left(z\right)$) 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 $\mathrm{\Delta }{\varphi }_T\left(z\right)$ produced by any scatterer within the sample is presented as

2.4. 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.

 

Fig. 3 Numerical 3D solid elastic media of 16 × 16 × 10 mm, showing boundary conditions (BC) of displacement and the ROI for measurements. (b) Four cases of material distribution defined in the simulation. Details on Young’s modulus for 3% and 5% material are specified in Table 1.

Download Full Size | PDF

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%, $\left[L3\%-R5\%\right]$ horizontally-distributed halves, and $\left[T3\%-B5\%\right]$ 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: $\left\{4,5,6,7,8,9\right\}$ mm in order to evaluate its effect on the generation of LSW wavefronts. Simulated particle velocity $v_z^{FEM}\left(x,y,z;t\right)$ 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_z^{FEM}\left(x,z;t\right)$ for the calculation of LSW speed.

Table 1. Elastic material parameters of 3% and 5% materials defined in Abaqus/CAE version 6.14-1.

View Table

2.5. 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

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_z^{Cplx}\left(x,z\right)$ obtained by: (1) Calculating the Fourier transform of the time signal at a given $\left(x_0,z_0\right)$ position: ${\Im }_{t_z}\left\{v_z\left(x_0,z_0;t\right)\right\}\left(x_0,z_0;\omega \right)$, and (2) evaluating the transformed signal at the angular frequency ${\omega }_0=2\pi f_0$, where f₀ is the excitation frequency: $v_z^{Cplx}\left(x_0,z_0\right)={\Im }_{t_z}\left\{v_z\left(x_0,z_0;t\right)\right\}\left(x_0,z_0;{\omega }_0\right)$. Defining the z-axis as the LSW propagation direction, for a given lateral position x₀, the first lag of the auto-correlation function R of $v_z^{Cplx}$ using a kernel window of size L along depth is:

3. Results

3.1. 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 $\left\{1.2,1.8,2.4\right\}$ 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 Fig. 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 Fig. 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.

 

Fig. 4 Sequence of motion frames showing LSW propagation when excited by a disk of 9 mm (top row), 6 mm (middle row), and 4 mm (bottom row) diameter, and one cycle of a 1 kHz signal in simulations for a uniform 5% material. Color bar represents displacement in μm. The dashed border box shows how the inward propagating shear wave reaches the center of the LSW faster when a smaller diameter disk is used.

Download Full Size | PDF

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) when excited by plate diameters of 4, 5, 6, 7, 8, and 9 mm, and shown in Fig. 5. In Figs. 5(a-c), we show how 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.

 

Fig. 5 Space-time map representation of LSW propagation taken along depth and time at position x₁ = 8 mm in Fig. 4, when a disk diameter of 9 mm (a), 6 mm (b), and 4 mm (c) is used as displacement loading, and excited using one cycle of a 1 kHz signal. Solid line indicates LSW propagation and the dashed border box shows the inward propagating shear wave distorting the LSW speed and amplitude. Color bar represents displacement in μm. (d) Displacement produced by the main LSW wavefront peak, while propagating in space and time, produced by disk-distribution displacement of different diameters. (e) Wave speed measured by tracking the main peak of the LSW wavefront produced by disk-distribution displacement of different diameters.Ground truth (GT) shear wave speed in the simulation is 1.943 m/s.

Download Full Size | PDF

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 $\left[L3\%-R5\%\right]$ characterization of the horizontally-distributed materials is explored. Figure 6(a) shows a sequence of motion frames at instants $\left\{1.2,1.8,2.4\right\}$ 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. Figs. 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. 6 Sequence of motion frames showing LSW propagation when excited by a 9 mm diameter disk at one cycle of a 1 kHz signal in simulations for the [L3% − R5%] material. The color bar represents displacement in μm. Space-time representations of LSW propagation are shown in (b) and (c) when measured along depth and time at lateral positions $x_1=6$ mm and $x_2=10$ mm, respectively, in (a). Depth-time tracking of LSW wavefront peak and speed calculation measured from (b) and (c).

Download Full Size | PDF

In Fig. 7, the potential of utilizing LSW for the characterization of the layered $\left[T5\%-B3\%\right]$ vertically distributed media is explored. Figure 7(a) shows a sequence of motion frames at instants $\left\{1.2,1.8,2.4\right\}$ 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

 

Fig. 7 (a) Sequence of motion frames showing LSW propagation when excited by a 9 mm diameter disk at one cycle of a 1 kHz signal in simulations for the [T5% − B3%] material. The color bar represents displacement in μm. (b) Space-time representation of LSW propagation measured along depth and time at lateral position x₁ = 8 mm in (a). (c) Depth-time tracking of LSW wavefront peak and speed calculations measured from (b). (d) Depth-dependent LSW speeds measured from (c) and compared to ground truth (GT) shear wave speed values of each layer material obtained in simulations. A sigmoid fitting is shown using a dashed line for resolution characterization.

Download Full Size | PDF

3.2. LSW propagation in tissue-mimicking phantoms

In Figure 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(b). Average LSW speeds of $0.80\pm 0.04$ m/s and $1.67\pm 0.09$ m/s are reported for the uniform 3% and 5% materials, respectively, considering all excitation frequencies and disk diameters. Experimental results have the same tendency as mechanical measurement results: $c_s^{3\%}=1.15\pm 0.03$ m/s, and $c_s^{5\%}=1.88\pm 0.06$m/s, for 3% and 5% concentrations, respectively. A bias error between experimental and mechanical measurement results is reported to be approximately 0.33 m/s.

 

Fig. 8 Experimental results of LSW propagation in uniform phantom materials 3% (top row), and 5% (bottom row) when excited by a 9 mm diameter disk with a continuous 1 kHz harmonic signal. (a) Motion frames recorded at time intervals of 1 ms. The colorbar is in arbitrary units representing particle velocity. (b) 2D SWSM representing local shear wave speeds using Eq. (6). The color bar is in m/s.

Download Full Size | PDF

Results for the $\left[L3\%-R5\%\right]$ horizontally-distributed phantoms are reported in Fig. 9 as a glass coverslip of 4 mm or 9 mm diameter is excited at 1 kHz and 2 kHz. In 9a, 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. 9b 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\pm 0.11$ m/s and $1.68\pm 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₂₀₈₀ 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. 9 Experimental results of LSW propagation in the [L3% − R5%] phantom material when excited by 9 mm (top two rows) and 4 mm (bottom two rows) diameter disks with a continuous harmonic signal of 1 kHz and 2 kHz. (a) Motion frames recorded at time intervals of 1 ms. The color bar is in arbitrary units representing particle velocity. (b) 2D SWSM representing local shear wave speed using Eq. 6. The color bar is in m/s.

Download Full Size | PDF

Results in $\left[T3\%-B5\%\right]$ and $\left[T5\%-B3\%\right]$ layer media are reported in 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 Fig. 10(a). 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.

 

Fig. 10 Experimental results of 2 kHz harmonic LSW propagation in [T3% − B5%] (top two rows) and [T5% − B3%] (bottom two rows) phantom materials when excited by 9 mm and 4 mm diameter disks. (a) Motion frames recorded at time intervals of 1 ms. The color bar is in arbitrary units representing particle velocity. (b) 2D SWSM representing local shear wave speed using Eq. (6). The color bar is in m/s.

Download Full Size | PDF

 

Fig. 11 Experimental LSW speed report in all phantom configurations: uniform, [L3% − R5%], [T3% − B5%], and [T5% − B3%] for different glass coverslip diameters and excitation frequencies. Bar heights and whisker size represent average and standard deviation of LSW speed in each case, respectively; and they are compared with ground truth (GT) shear wave speed values of 1.184 m/s for the 3%, and 1.943 m/s for the 5% layer material obtained in simulations.

Download Full Size | PDF

3.3. 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 $\left\{1.2,1.8,2.4\right\}$ ms.

 

Fig. 12 Experimental results in ex vivo mouse brain tissue from a 9 mm diameter glass coverslip excited continuously at 2 kHz. (a) 3D surface representation of the imaged brain area including brain components (cerebral cortex, midbrain, and cerebellum) and locations of cross-section analysis. (b) En face structural image (or transverse cut) of brain taken at a depth of 1.2 mm. (c) Motion snapshot of (b) at 1.8 ms instant showing LSW and Rayleigh wave propagation. The color bar represents particle velocity in arbitrary units. (d) LSW propagation analyzed at xz-plane cross-section (y₁ = 6.2 mm). The color code is the same as in (c). (e) B-mode (top) and 2D LSW elastogram (bottom) obtained after processing (d). 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.7 mm) and yz-plane (or sagittal cut) (x₁ = 3.1 mm), respectively. The color bar represents LSW speed in m/s.

Download Full Size | PDF

The correlation approach is used for LSW speed estimation in the selected cross-section (Fig. 12(e)-top) and a 2D speed elastogram is reported in Fig. 12(e)-bottom. Figure 12(e)-bottom shows that in a frontal slice, LSW in the midbrain region travels at LSW speed of approximately $1.18\pm 0.31$ m/s (averaged within a box shown in the same figure). Figures 12(f)-(g) show cross-sectional (frontal and sagittal cuts) B-mode images (top) and 2D elastograms (bottom) taken at the xz-plane ($y_2=2.7$ mm) and the yz-plane ($x_1=3.1$ mm). Figure 12(f)-bottom shows LSW speeds in the cerebral cortex of $2.91\pm 0.86$ m/s (averaged within a box shown in the same figure), demonstrating elastic heterogeneity. Figure 12(g)-bottom indicates that LSW speeds in the cerebral cortex are approximately $3.21\pm 0.51$ m/s and gradually decrease to approximately $0.94\pm 0.36$ m/s in the midbrain and cerebellum. Reported speed values were averaged within the box regions shown in Fig. 12(g)-bottom.

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 $\pm 0.1$ mm from each intersection: $1.09\pm 0.29$ m/s at x₁ in Fig. 12(e)-bottom, $2.99\pm 0.80$ m/s at x₁ in Fig. 12(f)-bottom, $0.81\pm 0.45$ m/s at y₁ Fig. 12(g)-bottom, and $2.91\pm 0.22$ m/s at y₂ 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.

4. 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 $\left[L3\%-R5\%\right]$ compared to the $\left[T3\%-B5\%\right]$ or $\left[T5\%-B3\%\right]$ cases since reflections created in the former case tend to propagate perpendicular to the boundary plane. Such reflections in the $\left[L3\%-R5\%\right]$ case produce strong artifacts in the 2D LSW speed maps and, therefore, discrepancies between 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₂₀₈₀ 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–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.

5. 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 (horizontally-distributed 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 neuro medicine studies and clinical applications.