Two-dimensional (2D) dynamic vibration optical coherence elastography (DV-OCE) for evaluating mechanical properties: a potential application in tissue engineering

: Mechanical properties in tissues are an important indicator because they are associated with disease states. One of the well-known excitation sources in optical coherence elastography (OCE) to determine mechanical properties is acoustic radiation force (ARF); however, a complicated focusing alignment cannot be avoided. Another excitation source is a piezoelectric (PZT) stack to obtain strain images via compression, which can affect the intrinsic mechanical properties of tissues in tissue engineering. In this study, we report a new technique called two-dimensional (2D) dynamic vibration OCE (DV-OCE) to evaluate 2D wave velocities without tedious focusing alignment procedures and is a non-contact method with respect to the samples. The three-dimensional (3D) Fourier transform was utilized to transfer the traveling waves ( x, y, t ) into 3D k -space ( k x , k y , f ). A spatial 2D wavenumber filter and multi-angle directional filter were employed to decompose the waves with omni-directional components into four individual traveling directions. The 2D local wave velocity algorithm was used to calculate a 2D wave velocity map. Six materials, two homogeneous phantoms with 10 mm thickness, two homogeneous phantoms with 2 mm thickness, one heterogeneous phantom with 2 mm diameter inclusion and an ex vivo porcine kidney, were examined in this study. In addition, the ARF-OCE was used to evaluate wave velocities for comparison. Numerical simulations were performed to validate the proposed 2D dynamic vibration OCE technique. We demonstrate that the experimental results were in a good agreement with the results from ARF-OCE (transient OCE) and numerical simulations. Our proposed 2D dynamic vibration OCE could potentially pave the way for mechanical evaluation in tissue engineering and for laboratory translation with easy-to-setup and contactless advantages. interacted waves in 3D k -space. A spatial 2D wavenumber filter and multi-angle directional filter were employed to decompose the waves into four individual directions. The 2D local wave velocity algorithm was used to evaluate a 2D wave velocity map. The ARF-OCE was used to evaluate wave velocities for comparison. Numerical simulations were performed to validate the proposed method. The experimental results were in a good agreement with the results from ARF-OCE and numerical simulations.


