Characterization of natural frequencies from nanoscale tissue oscillations using dynamic optical coherence elastography

: We demonstrate the use of OCT-based elastography for soft-tissue characterization using natural frequency oscillations. Sub-micrometer to sub-nanometer oscillations were induced in tissue phantoms and human cornea in vivo by perpendicular air-pulse stimulation and observed by common-path OCT imaging (sensitivity: 0.24 nm). Natural frequency and damping ratio were acquired in temporal and frequency domains using a single degree of freedom method. The dominant natural frequency was constant for diﬀerent stimulation pressures (4-32 Pa) and measured distances (0.3-5.3 mm), and decreased as the sample thickness increased. The dominant natural frequencies of 0.75-2% agar phantoms were 127-774 Hz (mean coeﬃcient of variation [CV]: 0.9%), and correlated with the square root of Young’s moduli (16.5-117.8 kPa, mean CV: 5.8%). These preliminary studies show repeatable in vivo corneal natural frequency measurements (259 Hz, CV: 1.9%). This novel OCE approach can distinguish tissues and materials with diﬀerent mechanical properties using the small-amplitude tissue oscillation features, and is suitable for characterizing delicate tissues in vivo such as the eye.


Introduction
Soft tissue biomechanical properties (e.g. stiffness or Young's modulus) are related to tissue health, and disease progression often changes the biomechanical properties of the affected tissues [1,2]. In a clinical setting, physical palpation is used to diagnose and locate diseases by feeling the stiffness changes in tissue. Analogous to palpation, elastography methods use mechanical tissue stimulation to assess subtle stiffness changes in soft tissues that may be caused by disease, such as the presence of a tumor [1][2][3]. Unlike the tactile nature of palpation, elastography relies on non-invasive high-resolution imaging techniques and can provide results that are both objective and quantitative.
Optical coherence elastography (OCE) is a recently developed technique that combines a mechanical loading system to induce tissue displacement and an optical coherence tomography (OCT) imaging system to detect the resulting tissue deformations. Tissue elastic properties can be reconstructed based on the relation between the tissue response and stimulation [4]. Recent advances in elastography methods include a variety of dynamic stimulation methods [5][6][7][8][9][10][11][12][13][14][15][16] and more sensitive phase detection methods [17][18][19][20] that have enabled the visualization and analysis of mechanical wave propagation in dynamic OCE [7,21]. The computational methods required to reconstruct tissue elasticity remains an active research area. Most of the current OCE analytical models are inherited from ultrasound [22][23][24] and MRI-based [25,26] elastography, including the commonly used shear-wave model to estimate Young's modulus based on transverse wave propagation velocity [27]. However, due to the differences in the detection field and penetration depth, the mechanical models used in ultrasound and MRI methods may not be appropriate for OCE applications in all tissues [28]. In tissues with complex geometries and multiple layers, such as cornea and skin, mechanical waves traveling along the surface contain multiple highly dispersive Rayleigh-Lamb components and become very complex compared to simple Rayleigh waves [29,30]. In this case, translation of the measured wave propagation speed into the shear wave model could lead to inaccurate estimation of tissue Young's modulus [28]. We recently proposed a modified Rayleigh-Lamb wave model to quantitatively assess the corneal viscoelasticity [30,31]. This method is limited in a first-order assumption that the cornea is isotropic, homogenous, and has a flat curvature. The development of robust computational methods and tissue modeling techniques is important to provide more robust tissue elasticity estimation from dynamic OCE [21].
Tissue natural frequency is an intrinsic property, and is defined as the frequency at which tissue tends to oscillate when disturbed [32]. Natural frequency oscillation in response to the excitation force is closely related to tissue elastic properties. Resonant ultrasound spectroscopy has been employed for decades to measure the resonant frequencies of samples with known size and mass, using oscillatory acoustic radiation force [33,34]. They are of considerable interest, but still limited in detection resolution [35]. A phase-sensitive OCE approach has been recently applied to detect and analyze the vibrational or resonant responses from samples and tissue with greater resolution, using a number of different tissue modulating forces including: acoustic radiation from ultrasound transducers [35], piezoelectric actuators [36] or mechanical wave drivers [37], magnetic force from embedded nanoparticle transducers [38][39][40], and sound waves from a speaker [8] etc. These dynamic OCE methods have demonstrated enhanced frequency-based B-scan contrast and volumetric imaging at certain excitation frequencies for ex vivo tissues and phantoms [35][36][37]40], and high-resolution measurement of resonant natural frequencies by sweeping the driving frequencies in step [35,36,[38][39][40]. Previous studies have demonstrated that the natural frequency is linearly related to the square root of Young's modulus in a simple elastic model [35,38].
However, the quantification of natural frequency for sensitive tissues (e.g. eye) using OCE is still constrained by the stimulation method and the sensitivity of OCT detection. The stimulation methods using any mechanical contact [36,37] or label agents [38][39][40] are largely unsafe and unsuitable for in vivo ocular OCE measurement. Most of the OCE methods sweep sinusoidal excitations over a defined frequency range to achieve the spectroscopy-based response [35][36][37]39,40]. Using square-wave tissue modulation [38], or impulse stimulation functions can provide a wide range of stimulus frequencies simultaneously. The requirement of frequency sweeping usually takes longer time, and could cause discomfort or harm during in vivo measurements for the human eye. The sound-induced OCE system was implemented to observe dominant response frequencies from bovine eyes, ex vivo [8]. These large-scale tissue vibrations (in the millimeter scale) require a large stimulation force that is potentially hazardous for ocular tissues in vivo. Instead of the sinusoidal excitation method, a single impulse stimulation force can provide a wide range of frequency excitation (e.g. 0 to kHz) simultaneously. A micro-scale air-pulse stimulation method was initially designed for ocular OCE applications [14], and was recently verified in our in vivo corneal elastography measurement [41]. The use of this micro-air-pulse stimulator can provide short duration (∼ 1 ms), micrometer-scale tissue displacements and a period of tissue damping oscillatory motion that ranges from sub-micrometer to nanometer-scale [42]. In previous work, we developed a relaxation model [12,13,43] to estimate tissue viscoelasticity by fitting exponential curves to model the tissue's primary deformation recovery response. This damped frequency response was limited to the primary deformation response only. Thus, this method was not a direct measurement of natural frequency, and is subject to assumptions and approximations.
We have recently introduced an OCE approach based on a higher resolution OCT technique and a perpendicular air-pulse stimulation method [42]. The newly developed common-path OCT detection method provides enhanced optical phase stability and detection sensitivity (0.24 ± 0.07 nm), and provides automatic compensation for polarization and dispersion. The improved lowforce (tens of Pascals) air-pulse tissue stimulation system was developed from an earlier oblique stimulation geometry [14,43] that now has a stimulation geometry normal to the surface. Loading normal to the surface improves tissue excitation efficiency, provides better wave propagation uniformity in radial directions, and simplifies modeling methods which can be used to derive the mechanical properties from the observed sample response [44]. Displacement amplitudes generated by this micro-force stimulation can be limited to sub-micrometer or a nanometer scale. This common-path OCT detection technique has shown better visualization and quantification of small-magnitude oscillations than conventional OCT [42]. The improved detection sensitivity of common-path OCT and the perpendicular micro-scale stimulation now enable direct observation of these small-magnitude damped oscillations [42].
Here we describe a natural frequency quantification method that employs small-magnitude damped oscillations by using the combination of perpendicular micro-air-pulse stimulation, high sensitivity common-path OCT detection, and a single degree of freedom (SDOF) model [36,45]. The oscillation features including the dominant natural frequencies, decay coefficients, and the damping ratios, can be analyzed using the OCE measurement and the SDOF quantification method. We verify the natural frequency concept using OCE measurements on agar phantoms at different stimulation pressures, at different measurement positions, and for various phantom concentrations and different thicknesses. We compare the measurement of natural frequencies and Young's moduli using the SDOF method and the elastic wave propagation method. We also report results form a pilot study using this novel OCE approach for in vivo measurements of human corneal natural frequency and damping ratios. We show that the small-amplitude oscillation measurements and the SDOF method can provide robust and precise quantification of the natural frequencies in biological tissues. It has the potential to be used in the further clinical applications, such as early disease detection and treatment evaluation.

