Open Access
Add to BibSonomy
Abstract
Microliter air-pulse optical coherence elastography (OCE) has recently been proposed for the characterization of soft-tissue biomechanics using transient, sub-nanometer to micrometer-scale natural frequency oscillations. However, previous studies have not been able to provide real-time air-pulse monitoring during OCE natural frequency measurement, which could lead to inaccurate measurement results due to the unknown excitation spectrum. To address this issue, we introduce a dual-channel air-pulse OCE method, with one channel stimulating the sample and the other being simultaneously measured with a pressure sensor. This allows for more accurate natural frequency characterization using the frequency response function, as proven by a comprehensive comparison under different conditions with a diverse range of excitation spectra (from broad to narrow, clean to noisy) as well as a diverse set of sample response spectra. We also demonstrate the capability of the frequency-response analysis in distinguishing samples with different stiffness levels: the dominant natural frequencies increased with agar concentrations (181–359 Hz, concentrations: 1–2%, and maximum displacements: 0.12–0.47 µm) and intraocular pressures (IOPs) for the silicone cornea (333–412 Hz, IOP: 5–40 mmHg, and maximum displacements: 0.41–0.52 µm) under a 200 Pa stimulation pressure. These frequencies remained consistent across different air-pulse durations (3 ms to 35 ms). The dual-channel OCE approach that uses transient, low-pressure stimulation and high-resolution imaging holds the potential to advance our understanding of sample frequency responses, especially when investigating delicate tissues such as the human cornea in vivo.
© 2024 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
The biomechanical properties, such as stiffness, elasticity, and viscosity, of soft tissue are closely related to tissue structure, health, age, and function [1]. Changes in these properties at organ, tissue, or cellular levels are often associated with pathological changes, such as cancer, swelling, or inflammation [2–4]. Alterations in tissue biomechanical properties have emerged as a characteristic feature in diagnosing various diseases. For instance, liver fibrosis is closely linked to cirrhosis [5]; arterial wall stiffness is associated with systemic inflammation and hypertension-related atherosclerosis [6]; and breast fibrosis or tumor formations cause stiffness changes in breast tissue that can be differentiated from normal or benign tissue formations [7].
Optical coherence elastography (OCE) is a rapidly evolving elastic imaging technique used to assess the biomechanical properties of soft tissue [8]. OCE integrates a mechanical loading system and optical coherence tomography (OCT) imaging to induce tissue deformations, measure the deformations, and calculate the tissue’s elastic properties based on its response to stimulation. Advances in phase-sensitive OCT imaging [9–11] enable visualization and quantification of tissue dynamics at a sub-nanometer scale [12]. The shear-wave method, adapted from ultrasound and magnetic resonance imaging (MRI)-based elastography [13], is the most commonly used OCE technique to estimate a tissue’s Young's modulus from its transverse mechanical wave propagation speed [14]. However, it is important to note that this method may not be universally suitable for all OCE applications across all tissue types [15]. For instance, tissues with complex geometries and multiple thin layers, such as the cornea and skin, are known to produce highly dispersive Rayleigh-Lamb components, resulting in mechanical waves that travel along the surface and manifest as more complex than simple Rayleigh waves [16]. Moreover, the presence of dispersion in the elastic waves within complex tissues, including the cornea (e.g., the symmetric (S0) and anti-symmetric (A0) waves), further complicates the interpretation of their biomechanical properties [17–20]. As a result, applying shear-wave methods in tissues with complex geometries and thin layers could lead to inaccurate estimations of tissue properties [15,16].
An alternative and complementary OCE approach has been developed to quantify soft tissue natural frequencies based on the induced vibrational responses [21–31]. Natural frequency, defined as the frequency at which a system tends to oscillate when disturbed, has been shown to be linearly related to the square root of Young's modulus in a single degree of freedom (DOF) model [22,25,27]. There are two main approaches to measure natural frequency: the forced-response method and the free-response method. In the forced-response scenario, both the input (stimulation) and the output (response) of the system are measured in the frequency domain to form the frequency response function (FRF) for natural frequency characterization. In the free-response scenario, only the output (response) of the system is utilized to calculate the natural frequency, with the assumption that the input (stimulation) has no effect on the output. In our recent work, we used a microliter air-pulse OCE system and the free-response quantification method to generate, observe, and quantify the natural frequency oscillations of samples with resonant magnitudes from sub-nanometers to micrometers [27]. To closely replicate the free-response scenario, we selected a specific portion of the response data to analyze, excluding sample displacements directly driven by the air-pulse and focusing on the sample’s subsequent oscillation around its equilibrium position following the cessation of the air-pulse excitation [27,29,32]. Our results show significantly higher repeatability and reproducibility in measuring the natural frequencies of human corneas in vivo (average coefficient of variation [CV]: 3.2%, 234–277 Hz for 20 eyes [29]) compared to the method that uses shear-wave propagation speeds (average CV: 19.3%, 2.4–4.2 m/s for 18 eyes [33]). With its noninvasive, low-force, and transient stimulation approach, along with high-resolution imaging and reliable natural frequency measurement, this method demonstrates strong potential for further clinical translation.
However, the free-response method is based on an approximation that is not always valid. The input (stimulation) does have an effect on the output (response), especially when the input is not uniform across all frequencies. Previous studies could not monitor the air-pulse pressures in real time during natural frequency OCE measurement. This limitation prevents the effective characterization of the excitation spectrum—the energy distribution across different frequencies in the input signal. It also hinders the ability to separate signals from noise within an unknown excitation spectrum and potentially prevents a more reliable identification of the natural frequency components from the oscillation response. Moreover, it may even cause erroneous measurement results for the natural frequencies, especially in cases where the excitation spectrum is irregular and narrow-banded. To address this issue, we developed a dual-channel air-pulse OCE method to enable simultaneous observations of both the stimulation and response signals in the time and frequency domains. Our proposed method demonstrates the capability to conduct frequency-response analysis in a forced-response scenario, considering both the external force and the inherent dynamic properties of the sample to quantify its natural frequency. To create a comprehensive comparison between the free- and forced-response scenarios on natural frequency measurement, we intentionally stimulated the sample with a range of pulse duration parameters (3–35 ms; 200 Pa). This choice allowed us to create a diverse range of excitation spectra, spanning from broad to narrow, as well as a diverse set of sample response spectra, whether clean or noisy, across different air-pulse durations and measurement locations. Consequently, we can compare the natural frequency measurement performance aspects such as the measurement reliability for these two scenarios under different test conditions. Finally, we performed frequency-response analysis on a set of samples with varying stiffness levels in the forced-response scenario to evaluate the method’s capability to differentiate between samples based on their dominant natural frequency values.
2. Theory
2.1 Dynamic equations in a multi-degree-of-freedom system
In a multi-DOF system with DOFs, the relationship between the forces applied to the system and the resulting motion is described by the dynamics equation. In a physical coordinate system, this equation is commonly represented as follows:
A special solution for Eq. (1) is ${\boldsymbol x} = {\boldsymbol \varphi }{{\boldsymbol e}^{{\mathbf \Lambda }t}}$, when ${\boldsymbol f}(t )= 0.$ The ${i^{th}}$ order eigenvalue ${\mathrm{\Lambda }_i}$ and mass-normalized eigenvectors ${{\boldsymbol \varphi }_{\boldsymbol i}} = {[{{\varphi_{1i\; }},\; {\varphi_{2i}} \ldots \; {\varphi_{ni}}} ]^\textrm{T}}$ can be derived accordingly (i = 1, …, n). ${\mathrm{\Lambda }_i}$ is a complex parameter that encompasses the influences of both the ${i^{th}}$ order natural frequency and the damping parameter, and its magnitude is equal to the natural frequency. Each normalized eigenvector corresponds to a specific natural frequency and mode shape, which describes the pattern of motion exhibited by the system at that frequency.
By defining ${\boldsymbol x} = {\boldsymbol \varPhi y}$, where ${\boldsymbol \varPhi } = [{{{\boldsymbol \varphi }_1},{\boldsymbol \; } \ldots {\boldsymbol \; }{{\boldsymbol \varphi }_{\boldsymbol i}}, \ldots {{\boldsymbol \varphi }_{\boldsymbol n}}} ]$, and ${\boldsymbol y}$ is a vector of the displacements of each DOF in modal coordinates, Eq. (1) can be converted into modal coordinates as follows:
2.2 Free- and forced-response scenarios for natural frequency estimation
By using Eq. (3), natural frequency (${w_{0i}}$ or ${f_{0i}}$) can be estimated for both free-response and forced-response vibrations. Free-response vibrations are the natural oscillations of a system when it is disturbed from its equilibrium position and allowed to vibrate without external forces ($f(t )= 0$); forced-response vibrations occur when an external force is applied to the system ($f(t )\ne 0$).
In an under-damped, free-response scenario ($f(t )= 0, $ 0 ≤ ${\zeta _{0i}}$ < 1), a system’s oscillation amplitudes and frequencies are determined by its initial conditions, such as the initial displacements and velocities. Over time, the energy of the system dissipates due to damping, causing the oscillations to gradually decrease in amplitude until the system returns to its equilibrium position. The relationship between the vibration frequency (${\omega _i}$) and natural frequency in a free-response scenario is described as follows [27]:
In a forced-response scenario ($f(t )\ne 0$), both the external force and the system's inherent dynamic properties influence the system’s mechanical behavior. Resonance occurs when the frequency of the external force matches the natural frequency of the system. Forced-response vibrations can be analyzed using the frequency-response function (FRF), which is a transfer function defined as the ratio of response to stimulation in the frequency domain:
$FRF({j\omega } )$ is a complex function consisting of both an amplitude (the response-force ratio) and a phase (the in- or out-of-phase relationship between the response and input). In this study, to evaluate the peak resonant frequencies, the amplitude value $|{FRF({j\omega } )} |$ was used, signifying the presence of the sample's natural frequencies. The FRF amplitudes can be represented as follows:
When the ${i^{th}}$ order FRF ${H_i}$ achieves its maximum magnitude, the relationship between the vibration frequency and natural frequency in a forced-response situation is described as follows:
It was noted that a slight shift between the observed resonant frequencies and natural frequencies occurred in both the free-response (Eq. (4)) and forced-response (Eq. (7)) scenarios. Previous studies have demonstrated that for natural frequencies (${f_{0i}}$) within the 50–1000 Hz range and damping ratios (${\zeta _{oi}}$) below 0.3, disparities between the resonant frequency and natural frequency are minimal (e.g., < 3 Hz) and can be ignored [27]. Consequently, it is generally acceptable to use the resonant peak frequencies as a representation of the natural frequencies (i.e., ${\omega _i} \approx {\omega _{0i}}$) in both cases. In a single-DOF system when $n = 1$, Eqs. (1–7) can be significantly simplified, and the dominant natural frequency can be used to represent the dominant resonance behaviors.
3. Materials and methods
3.1 Dual-channel air-pulse OCE system setup
We developed a dual-channel air-pulse OCE system based on the previously reported OCE system [34], which enabled simultaneous measurement of the air-pulse profiles and the induced sample’s dynamic response at the same data acquisition frequency of 70 kHz. This system comprises an air-pulse sub-system and a synchronized phase-sensitive spectral domain OCT sub-system, as illustrated in Fig. 1(a). The air-pulse sub-system included an air tank with 99.99% nitrogen, a pulse signal generator (BNC575, Berkeley Nucleonics Corp., San Rafael, CA, USA), air tubing, a high-speed solenoid valve (V114-6M-M5, SMC Corp., Tokyo, Japan), and dual-channel microbore cannulas (inner diameter: 0.3 mm, outer diameter: 0.6 mm, and length: 50 mm). The air pressure was regulated using the air tank valve, and the pulse signal generator was used to control the duration of the air pulses in both channels by opening and closing the solenoid valve. In Channel 1, the air cannula was bent to facilitate perpendicular sample stimulation. In Channel 2, the air cannula remained straight, and the emitted air was measured with a pressure sensor (MPXV7002DP, NXP Semiconductors, Eindhoven, Netherlands). The air-pulse pressure was converted into an analog voltage signal. This signal was subsequently transformed into a digital voltage signal using an analog-to-digital converter and was recorded at a frequency of 70 kHz.
Fig. 1. Detection system and sample. (a) Schematic of the dual-channel air-pulse OCE system. SLD: superluminescent laser diode (waveband: 1290 ± 40 nm). (b) An artificial eye model with a silicone cornea and adjustable water pressure (equivalent to intraocular pressure, IOP).
The details of the OCT platform have been previously documented [35]. In brief, the system had a wavelength of 1290 ± 40 nm, an A-scan rate of 70 kHz, a maximum imaging depth of 6.94 mm, a maximum sensitivity of 99.3 dB, a −6 dB fall-off range of approximately 0–3 mm, a maximum sensitivity fall-off of −28.6 dB at full depth, and an axial resolution of approximately 15 µm (calibrated in air). A common-path configuration [12] was used to integrate both the sample and reference arms into the same light path. This was achieved using a 4 mm thick acrylic reference plate with a hole in the center to allow air stimulation from the air cannula in Channel 1 and to protect the sample from touching the air cannula. The maximum output power at the sample surface was 1.8 mW. Notably, the stimulation was performed at one point of the sample, while the OCT beam was scanned in the adjacent area to record the sample’s response (Fig. 1(a)). Since the mechanical response was measured slightly off the excitation point, it can avoid the near-field effect, where the induced sample’s vibration and mechanical waves are significantly affected by the excitation source geometry, reflections, and scattering within the sample at close distances from the stimulation point (i.e., in the near field) [27].
3.2 Phase-sensitive detection methods
In phase-sensitive OCT, the induced sample’s motion signal aspects, such as displacement, vibration, or strains, can be reflected by the phase variations. The phase detection sensitivity is determined by the signal-to-noise ratio (SNR), as ${\sigma _{\Delta \varphi }} \approx \frac{1}{{\sqrt {SNR} }},\; ({SNR \gg 1} )$ [36], as well as the phase stability. Typically, the interface between the air and the sample yields a bright intensity signal in the A-scan. Consequently, we calculate the phase variation in a time $(t )$ series from repeated A-scans at the same interface position. This information is then utilized to track the sub-pixel tissue surface displacements $\Delta y(t )$ in air as
where ${\lambda _0}$ is the center wavelength and $\Delta \varphi (t )$ represents the phase change after unwrapping. The common-path configuration enabled super stable phase signals (4.3 ± 1.4 milliradians over 60 ms, repeated 10 times, sensitivity amplitude: 92.3 dB) by reducing the vibrational turbulence between two arms, which corresponded to displacements of 0.44 ± 0.14 nm [32].3.3 Samples
To evaluate the natural frequency measurement performance between the free- and forced-response scenarios, we used two types of samples for the replication of tissues with varying degrees of stiffness or under different measurement conditions. The first type was agar phantoms with concentrations of 1%, 1.5%, 2%, and 2.5%. These phantoms were made by mixing Biowest agarose 111860 powder and water according to a procedure described previously [27]. The agar phantoms were placed in a Petri dish weighing 7.55 g, with outer and inner diameters of 35 mm and 32 mm, respectively. The pure agar phantoms weighed 10 g and had thicknesses of 12.11 mm (1%), 12.09 mm (1.5%), 12.06 mm (2%), and 12.03 mm (2.5%). Phantom densities were 1.027 kg/m3, 1.029 kg/m3, 1.031 kg/m3, and 1.034 kg/m3.
The second sample was an artificial eye model with a silicone cornea and intraocular pressure (IOP) that could be precisely controlled from 5 mmHg to 40 mmHg, as illustrated in Fig. 1(b). The silicone cornea was produced by mixing two types of self-defoaming silica gels (A and B silicones from Sanjing Xinde Technology Inc., Beijing, China). After pouring the mixture into a corneal mold, it was refrigerated for 1 day to solidify into the desired shape. The silicone cornea had a 12 mm central diameter, mimicking the human cornea, and a flat surrounding edge (diameter: 12–20 mm; thickness: 1.5 mm) that was attached to the rest of the eye model. The central portion had molded anterior and posterior radii of 7.79 mm and 6.38 mm, respectively, with a central corneal thickness of 0.55 mm. The remaining components of the artificial eye model, including a 2.3 mm thick cover (inner diameter: 12.4 mm; outer diameter: 36.3 mm) and a chamber (height: 20 mm; inner and outer diameters: 20.0 mm and 36.3 mm, respectively), which enclosed the flat edge of the silicone cornea, were 3D printed. The pressure within the eye model was controlled by injecting or removing water from the silicone cornea and chamber, and a pressure gauge connected to the chamber allowed the IOP to be monitored. The rise in the IOP from 5 mmHg to 40 mmHg coincided with changes in the silicone cornea's geometry. OCT measurements revealed a decrease in the anterior radius of curvature from 7.57 mm to 6.89 mm, as indicated by the linear fitting equation $y ={-} 0.020x + 7.678,\; ({{R^2} = 0.996} )$, accompanied by a 0.50 mm increase in the center height ($y = 0.014x-{-}0.07,\; {R^2} = 0.999$) [34].
3.4 Measurement protocols
We conducted two sets of measurements for natural frequency quantification. In the first set, our goal was to compare the performance between free-response and forced-response methods. To accomplish this, we intentionally stimulated the sample with a range of pulse parameters (pressure: 200 Pa, duration: 3–35 ms) to create a diverse range of excitation spectra, spanning from broad to narrow, as well as a diverse set of sample response spectra, from clean to noisy. The air pulse was delivered at the center of the sample, and radial-scan OCT measurements [34] were performed in the 45°–315° directions (18° increments). Each radial direction had 17 sampling points, covering a distance of 0.93–2.80 mm. The measurement employed an M-B mode scan protocol to measure the surface displacement as a function of time (Eq. (8)) and to spatially map the displacement profiles. The air pulse was stimulated every 0.5 s, and the OCE measurement was synchronized to the stimulation, capturing the temporal displacement profile of each point for 85.7 ms. Subsequently, we generated the frequency response spectra for both the free-response and forced-response scenarios and compared their natural frequency estimation performance aspects, such as measurement reliability, under different test conditions. In the second set of measurements, our goal was to perform a fast natural frequency measurement to differentiate agar phantoms with different concentrations (1%–2%) or the silicone cornea for different IOPs (5–40 mmHg). We used a narrow air pulse (duration: 3 ms, maximum pressure: 200 Pa) to stimulate the center of the sample and we measured the displacement (or oscillation frequencies) at a single point of the sample located 1 mm away from the stimulation point in the lateral direction.
4. Results
4.1 Dual-channel air-pulse calibration
We first measured the air-pulse pressures from both channels by placing two identical pressure sensors 1 mm from the ends of the air cannulas. The pulse durations ranged from 3 to 35 ms, and the tested pressures ranged from 60 Pa to 2 kPa. As shown in Fig. 2(a), the air pulses in Channel 1 had slightly larger magnitudes than those in Channel 2 (mean pressure: ∼500 Pa), whereas their profiles remained remarkably similar in the temporal domain. We calibrated the pressure differences that occurred during a relatively flat period of ∼30 ms for 35 ms duration air pulses. The absolute pressure differences ranged from 14 ± 7 Pa (stimulation: 62 Pa) to 187 ± 37 Pa (stimulation: 2 kPa), with most of the relative pressure differences within 10% (344–2004Pa), as shown in Fig. 2(b).
Fig. 2. Dual-channel air-pulse OCE calibration. (a) Air pressure profiles in both channels (duration: 3–35 ms; magnitude: ∼500 Pa). (b) Pressure difference for 35 ms pulses (∼60 Pa to 2 kPa). (c) Frequency analysis for (a), presented in log scale.
To characterize the frequency feature of the air pulse (i.e., excitation spectrum), we performed a fast Fourier transform (FFT) analysis on the pulses with a duration of 3–35 ms and a magnitude of ∼500 Pa, as shown in Fig. 2(c). Both channels displayed nearly identical frequency profiles, with shorter duration air pulses covering a wider excitation spectrum (e.g., a 3 ms duration air pulse yielded an excitation spectrum of ∼0–1 kHz). Notably, even longer duration pulses (10–35 ms) could potentially span a broader excitation spectrum, as evidenced by the oscillation of spectrum magnitudes in a higher frequency domain. The fluctuating and expanded excitation spectrum in the frequency domain was a result of signal jitters in the temporal domain. This phenomenon might have potentially generated both signal and noise in the sample’s response spectrum during the OCE stimulation and measurement procedure.
We also compared the OCE-measured maximum displacements (Channel 1) in response to different pressures (Channel 2: 103–2374 Pa) in the agar phantoms and the silicone cornea of the artificial eye model, as shown in Fig. 3. As expected, increasing the pressure led to a clear increase in the maximum displacements of all samples. In addition, the maximum displacements decreased with the increase in the concentration of the agar phantoms (1–2.5%) and with the elevation of IOPs (5–40 mmHg) in the silicone cornea under the same pressures. In the following sections, we report OCE natural frequency measurements acquired using an air-pulse pressure of 200 Pa, which generated sub-micrometer maximum deformations and even smaller oscillation magnitudes on these two samples.
Fig. 3. Maximum induced displacements by Channel 1 in 1–2.5% agar phantoms and a silicone cornea (IOP: 5–40 mmHg) in response to air-pulse pressures, measured at Channel 2.
4.2 Comparison between free- and forced-response natural frequency quantification
We utilized Fig. 4 to illustrate the natural frequency estimation processes using both the free-response and the forced-response methods. OCE measurement was conducted at a single point of the 1% agar phantom sample, positioned 1 mm laterally from the stimulation point. The applied air-pulse pressure was 200 Pa for a duration of 3 ms, resulting in a maximum displacement (primary deformation) of −0.62 µm and oscillation magnitudes of approximately ±0.18 µm.
Fig. 4. Demonstration of natural frequency estimation process in free-response ($f(t )= 0$)) and forced-response ($f(t )\ne 0$) scenarios. Sample: 1% agar phantom; air pressure: 200 Pa; and air-pulse duration: 3 ms. (a) Temporal profiles of the stimulation and response. (b) Fast Fourier transform (FFT) applied to the oscillation period of the response (sub-window in panel (a)) to estimate the natural frequencies in the free-response scenario. (c, d) Forced-response scenario, with (c) representing the FFT results for both the stimulation and the response, and (d) illustrating the natural frequency analysis using the frequency-response function (FRF), defined as the ratio of response to stimulation in the frequency domain.
The temporal profiles of the stimulation and response are displayed in Fig. 4(a). It was observed that the sample’s response included three major periods: a period before sample excitation ($f(t )= 0$), a primary deformation period directly driven by the air pressure ($f(t )\ne 0$), and a post-stimulation sample oscillation period ($f(t )= 0$). Since the key difference between free-response and forced-response scenarios lies in the presence of an external force, we could consider the entire displacement profile to represent forced-response scenarios. In contrast, during the oscillation period (shown in a sub-window of Fig. 4(a)), the stimulation force was zero, allowing the sample to oscillate around its equilibrium position from an initialized displacement. Therefore, this period could be approximated as a free-response scenario, and the response oscillation could be used to calculate the natural frequencies via FFT, as depicted in Fig. 4(b). Figure 4(c) shows the frequency spectra for both the stimulation and response. Notably, the excitation spectrum was smooth and had more energy in the low-frequency range than in the high-frequency range. We could account for the effect of the excitation spectrum, by considering both the stimulation and response in FRF, as illustrated in Fig. 4(d).
By comparing Figs. 4(b) and (d), it was evident that the free- and forced-response methods both effectively identified the dominant natural frequency (162 Hz) for the 1% agar phantom sample. However, it was more difficult to distinguish higher-order natural frequencies, particularly those surpassing 500 Hz (indicated by the red arrows), using the free-response method. The air-pulse stimulation resulted in an excitation spectrum with more pronounced energy in the low-frequency range than in the high-frequency range (Fig. 4(c)), which implied that the response also had to exhibit higher energy in the low-frequency range than in the high-frequency range. Even though we did not account for the force-driven primary deformation period and assumed that the oscillation period was unaffected by any external driving force, the response spectrum of the sample’s oscillation period proved to be susceptible to the characteristics of the excitation spectrum. Consequently, in the free-response case, we obtained a frequency spectrum that accentuated the low-frequency components and underestimated the high-frequency components (Fig. 4(b)). In contrast, the forced-response approach, utilizing frequency-response analysis, mitigated the influence of excitation and provided a more effective method for identifying the characteristics of higher-order natural frequencies (Fig. 4(d)).
To more comprehensively compare the performance between these two scenarios, we stimulated the sample with a range of pulse parameters (pressure: 200 Pa, duration: 3–35 ms) to create a diverse range of excitation spectra, spanning from broad to narrow, as well as a diverse set of sample response spectra, from clean to noisy. We mapped the spatial-temporal responses using radial-scan measurements (radial direction: 45°–315°, radial distance: 0.93–2.80 mm), utilized FFT and FRF analyses to capture the spatial-frequency characteristics for both the free-response and forced-response scenarios, and evaluated the natural frequency estimation performance in the spectrum of 0–700 Hz. Figure 5(a–d) illustrates several examples of spatial-frequency analysis with durations of 3 ms, 10 ms, 20 ms, and 35 ms. In each panel, the spatial-frequency maps resulting from free-response and forced-response analyses are represented by the left and middle cylinders, while the normalized frequency spectra (mean ± standard deviation [SD]) and excitation spectrum are presented in the right chart. The excitation spectrum is plotted on a logarithmic scale to emphasize its oscillation effect on the response spectrum. As shown in Fig. 5(a), the transient air pulse (i.e., 3 ms duration) produced a smoothly degenerated excitation spectrum, resulting in relatively high-quality spatial-frequency maps and similar estimated dominant natural frequencies in both scenarios. As the duration of the air pulse increased (Figs. 5(b–d)), the quality of the excitation spectra deteriorated, with noticeable oscillations of magnitudes in the excitation spectra. This oscillation became progressively more pronounced as the duration increased. Thus, the response spectra and the spatial-frequency maps deteriorated in both scenarios. In the free-response scenario, we observed turbulence noises, mainly in the low-frequency spectra (e.g., ∼0–400 Hz), caused by unmitigated fluctuations in the excitation spectra, as shown in the windows within the spectra charts in Figs. 5(c) and (d). In the forced-response scenario, we noticed frequency noises in the high-frequency spectra (∼500–700 Hz) during the long-duration air-pulse instances, particularly at 35 ms, as indicated by the arrows in Fig. 5(d). These noises originated from the formation of FRF (Eq. (5)), where the excitation levels in the high-frequency spectra were insufficient. As a result, the division of the response by the stimulation in FRF could lead to localized peaks of the noises in the high-frequency range.
Fig. 5. Spatial characterization of the natural frequency response. Sample: 2.5% agar phantom. Air-pulse was delivered at the sample center (duration: 3–35 ms; pressure: 200 Pa), and OCT measurements were performed in the radial directions (45°–315°; distance: 0.93–2.80 mm). (a–d) Comparison examples at different air-pulse durations of 3 ms, 10 ms, 20 ms, and 35 ms. Left and middle panels: Spatial-frequency maps for the resonant frequency analysis using free-response and forced-response analysis, respectively. Right charts: frequency-magnitude plots for the free-response and forced-response scenarios (linear scale) and the stimulation curves (logarithmic scale).
In the optimal case of 3-ms air-pulse stimulation, the dominant natural frequencies (mean ± SD) were 449 ± 1 Hz for the free-response scenario and 448 ± 1 Hz for the forced-response scenario. The slight difference was attributed to the assumption of ${\omega _i} \approx {\omega _{0i}}$, despite a different description of ${\omega _i}$ in both scenarios (see Eqs. (4) and (7)). The 3-ms transient air-pulse generated a clean excitation spectrum as well as a clean response spectrum of the sample, making it suitable as a reference for all spatial measurements spanning pulse durations from 3 ms to 35 ms. We performed a count of dominant frequencies within 10 Hz intervals, ranging from 203 Hz to 593 Hz, and presented the quantitative analysis results in Fig. 6. Figures 6(a) and (b) show the histograms for the dominant natural frequency counts for these two scenarios, and Fig. 6(c) illustrates the measurement accuracy assessment within the count range of 10 Hz (443–453 Hz) and 30 Hz (433–463 Hz). Upon comparison of Figs. 6(a) and (b), it was observed that in all cases, 92.4% of the dominant natural frequencies fell within the 433–463 Hz range in the forced-response scenario using FRF analysis, significantly higher than the 76.6% observed in the free-response scenario. The decreased accuracy in the free-response scenario, particularly in longer air-pulse duration cases, was primarily due to overestimation of the lower-frequency components, with an estimation of 20.0% of falsely predicted dominant natural frequencies falling within the 353–383 Hz range. When comparing the measurement accuracy for 5–35 ms air-pulse duration cases, the forced-response showed 17.5 ± 7.1% and 16.7 ± 3.4% enhancement over the free-response scenario in the 10 Hz and 30 Hz counting ranges, respectively. Therefore, by considering both the stimulation and response characteristics in the frequency domain, the use of FRF could define the excitation spectrum, predict the potential signal and noise for the frequency peaks, and improve the natural frequency measurement accuracy in air-pulse OCE measurement.
Fig. 6. Quantitative analysis of the dominant natural frequency measured spatially. Sample: 2.5% agar phantom. Air-pulse was delivered at the sample center (duration: 3–35 ms; pressure: 200 Pa), and OCT measurements were performed in each radial direction from 45° to 315° (increment: 18°). Each radial direction contained 17 sampling points to cover a distance from 0.93 to 2.80 mm. (a, b) Histograms for dominant natural frequency quantification in free-response and forced-response scenarios. (c) Measurement accuracy assessment.
4.3 Frequency-response analysis for agar phantoms and artificial cornea
To demonstrate the capability of the frequency-response analysis in distinguishing samples with different stiffness levels based on their dominant natural frequencies, we performed dual-channel air-pulse OCE measurements and frequency-response analysis on the agar phantoms and the artificial cornea. The air-pulse durations were 5–35 ms and the pressure was 200 Pa. The cannula was positioned 1 mm away from the sample surface, and the vibrational response from a single point of the sample was measured 1 mm away from the stimulation point in the lateral direction.
Figure 7 shows the frequency-response analysis of the agar phantoms with 1%, 1.5%, and 2.0% concentrations. The maximum displacements were 0.14–0.47 µm for 1%–2% agar phantoms. Figures 7(a, c, e) show the temporal profiles of stimulation and response, and Figs. 7(b, d, f) display the maps of the FRF magnitudes within the frequency range of 0–700 Hz. As shown in Fig. 7(g), the dominant natural frequencies of the 1%, 1.5%, and 2.0% agar phantoms were 161 ± 2 Hz, 277 ± 2 Hz, and 359 ± 4 Hz, respectively, based on the mean ± SD of 10 repeated measurements. These findings indicate an increase in sample stiffness as the concentration of agar increases. Notably, the measured dominant natural frequencies exhibited consistency for the same agar phantom across different stimulus durations. However, the CV, a metric used to measure dispersion in measurements (CV = SD/mean), increased from 0% to 1.3% as the duration progressed from 3 ms to 35 ms (Fig. 7(h)).
Fig. 7. Frequency-response analysis for agar phantoms. Air-pulse duration: 3–35 ms. (a, c, e) Stimulation and response temporal profiles for 1, 1.5, and 2.0% agar phantoms, respectively. (b, d, f) Magnitude maps of the FRF with the dominant natural frequencies indicated by arrows (10 repetitions). These frequencies are quantified in (g) using bar plots (mean ± SD) and in (h) using the coefficients of variation (CVs).
Figure 8 displays the frequency-response analysis of the silicone cornea in an artificial eye model, covering an IOP range from 5 mmHg to 40 mmHg and air-pulse durations from 3 to 35 ms. The maximum displacements were 0.41–0.52 µm under a 200 Pa pressure. Figures 8(a, c, e, g) present the temporal profiles for the silicone cornea under different IOP conditions (5, 15, 30, and 40 mmHg). Correspondingly, Figs. 8(b, d, f, h) show the maps of FRF magnitudes, with each measurement repeated 10 times. Notably, the silicone cornea's dominant natural frequencies were found to be correlated with the IOP levels, independent of the stimulus durations. In Fig. 8(i), the dominant natural frequencies were 333 ± 7 Hz (5 mmHg), 352 ± 9 Hz (10 mmHg), 363 ± 7 Hz (15 mmHg), 375 ± 9 Hz (20 mmHg), 385 ± 7 Hz (25 mmHg), 395 ± 8 Hz (30 mmHg), 401 ± 9 Hz (35 mmHg), and 412 ± 10 Hz (40 mmHg). The positive correlation indicates increases in the elasticity and stiffness of the silicone cornea with rising IOP. Figure 8(j) demonstrates the measurement CV for the dominant natural frequencies, which has a general increasing trend with longer duration times, consistent with the observations in Fig. 7(h).
Fig. 8. Frequency-response analysis for the silicone cornea in an artificial eye model (IOP: 5–40 mmHg; air-pulse duration: 3–35 ms). (a, c, e, g) Temporal profiles of the stimulation and response at IOPs of 5 mmHg, 15 mmHg, 30 mmHg, and 40 mmHg, respectively. (b, d, f, h) Corresponding magnitude maps of the FRF with the dominant natural frequencies indicated by arrows (10 repeated measurements). (i) Bar plot (mean ± SD) of the dominant natural frequencies corresponding to IOP (5–40 mmHg, increment: 5 mmHg). (j) Coefficient of variation (CV) corresponding to air-pulse durations.
5. Discussion
In this study, we developed a dual-channel air-pulse OCE system that enables simultaneous recording of microliter air-pulse signals and micro- to sub-nanometer-scale tissue displacement profiles in both the temporal and frequency domains. With this system, frequency-response analysis was conducted by measuring both the stimulation and response signals, allowing us to investigate the forced-response method alongside the previously reported free-response method for natural frequency quantification. As expected, the use of a transient air-pulse stimulation resulted in a wider excitation spectrum (e.g., a 3 ms pulse produced ∼0–1000 Hz frequency spectrum). Interestingly, even longer duration air pulses (10–35 ms) had the potential to cover a broader excitation spectrum. This was evident from the fluctuations in magnitudes in the higher frequency range, primarily caused by temporal domain jitter signals (Fig. 2(c)). Although a pressure amplitude difference of ∼10% was observed between the two air-pulse channels (as shown in Figs. 2(a, b)), this discrepancy did not lead to measurement errors when employing the frequency response analysis for natural frequency quantification. Because the variance in magnitude did not influence the normalized frequency profiles, as illustrated in Fig. 2(c), we only used the normalized spectrum for both the stimulation and response to calculate the FRF.
Figures 4–6 confirm that the forced-response method was superior to the free-response method by accounting for the effect of the input (stimulation) on the output (response) and providing a more accurate and reliable frequency spectrum that reflects the natural frequency features. We compared transient (e.g., 3–5 ms durations) and long-duration (10–35 ms) air-pulse cases. In the transient air-pulse case (e.g., 3–5 ms durations), the excitation spectrum was smooth and had more energy in the low-frequency range than in the high-frequency range. This meant that the output (response) also had to have more energy in the low-frequency range than in the high-frequency range. As a result, directly applying FFT for the oscillation period in the free-response scenario leads to a frequency spectrum that overestimates the low-frequency components and underestimates the high-frequency components. Consequently, it becomes more difficult to distinguish the higher-order natural frequencies for the 1% agar phantom samples in the free-response scenario than in the forced-response scenario (Figs. 4(b) and (d)). In the case of long-duration air-pulse stimulation (e.g., 10–35 ms), the excitation spectrum displayed non-uniformity, narrowness, noise, and variations in energy across the frequency spectrum. This contributed to noise and errors in the response spectrum in both the forced-response and free-response scenarios, with the free-response scenario being more severely affected. The variations in the excitation spectrum could cause false predictions of the dominant natural frequency values in the free-response scenario, as shown by the quantitative results in Fig. 6. Since the dual-channel OCE system is capable of providing air-pulse profiles in the time-frequency domains, we can further reduce the noise and errors by selecting a smoother portion of the excitation spectrum for the FRF calculation. Instead of using the entire spectrum for all cases, this approach allows us to achieve a more precise measurement of natural frequency using the FRF method. Consequently, employing the dual-channel OCE system could also help establish an effective excitation spectrum for differentiating potential signals and noise in the response spectrum, thus resulting in a more accurate measurement of natural frequency through frequency-response analysis.
The use of low-pressure, transient, and measurable air-pulse profiles in the dual-channel OCE approach made this approach particularly suitable for characterizing delicate tissues such as the human eye in vivo. While corneal natural frequency vibrations are also observable using clinical instruments, such as CorVis ST (Oculus Optikgeräte GmbH, Wetzlar, Germany) [37], the duration of the clinical air-puff is typically 10–35 ms, potentially leading to a noisy excitation spectrum for natural frequency analysis, similar to what is depicted in Figs. 5 and 6. Figure 8 shows a positive correlation between the natural frequency and elevated IOP (333–412 Hz, IOP: 5–40 mmHg), which served as the basis for potential use in detecting ocular tissue changes caused by glaucoma. However, it should be noted that the ocular natural frequency is not only determined by tissue biomechanics (such as Young’s modulus), but also influenced by various other factors, including the mass, shape, and boundary conditions [29]. Unfortunately, we currently lack an accurate analytical model to analyze the corneal biomechanical property from the natural frequency value. Further in-depth investigations into corneal natural frequencies and their correlation with other ocular structural parameters (e.g., corneal center thickness, IOP, and corneal topography) are necessary to gain a comprehensive understanding of the limitations and potential clinical applications of these natural frequency measurements. Figure 4(d) displays the multiple frequency-response peaks observed, allowing the use of this novel OCE approach to facilitate multi-DOF analysis. In addition, by spatially quantifying the magnitudes and phases of the FRFs, this approach could be further used to identify all the modal parameters of the sample, including the natural frequencies, damping ratios, and mode shapes. Spatial variations in the modal parameters of the human cornea reflect changes in oscillation phases and amplitudes at different measurement positions, which may arise due to subtle changes in local stiffness or corneal geometry. Therefore, the detection of subtle changes in modal parameters might be a useful alternative to the mechanical wave OCE approach for diagnosing subtle morphological changes in the cornea caused by corneal disease (e.g., keratoconus or ectasia), wound healing, or treatments (e.g., corneal cross-linking or refractive surgery) [38–40]. We plan to develop novel analytical methods and finite element models incorporating shape, multi-layered structure, IOP, boundary conditions, tissue biomechanics, and morphological or surgical changes to investigate relationships between modal parameters (natural frequency in particular) and these factors. Our proposed dual-channel air-pulse OCE approach has the potential to provide additional information on corneal properties, thus complementing current clinical assessments. Furthermore, it could be a useful tool for clinical detection of ocular disease and evaluation of medical or surgical treatment outcomes.
Funding
National Natural Science Foundation of China (61975030); Guangdong Basic and Applied Basic Research Foundation (2021A1515011981); Guangdong Provincial Pearl River Talents Program (2019ZT08Y105); Guangdong Eye Intelligent Medical Imaging Equipment Engineering Technology Research Center (2022E076); Guangdong-Hong Kong-Macao Intelligent Micro-Nano Optoelectronic Technology Joint Laboratory (2020B1212030010); National Eye Institute (P30EY07551, R01EY022362, R01EY022362).
Disclosures
GL, JX and YH: Weiren Meditech Co., Ltd. (C); JQ. and LA: Weiren Meditech Co., Ltd. (E).
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. P. H. Wu, D. R.-B. Aroush, A. Asnacios, et al., “A comparison of methods to assess cell mechanical properties,” Nat. Methods 15(7), 491–498 (2018). [CrossRef]
2. T. R. Cox, “The matrix in cancer,” Nat. Rev. Cancer 21(4), 217–238 (2021). [CrossRef]
3. C. Guillot and T. Lecuit, “Mechanics of epithelial tissue homeostasis and morphogenesis,” Science 340(6137), 1185–1189 (2013). [CrossRef]
4. D. T. Butcher, T. Alliston, V. M. Weaver, et al., “A tense situation: forcing tumour progression,” Nat. Rev. Cancer 9(2), 108–122 (2009). [CrossRef]
5. V. De Lédinghen, J. Vergniol, C. Barthe, et al., “Non-invasive tests for fibrosis and liver stiffness predict 5-year survival of patients chronically infected with hepatitis B virus,” Aliment. Pharmacol. Ther. 37(10), 979–988 (2013). [CrossRef]
6. A. Mahmud and J. Feely, “Arterial stiffness is related to systemic inflammation in essential hypertension,” Hypertension 46(5), 1118–1122 (2005). [CrossRef]
7. A. Samani, J. Zubovits, D. Plewes, et al., “Elastic moduli of normal and pathological human breast tissues: an inversion-technique-based investigation of 169 samples,” Phys. Med. Biol. 52(6), 1565–1576 (2007). [CrossRef]
8. J. Schmitt, “OCT elastography: imaging microscopic deformation and strain of tissue,” Opt. Express 3(6), 199–211 (1998). [CrossRef]
9. Y. Zhao, Z. Chen, C. Saxer, et al., “Phase-resolved optical coherence tomography and optical Doppler tomography for imaging blood flow in human skin with fast scanning speed and high velocity sensitivity,” Opt. Lett. 25(2), 114–116 (2000). [CrossRef]
10. M. Sticker, C. K. Hitzenberger, R. Leitgeb, et al., “Quantitative differential phase measurement and imaging in transparent and turbid media by optical coherence tomography,” Opt. Lett. 26(8), 518–520 (2001). [CrossRef]
11. S. J. Kirkpatrick, R. K. Wang, D. D. Duncan, et al., “OCT-based elastography for large and small deformations,” Opt. Express 14(24), 11585–11597 (2006). [CrossRef]
12. G. Lan, M. Singh, K. V. Larin, et al., “Common-path phase-sensitive optical coherence tomography provides enhanced phase stability and detection sensitivity for dynamic elastography,” Biomed. Opt. Express 8(11), 5253–5266 (2017). [CrossRef]
13. E. L. Carstensen and K. J. Parker, “Physical models of tissue in shear fields,” Ultrasound Med. Biol. 40(4), 655–674 (2014). [CrossRef]
14. F. Zvietcovich and K. V. Larin, “Wave-based optical coherence elastography: the 10-year perspective,” Prog. Biomed. Eng. 4(1), 012007 (2022). [CrossRef]
15. I. Pelivanov, L. Gao, J. Pitre, et al., “Does group velocity always reflect elastic modulus in shear wave elastography?” J. Biomed. Opt. 24(07), 1 (2019). [CrossRef]
16. J. J. Pitre Jr, M. A. Kirby, L. Gao, et al., “Super-shear evanescent waves for non-contact elastography of soft tissues,” Appl. Phys. Lett. 115(8), 083701 (2019). [CrossRef]
17. A. Ramier, A. M. Eltony, Y. Chen, et al., “In vivo measurement of shear modulus of the human cornea using optical coherence elastography,” Sci. Rep. 10(1), 17366 (2020). [CrossRef]
18. G. Li, X. Feng, S. Yun, et al., “Simultaneous tensile and shear measurement of the human cornea in vivo using S0- and A0-wave optical coherence elastography,” Acta Biomater. 175, 114–122 (2024). [CrossRef]
19. X. Feng, G.-Y. Li, S.-H. Yun, et al., “Ultra-wideband optical coherence elastography from acoustic to ultrasonic frequencies,” Nat. Commun. 14(1), 4949 (2023). [CrossRef]
20. M. G. Sun, T. Son, J. Crutison, et al., “Optical coherence elastography for assessing the influence of intraocular pressure on elastic wave dispersion in the cornea,” J. Mech. Behav. Biomed. 128, 105100 (2022). [CrossRef]
21. X. Liang, A. L. Oldenburg, V. Crecea, et al., “Optical micro-scale mapping of dynamic biomechanical tissue properties,” Opt. Express 16(15), 11052–11065 (2008). [CrossRef]
22. V. Crecea, . L. Oldenburg, X. Liang, et al., “Magnetomotive nanoparticle transducers for optical rheology of viscoelastic materials,” Opt. Express 17(25), 23114–23122 (2009). [CrossRef]
23. S. G. Adie, X. Liang, B. F. Kennedy, et al., “Spectroscopic optical coherence elastography,” Opt. Express 18(25), 25519–25534 (2010). [CrossRef]
24. T. H. Fischer, “Imaging and elastometry of blood clots using magnetomotive optical coherence tomography and labeled platelets,” IEEE J. Sel. Top. Quantum Electron. 18(3), 1100–1109 (2012). [CrossRef]
25. W. Qi, R. Li, T. Ma, et al., “Resonant acoustic radiation force optical coherence elastography,” Appl. Phys. Lett. 103(10), 75–96 (2013). [CrossRef]
26. B. I. Akca, E. W. Chang, S. Kling, et al., “Observation of sound-induced corneal vibrational modes by optical coherence tomography,” Biomed. Opt. Express 6(9), 3313–3319 (2015). [CrossRef]
27. G. Lan, K. V. Larin, S. Aglyamov, et al., “Characterization of natural frequencies from nanoscale tissue oscillations using dynamic optical coherence elastography,” Biomed. Opt. Express 11(6), 3301–3318 (2020). [CrossRef]
28. F. H. Silver, N. Kelkar, T. Deshmukh, et al., “Mechano-vibrational spectroscopy of tissues and materials using vibrational optical coherence tomography: a new non-invasive and non-destructive technique,” Recent Prog. Mater. 2, 1–16 (2020). [CrossRef]
29. G. Lan, S. Aglyamov, K. V. Larin, et al., “In vivo human corneal natural frequency quantification using dynamic optical coherence elastography: repeatability and reproducibility,” J. Biomech. 121, 110427 (2021). [CrossRef]
30. M. A. Crespo, H. J. Jimenez, T. Deshmukh, et al., “In vivo determination of the human corneal elastic modulus using vibrational optical coherence tomography,” Transl. Vis. Sci. Technol. 11(7), 11 (2022). [CrossRef]
31. R. McAuley, A. Nolan, A. Curatolo, et al., “Co-axial acoustic-based optical coherence vibrometry probe for the quantification of resonance frequency modes in ocular tissue,” Sci. Rep. 12(1), 18834 (2022). [CrossRef]
32. G. Lan, Q. Shi, Y. Wang, et al., “Spatial assessment of heterogeneous tissue natural frequency using micro-force optical coherence elastography,” Front. Bioeng. Biotechnol. 382, 1 (2022). [CrossRef]
33. G. Lan, S. R. Aglyamov, K. V. Larin, et al., “In vivo human corneal shear-wave optical coherence elastography,” Optom. Vis. Sci. 98(1), 58–63 (2021). [CrossRef]
34. G. Ma, J. Cai, R. Zhong, et al., “Corneal surface wave propagation associated with intraocular pressures: OCT elastography assessment in a simplified eye model,” Bioengineering 10(7), 754 (2023). [CrossRef]
35. G. Lan, J. Xu, Z. Xu, et al., “Design of 1300-nm spectral domain optical coherence tomography angiography system for iris microvascular imaging,” J. Phys. D: Appl. Phys. 54(26), 264002 (2021). [CrossRef]
36. B. Hyle Park, M. C. Pierce, B. Cense, et al., “Real-time fiber-based multi-functional spectral-domain optical coherence tomography at 1.3 µm,” Opt. Express 13(11), 3931–3944 (2005). [CrossRef]
37. Z. Han, C. Tao, D. Zhou, et al., “Air puff induced corneal vibrations: theoretical simulations and clinical observations,” J. Refract. Surg. 30(3), 208–213 (2014). [CrossRef]
38. G. Regnault, M. A. Kirby, R. K. Wang, et al., “Possible depth-resolved reconstruction of shear moduli in the cornea following collagen crosslinking (CXL) with optical coherence tomography and elastography,” Biomed. Opt. Express 14(9), 5005–5021 (2023). [CrossRef]
39. X. Han, Y. Zhang, G. Shi, et al., “Quantitative assessment of corneal elasticity distribution after FS-LASIK using optical coherence elastography,” J. Biophotonics 17, e202300441 (2024). [CrossRef]
40. V. S. De Stefano, M. R. Ford, I. Seven, et al., “Depth-dependent corneal biomechanical properties in normal and keratoconic subjects by optical coherence elastography,” Transl. Vis. Sci. Technol. 9(7), 4 (2020). [CrossRef]