Introduction
Tissue mechanics have been comprehensively investigated for the past two decades with elastography techniques because tissue mechanical properties are significantly associated with disease states. For example, fibrosis and cancer cause prominent changes of mechanical properties in tissues [1]. The changes of tissue stiffness correlate with extracellular matrix (ECM) fibers interacting with fibroblast proliferation and differentiation [2]. In addition, mechanical properties of tissues can be altered by cells themselves. Those cells dominate the machinery of tissues by their mechanical forces such as stiffness of cytoskeleton and their responses to surrounding microenvironment. In the case of cancer, solid stress, caused by growing tumors in a constrained physical volume, will increase abnormal cell proliferation. The abnormal proliferation compresses healthy surrounding tissue and eventually increases the stiffness of malignant lesions [3]. Therefore, mechanical properties of tissues are an important indicator to better understand various disease states and physiological conditions [4].
Atomic force microscopy (AFM) is a powerful tool for the quantification of mechanical properties of soft tissues and cells; however, it includes the risk that samples are regionally damaged by the compression from a scanning cantilever, and has high cost and involves timeconsuming operation [5,6]. Wave propagation has been extensively used to evaluate biomechanical properties of soft tissues by using ultrasound shear wave elastography (SWE) [7][8][9][10][11] and magnetic resonance elastography (MRE) [12,13] in clinical applications. Generally, the spatial scales of SWE and MRE images are in the ranges of millimeters to centimeters such that it constrains investigations to macroscopic levels with an organ-sized field-of-view (FOV) [4,14]. Optical coherence tomography (OCT) is an alternative modality used to image tissue displacements on a micron-scale [14]. A number of advantages of OCT imaging include it being noninvasive, noncontact, fast, having high spatial resolution, and being sensitive to the topology of a surface. OCT-based elastography, named optical coherence elastography (OCE), was first proposed by J. Schmitt in 1998 to assess mechanical properties of the muscle and fat layers in a porcine meat sample and skin tissue by using compressive stress [15]. Numerous reports have shown the capability of OCE equipped with various excitation tools to generate shear wave or surface wave propagation in two-dimensional (2D) space for characterizing mechanical properties of tissues [16]. So far, it has been widely used for investigation of mechanical properties in biological tissues and biological fluids during the last two decades [16][17][18][19]. In wave-based OCE techniques, the spatial resolution closely depends on excitation frequencies, pulse duration for transient excitation and mechanical properties of samples (e.g., stiff or soft).
Many OCE experiments involve some form of contact method or are ex-vivo studies, i.e., taking a non-living tissue sample for measurement. The excitation tools either directly contact samples or contact containers used to hold samples, like a Petri dish or glass plate. For transient excitations with a contact method, Li, et al., used a piezo-vibrator contacting at approximately a 45°angle over a sample surface to measure surface waves and determine the elastic properties in skin, chicken breast and heterogeneous phantoms [20,21]. In a non-contact method, Ambroziński, et al., utilized a home-made air-coupled piezoelectric transducer to generate an air puff-based acoustic radiation force on an ex vivo porcine cornea for evaluating wave propagation and phase velocities [22]. A similar study was also conducted by S. Wang and K. V. Larin to assess the biomechanics of ex vivo rabbit cornea by using air-puffs generated via an air-gate controller [23]. Z. Jin et al. also used air-coupled ultrasound to measure 3D corneal elasticity in-vivo in a 3D manner [24]. A needle tip as a contact-based loading was reported for breast tissue characterizations [25,26]. Acoustic radiation force (ARF) is a reliable excitation method broadly used to evaluate mechanical properties of biological tissues and biological fluids such as carotid artery [27], porcine blood and plasma [18], cornea [28][29][30], retinal tissue [31] and crystalline lens [32][33][34].
However, using transient excitation methods in OCE can increase the difficulty in a translational setting for tissue engineering applications due to the need for complicated alignment and smaller samples. Besides, transient excitations can result in generating undesired reflected waves from rigid boundaries, and could increase the difficulty to evaluate correctly estimate wave velocities and affect the evaluation of mechanical properties [35]. Harmonic OCE has been reported that utilize waves of a fixed frequency in confined media. F. Zvietcovich, et al., used two vibrators with 400 Hz and 400.4 Hz harmonic frequencies respectively as excitation signals for generating waves in a heterogeneous phantom to evaluate shear wave velocities and mechanical measurements [36]. P. Meemon, et al., also used a similar system to distinguish two components of a phantom based on crawling waves [37]. This approach is limited by the sample size needed to allow waves to interfere, which may be relatively large compared with samples used for tissue engineering. Another approach took advantage of reverberant wave behavior to estimate the elasticity of individual corneal layers [38]; however, a customized multi-pronged excitation ring is required and contact with tissue cannot be avoided. On the other hand, those method could be difficultly implemented for tissue engineering applications due to samples usually consisted of hydrogel-tissue hybrid placed in a limited culture space (35 mm diameter wells as an example).
In this study, we report a new concept called 2D dynamic vibration OCE (DV-OCE) technique to evaluate mechanical properties of the samples in a confined environment to mimic the situation for tissue engineering applications. This proposed method is non-contact with respect to the sample and avoids a tedious excitation alignment process. This method could potentially be used for tissue engineering applications to measure mechanical properties of tissues growing in culture. A coil-based vibrator with harmonic excitation was used to generate traveling waves in two homogeneous materials of different thicknesses (7% and 14% gelatin, 10 and 2 mm thickness), a heterogeneous material and ex vivo porcine kidney tissue. The samples in 35 mm diameter Petri dishes were simply placed on the surface of vibrator without complicated excitation focusing alignments. A customized 3D OCT scan pattern was used to create 4D wave propagations (x, y, z, t). In addition, 2D wavenumber filter, multi-angle directional filter and 2D local wave velocity algorithm were implemented to determine wave velocities. Transient ARF-OCE was utilized to compare with our new 2D dynamic vibration OCE technique. Numerical simulations in the two homogeneous materials were performed to validate the experimental results. Our proposed 2D DV-OCE technique is robust in determining the wave velocities in the homogenous phantoms, the heterogenous phantom and the biological tissue without tedious focus alignment and contactless. The proposed method is a promising technique that could be applied for tissue engineering applications and laboratory translation under easy-to-setup procedures [1].