Human subjects
A pilot study was performed on the left eye of a healthy subject (35 years). He had no ocular disease or surgical history, except myopia (left eye: -3 D). The intraocular pressure was measured as 13.1 mmHg using a Goldmann tonometer. The research protocol was reviewed and approved by the institutional review board (IRB) of the University of Alabama at Birmingham and adhered to the tenets of the Declaration of Helsinki.

Common-path PhS-OCE instrumentation
A common-path OCE system ( Fig. 1) was described previously [42]. Briefly, a low-force air-pulse stimulator [14] was set perpendicular to the tissue or sample surface to provide short duration (≤ 1 ms), localized (diameter of 150 µm), and low-pressure (0-60 Pa) stimulation. Each air pulse provided a range of ∼ 0-1.5 kHz excitation frequency. A common-path phase-sensitive OCT system was synchronized to the stimulation to detect the subtle displacements in response to the applied force. The light source of the OCT system was a superluminescent laser diode (SLD, D-855, Superlum Diodes Ltd.) with a central wavelength of 845 nm, and a waveband of 100 nm. The common-path design used a shared common optical path for the sample and reference arms with a reference plane defined as the optical surface of a 5-mm thick acrylic plate, kept proximal to the sample. A telecentric scan lens (LSM04-BB, Thorlabsc Inc. New Jersey, USA) was inserted between two-dimensional Galvo scanners (for simplicity, only one scanner is shown in Fig. 1) and the reference plane. The lens enabled illumination that was parallel to the optical axis and ensured uniform illumination to the sample during lateral scans. The structural resolution, as calibrated in air, was ∼3.3 µm in the axial direction and ∼7.8 µm in the lateral direction, and the maximum imaging depth was ∼6.76 mm. The detection sensitivity for the dynamic tissue displacements was dependent on phase stability, which was calibrated as ∼0.24 nm in the depth of 0.33 mm to 6.66 mm using a mirror (signal sensitivity: 102.4 dB to 66.4 dB) [42]. Common-path PhS-OCE combines an air pulse stimulator to induce micrometer-scale tissue deformation and a high-sensitivity common-path PhS-OCT to detect tissue response (sensitivity 0.24 nm) [42]. In common-path OCT, the interference signal is produced by combing returned light from the sample and a reference plane adjacent to the sample.

