Elasticity and Viscoelasticity Imaging Using Small Particles Exposed to an External Force

1. Introduction

There are various imaging methods such as X-ray imaging, computed tomography (CT), magnetic resonance imaging (MRI), ultrasound imaging, and optical imaging to detect abnormalities [1]. X-ray imaging can provide detection of lesions when they absorb X-ray radiation differently from normal tissue. MRI produces multiplanar imaging and clear soft tissue contrast by distinguishing structures based on their water content, etc. Sonography can be used to detect abnormalities due to their hypoechoic or hyperechoic appearance on the B-scan and/or the acoustic shadowing effect they cause. However, these medical imaging techniques display tissue properties that are generally not related to the distinct mechanical properties of tissue, such as shear modulus and viscosity. The diagnostic value of characterizing mechanical properties of tissues resulted in development of methods such as magnetic resonance elastography (MRE) and ultrasound elastography for imaging tissue mechanical properties [2]. In these methods, a static, quasi-static or dynamic force is applied to the tissue and the deformation due to the force is observed or measured using a conventional imaging technique, such as MRI and ultrasound imaging. Many attempts have been made to develop elasticity imaging methods with distinct capabilities, especially using ultrasound technology [3,4]. Some of these methods include vibro-acoustography [5], shear wave elasticity imaging (SWEI) [6], acoustic radiation force impulse imaging (ARFI) [7], harmonic motion imaging (HMI) [8], supersonic shear imaging (SSI) [9], and shear wave spectroscopy (SWS) [10]. Readers may refer to the review papers for further information about ultrasound-based elasticity imaging methods [11,12] and clinical use of ultrasound elastography [13,14,15,16,17].

In addition to the aforementioned methods, acoustic radiation force based elasticity methods that use gas bubbles or rigid spheres located within the tissue have been proposed to identify the mechanical properties of the surrounding medium [18,19,20,21,22]. However, these techniques use a high-powered laser to generate the bubble inside the medium which limits its application to shallow targets and requires local destruction of the material [20]. Additionally, there is a need to locate a bubble or rigid sphere within the tissue in a way which may not be permissible in human tissue. A more practical method based on the use of particles located at the tissue interface has recently been proposed for the identification of tissue viscoelastic properties [23,24]. It should be remembered that the acoustic radiation force-based elasticity imaging methods that remotely palpate the tissue produce a complicated displacement field [20]. The displacement field depends on the transducer focusing characteristics and target geometry. For these absorption-based techniques, the focal volume affects elasticity reconstruction. On the other hand, for elasticity imaging based on the use of small particles, particle displacement is independent of the focal zone of the transducer and is affected only by the viscoelastic properties of the local tissue. As the object in the elasticity imaging techniques based on the use of small particles can be considered as almost a discrete target, simple and accurate elasticity reconstruction can be performed. Thanks to the analytical models for the response of the small particles located within the medium or at the tissue interface, in addition to elastic properties of tissue, the density and viscosity of tissue can be identified. It should be noted that the paramaters such as viscosity can be used as additional biomarkers for the evaluation of tissues [25].

Elasticity imaging based on the absorption mechanism and the use of small particles is illustrated in Figure 1. Displacements change widely over the region exposed to the acoustic radiation force due to differences in the applied force over the same region for the elasticity imaging methods based on the absorption mechanism (Figure 1a). On the other hand, for the elasticity imaging techniques based on the use of small particles, only the particle displacement is measured, which can be achieved to align with the focus of the applied force (Figure 1b). Because the radiation force applied to the particle located within the tissue or at the tissue interface is generally much greater than that arising in purely absorbing media, lower acoustic power can be used in moderately attenuating media for the elasticity imaging techniques based on the use of small particles. Similar displacements can be achieved with far lower acoustic intensities in particle-based methods. The elasticity imaging methods based on the use of small particles located within the tissue and at the tissue interface have potential to determine local tissue properties with improved measurement resolution at variable depths.

Elasticity imaging based on the use of small particles is illustrated in Figure 2. As stated before, the particle can be located within the tissue or at the tissue interface for elasticity imaging based on the use of small particles. Although the particle is pushed mostly using ultrasound, the excitation force ($f_{\mathrm{e}}$) can be applied using any other method such as magnetic or mechanical excitation. Monitoring the displacement of the particle ($u$) using ultrasound is practical, yet particle imaging can be performed using other methods, such as optical imaging or magnetic resonance imaging. Although there is a need to use a coupling material such as gel or water between the tissue and transducer for ultrasonic excitation, there is no need for a coupling medium (or only air is around the medium) for other excitation methods such as a magnetic force. As stated before, thanks to the analytical models for the small particles located within the tissue and at the tissue interface, in addition to elastic properties of tissue, the density, Poisson’s ratio and viscosity of tissue can be identified in the elasticity imaging based on the use of small particles. Therefore, instead of the use of the term elasticity imaging, we suggest the use of the term viscoelasticity imaging.

Although there have been a huge number of research studies (e.g., [26,27]) and there are some review papers (e.g., [11,12]) on the elasticity imaging methods based on the absorption mechanism, there have been a limited number of research studies on the elasticity imaging methods based on the use of small particles. To the best knowledge of authors, there is currently no review study focusing on the elasticity imaging techniques based on the use of small particles located within the tissue and at the tissue interface. We believe that elasticity imaging based on the use of small particles have potential to change the field, as they can provide the determination of accurate and local tissue elastic properties as well as density and viscosity. Here, we examined the mathematical models for small particles located within the tissue and at the tissue interface. In addition, we compared the mathematical models for small particles located within the tissue and at the tissue interface and provided elasticity/viscoelasticity imaging results based on these models and experiments. This paper is the first review study on the elasticity and viscoelasticity imaging based on the use of small (nanometer-, micrometer-, or millimeter-sized) particles (e.g., bubbles and rigid and soft spheres) located within the tissue and at the tissue interface exposed to static and dynamic external forces. Elasticity/viscoelasticity imaging methods based on the use of particles (bubbles and rigid spheres) located within the tissue are presented in Section 2. Elasticity/viscoelasticity imaging techniques based on the use of particles (bubbles, rigid and soft spheres) located at the tissue interface are covered in Section 3. It should be noted that, the viscoelastic medium was mostly modelled using the Kelvin-Voigt model in the models available in the literature. In this paper, we extended the responses of particles located within the tissue and at the tissue interface and generalized the relations so that, in addition to the Kelvin-Voigt model, any viscoelastic rheological model such as the Maxwell, generalized Maxwell, standard linear solid and Kelvin-Voigt fractional derivative models can be used to simulate the medium. Furthermore, based on the new advances, in this paper, we updated some of the models for the responses of the particles located within the tissue and at the tissue interface available in the literature. The mathematical models for the small particles located within the tissue and at the tissue interface are compared and the elasticity and viscoelasticity imaging methods based on the use of the particles are evaluated in Section 4. Some concluding remarks are given in Section 5.

2. Elasticity/Viscoelasticity Imaging Using Particles Located within Tissue

2.1. Bubble Located within Tissue