Fabrication of homogeneous phantoms and a heterogeneous phantom
A total of five phantoms were fabricated for this study. Four phantoms were homogeneous phantoms and one is a heterogeneous phantom. The two gelatin concentrations, 7% v/v and 14% v/v, with 10 mm thickness were fabricated by using gelatin powder (gel strength 300 type A, G2500-1KG, Sigma-Aldrich, St. Louis, MO, USA) to make the homogeneous phantoms and 1 g of titanium dioxide (TiO 2 , ReagentPlus ≥ 99%, Sigma-Aldrich, St. Louis, MO, USA) was used to provide optical scatterers. To provide cases that may be similar to those in tissue engineering applications, we created 7% v/v and 14% v/v homogeneous phantoms with 2 mm thickness as the advanced cases. A total volume of 100 mL tap water in a 500 mL beaker was heated to approximate 70°C. The 14% v/v gelatin powders and 1 g TiO 2 were added with stirring to the beaker for approximately five minutes to homogenize the solution. The mixed solution was placed in a de-gassing chamber to remove small bubbles in the fluid.
After de-gassing, the mixed gelatin solution was poured into two regular 35 mm diameter × 10 mm height Petri dishes (for 10 mm and 2 mm thickness) for DV-OCE experiments and a custom 85 mm diameter × 10 mm height Petri dish with Mylar film bottom for ARF-OCE experiments to obtain the comparison results. The ARF has a small attenuation due to only 100 µm thickness of the Mylar film, which can be considered acoustically transparent [5,39,40]. The Petri dishes were transferred to a 4°C refrigerator for congealing. The above procedure was repeated for the 7% v/v homogeneous gelatin phantom. For the fabrication of the heterogeneous phantom, the same 14% v/v gelatin solution was poured into another regular 35 mm diameter × 10 mm height Petri dish with a 2 mm diameter pillar placed in the middle of the Petri dish to create a small circular well for a 2 mm diameter inclusion. The pillar was removed once the 14% v/v gelatin solution was completely congealed, and the 7% v/v gelatin solution, following the above fabrication description, was poured into the small circular well. The heterogeneous phantom had an approximate 10 mm thickness.

Porcine kidney tissue preparation
A porcine kidney was excised immediately after sacrifice and frozen. The kidney was thawed and placed in saline solution at room temperature. The kidney was sliced in half in the long dimension and reflected. A slice of renal cortex with the approximate dimensions of 20 × 20 × 4 mm was cut for the experiment.

Three-dimensional (3D) numerical models
In order to validate our method, numerical simulations were performed using a finite element method. An explicit solver was used to solve the differential elastic wave equations in OnScale (OnScale, Redwood City, California). Two three-dimensional models with elastic material properties were examined. The models, consisting of 2,899,232 grid elements and 2,962,971 grid nodes, were used to create simulated phantoms with 35 mm diameter × 10 mm thickness. To reduce the computational burden and solution time, the models were modelled as axisymmetric models in two in-plane axes. Figure 1(a) presents the top and front views of the simulated phantom with specified the type and locations of the boundary conditions. The axes of x and y are symmetry. By assuming that the Petri dish is infinitely stiffer than the tested samples, we can use a fixed displacement boundary condition applied at the bottom and sides of the numerical phantoms. The mechanical transient analyses were performed using an explicit solver. The simulation time step is 0.0523 µs. Due to ample research demonstrating that experimental results from mechanical testing are highly consistent with ARF-based shear wave elastography results in homogeneous samples, the results from the ARF-OCE measurements were used as a reference result for this study [41,42]. The parameters of the numerical simulations were based on the results from the ARF-OCE experiments. The shear wave speed was set to be 0.94 m/s for 7%, and 1.66 m/s for 14% gelatin concentrations, respectively. The Young's modulus E can be simply estimated by 3ρ(1.05C) 2 , where C is the group velocity of surface waves, and the 1.05 factor is a correction factor for surface waves on an air-solid interface [43]. The compressional wave speed was 1530 m/s, and density ρ was 1000 kg/m 3 . The hexahedron element size used for numerical spatial discretization was equal to 15 elements per wavelength for the 14% gelatin material, which corresponded to 0.102 mm. A harmonic force excitation using a 300 Hz sinusoidal function was applied to the bottom surface of the model in an out-of-plane direction. Particle velocity motion data were measured from the top surface and further processed.