Dynamics of single degree of freedom spring-mass-damper system
We used a single degree of freedom (SDOF) model [36] to quantify the tissue oscillation dynamics. Figure 2(a) demonstrates an ideal SDOF spring-mass-damper system [45], where k is the spring stiffness coefficient, c is the viscous damping coefficient, and m is a mass.
Tissue natural frequency (f n ) is an intrinsic property that is determined by the tissue stiffness, mass, boundary conditions, thickness, shapes etc. In this SDOF model, the natural frequency f n can be calculated based on the spring stiffness k and the mass m, as The damping ratio is defined as Based on the values of ε, the response oscillation can be described as three different oscillation regimes: critical-damping (damping ration ε = 1), under-damping (0 ≤ ε < 1), and over-damping (ε > 1) [12,13,43]. Figure 2(b) demonstrates a decaying SDOF response of the spring-massdamper system in an under-damped condition (ε <1). The equation of motion to describe the free response of a SDOF system is: where y A (t) is a displacement of the center of mass. The solution of this equation is: where y E (t) is the envelope function, A is the maximum amplitude, B is the decay coefficient, and φ is a phase value: The natural frequency f n can be deduced based on the damped natural frequency f d and the damping ratio ε as When the damping ratio ε is small, the damped natural frequency f d is nearly equal to the undamped natural frequency f n . Figures 3(a) and (b) show the estimated damping ratio ε and the difference between the f n and f d based on Eqs. (6), (7). When the natural frequency f n is in the range of 50-1000 Hz, and the decay coefficient B is in the range from -100 to -10 s −1 , the damping ratio is smaller than 0.3 and the difference between the f n and f d is smaller than 3 Hz. In this situation, we can assume that the acquired damped natural frequency f d equal to the natural frequency f n .

Oscillation characterization using SDOF method
In PhS-OCE, the axial surface tissue displacement z(t J -t 0 ) of a point on an air/sample interface at time t j , relative to time t 0 is given by [48]: where ϕ z (t J -t 0 ) is the phase change, λ 0 is the center wavelength and n is the refractive index (n = 1 in air). A typical air-pulse induced displacement on the sample surface is demonstrated in Fig. 4(a), which was acquired from a 2% agar phantom (weight: 6.7 g, thickness: 7.1 mm) by common-path OCE with an applied pressure of 4 Pa. The figure shows a baseline period (0 ∼ 11.5 ms) followed by an initial surface displacement that is driven by the excitation force (primary deformation with amplitude A 0 ; 11.5 ms ∼ 13.1 ms). Then follows a recovery response period where the displacement amplitude returns from A 0 to zero for the first time (13.1 ms ∼ 14.6 ms), and a period of damped oscillations (14.6 ms ∼ 90 ms). The decay feature and the oscillation frequencies can be calculated from the damped oscillations depicted in the red-dashed window area (20 ms to 90 ms). Based on Eq. (5), the decay envelope can be described as y E (t) = A 1 e B(t -t 1 ) . A 1 is the decay amplitude, defined as the maximum negative displacement amplitude in the red-dashed window area, t 1 is the time when the maximum negative displacement A 1 occurs, and B is the decay coefficient (defined in Eq. (6)). A 1 corresponds to the general scale of the damping amplitudes and B corresponds to the damping speed of the oscillation amplitudes. In the representative oscillations shown in Fig. 4(a), the values for the exponential decay curve were A 1 = −0.123 µm, B = −37.5 s −1 , and R 2 = 0.98.
We used zero padding method to expand the displacement data from 90 ms detection period to 0.5 s, the stimulation interval between adjacent measurements. We then employed fast Fourier transform (FFT) method to analyze the oscillation frequencies with a frequency resolution of 2 Hz. In the frequency components ( Fig. 4(b)), the ∼ 20 Hz low frequency was identified previously as the phase noise caused by environmental factors such as vibration [42], the damped frequencies were in the range of 340-1230 Hz, and the dominant damped frequency f d was 776 Hz. Based on this dominant f d (776 Hz), and the fitted A 1 (−0.123 µm) and B (−37.5 s −1 ) values, we estimated the damping ratio ε as 0.008 based on Eqs. (6), (7). Since ε is very small, the natural frequency f n equals to the damped natural frequency f d . In Fig. 4(c) the original damping oscillation data is compared with the SDOF fitting data (Eq. (2)). Both the original data and the SDOF fitting data had the similar oscillation frequency and decay trend (R 2 = 0.81). The residual mismatch (root mean squared error: 0.01 µm) between the original data and the fitting data was because SDOF was a simplified method that discarded other frequency components. The residual errors can be reduced if more frequency components were used to describe the damping oscillation data using a multi degree of freedom (MDOF) method [38,45].