There many mathematical models for radial oscillations of a bubble inside the elastic and viscoelastic media exposed to sound (e.g., [28,29,30]). However, these models generally are the extensions of the well-known Rayleigh-Plesset equation (e.g., [31]) and were proposed for studying the dynamics of bubbles inside the tissue, rather than for using them for elasticity imaging. The first important study for the purpose of elasticity imaging using a bubble located within the medium was performed by Ilinskii et al. [19] in 2005. In this model, the radiation pressure on the target is estimated by integrating the acoustic field in the surrounding medium with respect to the target surface (i.e., a technique similar to the one applied to estimate Stockes’ drag of a solid sphere inside viscous incompressible liquid was used [32]). Hence, by assuming the elastic medium as isotropic, homogeneous, incompressible, and non-viscous, Ilinskii et al. [19] proposed models for the static and dynamic displacements of a bubble located within the soft medium in response to acoustic radiation force. The model for an external static force $f_{\mathrm{e}}=f_0$ and elastic medium is given by [19]:

$$u=\frac{f_0}{4\pi GR}$$

(1)

where $G$ is the shear modulus of the medium, and $R$ is the radius of the bubble. It should be noted that the Young’s modulus is related to the shear modulus by $E=2G\left(1+\nu \right)$ for homogeneous isotropic materials. Ilinskii et al. [19] obtained the following equation of motion for the dynamic displacement of a bubble located within the soft medium in response to acoustic radiation force:

$$\underset{\_}{u}+\underset{\_}{\dot{u}}+\frac{1}{6}\underset{\_}{\ddot{u}}+\frac{1}{18}\underset{\_}{\stackrel{⃛}{u}}={\underset{\_}{f}}_{\mathrm{e}}$$

(2a)

where

$$\underset{\_}{u}=\frac{u\left(t\right)}{R}$$

(2b)

$${\underset{\_}{f}}_{\mathrm{e}}=\frac{f_{\mathrm{e}}}{4\pi G{R}^2}$$

(2c)

$$\underset{\_}{t}=\frac{t\sqrt{G/\rho }}{R}$$

(2d)

Here, $\rho$ is the density of the medium, and $t$ shows the time. The solution of Eq. (2) provides the dynamic force-displacement relation for a bubble located within the soft elastic medium. The solution of Eq. (2) can be determined analytically (e.g., using the Laplace transformation) or numerically. Readers can refer to the reference [19] to see the details for the derivations of the models (i.e., Eqs. (1) and (2)). By assuming the medium as isotropic, homogeneous, and incompressible, a mathematical model for the dynamic response of a bubble located within the soft viscoelastic medium in response to acoustic radiation force was proposed. The frequency-domain equation that relates the displacement of the bubble located within the soft viscoelastic medium for a dynamic force is given by [19]:

$$4\pi \left(G-j\omega \eta \right)RU\left(1-jkR-\frac{1}{6}{k}^2{R}^2+\frac{1}{18}j{k}^3{R}^3\right)=F_{\mathrm{e}}$$

(3)

where $\eta$ is the viscosity of the medium, $k=\frac{\omega }{\sqrt{\left(G-j\omega \eta \right)/\rho }}$ is the complex wave number of the shear wave with frequency $\omega$, $j=\sqrt{-1}$, $F_{\mathrm{e}}$ is the Fourier transform of the external force $f_{\mathrm{e}}\left(t\right)$, and $U$ is a spectral component of the bubble displacement $u\left(t\right)$. Readers can refer to the reference [19] to see the details for the derivation of the model (i.e., Eq. (3)). The response to a rectangular pulse with an amplitude of $f_0$ and a duration of $\tau$ represents the impulse response for small $\tau$ values and the step response for large $\tau$ values with $0\le t\le \tau$. Therefore, the external force is assumed as a rectangular pulse here. It is noted that the Fourier transform of the rectangular pulse is $F_{\mathrm{e}}=-\frac{j{f}_0}{\omega }\left(e^{j\omega \tau }-1\right)$. In conclusion, the dynamic displacement of the bubble located within the soft viscoelastic medium can be calculated using the inverse Fourier transformation by [33]:

$$u\left(t\right)=-\frac{j{f}_0}{24{\pi }^{2}R}\underset{-\infty }{\overset{\infty }{\int }}\frac{\left(e^{j\omega \tau }-1\right)\left(3-jkR\right)e^{-j\omega t}}{\omega \left(G-j\omega \eta \right)\left(1-jkR-\frac{1}{6}{k}^2{R}^2+\frac{1}{18}j{k}^3{R}^3\right)}\mathrm{d}\omega$$

(4)

By dividing the time duration into N points (e.g., N = 1000), and the calculations can be repeated over the entire period of interest to find the dynamic response of the bubble located within the medium numerically using Eq. (4).

Exploiting the static model (i.e., Eq. (1)), Erpelding et al. [20] identified the Young’s moduli of gelatin phantoms. For this purpose, acoustic radiation force was applied to individual bubbles generated by laser-induced optical breakdown. A two-element confocal ultrasonic transducer generated acoustic radiation force on individual bubbles while monitoring their displacement within the medium. Tone burst pushes of varying duration were applied by the outer element at 1.5 MHz. The inner element received pulse-echo recordings at 7.44 MHz. Maximum bubble displacements were used to the determine the Young’s moduli for different gel phantoms; average Young’s moduli varied from 1.7 kPa to 4.8 kPa as the gel concentration was increased from 5 to 9%. Erpelding et al. [34] used the same technique to map age-related elasticity changes in porcine lenses. Hollman et al. [35] used the same technique to map elasticity in human lenses; for middle-age lenses, Young’s modulus ranged from 5.2 kPa in the center to 1.1 kPa on the periphery and for old-age lenses, Young’s modulus ranged from 10.6 kPa in the center to 1.4 kPa on the periphery. Similarly, Mikula et al. [36] exploited this simple model to determine corneal elasticity. For this purpose, a low-frequency, high-intensity acoustic force is used to displace a femtosecond laser-generated microbubble, while high-frequency, low-intensity ultrasound is used to monitor the position of the microbubble within the ex vivo human cornea. They determined the Young’s modulus of E = 1.39 ± 0.28 kPa for the central anterior cornea and E = 0.71 ± 0.21 kPa for the central posterior cornea. In another study, using the same technique, Mikula et al. [37] determined elasticity map across the ex vivo human cornea.

Using the dynamic model presented above, Yoon et al. [38] estimated the mechanical properties of gelatin phantoms using a laser-induced microbubble interrogated by the acoustic radiation force. For this purpose, A 3.5 MHz ultrasound transducer was used to generate acoustic radiation force to excite a laser-induced microbubble and the motion of the microbubble was tracked using a 25 MHz imaging transducer. By using both maximum displacements and the times of maximum displacement of bubbles, they determined the Young’s modulus of the 3% gelatin phantom as E = 4.37 ± 0.26 kPa. In another study, exploiting the same model, the Young’s moduli and viscosities of ex vivo bovine and porcine crystalline lenses were identified using the time characteristics of the microbubble dynamics [39].

In the model presented above (i.e., Eqs. (3) and (4)), the viscoelastic medium was modelled using the Kelvin-Voigt model. In this paper, we extended the response of a bubble located within the soft viscoelastic medium and generalized the relations so that, in addition to the Kelvin-Voigt model, any viscoelastic rheological model such as the Maxwell, generalized Maxwell, standard linear solid and Kelvin-Voigt fractional derivative models can be used to simulate the medium. In conclusion, the generalized model for the response of a bubble located within the soft viscoelastic medium can be represented as:

$$u\left(t\right)=-\frac{j{f}_0}{24{\pi }^{2}R}\underset{-\infty }{\overset{\infty }{\int }}\frac{\left(e^{j\omega \tau }-1\right)\left(3-jkR\right)e^{-j\omega t}}{\omega G^{*}\left(1-jkR-\frac{1}{6}{k}^2{R}^2+\frac{1}{18}j{k}^3{R}^3\right)}\mathrm{d}\omega$$

(5)

where $G^{*}\left(\omega \right)=G^{\text{'}}\left(\omega \right)-j{G}^{\text{'}\text{'}}\left(\omega \right)$ is the complex modulus of the medium material; $G^{\text{'}}\left(\omega \right)$ and $G^{\text{'}\text{'}}\left(\omega \right)$ representing the storage and loss moduli, respectively. Readers may refer to the reference [22] for the expressions of $G^{\text{'}}\left(\omega \right)$ and $G^{\text{'}\text{'}}\left(\omega \right)$ for different rheological models.

2.2. Sphere Located within Tissue

The response of a sphere located within the medium has been studied for over the past 70 years. Oestreicher [40] obtained mechanical impedance of an oscillating sphere inside a viscoelastic medium in frequency domain by modelling the medium using the Kelvin-Voigt model. Chen et al. [18] used the measured response of a vibrating sphere inside gelatin phantom excited by modulated ultrasound and the model in the reference [40] for the determination of the material viscoelastic properties. However, the first important study for the purpose of elasticity imaging using a sphere located within the medium was performed by Ilinskii et al. [19] in 2005. By assuming the elastic medium as isotropic, homogeneous, incompressible, and non-viscous, Ilinskii et al. [19] proposed models for the static and dynamic displacements of a rigid sphere located within the medium in response to the acoustic radiation force. The model for a static force $f_0$ and elastic medium is given by [19]:

$$u=\frac{f_0}{6\pi GR}$$

(6)

For a dynamic force $f_{\mathrm{e}}$ and an elastic medium, they obtained the following equation of motion [19]:

$$\underset{\_}{u}+\underset{\_}{\dot{u}}+\frac{1}{9}\left(1+2\frac{{\rho }_{\mathrm{s}}}{\rho }\right)\underset{\_}{\ddot{u}}={\underset{\_}{f}}_{\mathrm{e}}$$

(7a)

where

$$\underset{\_}{u}+\underset{\_}{\dot{u}}+\frac{1}{9}\left(1+2\frac{{\rho }_{\mathrm{s}}}{\rho }\right)\underset{\_}{\ddot{u}}={\underset{\_}{f}}_{\mathrm{e}}$$

(7a)

$$\underset{\_}{u}=\frac{u\left(t\right)}{R}$$

(7b)

$${\underset{\_}{f}}_{\mathrm{e}}=\frac{f_{\mathrm{e}}}{6\pi G{R}^2}$$

(7c)

$$\underset{\_}{t}=\frac{t\sqrt{G/\rho }}{R}$$

(7d)

here, ${\rho }_{\mathrm{s}}$ is the density of the sphere. The solution of Eq. (7) provides the dynamic force-displacement relation for a sphere located within the elastic medium. The solution of Eq. (7) can be determined analytically (e.g., using the Laplace transformation) or numerically. Readers can refer to the reference [19] to see the details for the derivations of the models (i.e., Eqs. (6) and (7)). By assuming the medium as isotropic, homogeneous, and incompressible, Aglyamov et al. [21] developed a mathematical model for the dynamic response of a sphere located within the viscoelastic medium in response to acoustic radiation force in 2007. The frequency-domain equation that relates the displacement of the sphere located within the viscoelastic medium for a dynamic force is given by [21]:

$$-m_{\mathrm{s}}{\omega }^{2}U+6\pi \left(G-j\omega \eta \right)RU\left(1-jkR-\frac{1}{9}{k}^2{R}^2\right)=F_{\mathrm{e}}$$

(8)

where $m_{\mathrm{s}}=\frac{4}{3}\pi R^3{\rho }_{\mathrm{s}}$ is the mass of the sphere. Readers can refer to the reference [21] to see the details for the derivation of the model (i.e., Eq. (8)). By performing some operations, the dynamic displacement of the sphere located within the viscoelastic medium can be found by [33]:

$$u\left(t\right)=-\frac{j{f}_0}{12{\pi }^{2}R}\underset{-\infty }{\overset{\infty }{\int }}\frac{\left(e^{j\omega \tau }-1\right)e^{-j\omega t}}{\omega \left[-\frac{2}{9}{R}^{2{\rho }_{\mathrm{s}}}{\omega }^2+\left(G-j\omega \eta \right)\left(1-jkR-\frac{1}{9}{k}^2{R}^2\right)\right]}\mathrm{d}\omega$$

(9)

By dividing the time duration into N points (e.g., N = 1000), and the calculations can be repeated over the entire period of interest to find the dynamic response of the sphere located within the viscoelastic medium numerically using Eq. (9).

Chen et al. [18], by using the measured response of a vibrating sphere inside the gelatin phantom excited by modulated ultrasound and the expression for the mechanical impedance of a vibrating sphere inside the gelation phantom, determined the viscoelastic material properties of gelatin. Karpiouk et al. [41] used the dynamic model presented above for a sphere inside the tissue to determine the shear modulus of the gelatin phantom. For this purpose, a short pulse produced by a 1.5-MHz highly focused ultrasound transducer was used to initiate the motion of a rigid sphere inside the elastic medium and another 25 MHz focused ultrasound transducer operating in pulse-echo mode was used to track the displacement of the sphere. They determined the shear modulus of the gelatin phantom using maximum displacement of the sphere and the time required for the rigid sphere to reach the maximum displacement (no need for knowing the applied force). Using a similar test setup, Shih et al. [42] exploited the spatiotemporal behavior of the displacement of the sphere and determined the viscoelastic properties of porcine clot (with different hematocrits); the value of shear modulus of blood clot decreased from 585±127 Pa at plasma to 168±26 Pa at 40% hematocrit and the viscosity of blood clot increased from 0.125±0.025 Pa s at plasma to 0.28±0.12 Pa s at 40% hematocrit. Similarly, Huang et al. [43], by using the dynamic model, identified the shear modulus and viscosity of porcine blood clots based on the spatiotemporal behavior of the sphere displacement.

It should be noted that the medium was modelled using the Kelvin-Voigt model in the model presented above (i.e., Eqs. (8) and (9)). Urban et al. [22] extended the vibration response of a sphere located within the viscoelastic medium and generalized the relations so that, in addition to the Kelvin-Voigt model, any viscoelastic rheological model such as the Maxwell, generalized Maxwell, standard linear solid and Kelvin-Voigt fractional derivative models can be used to simulate the medium. In conclusion, the generalized model for the vibration response of a sphere located within the viscoelastic medium can be written as:

$$u\left(t\right)=-\frac{j{f}_0}{12{\pi }^{2}R}\underset{-\infty }{\overset{\infty }{\int }}\frac{\left(e^{j\omega \tau }-1\right)e^{-j\omega t}}{\omega \left[-\frac{2}{9}{R}^{2{\rho }_{\mathrm{s}}}{\omega }^2+G^{*}\left(1-jkR-\frac{1}{9}{k}^2{R}^2\right)\right]}\mathrm{d}\omega$$

(10)

Readers may refer to the reference [22] for the expressions of $G^{*}\left(\omega \right)=G^{\text{'}}\left(\omega \right)-j{G}^{\text{'}\text{'}}\left(\omega \right)$ for different rheological models. Urban et al. [22] determined the viscoelastic properties of gelatin and rubber phantoms by matching the measured and theoretically predicted responses of the spheres inside phantoms using different rheological models for the medium.

It should be noted that the particle located within the tissue has been mostly pushed using acoustic excitation (e.g., [22,41]). However, the particle can be pushed using any proper excitation method such as a magnetic force. Recently, Levy and Oldenburg [44] used magnetic force to excite a sphere inside the medium for the elasticity imaging purpose. They formulated the motion of a magnetic particle exposed to the special magnetomotive ultrasound driving force as in the frequency domain and then obtained an analytical expression by using the inverse Fourier transform as following:

$$\begin{gathered} u\left(t\right)=\frac{f_0}{6\pi GR}\left[1-\frac{1}{\sqrt{{\chi }^2+{\psi }^2}}\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{0}t-{\varphi }_{\mathrm{s}}\right) \\ +{\omega }_0^2{\mathrm{e}}^{-{\omega }_{\mathrm{b}}t}\frac{\sqrt{Z^2+1}}{\sqrt{X^2+Y^2}}\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{\mathrm{a}}t-{\varphi }_{\mathrm{t}}\right)\right] \end{gathered}$$

(11)

where $f_0$ and ${\omega }_0$ are the amplitude and frequency of the driving force, respectively. Readers may refer to the reference [44] for the explanation and expressions of the parameters ${\omega }_{\mathrm{a}}$, ${\omega }_{\mathrm{b}}$, $\chi$, $\psi$, $X$, $Y$, ${\varphi }_{\mathrm{s}}$ and ${\varphi }_{\mathrm{t}}$. Cebrecos et al. [45] evaluated the complex shear modulus of gelatin based on the application of a magnetic force to a sphere inside the medium. As obtaining an ideal rectangular waveform for the current is not possible due to the inductance of the coil, they assume that the temporal evolution of the magnetic force is given by a trapezoidal pulse with a finite rise time $t_{\mathrm{r}}$. The Fourier transform of the trapezoidal pulsed force is $F_{\mathrm{e}}=-\frac{j{f}_0}{t_{\mathrm{r}}{\omega }^2}\left(e^{j\omega \tau }-1\right)\left(e^{j\omega t_{\mathrm{r}}}-1\right)$. Hence, the displacement of the sphere located within the medium can be found using:

$$u\left(t\right)=-\frac{j{f}_0}{12{\pi }^{2}R{t}_{\mathrm{r}}}\underset{-\infty }{\overset{\infty }{\int }}\frac{\left(e^{j\omega \tau }-1\right)\left(e^{j\omega t_{\mathrm{r}}}-1\right)e^{-j\omega t}}{{\omega }^2\left[-\frac{2}{9}{R}^{2{\rho }_{\mathrm{s}}}{\omega }^2+\left(G-j\omega \eta \right)\left(1-jkR-\frac{1}{9}{k}^2{R}^2\right)\right]}\mathrm{d}\omega$$

(12)

A phased-array system using 3 MHz ultrasound probe operating in pulse-echo mode was used to track the displacement of the sphere. Using the measured response and the theory (i.e., Eq. (12)), the shear moduli and viscosities of two different gelatin phantoms were determined as 1.86 ± 0.005 and 4.96 ± 0.025 kPa and 0.61 ± 0.01 and 1.05 ± 0.03 Pa s, respectively.

3. Elasticity/Viscoelasticity Imaging Using Particles Located at Tissue Interface

3.1. Bubble Located at Tissue Interface

Although there many mathematical models for the radial oscillations of a bubble inside the elastic and viscoelastic media exposed to sound (e.g., [28,29,30,31]) and some studies on the elasticity imaging based on the use of a bubble located within the elastic and viscoelastic media (e.g., [19,20,34,35]), there are a few studies on the responses of the bubbles located at the medium interfaces [24,46,47,48]. By assuming the elastic medium as isotropic, homogeneous, incompressible, and non-viscous, Koruk and Choi [46] proposed models for the static and dynamic displacements of a bubble located at the soft medium interface in response to an external force in 2018. The model for an external static force $f_{\mathrm{e}}=f_0$ and a soft elastic medium is given by [46]:

$$u=\frac{f_0}{2\pi GR\left[1-{\left(1-\frac{u}{R}\right)}^3\right]}$$

(13)

For a dynamic force $f_{\mathrm{e}}$ and a soft elastic medium, they obtained the following equation of motion [46]:

$$\underset{\_}{u}+\underset{\_}{\dot{u}}+\frac{1}{6}\underset{\_}{\ddot{u}}+\frac{1}{18}\underset{\_}{\stackrel{⃛}{u}}={\underset{\_}{f}}_{\mathrm{e}}$$

(14a)

where

$$\underset{\_}{u}=\frac{u\left(t\right)}{R}$$

(14b)

$$\underset{\_}{u}+\underset{\_}{\dot{u}}+\frac{1}{6}\underset{\_}{\ddot{u}}+\frac{1}{18}\underset{\_}{\stackrel{⃛}{u}}={\underset{\_}{f}}_{\mathrm{e}}$$

(14c)

$$\underset{\_}{t}=\frac{t\sqrt{G/\rho }}{R}$$

(14d)

where ${\theta }_{\mathrm{u}}$ is the angle corresponding to the displacement $u$ (see Figure 2b). The solution of Eq. (14) provides the dynamic force-displacement relation for a bubble located at the soft elastic medium interface. The solution of Eq. (14) can be determined numerically. Readers can refer to the reference [46] to see the details for the derivations of the models (i.e., Eqs. (13) and (14)). Later, Koruk and Choi [47], by assuming the medium as isotropic, homogeneous, and incompressible, proposed a mathematical model for the dynamic response of a bubble located at the soft viscoelastic medium interface in response to an external force in 2019. The frequency-domain equation that relates the displacement of the bubble located at the soft viscoelastic medium interface for a dynamic force can be expressed as:

$$2\pi \left(G-j\omega \eta \right)RU\left(1-{cos}^3{\theta }_{\mathrm{u}}\right)\left(1-jkR-\frac{1}{6}{k}^2{R}^2+\frac{1}{18}j{k}^3{R}^3\right)=F_{\mathrm{e}}$$

(15)

In conclusion, the dynamic displacement of the bubble located at the soft viscoelastic medium interface in response to a rectangular pulse input can be calculated using the inverse Fourier transformation [47,49]:

$$u\left(t\right)=-\frac{j{f}_0}{12{\pi }^{2}R\left\{1-{\left[1-\frac{u\left(t\right)}{R}\right]}^3\right\}}\underset{-\infty }{\overset{\infty }{\int }}\frac{\left(e^{j\omega \tau }-1\right){\left(3-jkR\right)e}^{-j\omega t}}{\omega \left(G-j\omega \eta \right)\left(1-jkR-\frac{1}{6}{k}^2{R}^2+\frac{1}{18}j{k}^3{R}^3\right)}\mathrm{d}\omega$$

(16)

Readers can refer to the references [46,47] to see the details for the derivation of the model (i.e., Eqs. (15) and (16)). By dividing the time duration into N points (e.g., N = 1000), and the calculations can be repeated over the entire period of interest to find the dynamic response of the bubble numerically using Eq. (16).

Koruk et al. [23] emitted ultrasound pulses across a tunnel of a gelatin material filled with a cloud of microbubbles. Ultrasound stimulated the microbubble cloud to move in the direction of wave propagation, press against the distal surface, and cause deformations. Shear waves propagated away from the palpation site with a velocity that was used to determine the Young's moduli of the gelatin phantoms ($v_{\mathrm{s}}=\frac{E}{3\rho }$ where $\rho$ is the gel density). Hence, they determined E = 0.46 ± 0.06 kPa for a 2.5% gelatin phantom and E = 1.54 ± 0.32 kPa for a 5% gelation phantom. Bezer et al. [24] administered individual microbubbles into a wall-less hydrogel channel (diameter: 25–100 µm) and exposed to an acoustic pulse (centre frequency: 1 MHz, pulse length: 10 ms, peak-rarefactional pressures: 0.6–1.0 MPa). Using high-speed microscopy, each microbubble was tracked as it pushed against the hydrogel wall. Using the mathematical model explained above (i.e., Eq. 16) and the measured response of each microbubble located at the gel interface, they evaluated the Young’s moduli (2-8.7 kPa) and the viscosities (0.12-0.20 Pa s) of the hydrogels.

In the model presented above (i.e., Eqs. (15) and (16)), the soft viscoelastic medium was modelled using the Kelvin-Voigt model. In this paper, we extended the response of a bubble located at the soft viscoelastic medium interface and generalized the relations so that, in addition to the Kelvin-Voigt model, any viscoelastic rheological model such as the Maxwell, generalized Maxwell, standard linear solid and Kelvin-Voigt fractional derivative models can be used to simulate the medium. In conclusion, the generalized model for the response of a bubble located at the soft viscoelastic medium interface can be represented as:

$$u\left(t\right)=-\frac{j{f}_0}{12{\pi }^{2}R\left\{1-{\left[1-\frac{u\left(t\right)}{R}\right]}^3\right\}}\underset{-\infty }{\overset{\infty }{\int }}\frac{\left(e^{j\omega \tau }-1\right){\left(3-jkR\right)e}^{-j\omega t}}{\omega G^{*}\left(1-jkR-\frac{1}{6}{k}^2{R}^2+\frac{1}{18}j{k}^3{R}^3\right)}\mathrm{d}\omega$$

(17)

Readers may refer to the reference [22] for the expressions of $G^{*}\left(\omega \right)=G^{\text{'}}\left(\omega \right)-j{G}^{\text{'}\text{'}}\left(\omega \right)$ for different rheological models.

3.2. Sphere Located at Tissue Interface

The displacement of a rigid sphere located at the tissue interface exposed to a static force can be calculated using the well-known Hertz Model [49]. However, the Hertz model can produce accurate estimates when the ratio of the sphere displacement to the sphere radius is less than 0.1 [50]. Some modified Hertz models have been proposed for more accurate estimates for static loading (e.g., [51,52,53,54]). Koruk [55] proposed a Hertz model for small and large sphere displacements and a Poisson’s ratio of v = 0.45 given by:

$$f_0={\left[\frac{1-0.1\frac{u}{R}}{1+0.5\left(1-0.1\frac{u}{R}\right)}\right]}^{3/2}\frac{4E\sqrt{R}}{3\left(1-v^2\right)}{u}^{3/2} \mathrm{f}\mathrm{o}\mathrm{r} v=0.45$$

(2)

Koruk [56] recently proposed an improved Hertz model based on the finite element analyses for small and large sphere displacements and different values of the Poisson’s ratio of the medium given below:

$$f_0={\left\{\frac{1-\left(v-0.35\right)\frac{u}{R}}{1+0.5\left[1-\left(v-0.35\right)\frac{u}{R}\right]}\right\}}^{3/2}\frac{4E\sqrt{R}}{3\left(1-v^2\right)}{u}^{3/2}$$

(3)

Among the many investigated models, this model was shown to be more accurate model to accurately estimate the sphere displacement for small and large sphere displacements (i.e., $\frac{u}{R}=0-0.6$) and practical values of the medium Poisson’s ratio (i.e., v = 0.30 – 0.49).

Using the force-displacement relation in Eq. (18) and considering the damping of the oscillations of the sphere due to the radiation of shear waves, Koruk [55] obtained an analytical model for the dynamic response of a rigid sphere located at the elastic medium interface. However, as stated before, the force-displacement relation in Eq. (19) was shown to produce more accurate estimates. Therefore, the analytical model for the dynamic response of a rigid sphere located at the elastic medium interface in response to a rectangular pulse input in the reference [55] was updated using the force-displacement relation in Eq. (19) and the solution is presented in this paper. Hence, the solution is obtained as:

$$u\left(t\right)=\frac{f}{k}-\frac{f}{k\sqrt{1-{\zeta }^2}}{e}^{-\zeta {\omega }_{\mathrm{n}}t}\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{\mathrm{d}}t-\phi \right) \mathrm{for} 0\le t\le \tau$$

(20a)

$$u\left(t\right)=\frac{f{e}^{-\zeta {\omega }_{\mathrm{n}}t}}{k\sqrt{1-{\zeta }^2}}\left\{e^{\zeta {\omega }_{\mathrm{n}}\tau }\mathrm{c}\mathrm{o}\mathrm{s}\left[{\omega }_{\mathrm{d}}\left(t-\tau \right)-\phi \right]-\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{\mathrm{d}}t-\phi \right)\right\} \mathrm{for} t>\tau$$

(20b)

where

$$k=\left[1-\left(v-0.35\right)\frac{\delta }{R}\right]1.5{\left(\frac{4E^{*}\sqrt{R}}{3}\right)}^{2/3}{f}_0^{1/3}$$

(20c)

$$c=\frac{1}{2}\left(0.5+\frac{\delta }{R}\right)\left(\sqrt{\frac{\rho }{G}}R\right)\left[1-\left(v-0.35\right)\frac{\delta }{R}\right]1.5f_0^{1/3}{\left(\frac{4E^{*}\sqrt{R}}{3}\right)}^{2/3}$$

(20d)

$$m=\frac{1}{3}\pi R^3\left(4{\rho }_{\mathrm{s}}+\frac{\delta }{R}\rho \right)$$

(20e)

$$f=\left\{1+0.5\left[1-\left(v-0.35\right)\frac{\delta }{R}\right]\right\}{f}_0$$

(20f)

$${\omega }_{\mathrm{n}}=\sqrt{\frac{k}{m}}$$

(20g)

$${\omega }_{\mathrm{d}}={\omega }_{\mathrm{n}}\sqrt{1-{\zeta }^2}$$

(20h)

$$\zeta =\frac{c}{c_{\mathrm{c}\mathrm{r}}}=\frac{c}{2\sqrt{km}}$$

(20i)

$$\phi ={\mathrm{t}\mathrm{a}\mathrm{n}}^{-1}\frac{\zeta }{\sqrt{1-{\zeta }^2}}$$

(20j)

$$\delta ={\left(\frac{3f_0}{4E^{*}\sqrt{R}}\right)}^{2/3}$$

(20k)

here, $E^{\mathrm{*}}$ is the reduced Young’s modulus for a system consisting of sphere located at a flat medium interface with a Young’s modulus $E$ and Poisson’s ratio $\nu$. For a rigid sphere, the reduced Young’s modulus becomes $E^{\mathrm{*}}=E/\left(1-{\nu }^2\right)$. Koruk et al. [57,58] pushed a sphere located at the gelatin-phantom interface using a magnetic force. Using high-speed camera, the sphere was tracked as it pushed against the gelatin phantom. For example, using the mathematical model explained above (i.e., Eq. (20)) and the measured response of the sphere located at the gelatin-phantom interface, the shear moduli of some gelatin phantoms were determined to be 1.5-2.8 kPa [57]. In addition, by using an equivalent viscous damping ratio for the soft sample in the model (i.e., Eq. (20)), the equivalent viscous damping ratios of the gelatin phantoms were determined to be 0.12-0.38 [57].

Later, Koruk [56], by using the force-displacement relation in Eq. (19) and considering the medium viscosity via the Kelvin-Voigt model, updated the model for the dynamic response of a rigid sphere located at the medium interface in the reference [55]. Hence, the dynamic response of the sphere located at the viscoelastic medium interface in response to a rectangular pulse input in the frequency domain was obtained as:

$$U=\frac{\left\{1+0.5\left[1-\left(v-0.35\right)\frac{\delta }{R}\right]\right\}\left[-\frac{j{f}_0}{\omega }\left(e^{j\omega \tau }-1\right)\right]}{\frac{1}{3}\pi R^3\left(4{\rho }_{\mathrm{s}}+\frac{\delta }{R}\rho \right)\left(-{\omega }^2\right)+\left[\frac{1}{2}\left(0.5+\frac{\delta }{R}\right)\left(\sqrt{\frac{\rho }{G-j\omega \eta }}R\right)\left(-j\omega \right)+1\right]\left[1-\left(v-0.35\right)\frac{\delta }{R}\right]1.5f_0^{1/3}{\left(\frac{4{\tilde{E}}^{*}\sqrt{R}}{3}\right)}^{2/3}}$$

(21)

where ${\tilde{E}}^{*}=2\left(G-j\omega \eta \right)\left(1+\nu \right)/\left(1-{\nu }^2\right)$ for a homogeneous isotropic material and a non-deformable sphere. The time-domain response of the sphere located at the viscoelastic interface was obtained using the inverse Fourier transform as following [56]:

$$u\left(t\right)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\frac{\left\{1+0.5\left[1-\left(v-0.35\right)\frac{\delta }{R}\right]\right\}\left[-\frac{j{f}_0}{\omega }\left(e^{j\omega \tau }-1\right)\right]e^{-j\omega t}}{\frac{1}{3}\pi R^3\left(4{\rho }_{\mathrm{s}}+\frac{\delta }{R}\rho \right)\left(-{\omega }^2\right)+\left[\frac{1}{2}\left(0.5+\frac{\delta }{R}\right)\left(\sqrt{\frac{\rho }{G-j\omega \eta }}R\right)\left(-j\omega \right)+1\right]\left[1-\left(v-0.35\right)\frac{\delta }{R}\right]1.5f_0^{1/3}{\left(\frac{4{\tilde{E}}^{*}\sqrt{R}}{3}\right)}^{2/3}}\mathrm{d}\omega$$

(22)

In addition to the elastic properties of the medium and size of the sphere, this model includes the corrected models for the inertia force due to the medium involved in motion, the inertia force of the sphere, the damping due to the oscillations of the sphere due to the radiation of shear waves, and force-displacement relation that is valid for small and large sphere displacements and practical values of the medium Poisson’s ratio. Therefore, the model given above (i.e., Eqs. (21) and (22)) is considered as the most comprehensive mathematical model for the response of a sphere located at the medium interface in the literature. It should be noted that the model presented above (i.e., Eqs. (21) and (22)) can be extended for the response of a soft sphere located at the medium interface by using the reduced Young’s modulus ${\tilde{E}}^{*}={\left[\frac{1-v_1^2}{E_1}+\frac{1-v_2^2}{E_2}\right]}^{-1}$ where $E_1$ and $E_2$ and $v_1$ and $v_2$ show the Young’s moduli and Poisson’s ratios of the medium and deformable sphere, respectively. Furthermore, the model can be extended for a sphere located at a spherical interface by replacing $R$ with the relative radius $R^{*}={\left[\frac{1}{R_1}+\frac{1}{R_2}\right]}^{-1}$where $R_1$ and $R_2$ show the radii of the spherical interface and sphere, respectively.

It should be noted that the particles such as nondeformable spherical objects are commonly used in conventional indentation tests [59,60] and AFM or atomic force microcopy [61] to identify the material properties. However, a sphere is attached to a rigid holder in conventional indentation tests and to a flexible beam in AFM. Usually, only the Young’s modulus of material can be determined via conventional indentation tests while an AFM system is quite complicated and expensive. Furthermore, these systems are used for material identification in laboratory and can not be used for in vivo elasticity imaging. For example, Huth et al. [62] determined the Young’s moduli of some hydrogels using AFM indentation tests and Hertz model. Orikasa et al. [63] developed an indentation-based testing platform for measuring ultra-soft materials’ stiffness at meso/macro-scales. However, the model proposed by Koruk [56] can be used to identify the shear or Young’s modulus and viscosity of the tissue. Even, this model can be used to evaluate the density and Poisson’s ratio of the tissue.

Koruk et al. [64] has recently been updated the mathematical model in [56] for a trapezoidal pulsed force. Furthermore, Koruk and Pouliopoulos [65] have recently extended the mathematical model for a spherical object located at elastic and viscoelastic medium interface in [55,56] for different practical loading cases. The mathematical models in [65] simulate the change in the frequency of oscillations of the spherical object located at medium interface with the applied load. It was shown that the frequencies of oscillations were different in the loading and unloading phases, and the frequency of oscillations decreased with decreasing external load [65]. In this study, the shear modulus and viscosity of a gelatin sample with a density of 1133 kg/m³ were determined to be 6800 Pa and 6.0 Pa⋅s, respectively [65].

4. Evaluation of Elasticity/Viscoelasticity Imaging Using Particles Located within Tissue and at Tissue Interface

4.1. Evaluation of Mathematical Models

Here, we compared the most updated mathematical models for small particles located within the tissue and at the tissue interface and evaluated the models. The radius of the bubbles used in practical ultrasound applications are a few micrometers (e.g., [48,66]). The radius of the bubbles and spheres used for material characterization in laboratory mostly range from tens of micrometers to a few millimeters (e.g., [21,38]). Therefore, here the results for spherical particles with radii of $R=3$ and 50 $\mu$m were presented. The shear modulus of tissue mostly ranges from a few hundred Pascal to ten thousand of Pascal, while it is a few thousand Pascal for most tissues (e.g., [36,67]). For example, the shear modulus, density, and viscosity of the liver are around 2000 Pa [67], 1000 kg/m³ [68], and 0.5 Pa$\bullet$s [69], respectively. Here, the analyses were performed for physiologically relevant materials (i.e., G = 2000–6000 Pa, ρ = 1000 kg/m³, and η= 0.01–1.6  Pa$\bullet$s). It should be noted that, for all the models presented above, scripts were written using Matlab (MathWorks, Natick, MA, USA) to find the responses of the small particles located within the tissue and at the tissue interface.

The static force required for a specific normalized displacement (i.e., $u/R$) for a bubble inside the medium (Eq. (1)), a sphere inside the medium (Eq. (6)), a bubble located at the medium interface (Eq. (13)) and a sphere located at the medium interface (Eq. (19)) are presented in Figure 3a ($R=3$ $\mu$m) and 3c ($R=50$ $\mu$m). The stiffnesses for the four cases are plotted as a function of the normalized displacement in Figure 3b ($R=3$ $\mu$m) and 3d ($R=50$ $\mu$m). It is seen that the stiffness of the bubble and sphere inside the medium are constant, while the stiffness changes nonlinearly with the displacement for the particles located at the medium interface. This stiffness for a normalized displacement of $\frac{u}{R}=0.5$ is 0.05, 0.08, 0.19 and 0.28 N/m when $R=3$ $\mu$m and 0.81, 1.37, 3.14 and 4.71 N/m when $R=50$ $\mu$m for the sphere located at the medium interface, bubble located at the medium interface, bubble inside the medium and sphere inside the medium, respectively. As the stiffness of the system increases, the force required for a specific displacement increases. It is seen that the value of the force for a normalized displacement of $\frac{u}{R}=0.5$ is 0.07, 0.12, 0.28 and 0.42 $\mu$N when $R=3$ $\mu$m and 20.35, 34.36, 78.54 and 117.8 $\mu$N when $R=50$ $\mu$m for the sphere located at the medium interface, bubble located at the medium interface, bubble inside the medium and sphere inside the medium, respectively. This clearly shows the order of the external force needed to be applied for a specific particle displacement for different systems.

The dynamic responses for a bubble inside the medium (Eq. (4)), a sphere inside the medium (Eq. (9)), a bubble located at the medium interface (Eq. (16)) and a sphere located at the medium interface (Eq. (22)) for $R=3$ $\mu$m (a and b) and $R=50$ $\mu$m (c and d) and two different medium viscosities ($\eta =0.02$ and 0.2 Pa$\bullet$s) are presented in Figure 4. Among the four small particles with $R=3$ $\mu$m, only the sphere located at the medium interface oscillates when $\eta =0.02$ Pa$\bullet$s. However, this small sphere does not have any oscillations when the viscosity is increased to $\eta =$0.2 Pa$\bullet$s. On the other hand, when the particle size is increased (i.e., $R=50$ $\mu$m), both the sphere inside the medium and the sphere located at the medium interface oscillate for the viscosity $\eta =0.02$ Pa$\bullet$s. It is seen that this larger sphere located at the medium interface still oscillates for the viscosity $\eta =0.2$ Pa$\bullet$s. As expected, the period of oscillations of the particle increases (or the frequency of oscillation decreases) as the size of the particle increases. For example, the frequencies of oscillations are 38168 and 2825 Hz for the $3$ and $50$ $\mu$m spheres, respectively, when $\eta =0.02$ Pa$\bullet$s. The viscosity has small effect on the frequency of oscillation. For example, the frequency of oscillations decreases from 2825 to 2674 Hz for the $50$ $\mu$m sphere located at the medium interface when the medium viscosity increased from $\eta =0.02$ to 0.2 Pa$\bullet$s. The time needed for the particle to reach the steady state increases as medium viscosity increases. However, the viscosity does not change the steady-state displacements of the particles.

The effects of medium shear modulus and viscosity on the dynamic responses of a bubble inside the medium (Eq. (4)), a sphere inside the medium (Eq. (9)), a bubble located at the medium interface (Eq. (16)) and a sphere located at the medium interface (Eq. (22)) for $R=50$ $\mu$m are presented in Figure 5 and Figure 6, respectively. It is seen that the displacements of the particles decrease as the medium shear modulus increases. The particles react faster and the time to reach the steady-state decreases as the medium shear modulus increases and medium viscosity decreases. The period of oscillations decreases (or the frequency of oscillation increases) as the medium shear modulus increases. Although, the displacements of the particles at a specific time before steady-state decreases with increasing medium viscosity, the viscosity does not change the steady-state displacements of the particles. It is seen that the mathematical models for all the particles can properly simulate the effects of medium shear modulus and viscosity on the dynamic responses of the particles.

In addition to the time-domain data, the frequency-domain data (spectra) can be used to identify the tissue properties. The FFTs or Fast Fourier Transforms of the response of the sphere located at the medium interface for different medium shear moduli, densities and viscosities are presented in Figure 7. The spectra clearly show that the displacement of the particle decreases and its frequency of oscillations increases with increasing medium shear modulus. The displacement of the particle and its frequency of oscillations decrease with increasing medium density. The displacement of the particle considerably decreases, and its frequency of oscillations slightly decreases with increasing medium viscosity.

In addition to the part of the curve during loading, the part of the curve corresponding to unloading can be used to identify the tissue properties. The external force function and the dynamic responses for a bubble inside the medium (Eq. (4)), a sphere inside the medium (Eq. (9)), a bubble located at the medium interface (Eq. (16)) and a sphere located at the medium interface (Eq. (22)) for $R=3$ $\mu$m (a and c) and $R=50$ $\mu$m (b and d) and two different medium viscosities ($\eta =0.8$ and 1.6 Pa$\bullet$s) are presented in Figure 8. It is observed that as the viscosity of the medium increases, the response of the particle slows down during loading and unloading. The part of the curve corresponding to loading or unloading can be used to identify medium viscosity. The response rate is lowest for the bubble located at the medium interface and the response rate is quite high for the sphere located at the medium interface.

The displacement $u$ presented so far shows the displacement of the tip point (the tip contact point between the particle and medium). It should be noted the displacement $u$ is mostly used for elasticity and viscoelastic imaging purpose. However, if it is needed, the radial $u_r\left(r,\theta \right)$ and polar $u_{\theta }\left(r,\theta \right)$ displacement components and radial ${\sigma }_r\left(r,\theta \right)$ and polar ${\sigma }_{\theta }\left(r,\theta \right)$ stress components at every point in the medium can be determined (see Figure 2b). Readers may refer to the references [19,33] for the corresponding displacement and stress expressions for a bubble and sphere inside the medium. Readers may refer to the references [33,46] for a bubble located at the medium interface and to the reference [33] for the sphere located at the medium interface for the corresponding displacement and stress expressions.

4.2. Evaluation of Tissue Identification Systems

Some of the elasticity and viscoelasticity imaging studies based on measurements and the mathematical models presented above are summarized in Table 1. It is seen that ultrasonic excitation and monitoring is mostly used to identify material properties, though magnetic and mechanical excitation and optical and MRI imaging can be used for elasticity and viscoelastic imaging based on the use of small particles located within the tissue and at the tissue interface (see Figure 2a). It is seen that there are still no applications for identification of tissue properties based on the use of small particles located within the tissue and at the tissue interface for in vivo. One of the reasons for this is that, although the particles inside the tissue have been used for the last 70 years, this technique use a high-powered laser to create the bubble inside the medium which limits its application to shallow targets and requires local destruction of the material or it is needed to locate a bubble or rigid sphere inside the tissue in a way (e.g., [20,43]), which may not be permissible in human tissue. However, thanks to the mathematical models recently proposed for the particles located at the tissue interfaces [46,47,55,56,64,65], it is believed that tissue identification for in vivo can be possible in future. We believe that tissue identification based on the use of small particles located within the tissue and at the tissue interface have potential to change the field, as they can provide the determination of accurate and local tissue elastic properties as well as density, Poisson’s ratio and viscosity, thanks to the developed sophisticated mathematical models. As different from conventional elasticity imaging that map elastic properties of tissue, some other properties such as viscosity of the tissue can be identified, we suggest the use of the term viscoelasticity imaging.

As mentioned before, because there is a need to locate the particles inside the medium in a way for the use of mathematical models for the particles inside the medium, these models are even difficult to be used for the identification of tissue-mimicking materials or ex vivo in laboratory. However, as it is straightforward to locate a particle at the medium interface and this does not alternate material properties, the mathematical models for the particles located at the medium interface can be easily used for the identification of material properties of tissue-mimicking materials and ex vivo in laboratory. Because, in addition to the elastic properties of the medium and size of the sphere, the model for the sphere located at the medium interface includes the corrected models for the inertia force due to the medium involved in motion, the inertia force of the sphere, the damping due to the oscillations of the sphere due to the radiation of shear waves, and force-displacement relation that is valid for small and large sphere displacements and practical values of the medium Poisson’s ratio and its experimental setup is very straightforward, this system seems very promising for tissue identification purposes.

The mathematical models presented in this paper can be used to identify material properties in macroscopic or microscopic scales. If only the identification of the elasticity (or the Young’s modulus) or shear modulus of tissue is required (i.e., elasticity imaging), it is only needed to measure the displacement of the particle exposed to a static force or the steady-state displacement of the particle exposed to a dynamic force. Using the corresponding mathematical models for a static external force (i.e., Eqs. (1), (6), (13 or (19)), the elastic properties of the tissue can be identified. However, in addition to the elasticity or shear modulus, if the identification of the density and/or viscosity of tissue is required (i.e., viscoelasticity imaging), there is a need to measure the dynamic response of the particle and to use the corresponding mathematical models for dynamic loading (i.e., Eqs. (4), (9), (16) or (22)). By performing curve fitting using experimental data and the mathematical models, the elasticity modulus, density, and viscosity of tissue can be identified. Alternatively, the elastic properties can be determined from the steady-state displacement and the viscosity can be identified from the part of the measured curve corresponding to loading or unloading and the density of the tissue can be determined by matching the measured and theoretical oscillating frequency of the particle. Furthermore, as it may be difficult to know the magnitude of the applied force in practice, without needing the amplitude of the applied force, the elasticity or shear modulus and/or the density of tissue can be identified by matching the measured and theoretical oscillating frequency of the particle and the viscosity of tissue can be identified from the part of the measured curve corresponding to loading or unloading. In addition, the frequency of oscillation of the particle and damping of tissue can be determined using the spectrum of the time-domain data and the modal analysis methods such as half-power, circle-fit or line-fit (e.g., [70,71]).

5. Concluding Remarks

We investigated the elasticity/viscoelasticity imaging based on the use of small particles located within the tissue and at the tissue interface exposed to static and dynamic external loads in this paper. We compared the mathematical models for the small particles exposed to static and dynamic external loads and evaluated the elasticity and viscoelasticity imaging methods based on the use of small particles located within the tissue and at the tissue interface. Some concluding remarks are listed below.

• The acoustic radiation force-based elasticity imaging methods that remotely palpate tissue produce a complicated displacement field. The displacement field depends on the transducer focusing characteristics and target geometry. For these absorption-based techniques, the focal volume affects elasticity reconstruction. For elasticity imaging based on the use of small particles located within the tissue and at the tissue interface, particle displacement is independent of the focal zone of the transducer and is affected only by the viscoelastic properties of the local tissue and simple and accurate elasticity reconstruction can be performed. As sophisticated mathematical models for a particle located within the tissue and at the tissue interface can be developed, in addition to elastic properties, density and viscosity of tissue can be determined in the elasticity imaging based on the use of small particles.
• Because the radiation force applied to the particle located within the tissue and at the tissue interface is generally much greater than that arising in purely absorbing media, lower acoustic power can be used in moderately attenuating media for the elasticity imaging techniques based on the use of small particles.
• For the particles located inside the tissue, there is a need to use a high-powered laser to generate the bubble inside the medium which limits its application to shallow targets and requires local destruction of the material or it is needed to locate a bubble or rigid sphere inside the tissue in a way, which may not be permissible in human tissue. However, thanks to the models recently proposed for the particles located at the tissue interfaces, it is believed that tissue identification for in vivo can be possible in future. Furthermore, thanks to well-developed mathematical models and the ease of their use, the particles located at the medium interface can be conveniently used to identify tissue-mimicking materials and ex vivo in laboratory.

• The stiffness of the bubble and sphere inside the medium is constant, while the stiffness changes nonlinearly with the displacement for the particles located at the medium interfaces. The stiffness, hence, the force required for a specific displacement of the particle, increases from the sphere located at the medium interface to bubble located at the medium interface, bubble inside the medium and sphere inside the medium. The value of force needed for a ratio of 0.5 particle displacement to the radius of the special object and a medium with shear modulus of 5000 Pa is around 70-420 nano Newton and 20-120 micro-Newton for particles with radius of $3$ and $50$ $\mu$m, respectively.

• The period of oscillations of the particle significantly increases (or the frequency of oscillation significantly decreases) as the size of the particle increases. However, even the sphere located at the medium interface does not have any oscillations if the particle size is small and the medium viscosity is considerably high.
• The displacements of the particles decrease as the medium shear modulus increases. The particles react faster and the time to reach to the steady-state decreases as the medium shear modulus increases. The period of oscillations decreases (or the frequency of oscillation increases) as the medium shear modulus increases.
• The displacements of the particles at a specific time before steady-state decreases with increasing medium viscosity. The particles react faster, and the time needed for the particle to reach the steady state decreases as medium viscosity decreases. The viscosity does not change the steady-state displacements of the particles. The frequency of oscillation of the particle interacting with tissue slightly decreases as medium viscosity increases.
• The frequency of oscillations of the particle interacting with tissue deceases with increasing medium density.
• If only the identification of the elasticity (or Young’s modulus) or shear modulus of tissue is required (i.e., elasticity imaging), it is only needed to measure the displacement of the particle exposed to a static force or the steady-state displacement of the particle exposed to a dynamic force. Using the corresponding mathematical models for a static external force, the elastic properties of the tissue can be identified. However, in addition to the elasticity or shear modulus, if the identification of the density and/or viscosity of tissue is required (i.e., viscoelasticity imaging), there is a need to measure the dynamic response of the particle and to use the corresponding mathematical models for dynamic loading. By performing curve fitting using experimental data and the mathematical models, the Young’s or shear modulus, density and viscosity of tissue can be identified.
• Alternatively, the elastic properties can be determined from the steady-state displacement and the viscosity can be identified from the part of the measured curve corresponding to loading or unloading and the density of the tissue can be determined by matching the measured and theoretical oscillating frequency of the particle. Furthermore, as it may be difficult to know the magnitude of the applied force in practice, without needing the amplitude of the applied force, the Young’s modulus or shear modulus and/or the density of tissue can be identified by matching the measured and theoretical oscillating frequency of the particle and the viscosity of tissue can be identified from the part of the measured curve corresponding to loading or unloading.
• Although ultrasonic excitation and monitoring is mostly used to identify material properties based on the use of small particles, magnetic and mechanical excitation and optical and MRI imaging can be used for elasticity and viscoelastic imaging based on the use of small particles.

We believe that tissue identification based on the use of small particles located within the tissue and at the tissue interface have potential to change the field, as they can provide the determination of accurate and local tissue elastic properties as well as density, Poisson’s ratio, and viscosity, thanks to the developed sophisticated mathematical models. As different from conventional elasticity imaging that map elastic properties of tissue, some other properties such as density, Poisson’s ratio and viscosity of tissue can be identified using small particles located within the tissue and at the tissue interface, we suggest the use of the term viscoelasticity imaging.