Three-dimensional (3D) dynamic vibration optical coherence elastography (DV-OCE) system
The basic optical layout of the spectral-domain optical coherence tomography (SD-OCT) system is illustrated in Fig. 2. A low coherence broadband source is split into a reference beam directed toward a stationary reference mirror and a sample beam directed toward the samples. According to the basic principle of the Michelson interferometer, the back-reflected and back-scattered light from samples and retroreflected light from the reference mirror are recombined by a coupler. Eventually, a spectral interferogram is formed by spectrometers and data were collected by a frame grabber card. The fast Fourier transform (FFT) is utilized to form an A-scan from the receiver array in the SD-OCT system. In the study, the SD-OCT system (TEL320C1, Thorlabs Inc., Newton, NJ, USA) is equipped with a 1300 nm source and LK3 lens kit (Thorlabs Inc., Newton, NJ, USA) to provide a 10 mm × 10 mm FOV and the maximum imaging depth for the investigated phantoms and the kidney tissue were approximately 400-500 µm in this study. Two excitation methods, transient and harmonic excitation, were used in the study. For the transient excitation, a 7.5 MHz highly focused single-element transducer (ISO703HR, Valpey-Fisher, Hopkinton, MA, USA) with a 3,750 cycle burst (0.5 ms) was used to provide an ARF exerted on the surface of the samples to generate Rayleigh wave propagation. The focal distance of the transducer is 11.84 mm measured by a pulse-echo test and the f -number of 1.07 was obtained by the definition of the focal distance divided by the 11 mm aperture size. The sinusoid burst signal, provided by function generator 2 (33250A, Agilent, Santa Clara, CA, USA), was amplified 50 dB by a radiofrequency (RF) power amplifier (240L, Electronics and Innovation, LTD, Rochester, NY, USA) to drive the transducer. For the harmonic excitation, a coil-based vibrator was used to generate waves propagating in multiple directions. Function generator 2 was set to apply a 300 Hz signal with a 15-cycle burst (50 ms) and this signal was amplified by a 26 dB stereo amplifier (D150A, Crown, Northridge, CA, USA) to drive the vibrator. Function generators 3 and 1 (both 33500B, Keysight, Santa Rosa, CA, USA) were used to provide a trigger signal for the SD-OCT and a master trigger to synchronize the timing for the entire system, respectively.
A 10 kHz A-scan rate was selected to obtain a better penetration depth and a customized two-dimensional (2D) acquisition with 600 axial scans (M-B scan) was set to facilitate high frame rates to track dynamic processes in the space-time domain [44]. Each M-scan provides a 60 ms recording where the first 50 ms is for recording wave propagations and last 10 ms is a rest time for the vibrator to avoiding incorrect vibrations described by Faraday's law of induction at the beginning and end of the applied drive voltage signal. A total of 100 discrete positions (100 µm scanning step) in both the lateral, x, and elevational, y, directions were evaluated by using our customized three-dimensional (3D) OCT acquisition pattern to build a four-dimensional (4D) wave propagation map. Each 4D acquisition is composed of a dataset with dimensions (z, x, y, t), where z is imaging depth, x is the index for points in the B-scan, y is the index for elevational positions and t is the index for the M-scan acquisitions.
The 1-D autocorrelation was used to determine wave motion inside the phantoms and the kidney tissue [45,46]. Each pixel includes a real value and an imaginary value (in-phase/quadrature, IQ) from which magnitude and phase can be calculated. The phase unwrapping algorithm was done right after motion estimation by using the MATLAB 'unwrap' function (Mathworks, Natick, MA) when phase wraps were observed. A 3D wave propagation at a specific depth z was selected to create a C-plane (x, y, t) image. Figures 3(a)   To evaluate wave velocity in a small container, the omni-directional wave pattern needs to be decomposed by a multi-angle directional filter and the undesired frequency components should be filtered by a wavenumber filter. The 3D fast Fourier transform (3D-FFT) was applied to the 2D wave propagations to obtain a 3D k-space (k x , k y , f ). Before performing 3D-FFT, Tukey windows with a cosine taper parameter of 0.1 were applied to spatio-temporal signals to minimize the edge effects caused by the Fourier transform [47]. Figure 4(a) shows an example of a 2D k-space with maximum energy at the excitation frequency of 300 Hz, which derives from the wave propagations in 14% gelatin homogeneous phantom. In addition, a high cutoff wavenumber k h (large dash-line circle) and a low cutoff wavenumber k l (small dash-line circle) can be selected in appropriate spatial ranges. In general, a common range of shear wave velocity in soft tissues ranges from approximately 1-3.5 m/s, for pancreas, kidney or fibrotic liver [1,48,49]. Therefore, a 2D spatial wavenumber filter (Butterworth filter) [47] with a spatial frequency limit k l of 100 m −1 and k h of 400 m −1 were selected to limit upper speeds to 3 m/s and a lower speed limit as 0.75 m/s for the excitation frequency in the study. Furthermore, the filter is able to remove high frequency noise and compressional wave components [50]. The 6 th order wavenumber filter illustrated in Fig. 4(b) was used to avoid Gibbs ringing effect [47] and the filtering result is displayed in Fig. 4(c). Due to the multi-directional nature of the wave field due to reflections from the hard wall of a Petri dish or tissue boundary, the waves need to be decomposed by a multi-angle directional filter for the time-of-flight wave velocity estimation [47,51,52]. Figure 4(d) shows the multi-angle directional filter response in the 90°direction. For each direction, the wave propagating with the primary direction plus an angular range is able to pass but will be suppressed in other directions, presented in Fig. 4(e). After the two filters are applied, the 3D k-space was transformed to the spatio-temporal domain by using an inverse 3D-FFT. Figure 5 demonstrates the wave propagation at 17.5 ms in a 14% gelatin phantom, presented in Fig. 4(c), which were decomposed into 4 directions: south-to-north (SN, 90°, Fig. 5(a)), east-to-west (EW, 180°, Fig. 5(b)), west-to-east (WE, 0°, Fig. 5(d)) and north-to-south (NS, 270°, Fig. 5(e)), after applying the 2D wavenumber filter and 4-angle directional filter. The waves in different directions were separately processed by using a 2D local wave velocity algorithm, which is based on cross-correlation to estimate the time delays between temporal signals at different locations in the x and y directions with a certain distance between them to calculate the local wave velocity [53,54]. This is performed in two-dimensions (lateral and axial for ultrasound and within the two-dimensional plane of the surface when measured with DV-OCE). The 2D local wave velocity algorithm incorporates a sliding patch within a larger window to avoid measurement errors caused by small step size calculations and to increase the robustness of the algorithm by using a 2D processing window [52,55]. A small window size results in a higher possibility for error in time delay estimation but obtains a higher resolution result [47]. In the study, the window was set as 2 mm (20 pixels) and the patch was set to be 1 mm (10 pixels). The final wave velocity can be obtained by averaging the measured wave velocity in each direction. Weighting of the wave velocity maps was carried out based on the wave energy as described by Song et al. [52]. A flowchart presented in Fig. 6 illustrates the whole process of wave velocity estimation. Interested readers can refer to our previous literature for the details of 2D wavenumber filter, multi-angle directional filter and 2D local wave velocity algorithm [47,50,52,54]. The Poisson's ratio ν for incompressible materials and the density ρ for soft tissue are typically assumed to be 0.50 and 1000 kg/m 3 , respectively. The above procedure was programmed by using MATLAB R2019a software implemented in the desktop computer with Intel Core i5-8500 CPU at 3 GHz processor, 8 GB memory and 64 bit Windows 10 operating system.