Phantom measurement
Natural frequency is an intrinsic property of tissues, determined by factors such as spring stiffness and mass (Eq. (1)), and is not determined by the stimulation force. For an ideal homogenous material, the natural frequency should be the same, regardless of the stimulation force and the measurement location. We verified the natural frequency concept using OCE measurements on agar phantoms at different stimulation pressures (Section 3.1), at different measurement positions (Section 3.2), and for various phantom concentrations and different thicknesses (Section 3.3). We also compared the measurements of natural frequencies and Young's moduli using the SDOF method and the elastic wave propagation method in Section 3.4. The A-line measurement speed was 70 kHz. Displacement profiles were obtained from 6000 A-lines (90 ms). The stimulation interval between sequential measurements was 500 ms for all the phantoms tested.

Oscillation frequency for different stimulation forces
The measurements were performed at the same location on the surface of 2% agar phantom, 0.3 mm from the excitation point. The air-pulse pressure was increased from 4 Pa to 32 Pa in step of 4 Pa, and five measurements were made for each pressure in M-mode (repeated A-scan acquisitions over time at the same location). Figure 5(a) shows the surface displacements for each pressure as well as the oscillation features in response the range of stimulation pressures tested (4 Pa to 32 Pa). As the stimulation pressure was increased, the primary displacement amplitude A 0 increased, but the damping oscillation behavior remained similar. Figure 5(b) shows that the primary displacement amplitudes A 0 ranged from -0.2 µm to -4.0 µm (mean coefficient of variation (CV) 1.87%), and fit linearly in response to the applied force (y = −0.144x + 0.638, R 2 = 0.983). Figure 5(c) shows the FFT for each stimulation force. The dominant damped frequency f d for all of the measurements was 776 Hz. Figure 5(d) and 5(e) show the results of fitting the damped oscillations. The decay amplitudes A 1 ranged from -0.12 µm to -0.17 µm (mean CV 3.7%) and varied as the force changed. There was no observable relation between the decay features (A 1 and B) and the primary displacement amplitudes A 0 . The decay coefficients B and the damping ratio ε were not sensitive to the applied force. The mean value of the B and ε for all the measurements were -37.6 s −1 and 0.008 (CV: 7.7%), respectively. We assumed that the dominant natural frequency f n was close to the dominant damped natural frequency f d as 776 Hz since ε is small.

Oscillation frequency for different measurement positions
The measurements were designed to assess the effect of sample position on oscillatory motions. The measurements ranged between 0.3 mm to 5.3 mm from the stimulation point; the air-pulse stimulation was fixed at 20 Pa and the measurements were repeated 5 times for each position. Figure 6(a) shows the surface displacements for each position and the enlarged area shows oscillation features with 0.3 mm at the bottom and 5.3 mm at the top for the windowed area. As the measurement distance was increased, the primary displacement amplitude (A 0 ) decreased, the tissue oscillation behavior remained similar, and the time when displacements occurred were delayed. These observed time-shifts and the measurement positions were used to calculate the wave propagation speed and Young's modulus [27]. Figure 6(b) shows that as the measurement distance increased, A 0 decreased in absolute value from -1.86 µm to -0.04 µm (mean CV 5.4%). Employing a previously reported method [49], the decrease in primary deformation with position was fitted to an attenuation curve as y = ae b(x−0. 3) , where a is the amplitude, b is the damping coefficient, and the scales for x and y are millimeter and micrometer, respectively. Here, a = -1.86 µm, b = -1.05 mm −1 , and R 2 = 0.999. In Fig. 6(c), the dominant frequency was calculated to be 778 ± 1 Hz. Figure 6(d) and (e) show the decay fitting results for the damped oscillations. The decay amplitudes A 1 were from -0.08 µm to -0.19 µm (mean CV: 6.4%). There was no observable relation between the decay features (A 1 and B) and the primary displacement amplitudes A 0 . The decay coefficients B were from -31.2 s −1 to -63.3 s −1 and the damping ratios ε were from 0.006 to 0.013 (mean CV: 6.7% for all measurement positions). Figures 5 and 6 showed that the tissue phantom oscillated at the same dominant frequency regardless of the external applied forces (4 Pa to 32 Pa) or the distances between stimulation location and the points of measurements (0.3 mm to 5.3 mm). Because the measured damping ratio ε is small (less than 0.013), the damped natural frequency is nearly equal to the undamped natural frequency. Therefore, the SDOF model is an effective analytical model to determine the dominant natural frequency from the induced oscillation process.

