Open Access
Add to BibSonomy
Abstract
Skin-elasticity measurements can assist in the clinical diagnosis of skin diseases, which has important clinical significance. Accurately determining the depth-resolved elasticity of superficial biological tissue is an important research direction. This paper presents an optical coherence elastography technique that combines surface acoustic waves and shear waves to obtain the elasticity of multilayer tissue. First, the phase velocity of the high-frequency surface acoustic wave is calculated at the surface of the sample to obtain the Young's modulus of the top layer. Then, the shear wave velocities in the other layers are calculated to obtain their respective Young's moduli. In the bilayer phantom experiment, the maximum error in the elastic estimation of each layer was 2.2%. The results show that the proposed method can accurately evaluate the depth-resolved elasticity of layered tissue-mimicking phantoms, which can potentially expand the clinical applications of elastic wave elastography.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
The skin has a layered structure, and the elasticity of the epidermis is different from that of the dermis [1]. Variations in skin elasticity have been linked to pathological changes in the skin in some circumstances [2]. The elasticity of cutaneous lesions, such as systemic sclerosis [3] and melanoma [4], is different from that of healthy skin; therefore, the measurement of skin elasticity is useful for assessing skin health. Furthermore, the depth of a cutaneous lesion can indicate its prognosis [5]. For example, melanoma only affects the epidermis or more superficial layers of the skin at an early stage; thus, early diagnosis can be achieved by detecting elastic changes in shallow tissue, thereby facilitating an accurate prognosis [5,6]. However, if the diagnosis is made at a more advanced stage, the tumor may affect the dermis or deeper layers, and the 5-year survival rate will decrease to 10%–15% [5]. Thus, early diagnosis of skin diseases by accurately evaluating the depth-resolved elasticity of skin tissue to locate the depth of lesions is important.
Skin elasticity is typically qualitatively judged by dermatologists through palpation; however, this method is highly subjective. Currently, ultrasound elastography [7] and magnetic resonance elastography (MRE) [8] have been applied in clinical diagnosis owing to their advantageous non-invasiveness and high penetration depth. However, MRE has low spatial resolution [9]and requires long clinical waiting times. To improve image resolution, high-frequency ultrasound ($\ge $ 18 MHz) is applied in ultrasound elastography for assessing the stiffness distributions in small, superficially regions [10,11]. In contrast, optical coherence elastography (OCE) provides the advantages of micro-level spatial resolution of optical coherence tomography (OCT) [12–15] and a higher sensitivity to the mechanical deformation of tissue compared with ultrasound elastography [16]. Therefore, it enables the non-invasive detection of small elastic variations in the tissue at the early stage of the tissue lesions.
Compared with static and quasi-static OCE, the elastic wave detection methods used in dynamic OCE are easier to implement for quantifying elastic information. Surface acoustic waves (SAWs) are used to study the elasticity of biological tissue in superficial areas, such as the skin [1,17] and cornea [18]. Generating SAWs is easy when excitation is focused on the free sample surface. Additionally, most of the energy of SAWs remains on the free sample surface, which is easier to detect [19]. Recently, many studies have used the dispersion curve of SAWs to assess the depth-resolved elasticity of superficial biological tissues [1,20–22]. The forward model based on wave equations has high computing costs and a long inversion time. Li et al. [20]. proposed a technique based on the wavelength-depth relationship of SAWs to directly estimate the elastic properties of soft-tissue-mimicking materials. The Young's moduli of the samples at different depths were calculated using the velocities of different frequencies to distinguish the different layers in heterogeneous materials. However, the wavelength-depth inversion model did not demodulate the combined phase velocities at different depths; thus, the Young's modulus of the lower-layer material is biased. Zhou et al. [1]. proposed a weighted average phase-velocity inversion model. The model matched the measured phase-velocity dispersion curve to that generated from numerical simulations to assess the depth-resolved elasticity. However, the implementation of this model is limited because it requires accurate prior knowledge of the layer thickness to calculate the model inversion error, which varies with different top-layer thicknesses and top-bottom shear wave velocity ratios. Chen et al. [22]. proposed a total variance-regularization iterative algorithm for obtaining the elastic modulus distribution along the skin depth from the dispersion curve of SAWs. This method also requires accurate information on the layer thickness, and it is only applicable to layered structures with a positive dispersion. However, healthy skin has a layered structure with a negative dispersion, and the epidermis is harder than the dermis.
This study proposes a tractable method to detect the depth-resolved elasticity of layered tissue. According to the principle of SAW penetration depth, the elasticity of the different layers of the samples can be obtained by acquiring the SAW and shear waves that are generated simultaneously after the samples are excited. In this study, a focused ultrasound transducer is used as the excitation source and phase-sensitive OCT (PhS-OCT) using an M-B scan mode is used as the detector. The velocity of high-frequency SAW on the sample surface is obtained for calculating the Young's modulus of the top layer, and the shear wave velocities are obtained in the other layers to calculate the Young's modulus of each layer.
2. Theory
When a sample (tissue or tissue-mimicking phantom) is stimulated, the free surface of the sample produces SAWs. Consider a dual-layer model as an example. Assuming that each layer of the sample is isotropic and homogeneous, the longitudinal and shear wave equations can be obtained by substituting the Stokes–Helmholtz decomposition of the displacement field into Navier’s equilibrium equation with zero equivalent force [23]. This can be written according to the compression wave velocity ${C_p}$ and shear wave velocity ${C_s}$ [23] as follows:
Considering a plane wave propagation in the x-z plane, the displacement is independent of the y-axis; thus, the curl operator in the Cartesian coordinate system can be simplified to obtain the potential energy functions $\phi $ and $\psi $ [23] as follows:
The sample material parameters were substituted to solve the wave velocity equation. The parameters of common imitation materials were selected as follows: the upper-sample elasticity was 60 kPa, the upper-sample density was 1049 kg/m3, the upper-sample thickness h was 1.5 mm, the lower-sample elasticity was 20 kPa, the lower-sample density was 1000 kg/m3, and Poisson's ratio was 0.49. The obtained dispersion curve of the SAW is shown in Fig. 2. It shows that the different frequency components of the SAW have different phase velocities. The high-frequency phase velocity tends to stabilize and converge to the Rayleigh-wave velocity of the top layer. Similar SAW dispersion curves are obtained when the elasticity of the upper layer is greater than that of the lower layer, which can be explained by the relationship between the penetration depth and wavelength of the SAW in homogeneous samples.
Fig. 2. Wave dispersion curve of the bilayer sample surface. ${C_{R1}}$ is the Rayleigh wave velocity calculated by substituting the Young's modulus of layer 1 into Eq. (8).
On the surface of a planar, stress-free, and isotropic medium, the penetration depth l of the SAW is approximately equal to the wavelength $\lambda $ of the SAW, and is expressed as follows [24]:
where f and ${C_R}$ are the frequency and velocity of the SAW, respectively. The higher frequency component of a SAW that has a shorter wavelength penetrates the shallower region of the medium. In contrast, the lower frequency component of a SAW that has a longer wavelength penetrates the deeper region below the surface [25]. Therefore, the high-frequency phase velocity that tends to stabilize is only related to the elasticity of the top layer, and the low-frequency phase velocity can be assumed to be a weighted superposition of the phase velocity from the surface to its maximum penetration depth [21].When an ultrasonic transducer is used as an excitation source, the elastic wave generated in the sample is a combination of the SAW at the surface and shear waves below the surface. The phase-velocity dispersion curve is also obtained at a measured depth, ${z_0}$, within a multi-layer sample, other than at the top layer. Because the penetration depth, l, of the SAW at a frequency f ($f$>${f_3}$, ${f_3} \approx {{{C_R}} / {{z_0}}}$) is smaller than the measured depth ${z_0}$, the phase velocity at f obtained at ${z_0}$ is not affected by the SAW, which is the shear wave velocity ${C_s}$.
Thus, in this study, the phase velocity at a high frequency, which tends to stabilize in the phase-velocity curve of the sample surface, represents the Rayleigh-wave velocity, ${C_R}$, of the top layer when the homogeneous sample surface is in air. The Young's modulus E of the top layer is calculated using the relationship between the Rayleigh-wave velocity ${C_R}$ and the Young's modulus as follows:
Moreover, the high-frequency phase velocity, which tends to stabilize in the phase-velocity curves of the remaining layers, represents the shear wave velocity ${C_\textrm{s}}$ of each layer. The Young's modulus E of each of the remaining layers is calculated using the relationship between the shear wave velocity ${C_\textrm{s}}$ and Young's modulus E, as follows:
To further verify the above method, the commercial software COMSOL Multiphysics (COMSOL AB, Stockholm, Sweden) was used for a finite element simulation. A two-dimensional axisymmetric model of the bilayer tissue was created to simplify the calculation because the elastic wave generated by ultrasonic excitation propagates cylindrically. The excitation actuator was modeled as a 5 MHz focused ultrasound transducer, and the pulse excitation time was 200 $\mathrm{\mu s}$. The ultrasonic focus area is determined by the full width at half-maximum of the sound pressure field, which was approximately $0.58 \times 3.9$ mm (transverse diameter ${\times} $ axial width). The acoustic average grid quality was 0.95. A pure elastic model with dimensions size $50 \times 50$ mm was created. Table 1 lists specific parameters. The mechanical average mesh quality was 0.96, and the bottom surface was fixed. Because half of the cross-section was modeled, only the area with a width of 5 mm on the right side was set as a perfect matching layer to reduce the border effect. The area of the focal point was set just below the tissue-model surface to increase the energy of the shear wave in the second layer. Twenty-five probes were set along the depth direction from the surface of the tissue model, 2.5 mm away from the load center, at intervals of 0.25 mm. The longitudinal displacements of the tissue model at the probes were acquired for 10 ms.
Table 1. Material parameters of the simulation model
The surface phase-velocity curve was obtained by processing the displacement data. As shown in Fig. 3(b), the high-frequency phase velocity in the surface phase velocity curve tends to be stable, and it converges to the Rayleigh-wave velocity, ${C_{R1}}$, of the top layer. The specific data processing is presented in Section 3.3. The signals from 0.2 ms to 4 ms in the spatial-temporal maps were selected to obtain the wavenumber-frequency domain. The Young's modulus of the first layer, which was calculated using the high-frequency SAW, was $59.6 \pm 1.31$ kPa, and the error of the Young's modulus was 0.6%. The phase-velocity curve was obtained from the displacement data at the depth of the second layer, i.e., 1 mm from the interlayer boundary. As shown in Fig. 3(d), the high-frequency phase velocity also tends to be stable, and it converges to the shear wave velocity, ${C_{s2}}$, of the second layer. According to Eq. (9), the Young's modulus of the second layer was $19.9 \pm 0.31$ kPa with an error of 0.5%. The simulation results verify the feasibility of combining the SAW and shear waves to calculate the elasticity of multi-layer samples.
Fig. 3. The simulation results. (a), (b), and (c) are the spatial-temporal displacement gradient map, the wavenumber-frequency diagram, and phase-velocity curve extracted from the surface, respectively. (d), (e), and (f) are the spatial-temporal displacement gradient map, the wavenumber-frequency diagram, and phase-velocity curve extracted from the inside of the second layer, respectively. Mean value ${\pm}$ standard deviation. The black lines represent the wavenumber with the maximum intensity at each frequency between high and low cut-off frequencies. The dotted red lines represent the high-frequency average phase velocities calculated from the red regions.
3. Materials and methods
3.1 Phantom samples
In this study, the agar-agar bilayer phantoms were used to investigate the elasticity of bilayer materials. Agar (Marklin, China) phantoms with different concentrations (1%, 1.5%, and 2%) were used to represent materials with different elasticities. Only 0.3% titanium dioxide (Marklin, China) was added to the second-layer agar solution. This was done to increase the scattering of the layer to distinguish the top layer and the second layer of the bilayer phantom sample in the B-scan image. The liquid agar was poured into a rectangular container with a width of 20 mm and a height of 25 mm. To obtain a bilayer phantom, the agar solution of the top layer was poured onto the solidified bottom layer. Table 2 lists the detailed parameters of the bilayer phantoms. The thickness of the top layer was measured using the OCT structural image. Because the epidermis is harder than the dermis, the designed concentration of the top-layer agar was higher than that of the second-layer agar. Additionally, three types of homogeneous phantoms with concentrations of 1%, 1.5%, and 2% and a size of approximately 40 (L) ${\times} $ 40 (W) ${\times} $ 12 mm (H) were prepared. The phantoms were used for traditional SAW detection as a reference for the bilayer phantom experiments. The mass and volume of the phantoms were measured to calculate the density of those.
Table 2. Parameters of the experimental phantoms
3.2 Experimental setup
As shown in Fig. 4, the OCE system used in this study comprised an elastic wave excitation subsystem and a detection subsystem. The excitation subsystem was composed of a function generator (DG4062, RIGOL, China), power amplifier (BT00100-AlphaS-CW, TOMCO, USA), and 5 MHz single-element focused ultrasound transducer (A308S-SU-F-1.25-IN-PTF, Olympus, Japan) with a focal length of 31.8 mm. The function generator was used to generate a sine-wave signal with 1000 cycles per pulse at 5 MHz, which was amplified by the power amplifier to excite the elastic wave in the sample. The focused ultrasound transducer was attached to the bottom of a water-filled tank. The center of the focal point was set below the upper surface of the phantom, which was above the water level, to enhance the strength of the shear waves inside the phantom. The excitation trigger signal and the control signal for the galvanometer scanner were synchronized with the data acquisition.
Fig. 4. Schematic of the optical coherent elastography system. Red dashed box A is the elastic wave detection subsystem, purple dashed box B is the elastic wave excitation subsystem.
PhS-OCT was incorporated into the elastic wave detection subsystem, which included a swept-source laser (HSL-20, Santec, Japan) with a center wavelength of 1300 nm and a scan repetition rate of 50 kHz. Light from the laser passed through a 90/10 coupler; light from the 10% and 90% outputs entered the reference and sample arms, respectively. The light that returned from each arm interfered in a 50/50 coupler. Then, the interference signals were transformed into an electrical signal by a balanced detector; and detected using a data acquisition card. The axial displacement of the given location ${u_z}({z,x,t} )$ can be calculated as follows:
The PhS-OCT system was operated in the M-B scan mode to track the propagation of the elastic wave. Each M-scan was composed of 512 A-scans, which were acquired at the same location. At each location, an M-scan was acquired for each excitation. A total of 512 M-scans were horizontally repeated across the entire imaging plane to form a complete M-B scan. In this study, the trigger signal of the elastic wave excitation was delayed by 1 ms for each M-scan to prevent data loss caused by the system acquisition delay. The lateral distance of the imaging area was 4.12 mm with an axial resolution of 11.2 $\mathrm{\mu m}$ and lateral resolution of 14.5 $\mathrm{\mu m}$. The starting position of the OCT scanning was approximately 2.5 mm away from the ultrasonic excitation position to avoid the near-field effect.
3.3 Data processing
In this study, the OCT-structure image information was extracted from the M-B scanning data that was collected by the OCE system to identify the regions of different sample layers. To estimate the Young's moduli of different sample layers, the velocity of the SAW in the high-frequency region and the shear wave velocity were extracted from the phase-velocity dispersion curves on the surface and inside of each layer. Figure 5 shows the main data-processing steps for the bilayer samples.
Fig. 5. Main data processing steps for the elastic wave imaging of the bilayer sample. (I) is the B-scan image obtained from the M-B-scan data of the bilayer agar phantom. The size of the B-scan image for the phantom is 4.35 ${\times} $ 4.12 mm (depth ${\times} $ lateral distance). The red line represents the identified surface of the phantom, and the blue line is the position of the second layer of the phantom where the elastic wave dispersion curve is calculated. (a) and (b) are the spatial-temporal displacement gradient maps of the surface and the interior of the second layer, respectively. (c) and (d) are the wavenumber-frequency domains after Fourier transforms of the rectangular regions of (a) and (b), respectively. The black lines represent the wavenumber with the maximum intensity at each frequency between high and low cut-off frequencies. (e) and (f) are the phase-velocity dispersion curves of the surface and the interior of the second layer, respectively. The dotted red lines represent the high-frequency average phase velocities calculated from the red regions.
First, the unwrapping phase was extracted separately from the sample’s surface and any depth in the second layer. When the phase inside the sample was extracted, phase correction was required for eliminating the artifacts caused by surface motion, and refractive-index mismatch between air and the sample [26]. The phase-velocity dispersion curves were calculated using the phases obtained from the surface and the interior of each layer separately using the same procedure. The spatial-temporal slice of the axial displacement gradient was calculated using Eq. (10), as shown in Fig. 5(a-b). The regions of interest (ROIs) in the spatial-temporal slice were selected (the rectangular regions marked by the black rectangular in Fig. 5(a-b)), where the displacements significantly fluctuate with time and were affected by the SAW and shear waves. The sizes of the ROIs were set according to the reasonable velocity range of SAWs and shear waves in the sample. A two-dimensional Fourier transform was applied to the signal in the rectangular regions to obtain the wavenumber-frequency domain. Then, the threshold filter was applied to remove the signal with an intensity lower than −8 dB of the maximum intensity value in the wavenumber-frequency domain (Fig. 5(c-d)). The high cut-off frequency was determined by removing the signals with the low OCT signal-to-noise ratio. The low-frequency wave generated on the surface of the sample can be reflected from the sample boundaries to generate the complex wave field [27]. Thus, the low cut-off frequencies at 0.4kHz were determined according to the depth of the samples. Finally, the phase-velocity dispersion curve was obtained by selecting the wavenumber of the maximal intensity at each frequency and calculating the corresponding phase velocity ${C_{phase}}$ as follows:
where f is the frequency and k is the wavenumber. The phase-velocity curve, shown in Fig. 5(e-f), is relatively flat at high frequencies. In the case of the air-medium boundary, the high-frequency-limit velocity of the phase-velocity curve obtained from the surface is identified as the Rayleigh-wave velocity ${C_R}$. The Young's modulus E of the top layer was obtained using Eq. (8). The phase velocity in the high-frequency region obtained in the second layer at a depth of about 0.1–0.8 mm from the interlayer boundary converged, and its value corresponded to the shear wave velocity. The Young's modulus of the second layer was obtained using the shear wave velocity according to Eq. (9).If the sample is a multi-layer sample, the calculation method used to obtain the Young's modulus of the layers below the second layer is the same as that of the second layer. As observed in Fig. 5(a), the sample surface had a fast wave [28] with lower energy (marked by the red dashed line) compared to the SAW. Therefore, the window position that is far from the excitation distance in the spatial-temporal slice was selected, where the displacement fluctuations with time were primarily affected by the SAW.
4. Results
The elasticities of the different layers of the three bilayer phantoms were measured using the proposed method. The phase-velocity dispersion curves obtained from the surfaces and the interiors of the second layers of the three samples were plotted in Fig. 6. The dispersion curves of the elastic waves tended to be stable at high frequencies, which is in agreement with the simulations. The dotted red line in Fig. 6 represents the high-frequency average phase velocity. The Rayleigh-wave velocities (a-c), and shear wave velocities (d-f) were indicated in the figure. The measured densities of the 1%, 1.5%, and 2% agar phantoms were 964 kg/m3, 1024 kg/m3and 1040 kg/m3, respectively, and Poisson's ratio was 0.499. By substituting the Rayleigh-wave velocity into Eq. (8), the Young's moduli of the top layer of the three samples were $54.9 \pm 2.32$ kPa, $56.2 \pm 1.12$ kPa, and $115 \pm 3.28$ kPa, respectively. By substituting the shear wave velocity into Eq. (9), the Young's moduli of the second layer of the three samples were $22.9 \pm 1.84$ kPa, $23.0 \pm 1.25$ kPa, and $54.4 \pm 1.86$ kPa, respectively.
Fig. 6. Elasticity test results of different layers of samples. (a), (b), and (c) are the surface elastic wave dispersion curves of sample 1, sample 2, and sample 3 respectively. (d), (e), and (f) are the dispersion curves of elastic waves in the second layer of samples 1, 2, and 3, respectively. Mean value ${\pm}$ standard deviation.
Traditional SAW elasticity tests were performed on three types of homogeneous phantoms. The results were used as reference values for the experimental results obtained from the bilayer phantoms. As shown in Fig. 7(a–c), when calculating the Young's moduli of the homogeneous phantoms, the unwrapping phase of the phantom surface was selected to calculate the SAW phase-velocity curve. The average phase velocity above the center frequency was calculated as the Rayleigh-wave velocity [29]. The center frequency is defined as the frequency at the highest energy in the wavenumber-frequency domain [30]. Each sample was tested five times at different locations to obtain the average phase-velocity dispersion curve, as shown in Fig. 7(d–f). Table 3 summarizes the calculation results. The Rayleigh-wave velocities and the Young's moduli of the homogeneous phantoms are presented as mean values with standard error in the table.
Fig. 7. Algorithm processing and phase-velocity dispersion curves of homogeneous phantoms. (a) is the B-scan image, and the dotted red line represents the surface of the phantom. The size of the B-scan image for the phantom is 4.5 ${\times} $ 4.12 mm (depth ${\times} $ lateral distance). (b) The spatial-temporal displacement gradient image of the surface of the 1% agar phantom. (c) Wavenumber-frequency diagram after the Fourier transform from the rectangular region of (b). The red asterisk indicates the center frequency of the surface wave. The black lines represent the wavenumber with the maximum intensity at each frequency between high and low cut-off frequencies. (d), (e), and (f) are the 1%, 1.5%, and 2% agar phantom surface wave phase-velocity dispersion curves, respectively. SE is the standard error.
Table 3. Experimental results of homogeneous phantom
Uniaxial mechanical compressional tests (5943, Instron Corp., MA, USA) were performed on homogeneous agar phantoms at identical concentrations to evaluate the fidelity of the measurement using the SAW elastography for the homogeneous phantoms. The phantom used for mechanical compression was 40 mm in diameter and 16 mm in height. The uniaxial compression was performed using a re-load force of 0.1 N at 0.25 mm/s. The slope of the initial linear part of the force curve was selected to calculate the Young's modulus [31]. The Young's moduli of the 1%, 1.5%, and 2% agar phantom were 24.4 ± 0.9 kPa, 55.3 ± 1.7 kPa, 121.3 ± 5.0 kPa, respectively. The results of the mechanical compression tests and SAW elastography were in good agreement with the estimated Young's moduli. The slight difference between the OCE and mechanical test results may have occurred because the mechanical estimate is based on larger strain values. The mechanical response of the agar phantom based on high-strain values, where the stress-strain curve of the agar phantom is not completely linear, is different from that based on low-strain values [32].
The test results obtained using this method were compared with the elastic results obtained using the homogeneous phantoms, which are listed in Table 4. The Young's moduli of the homogeneous phantoms were calculated as the standard values of the bilayer phantom experiment using the traditional SAW model. The corresponding estimated error of each sample relative to the SAW test result of the homogeneous phantom is listed in Table 4, and the maximum estimation error is 2.2%. The experimental results show that the proposed method accurately estimated the elasticity of layered tissue-mimicking phantoms.
Table 4. Test results of the bilayer phantoms
5. Discussion
The frequency at which the phase velocity tends to stabilize in the dispersion curve of the elastic wave obtained from the surface is related to the penetration depth of the SAW. When the penetration depth of the SAW is less than the top-layer thickness, the phase velocity approaches and stabilizes at a constant value that is related only to the elasticity of the top layer. In the bilayer phantom experiment of this study, the top-layer thickness of sample 2 was less than that of sample 1. As shown in Fig. 6(a-b), the initial frequency of the phase-velocity curve at which the surface of sample 2 tends to stability is higher than that of sample 1. Therefore, the thinner the top layer, the higher the frequency of the stable SAW phase velocity. A phantom layer with a thickness consistent with the thickness of the epidermis (about 0.2mm) didn’t be created in the experiment, which didn’t affect the verification of the proposed method. In the case of the thinner top layer, wave velocity information at a higher frequency is required to obtain accurate top layer elasticity. The frequency of the Rayleigh wave is given by the following formula [33].
where ${C_R}$ is the Rayleigh-wave velocity, and d is the diameter of the push cross-section. In this study, a focused ultrasound transducer was used as the excitation source. The value of d of the focused ultrasonic transducer should be smaller than that of the mechanical vibrator. Thus, a focused ultrasound excitation can produce higher elastic wave frequency information than a mechanical vibrator. In addition, the diameter of the ultrasonic focal region can be reduced to obtain a higher frequency of SAW.It was observed that the frequency at which the phase velocity tends to stabilize, in the dispersion curve of the elastic wave obtained from the inside of the sample, is related to the penetration depth of the SAW and the strength of the shear wave. The proposed method emphasizes the extraction of the spatial domain and the analysis of the frequency domain. As shown in Fig. 5(d), the SAW and shear wave velocities are divided into two regions with increasing frequency. The dispersion curve of the elastic wave was obtained by selecting the wave velocity with the maximum intensity at each frequency. Therefore, if the center of the focal point of the ultrasound transducer is positioned inside the sample instead of on the surface, the shear wave intensity in the sample will be higher. Thus, the shear wave velocity can be obtained at a frequency that is lower than the frequency that corresponds to the penetration depth of the SAW.
In the case of plane wave propagation, the wave-velocity equation of the SAW in a dual-layer material was derived. The phase velocity of the wave generated by the ultrasonic excitation is not exactly equal to that of the plane wave. However, the slight discrepancy does not affect the trend analysis of the SAW dispersion curve. What’s more, a two-dimensional axisymmetric model of the bilayer tissue was created in a finite element simulation to meet the experiment case.
Compared with the wavelength-depth inversion model used in previous studies, the proposed method provides higher accuracy. The first layer significantly affects the second-layer elasticity calculated by the wavelength-depth inversion model, and only an approximate value of the second-layer elasticity can be obtained. In the proposed method, the shear wave in the second layer is selected, and no demodulation of the combined phase velocities at different depths is required. Moreover, the phase-velocity inversion model fitted using the SAW dispersion curve requires accurate layer-thickness information; however, using the proposed method, only an approximate region of the second layer needs to be identified.
The ultrasound transducer excitation module and OCT detection module used in the proposed method are placed on both sides of the sample, which limits the in-vivo application of this technology in clinical practice. In future research, a ring ultrasonic transducer will be used to realize the ipsilateral excitation and detection of OCE. In addition, in this study, it was assumed that each layer of tissue was isotropic and pure elastic. The anisotropy in the transverse range and the viscoelasticity of tissue will be considered in future research to ensure consistency with the characteristics of real skin.
6. Conclusion
This study presents a method to accurately obtain the depth-resolved elasticity of multi-layer tissue. Using a focused ultrasound transducer as the excitation source and PhS-OCT for detection, the depth-resolved elasticity of the sample was obtained by acquiring the SAW at a high frequency on the surface of the sample and shear waves inside the sample. The feasibility and accuracy of this method were proved via experiments on tissue-mimicking phantoms with different agar concentrations and thicknesses; a maximum estimating error of 2.2% was obtained. This study shows that the proposed method has great potential for wide application in the diagnosis and treatment monitoring of real skin diseases in the future.
Funding
National Natural Science Foundation of China (61971406, 81927801); Science and Technology Commission of Shanghai Municipality (19441910000); Youth Innovation Promotion Association of the Chinese Academy of Sciences.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (61971406 and 81927801), the Science and Technology Commission of Shanghai Municipality (19441910000), and the Youth Innovation Promotion Association CAS.
Disclosures
The authors declare no conflicts of interest.
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. K. Zhou, K. Feng, C. Li, and Z. Huang, “A Weighted Average Phase Velocity Inversion Model for Depth-Resolved Elasticity Evaluation in Human Skin In-Vivo,” IEEE Trans. Biomed. Eng. 68(6), 1969–1977 (2021). [CrossRef]
2. T. R. Tilleman, M. M. Tilleman, and M. H. A. Neumann, “The elastic properties of cancerous skin: Poisson's ratio and Young's modulus,” Isr. Med. Assoc. J. 6, 753–755 (2004).
3. X. M. Zhang, T. G. Osborn, M. R. Pittelkow, B. Qiang, R. R. Kinnick, and J. F. Greenleaf, “Quantitative assessment of scleroderma by surface wave technique,” Med. Eng. Phys. 33(1), 31–37 (2011). [CrossRef]
4. C. M. Botar-Jid, R. Cosgarea, S. D. Bolboaca, S. C. Senila, L. M. Lenghel, L. Rogojan, and S. M. Dudea, “Assessment of Cutaneous Melanoma by Use of Very-High-Frequency Ultrasound and Real-Time Elastography,” Am. J. Roentgenol. 206(4), 699–704 (2016). [CrossRef]
5. P. Ciarletta, L. Foret, and M. Ben Amar, “The radial growth phase of malignant melanoma: multi-phase modelling, numerical simulations and linear stability analysis,” J. R. Soc. Interface 8(56), 345–368 (2011). [CrossRef]
6. N. M. Fisher, J. V. Schaffer, M. Berwick, and J. L. Bolognia, “Breslow depth of cutaneous melanoma: Impact of factors related to surveillance of the skin, including prior skin biopsies and family history of melanoma,” J. Am. Acad. Dermatol. 53(3), 393–406 (2005). [CrossRef]
7. R. M. S. Sigrist, J. Liau, A. El Kaffas, M. C. Chammas, and J. K. Willmann, “Ultrasound Elastography: Review of Techniques and Clinical Applications,” Theranostics 7(5), 1303–1329 (2017). [CrossRef]
8. S. K. Venkatesh, M. Yin, and R. L. Ehman, “Magnetic resonance elastography of liver: Technique, analysis, and clinical applications,” J. Magn. Reson. Imaging 37(3), 544–555 (2013). [CrossRef]
9. B. Babaei, D. Fovargue, R. A. Lloyd, R. Miller, L. Juge, M. Kaplan, R. Sinkus, D. A. Nordsletten, and L. E. Bilston, “Magnetic Resonance Elastography Reconstruction for Anisotropic Tissues,” Med. Image Anal. 74, 102212 (2021). [CrossRef]
10. P. Sobolewski, M. Maslinska, J. Zakrzewski, L. Paluch, E. Szymanska, and I. Walecka, “Applicability of shear wave elastography for the evaluation of skin strain in systemic sclerosis,” Rheumatol. Int. 40(5), 737–745 (2020). [CrossRef]
11. V. A. Flower, S. L. Barratt, D. J. Hart, A. B. Mackenzie, J. A. Shipley, S. G. Ward, and J. D. Pauling, “High-frequency Ultrasound Assessment of Systemic Sclerosis Skin Involvement: Intraobserver Repeatability and Relationship With Clinician Assessment and Dermal Collagen Content,” J. Rheumatol. 48(6), 867–876 (2021). [CrossRef]
12. D. P. Popescu, L.-P. I. Choo-Smith, C. Flueraru, Y. Mao, S. Chang, J. Disano, S. Sherif, and M. G. Sowa, “Optical coherence tomography: fundamental principles, instrumental designs and biomedical applications,” Biophys. Rev. 3(3), 155–169 (2011). [CrossRef]
13. Z. Ibarra-Borja, C. Sevilla-Gutierrez, R. Ramirez-Alarcon, H. Cruz-Ramirez, and A. B. U’Ren, “Experimental demonstration of full-field quantum optical coherence tomography,” Photonics Res. 8(1), 51–56 (2020). [CrossRef]
14. D. Huang, F. Li, C. Shang, Z. Cheng, and P. K. A. Wai, “Reconfigurable time-stretched swept laser source with up to 100 MHz sweep rate, 100 nm bandwidth, and 100 mm OCT imaging range,” Photonics Res. 8(8), 1360–1367 (2020). [CrossRef]
15. J. Jerwick, Y. Huang, Z. Dong, A. Slaudades, A. J. Brucker, and C. Zhou, “Wide-field ophthalmic space-division multiplexing optical coherence tomography,” Photonics Res. 8(4), 539–547 (2020). [CrossRef]
16. C.-H. Liu, S. Assassi, S. Theodore, C. Smith, A. Schill, M. Singh, S. Aglyamov, C. Mohan, and K. V. Larin, “Translational optical coherence elastography for assessment of systemic sclerosis,” J. Biophotonics 12(12), e201900236 (2019). [CrossRef]
17. X. Liang and S. A. Boppart, “Biomechanical Properties of In Vivo Human Skin From Dynamic Optical Coherence Elastography,” IEEE Trans. Biomed. Eng. 57(4), 953–959 (2010). [CrossRef]
18. A. Ramier, A. M. Eltony, Y. Chen, F. Clouser, J. S. Birkenfeld, A. Watts, and S.-H. Yun, “In vivo measurement of shear modulus of the human cornea using optical coherence elastography,” Sci. Rep. 10(1), 17366 (2020). [CrossRef]
19. K. Zhou, C. Li, S. Chen, G. Nabi, and Z. Huang, “Feasibility study of using the dispersion of surface acoustic wave impulse for viscoelasticity characterization in tissue mimicking phantoms,” J. Biophotonics 12(1), e201800177 (2019). [CrossRef]
20. C. Li, G. Guan, X. Cheng, Z. Huang, and R. K. Wang, “Quantitative elastography provided by surface acoustic waves measured by phase-sensitive optical coherence tomography,” Opt. Lett. 37(4), 722–724 (2012). [CrossRef]
21. K. Zhou, N. Le, Z. Huang, and C. Li, “High-intensity-focused ultrasound and phase-sensitive optical coherence tomography for high resolution surface acoustic wave elastography,” J. Biophotonics 11(2), e201700051 (2018). [CrossRef]
22. J. Chen, Y. Hu, X. Lin, H. Liu, and C. Sun, “Optical coherence elastography of bilayer soft tissue based on harmonic surface wave spectroscopy,” Opt. Lasers Eng. 145, 106667 (2021). [CrossRef]
23. M. Bernal, I. Nenadic, M. W. Urban, and J. F. Greenleaf, “Material property estimation for tubes and arteries using ultrasound radiation force and analysis of propagating modes,” J. Acoust. Soc. Am. 129(3), 1344–1354 (2011). [CrossRef]
24. H.-C. Wang, S. Fleming, Y.-C. Lee, S. Law, M. Swain, and J. Xue, “Laser ultrasonic surface wave dispersion technique for non-destructive evaluation of human dental enamel,” Opt. Express 17(18), 15592–15607 (2009). [CrossRef]
25. C. Li, Z. Huang, and R. K. Wang, “Elastic properties of soft tissue-mimicking phantoms assessed by combined use of laser ultrasonics and low coherence interferometry,” Opt. Express 19(11), 10153–10163 (2011). [CrossRef]
26. S. Song, Z. Huang, and R. K. Wang, “Tracking mechanical wave propagation within tissue using phase-sensitive optical coherence tomography: motion artifact and its compensation,” J. Biomed. Opt. 18(12), 121505 (2013). [CrossRef]
27. X. Zhang, “Identification of the Rayleigh surface waves for estimation of viscoelasticity using the surface wave elastography technique (L),” J. Acous. Soc. Am. 140(5), 3619–3622 (2016). [CrossRef]
28. A. Ramier, B. Tavakol, and S.-H. Yun, “Measuring mechanical wave speed, dispersion, and viscoelastic modulus of the cornea using optical coherence elastography,” Opt. Express 27(12), 16635–16649 (2019). [CrossRef]
29. Z. Jin, Y. Zhou, M. Shen, Y. Wang, F. Lu, and D. Zhu, “Assessment of corneal viscoelasticity using elastic wave optical coherence elastography,” J. Biophotonics 13, e201960074 (2020). [CrossRef]
30. C. A. Carrascal, S. Chen, M. W. Urban, and J. F. Greenleaf, “Acoustic Radiation Force-Induced Creep-Recovery (ARFICR): A Noninvasive Method to Characterize Tissue Viscoelasticity,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 65(1), 3–13 (2018). [CrossRef]
31. U. Hamhaber, V. A. Grieshaber, J. H. Nagel, and U. Klose, “Comparison of quantitative shear wave MR-elastography with mechanical compression tests,” Magn. Reson. Med. 49(1), 71–77 (2003). [CrossRef]
32. Z. Han, J. Li, M. Singh, C. Wu, C.-H. Liu, S. Wang, R. Idugboe, R. Raghunathan, N. Sudheendran, S. R. Aglyamov, M. D. Twa, and K. V. Larin, “Quantitative methods for reconstructing tissue biomechanical properties in optical coherence elastography: a comparison study,” Phys. Med. Biol. 60(9), 3531–3547 (2015). [CrossRef]
33. Y. H. Sohn and S. Krishnaswamy, “Mass spring lattice modeling of the scanning laser source technique,” Ultrasonics 39(8), 543–551 (2002). [CrossRef]