Results and discussion
The transient ARF-OCE was used to evaluate the group velocities, k-space and dispersion analysis in 7% and 14% homogeneous gelatin phantoms with 10 mm thickness for comparison results. According to the experiments by using ARF-OCE, the group velocities, c g , in 7% and 14% were approximately 0.94 m/s and 1.66 m/s, respectively, presented in Fig. 7(a) and Fig. 7(d). The k-space (frequency-wavenumber) domain can be obtained by using a 2D Fourier transform to illustrate the energy distribution of the motion from the 7% and 14% gelatin phantoms, which is displayed in Fig. 7(b) and Fig. 7(e), respectively. Based on the k-space, the phase velocities associated with a range of frequency components in 7% and 14% can be calculated. With a range of the frequencies from 200-600 Hz, the mean phase velocity was 0.97 ± 0.003 m/s in 7% and 1.57 ± 0.022 m/s in 14% gelatin phantom, displayed in Fig. 6(c) and Fig. 6(f), respectively. We are also interested in the phase velocity, c p , at frequency of 300 Hz compared with the results from the 2D DV-OCE. The phase velocities in 7% and 14% at 300 Hz were 0.97 m/s and 1.57 m/s, respectively. We observed that the group velocities and the phase velocities were highly consistent and the fluctuation of the standard deviation (SD) in the phase velocity were tiny varying with frequencies, which implies the media are elastic and locally homogeneous. The 2D depiction of interacting waves (C-plane) at approximately 200 µm under the surface of the phantoms were illustrated in Fig. 8(a) for 7% and in Fig. 8(e) for 14% homogeneous phantom as examples. Visualization 1 and Visualization 2 display the 2D wave propagations excited by coil-based vibrator with 300 Hz in the 7% and 14% phantoms, respectively. Based on the comparison results by using ARF, the wave velocities in the numerical simulations were set as 0.94 m/s for 7% and 1.66 m/s for 14%, presented in Fig. 8(b) and (f). The simulation results showed similar wave propagation patterns with our experimental results. Visualization 3 and Visualization 4 exhibit the numerical simulations of 2D wave propagations in 7% and 14% phantom, respectively. In addition, the 4D wave propagation visualizations were constructed by using our customized 3D OCT acquisition pattern to provide a better demonstration of how waves propagate inside materials. Figures 8(c) and (g) illustrate the 4D wave propagation in the 7% and 14% phantoms, respectively, and the videos of 7% and 14% can be found in the Visualization 5 and Visualization 6, respectively. We also performed the numerical simulations of 4D wave interaction behaviors for 7% and 14% with the 35 mm diameter × 10 mm thickness cylindrical structure, presented in Fig. 8(d) (Visualization 7) and Fig. 8(h) (Visualization 8), respectively. The dashed-line square restricts the FOV into a 10 × 10 mm area, the same as the FOV in OCT. Based on our 2D and 4D numerical analysis, we report that the wave behaviors in simulation results were in a good agreement with those in the OCT experiments.  Figure 9 illustrates the 2D wave velocities in the homogeneous phantoms by using our proposed DV-OCE method after performing the 2D wavenumber filter, multi-angle directional filter and 2D local wave velocity algorithm. Based on the homogeneous phantom experiments, the 2D wave velocity was 1.09 ± 0.018 m/s in 7% and 1.53 ± 0.012 m/s in 14%, presented in Fig. 9(a) and Fig. 9(d), respectively. Furthermore, 2D wave propagation results from the numerical simulation were evaluated to validate the results for the phantom studies. For the wave propagations produced by the numerical simulation, the 2D wave velocity was 0.99 ± 0.002 m/s in 7% ( Fig. 9(b)) and 1.53 ± 0.005 m/s in 14% (Fig. 9(e)). To evaluate the level of artifacts, a profile of the wave velocity distribution in the 7% and 14% gelatin materials along the x direction at y = 5 mm was presented in Fig. 9(c) and Fig. 9(f), respectively. As the results described above, the results from the phantom experiments represent a good agreement with the results from numerical simulations and from the transient OCE. We can observe that the speed maps from simulation showed highly uniform speed distributions due to the perfect symmetric patterns of the simulated interacted waves. The results from the experiments have relatively more variability. Using ARF as an excitation source can obtain reliable measurements; however, an artifact effect occupies a certain amount of space at the excitation location and limits the wave velocity maps due to the absence of propagating waves in the excitation region and sacrifices the distance of wave traveling, which could result in incorrect evaluation of wave velocities [39]. Noticeably, no artifact effect is presented in the vibration results. Besides, the reflected waves also affect the accuracy in the evaluation of wave velocity by using ARF-OCE if samples are too small. In DV-OCE in a homogeneous medium, the boundary conditions cause wave reflections but do not generally affect the estimation of local wave velocity because the reflected waves can be decomposed into different components that can be evaluated utilizing the band-pass and directional filters. This assumes that the field-of-view has 1-2 wavelengths for accurate wave velocity estimation. In the study, the proposed 2D dynamic vibration OCE can provide a full-field 2D wave velocity map with no artifact effect from the excitation source. In general, samples from tissue engineering are often thin and usually embedded within a hydrogel to make the entire material semi-solid [35]. To achieve the possibility of applying this method for tissue engineering applications, our proposed 2D dynamic vibration OCE was applied to 7% and 14% gelatin with 35 mm diameters × 2 mm thickness phantoms as an advanced condition to mimic tissue engineering conditions. However, in phase-sensitive OCE, the measured particle velocities can be corrupted by phase wraps, especially in a thin layer structure [44,56,57]. The phase wrap will happen once scatterers travel such that the phase shift is greater than π radians. The interacted waves with a 2D phase wrap in 7% and 14% phantom was illustrated in Fig. 10(a) (Visualization 9) and Fig. 10(e) (Visualization 10) with the yellow arrows, and it also can be observed in spatio-temporal maps presented in Fig. 10(b) and Fig. 10(f) with the white rectangles. Fortunately, the phase wrap is a fundamental problem and had been addressed using phase-unwrapping algorithms or reduced excitation energy to make displacement amplitudes fall into nanometer ranges [44]. Figure 10(c) (Visualization 11) and Fig. 10(g) (Visualization 12) show the clean standing wave propagation with the phase unwrapped. Interested readers can refer to the previous literatures for phase unwrapping algorithm [44,56]. After phase unwrapping, the 2D wave velocities by using the 2D DV-OCE were determined as 1.07 ± 0.022 m/s for 7% (Fig. 10(d)) and 1.63 ± 0.011 m/s for 14% (Fig. 10(h)) gelatin phantom. These results corresponded well with the results from the numerical simulations and the transient OCE with 2 mm thickness, presented in Fig. 7(d) and (j). By using ARF-OCE, the mean phase velocities within the frequency range between 200 Hz and 600 Hz were determined as 1.04 ± 0.015 m/s for 7% with 2 mm thickness and 1.62 ± 0.105 m/s for 14% with 2 mm thickness gelatin phantom. Experiments based on a stricter condition for phantom thickness, less than 2 mm, will be conducted in our future studies.  Figure 11(a) and Fig. 11(b) illustrate a heterogeneous gelatin phantom with 2 mm diameter inclusion, and Fig. 11(c) displays a 2D-OCT image with the 2 mm diameter inclusion. The inclusion and background were 7% and 14% gelatin concentration, respectively. The 2D wave velocities evaluated using the 2D DV-OCE in the heterogeneous material was 1.10 ± 0.007 m/s for the inclusion and 1.55 ± 0.012 m/s for background, shown in Fig. 11(d), which is a good agreement with the results from ARF-OCE with 10 mm thickness, illustrated in Fig. 7(a) and (g). A profile of the wave velocity distribution along the x direction at y = 5 mm was presented for evaluating the presences of artifacts. Moreover, the inclusion size in the 2D speed map and in the 2D OCT image were highly consistent. The 2D wave velocities determined by transient ARF-OCE and by the proposed 2D vibration OCE were incorporated in the Table 1. Similar research was reported by X. Liang, et al., by using a piezoelectric (PZT) stack with multiple excitation frequencies to compress ex vivo rat tumor tissue against a glass window, and to highlight tissue components by strain images [49]. A frequency of 313 Hz was used as an excitation frequency to locate the tumor tissue region, which provided us with the idea that using an excitation frequency around 300 Hz for evaluating soft tissues would be appropriate. Another article reported strain imaging of prostate cancer using an electromagnetic actuator to compress the specimens with a contact method [58]. In this study, our proposed 2D dynamic vibration OCE could be a promising new method to evaluate 2D wave velocities and mechanical properties of small samples in a limited space, i.e., 35 mm Petri dish, with no direct contact with the sample and no complicated excitation alignments issues.  A biological tissue was tested in this study. Figure 12(a) shows the experiment of the proposed 2D DV-OCE applied to a porcine kidney cortex sample. A 2D wave propagation at 20 ms and a spatio-temporal map at y = 5 mm were illustrated in Fig. 12(b) (Visualization 13) and Fig. 12(c) as examples. The group wave velocity in the kidney tissue was 1.15 m/s measured by the transient ARF-OCE, presented in Fig. 12(e). The 2D wave velocity in the kidney tissue was determined as 1.11 ± 0.005 m/s, displayed in Fig. 12(d), by the proposed method, demonstrating good agreement between the methods. The dispersion curve evaluated by transient ARF-OCE is exhibited in Fig. 12(f) and shows the phase velocities associated with various frequencies. The phase velocity at the frequency of 300 Hz was 1.25 m/s, similar to that obtained with DV-OCE method. I. Grosum, et al. reported that the wave velocity in healthy human kidney is approximately 1.2 m/s [59]. In the study, our proposed 2D DV-OCE was capable of evaluating biological tissues with a satisfactory result. A limitation of this proposed method is that highly viscous materials or large bulk materials may not be appropriate due to wave amplitude dissipation, which may not yield interacted waves. There are numerous advantages of the proposed 2D DV-OCE. First, no artifact will be produced by using the vibration excitation source as wave velocities are evaluated. The artifact from ARF excitation occupied a sizable portion of the image, especially for a FOV of OCE, and this is eliminated with vibration [39]. Decreasing the ARF energy might alleviate the artifact effect; however, a weaker amplitude wave traveling inside samples results in incorrect evaluation of wave velocities, which would need to be considered. In addition, this is a contactless method, which would be appropriate for samples in culture where direct contact with the samples could change the microenvironment and intrinsic mechanical properties [60]. This is in contrast with a modality like AFM that requires direct contact with the sample. The proposed OCE method could allow for longitudinal studies for these types of samples, which may reduce costs for studies. Compared with mainstream excitation methods, coil-based vibrator is relatively economical without a RF amplifier, a PZT controller or an air-gate controller. Moreover, for the ARF excitation, acoustic energy will be dramatically lost as a sample is placed in a regular Petri dish. We have found that a Mylar film is needed to replace the bottom of a Petri dish, to avoid issues with reflected ultrasound waves [40]. In the proposed method, samples can be placed in a regular Petri dish and no complicated excitation focusing alignment is needed. According to the advantages above, the proposed 2D DV-OCE would be a promising method for mechanical evaluation in tissue engineering with easy-to-setup operations [1,35].
In this study, we examined lateral heterogeneity (Fig. 11). It is reasonable that the DV-OCE method could also be used to examine axial heterogeneity. The axial heterogeneity could take a couple of forms, either a layered medium with each layer having different mechanical properties or the example given with a stiff inclusion in a soft background, where the inclusion is at some depth. In theory, our DV-OCE system should have the ability to evaluate the wave velocity in layered samples. However, it could be a little bit challenging. First, our TEL320C1 OCT can usually penetrate ∼700-1500 µm depth depending on optical scatterers in the medium. Provided that wave motion is present and the motion can be reliably detected with the OCT signal, then wave velocity could be reconstructed. However, with the finite thicknesses of layers or samples, the wave behavior may be more complicated and may have some guided wave features that would need to be accounted for in the inversion. However, this is an interesting topic that may be explored in our future works, which can be potentially used for multilayered engineered tissue sheets for applications in tissue regeneration.

Conclusion
In the study, we reported a new technique, called two-dimensional (2D) dynamic vibration OCE (DV-OCE), to evaluate 2D wave velocities with no tedious focusing alignment procedures and contactless. The 3D-FFT was utilized to evaluate the interacted waves in 3D k-space. A spatial 2D wavenumber filter and multi-angle directional filter were employed to decompose the waves into four individual directions. The 2D local wave velocity algorithm was used to evaluate a 2D wave velocity map. The ARF-OCE was used to evaluate wave velocities for comparison. Numerical simulations were performed to validate the proposed method. The experimental results were in a good agreement with the results from ARF-OCE and numerical simulations. Using DV-OCE on thinner biological samples such as thin gelatin hydrogel embedded with tissues will be investigated in the future experiments. This may give rise to more complicated wave behavior and higher frequency excitations will be explored. On the other hand, if tissue size or the size of the field of view is small compared to excitation wavelength, the reconstruction accuracy could be affected. We chose the excitation frequency to make sure we had at least one wavelength in the field-of-view. The inclusion shape in theory may not be affected; however, the further experiments will be necessary to prove it. Our proposed 2D DV-OCE could potentially pave the way for mechanical evaluation in tissue engineering applications with easy-to-setup and contactless approaches.