Oscillation features for various sample concentrations and thicknesses
As demonstrated in Eq. (1), natural frequency is determined by spring stiffness and mass. We evaluated the oscillation features (f n , B, and ε) using SDOF method for agar phantoms with various concentrations (1%, 1.5%, and 2%) and different thicknesses (3 mm and 6 mm). The weights of the agar phantoms were 2.4-2.8 g for 3 mm thickness phantoms, and were 4.9-5.7 g for 6 mm thickness phantoms. The stimulation force was 20 Pa. The measurements were performed at five surface positions (0.3-1.3 mm from stimulation force, increment: 0.25 mm), and were repeated 6 times for each position. The data was selected based on the fitting R 2 for the decay coefficient B (≥ 0.9). Figure 7(a) shows the natural frequency values computed using the FFT method. For the 3 mm thickness phantoms, the measured natural frequencies and standard deviations (STDs) for 1%, 1.5% and 2% agar phantoms were 591 Hz ± 2 Hz (n = 14), 1166 Hz ± 6 Hz (n = 27), and 1328 Hz ± 8 Hz (n = 20). For the 6 mm thickness phantoms, the measurement natural frequencies and the STDs for 1%, 1.5% and 2% agar phantoms were 363 Hz ± 2 Hz (n = 19), 724 Hz ± 14 Hz (n = 9), and 907 Hz ± 4 Hz (n = 23). The natural frequencies quantified from FFT method (mean CV: 0.7%) were observed to be increased with phantom stiffness and decreased with phantom thickness. Figure 7(b) shows the decay coefficient computed using the curve fitting method and Eq. (5). For the 3 mm thickness phantoms, the calculated decay coefficients and the STDs for 1%, 1.5% and 2% agar phantoms were -62.0 s −1 ± 9.1 s −1 , -62.6 s −1 ± 7.3 s −1 , and -76.8 s −1 ± 16.7 s −1 . For the 6 mm thickness phantoms, the calculated decay coefficients and the STDs for 1%, 1.5% and 2% agar phantoms were -35.8 s −1 ± 3.7 s −1 , -56.4 s −1 ± 7.6 s −1 , and -68.1 s −1 ± 13.8 s −1 . The observed absolute values of the decay coefficient (mean CV: 15.4%) were also observed to increase with phantom concentrations and decrease with phantom thickness. However, the measurement errors were so large that the 1% and 1.5% agar phantoms with 3 mm thickness were not distinguishable. Figure 7(c) shows the calculated damping ratio ε based on Eq. (6). For the 3 mm thickness phantoms, the calculated damping ratios and the STDs for 1%, 1.5% and 2% agar phantoms were 0.017 ± 0.002, 0.009 ± 0.002, and 0.009 ± 0.002. For the 6 mm thickness phantoms, the calculated damping ratios and the STDs for 1%, 1.5% and 2% agar phantoms were 0.016 ± 0.002, 0.013 ± 0.002, and 0.011 ± 0.003. The calculated damping ratios decreased with phantom concentrations increased, and do not change obviously as the phantom thickness changes. Figures 7(d) and (e) summarize the relation among the natural frequencies, decay coefficients, and damping ratios of the agar phantoms. The decay coefficients correlate with the natural frequencies (y = -0.03x -34.7, R 2 = 0.32, p < 0.01), as predicted by Eq. (7). There was no obvious correlation between the measured natural frequencies and the damping ratios.

Relation between the natural frequency and the Young's modulus
Young's modulus is a measure of how easily a material deforms, and is defined as the ratio of the stress and strain [32]. Our previous studies have shown that Young's modulus can be estimated from the elastic surface wave propagation velocity [46,47,50]. Previous results have also shown good agreements between this OCE-based method and the gold standard-mechanical testing methods [46,47]. Here we compared the measurement of natural frequencies and Young's moduli using the SDOF method and the elastic wave propagation method. The measured agar phantoms had various concentrations from 0.75% to 2% (weight: 6.7g; thickness: from 7.1 mm to 8.5 mm). The stimulation pressure was 20 Pa.
Natural frequencies, decay coefficients, and damping ratios were measured using the SDOF method. The measurement positions were from 0.3 mm to 1.3 mm (step: 0.25 mm), relative to the stimulation points. Figures 8(a) and (b) show the representative normalized tissue damping oscillations for 0.75% -2% agar phantoms and the corresponding frequency components calculated using the FFT method. Figure 8(c) shows the natural frequency values, that ranged from 127 Hz to 774 Hz (mean CV: 0.9%) for the agar phantoms with a frequency resolution Fig. 8. Measurement of the natural frequencies and damping ratios using the SDOF method, and the Young's moduli from elastic wave propagation process on 0.75% to 2% agar phantoms (weight: 6.7g). The stimulation force was 20 Pa. (a) Representative normalized tissue damping oscillations for 0.75-2% agar phantoms. (b) Representative frequency components from the FFT analysis of tissue oscillations. The dominant frequencies (from 127 Hz to 774 Hz) were considered as the natural frequencies for the agar phantoms. Panels (c-d) show the quantifications (mean ± STD) of natural frequencies and damping ratios. Panel (e) shows the Young's moduli (mean ± STD). (f) The natural frequencies (mean ± frequency resolution) derived from the induced agar phantom tissue oscillations were linearly correlated to the square roots of the elastic modulus (mean ± standard deviation). The error bars in the x direction (mean STD: ± 2 Hz) were too small to be observed. of ± 2 Hz. Figure 8(d) shows the calculated damping ratios, that ranged from 0.022 to 0.008. The measurement CVs ranged from 53.9% (0.75% agar phantom) to 14.9% (2% agar phantom) with the mean CV of 27.9%.
In an isotropic homogeneous elastic material, Young's modulus E can be estimated from the speed c of elastic surface wave propagation, as follows [27]: where ρ is the material density and ν is the Poisson's ratio. The densities (ρ) were 820 kg/m 3 , 839 kg/m 3 , 871 kg/m 3 , 942 kg/m 3 , and 985 kg/m 3 , respectively, for the 0.75%, 1%, 1.25%, 1.5%, and 2% agar phantoms. The Poisson's ratio ν can be assumed as 0.5 [42,46]. The group velocity of the elastic surface wave is C g = d/t, where d and t are the distance and time delay of the primary surface displacement between measurement points [46]. The elastic waves were measured at 0.3 mm to 5.3 mm away from the stimulation point, in increment of 1 mm. The measurements were repeated at least five times at each location. The measured elastic wave propagation speeds in 0.75%, 1%, 1.25%, 1.5%, and 2% agar phantoms were 2.47 m/s, 2.78 m/s, 3.05 m/s, 3.93 m/s, and 6.02 m/s (mean CV: 2.9%), respectively; and the corresponding Young's moduli were 16.5 kPa, 21.4 kPa, 26.7 kPa, 48.0 kPa, and 117.8 kPa (mean CV: 5.8%), respectively, as shown in Fig. 8(e). The natural frequencies and the Young's moduli increased, while the damping ratios decreased as the agar concentration increased. The measurement repeatability for natural frequency (mean CV: 0.9%) and Young's modulus (mean CV: 5.8%) was high, but the measurement precision for the damping ratio was limited due to the low measurement precision of the decay coefficients (mean CV: 27.9%). Figure 8(f) summarizes the measurements of Young's moduli and natural frequencies of the agar phantoms. The natural frequencies closely correlated with the square root of the Young's moduli (y = 0.01x + 2.60, R 2 = 0.998, p < 0.01). The linear relationship between the natural frequency and the square root of the Young's modulus of a viscoelastic material was demonstrated previously in the experiments on silicone phantoms [35,38]. Our research was consistent with the results of previous studies, and confirmed that the natural frequency could be used to distinguish tissues and materials with different stiffness.

In vivo human cornea measurement
We performed a natural frequency measurement on the human cornea in vivo using OCE imaging and the SDOF method. During the corneal elastography imaging, the subject sat in a chair, placed his chin on a chin rest and forehead against a headband, and focused the eye on a fixation target. The stimulation pressure was 13 Pa, which was less than our previous work (20-60 Pa) [41]. The duration of the air force was 1 ms, and the time between two successive excitations was 100 ms. Corneal surface oscillations were recorded with an A-scan sampling rate of 20 kHz, and a total of 600 A-lines (30 ms) for each measurement position. The measurement points were scanned from 0.25 mm to 2.75 mm away from the stimulation point with the increment of 0.1 mm in the horizontal direction on corneal surface. Each set of measurements for the 2.5 mm distance took 2.6 s. Five measurements were repeated to assess measurement precision and repeatability. Figure 9(a) demonstrates selected corneal surface displacement profiles at 6 measurement locations. As wave-traveling distance increased, the primary displacement amplitude was decreased while the oscillation behavior (windowed area) remained similar. This was consistent with the phantom measurement results in Fig. 6. Corneal natural frequencies and damping ratios of each measurement position were acquired using the SDOF method for the windowed area. Figure 9(b) shows the fast Fourier Transform (FFT) analysis results in the frequency domain from 0-1200 Hz. Regardless of measurement positions, the corneal surface oscillated at similar damping frequencies (in the range of ∼ 0-600 Hz). Although we report only the damping ratios (ε) here, it is possible to derive the decay coefficients (B) using Eq. (6). The dominant damped frequencies (f d ) were converted to natural frequencies (f n ) based on Eq. (7). Figure 9(c) and Fig. 9(d) show the means and standard deviations (SDs) for the measured natural frequencies (f n ) and damping ratios (ε), respectively, for each measurement position (n = 5 repeat measurements).
The overall values for f n and ε were (mean ± SD) 259 ± 5 Hz and 0.084 ± 0.026, respectively.

Discussion
The improved detection sensitivity of common-path OCT and the perpendicular micro-scale air-pulse stimulation enable direct observation of small-magnitude damped oscillations [42] and a more precise natural frequency quantification that was possible to observe in vivo for the cornea. The low pressure (< 60 Pa), short duration (≤1 ms) perpendicular air-pulse stimulation can provide a range of ∼ 0-1.5 kHz excitation frequency, and induce oscillations that gradually decreases from sub-micrometer scale to sub-nanometer scale to zero. The common-path detection technique has enhanced optical phase stability (0.24 ± 0.07 nm) that make the small-magnitude oscillation behavior detectable [42]. In this study, we have characterized these small-magnitude oscillations, and computed the natural frequencies and damping ratios of agar phantoms and in vivo human cornea using a simple single degree of freedom (SDOF) method [36,45].
The dominant damped natural frequency f d and the decay coefficient B acquired from oscillation process are used in this SDOF method to calculate the natural frequency f n and the damping ratio ε (Eqs. (1)- (7)). We have demonstrated that the phantom oscillates at the same dominant frequency, regardless of the applied forces and stimulation distances (Figs. 5 and 6). Thereby, this SDOF method is effective to measure the dominant damped and undamped natural frequencies (f d and f n ) from the induced oscillation process. We have also demonstrated that the natural frequencies were larger in the agar phantoms with higher stiffness, and were lower in the agar phantoms with less weights/thicknesses (Fig. 7). This was consistent with the Eq. (1).
We measured the natural frequencies and Young's moduli for 0.75%, 1.0%, 1.25%, 1.5%, and 2.0% agar phantoms. The natural frequencies correlated to the square roots of the elastic modulus (Fig. 8). This result was consistent with previous publications on silicone phantoms [35,38]. We also confirmed that the natural frequency can be used to distinguish tissues and materials of different stiffness. The elastic wave propagation speeds in soft tissues are in the range of several meters per seconds. In this study, the propagation speeds were approximately 2.5-6.0 m/s for agar phantoms of concentrations between 0.75% and 2.0%. Most OCE systems are not fast enough to track the propagation of the elastic wave using one stimulation. The typical method to track the wave speed, also used in this study, is to use multiple stimulations and to measure data at different positions in the sample. Tissue motion, especially in the lateral direction, affects measurement repeatability and precision. Since natural frequency is not determined by the measurement position, it is not as sensitive to sample motions as the measurement of wave speed. Consequently, the measurement of natural frequency may be useful to determine the elastic properties of live tissues in vivo.
The damping ratios (ε) for all the measured agar phantoms (concentration: 0.75% to 2%) are very small, approximately from 0.01 to 0.04. Thereby, the damped natural frequency f d is close to the undamped natural frequency f n based on Eq. (7). We have shown that the damping ratio decreases as the agar concentration increases (Fig. 7(c) and Fig. 8(d)). In addition, the damping ratio is independent of phantom thickness (Fig. 7(c)), and there is no obvious correlation between the measured natural frequencies and the damping ratios (Fig. 7(e)).
In Fig. 8(d), the measurement CV of ε for 0.75% agar phantom is as large as 53.9%. The estimation of the damping ratio ε is limited by the measurement repeatability of the damping coefficient B (Eq. (6) and Eq. (7)). The absolute value of the coefficient B is larger in phantoms with higher agar concentration and smaller in thicker phantoms ( Fig. 7(b)), as well as, correlates with the natural frequencies (Fig. 7(c)). However, the measurement repeatability of B (CV: 15% ∼ 54%) are much worse than that of the natural frequency (CV: 0.3% ∼ 1.6%). The large measurement variation of B may limit its use for distinguishing tissues with different mechanical properties. The 1% and 1.5% agar phantoms with 3 mm thickness are not distinguishable based on the decay coefficients ( Fig. 7(b)) due to the large measurement error, but can be clearly distinguished based on the natural frequencies ( Fig. 7(a)). The large variation of the decay features in the homogeneous tissue phantoms might be caused by the limited fitting points for the decay envelop fitting (shown in Fig. 4(a)) as well as the superposition of waves propagating along the surface, inside the phantom, or reflected by phantom boundaries. This wave-superposition would be more complicated in biological tissues, where there are multiple layers or boundaries, or with illness that changes the local stiffness and induces boundaries. The analysis of the variation of the decay features or the wave-superposition patterns might be used to get information about the tissue boundaries. Future work will include development of new analytical models and use of finite element analysis method, to further explore the mechanisms of wave propagation and wave superposition in tissues with complex boundary conditions. Previous OCE methods have successfully shown the high-resolution resonant frequency measurement [35,36,[38][39][40] and have demonstrated that the natural frequency is linearly related to the square root of Young's modulus in a simple elastic model [35,38]. Our results are in agreement with these previous observations and have advantages over previous approaches. Here, we applied a non-contact air pulse stimulation method that provided a wide range of frequency excitation (e.g. 0 to kHz) simultaneously. We also used a high-sensitivity OCE system to detect the small-amplitude (sub-nanometer to sub-micrometer scale) damped oscillations with high resolution (0.24 nm). This air-pulse stimulation generated free oscillations after the primary deformation (Fig. 4). In this case, Eq. (3) described a SDOF method for the free oscillation response, which is slightly different from the forced oscillation response described in previous studies [36,37]. In Fig. 6, we have demonstrated that the free oscillation feature is independent of the measurement distance over a range of several millimeters. Therefore, there are potential benefits to use the tissue natural frequency to determine a global estimation of tissue properties. In addition, if the diseased tissue area is not easily accessible or it would be too invasive if measured directly, natural frequency measurements could provide an alternative way to assess tissue stiffness. In Fig. 5, we have demonstrated that the free oscillation feature is not dependent upon the stimulation force amplitude. Therefore, we can use very small stimulation forces to stimulate delicate tissues, such as the eye.
We have demonstrated in vivo measurements of human corneal natural frequency (f n ) and damping ratio (ε) in Fig. 9. In this pilot study, we applied a much lower air pressure (13 Pa) to stimulate, observe, and quantify the corneal oscillation features, compared to our previous work (20-60 Pa) [41]. We have shown good measurement precision and repeatability for the measurement of human corneal natural frequencies (mean ± SD: 259 ± 5 Hz, CV: 1.9%). The natural frequency value of the measured human cornea was very near to the 1.25% agar phantom, which was 261 ± 2 Hz (CV: 0.8%). The damping ratio (0.084 ± 0.026, CV: 31%) was much larger than that of 1.25% agar phantom (0.014 ± 0.005, CV: 34%).
It should be noted that the SDOF method used here is only an approximation and simplification method that describes the dominant oscillation frequencies and ignore other frequencies. As demonstrated in Fig. 4(b), the measured damped natural frequencies for a 2% agar phantom sample are in the range of 300-1500 Hz, but we only use the dominant damped natural frequency (776 Hz) in the SDOF model. This leads to some residual matching errors between the original damping oscillation data and the SDOF fitting data as shown in Fig. 4(c). The multi degree of freedom (MDOF) method can be used to describe a more complex motion system that the general vibration of the system consists of a sum of all the vibration modes and each vibration mode vibrates at its own frequency [45]. At this stage, we will focus on the use the dominant natural frequency in the SDOF model; in the near further we will also use the MDOF method for more detailed and more precise characterization of tissue biomechanics.
It should also be noted that the use of OCE in the measurement of natural frequency from the sub-micrometer to sub-nanometer tissue oscillations is only a preliminary study, and we still lack enough knowledge on tissue biomechanical property characterization using natural frequency values. First, it is not clear how the natural frequency spatially distributes in heterogeneous tissues of complex geometries with multiple tissue/liquid (or other) interfaces and thin-layers. Second, we do not know whether the natural frequency values can represent local variations of tissue biomechanical properties due to disease progression. Third, there is no direct model to determine Young's modulus from the observed natural frequency. In addition, we are still not sure whether this method is sensitive enough to distinguish changes in tissue biomechanical properties caused by tissue heterogeneities, as well as the resolution or repeatability of this approach in natural tissues. These are questions that should be investigated in future studies.

Conclusions
We have performed non-contact non-invasive natural frequency measurements on agar phantoms and in vivo human corneal imaging using an air-pulse based common-path OCE system and a single degree of freedom (SDOF) method. Small-amplitude (sub-nanometer to sub-micrometer scale) damped oscillations were induced by perpendicular air-pulse stimulation, and were directly observed using the common-path OCT with displacement resolution of 0.24. Tissue dominant natural frequency and damping ratio were obtained by a SDOF method in both the frequency and temporal domains using fast the Fourier transform (FFT) and curve-fitting methods. The tissue phantoms oscillated at the same dominant frequency, regardless of the applied forces and stimulation distances, showing that this oscillation frequency is the dominant natural frequency of the sample. By measuring the elastic properties of agar phantoms (0.75% to 2.0%), we also demonstrated that the dominant natural frequency correlated with the square root of Young's moduli and can be used to distinguish tissues and materials of different stiffness. Preliminary OCE imaging on the in vivo human cornea has shown good precision and repeatability for natural frequency measurement (259 Hz, CV: 1.9%) within a measurement distance of 2.5 mm from the point of stimulation on corneal surface.