3D Ultrasound Shear Wave Imaging Simulation of Heterogeneous Tissue

1. Introduction

Shear wave imaging has been used to quantitatively measure tissue stiffness noninvasively in vivo. Shear waves can be generated in soft tissues by either direct mechanical vibration [1] or transient ultrasound (US) radiation force excitation [2-4]. The shear wave speed in tissues is directly related to the underlying tissue shear modulus [3]. 3D shear wave imaging has been utilized to access tissue elasticity using two-dimensional matrix array US transducers [5, 6]. The shear wave propagations are captured from the radiation force excitation axis in three dimensions [5]. With continuous development of technology, shear wave imaging has been investigated for diagnosing some diseases, mainly for detecting breast cancer [4] and liver cirrhosis [7, 8]. To reconstruct the local tissue shear modulus, the 2D algebraic Helmholtz inversion (AHI) equation has been widely used [3, 8-11].

The 2D AHI equation reconstructs the shear modulus based on assumptions that planar shear wave propagates in lateral direction (across beam direction) and the shear wave diffraction in the elevational direction (out of plane direction) is not substantial [3]. While the 2D AHI equation can assess the elasticity of homogeneous tissues, its application is limited in assessing heterogeneous tissues. In heterogeneous tissues, high attenuation and refraction of shear waves at the boundary of different tissues have been reported [2], which reduces signal to noise ratio (SNR). The shear modulus reconstruction using 2D AHI equation happens to be sensitive to noise [2]. In addition, the shear wave diffraction in elevational direction is substantial in heterogeneous tissues, and should not be ignored in computations. The 2D AHI equation only uses shear wave propagation in lateral direction and this limits shear modulus reconstruction of heterogeneous tissues. Moreover, the shear wave is no longer planar in heterogeneous tissues due to wave refraction at the boundary of different tissues. Therefore, a 3D shear modulus reconstruction algorithm is necessary to increase accuracy of shear modulus calculations in heterogeneous tissues.

In this study, we developed a 3D AHI equation by expanding the 2D AHI equation. We evaluated the results from 2D and 3D AHI calculations by comparing the shear modulus in two types of heterogeneous tissues: early stage of cancer and advanced cancer. First, a 3D finite element (FE) tissue model, with a spherical inclusion (D: 8.0 mm, hard) embedded in a soft tissue block of 3.0×3.0×3.0 cm, was constructed. Second, shear wave propagations in the 3D FE tissue model were simulated. Finally, the shear modulus in the 3D FE tissue model was reconstructed using the 3D AHI equation. The shear modulus by the 3D AHI equation was compared with the shear modulus reconstructed by the 2D AHI in a cross-section of the FE model.

 

2. Materials and Methods

2.1 3D algebraic helmholtz inversion equation

The tissue elasticity map can be reconstructed using an axial displacements (along the US beam direction) field induced by the shear wave propagation [2, 3]. The tissue is assumed to be elastic, isotropic, and homogeneous. After application of a radiation force, the shear wave propagation equation is presented as:

In Eq. (1), μ is the shear modulus, u is the displacement vector, t is the time of shear wave propagation, and ρ is the density of the medium. In the US shear wave imaging, the shear wave propagation in lateral direction x is obtained from the displacements field of axial direction z. The shear wave propagation equation for the z-component is:

In the 2D AHI Eq. (2), it is assumed that the shear wave diffraction is not significant in the elevation direction y, and the shear modulus is computed from image area corresponding to a (x, z) plane. Thus, y-component is not included in the equation such as in the following:

As the displacement is a function of time, the local shear modulus can be written in the Fourier domain, in which F represents the Fourier transform in the time domain, and N is the number of frequencies taken into account:

In the 3D AHI equation, the diffraction of shear wave in the y direction is incorporated into the equation, and the shear modulus is computed using a (x, y, z) volume. The shear modulus can be computed as the following:

2.2 Finite element modeling

A 3D FE tissue model was developed using a commercially available software package (Comsol 3.5a Multiphysics software, COMSOL, MA, USA.). Fig. 1(a) shows the schematic of the 3D FE tissue model (3.0×3.0 ×3.0 cm) containing a hard spherical inclusion (D: 8.0 mm). The material properties of the surrounding tissue were set at 60 kPa shear modulus, isotropic, 1050 kg/m3 density. The Poisson’s ratio was input as ν=0.495 by assuming incompressible tissue. Note that the shear modulus was computed from input Young’s modulus and Poisson’s ratio. To simulate early stage of cancer, for the inclusion, the material properties were set at 150 kPa shear modulus, isotropic, 1050 kg/m3 density and Poisson’s ratio ν=0.495. To simulate advanced cancer, for the inclusion, the shear modulus of inclusion was set at 600 kPa with the same density and Poisson’s ratio.

Fig. 1.Schematic diagram of 3D finite element (FE) tissue modeling and simulation; (a) Schematic of the 3D FE tissue model containing a hard spherical inclusion; (b) A diagram of cross-section A located in the middle of the 3D FE model in Fig. 1(a).

The model was meshed using tetrahedral elements. Fig. 1(b) presents a diagram of cross-section A located in the middle of the 3D FE model in Fig. 1(a). A half section of the spherical inclusion was located at 1.5 cm from the bottom and 0.7 cm from the right boundary of the FE model. To simulate the shear wave propagation, body forces of a circular shape with a 1.0 mm diameter were applied downward through the entire depth of the FE model for 180 μs to create angled planar shear waves which are practically observed in US shear wave imaging. The bottom boundaries of the mesh were fixed to prevent any bulk motion from the local body force excitation for shear wave generation. The z displacements of the shear waves over time were extracted over the entire mesh of the FE model. The shear modulus of the 3D FE model was computed from the z displacements using the 3D AHI equation. For the shear modulus reconstruction, a 3D median filter, 3.0×3.0×3.0 mm in kernel size, was applied to smooth image.

The shear modulus of cross-section A reconstructed by the 3D AHI equation was compared with the 2D AHI equation. For the 2D AHI equation, we first extracted z displacements of the shear waves over time in cross-section A over the entire mesh. Then, the shear modulus of cross-section A was computed from the z displacements using the 2D AHI equation. For the shear modulus reconstruction, a 2D median filter, 3.0×3.0 mm in kernel size, was applied. Each shear modulus for the inclusion and the surrounding tissue in cross-section A was spatially averaged within the area of the inclusion size.

 

3. Results

Fig. 2(a) shows the axial displacement color map of the early stage of cancer model captured from the FE simulation. The shear wave propagation from the center of the FE model going outward was clearly visible on the top surface of the FE model. Fig. 2(b) presents the axial displacement color maps in cross-section A which were captured at 2.0, 2.6, and 3.2 ms after the shear waves were generated by body force. The red color represents the shear wave front, and the location of the hard inclusion is marked by a black dashed circle. The shear wave front advanced as soon as it was propagated into the inclusion area. Note that the dynamic range of the color map adjusted to clearly display the shear wave front.

Fig. 2.(a) The axial displacement color map of the FE tissue model captured from the FE simulation; (b) The axial displacement color maps in cross-section A which were captured at 2.0, 2.6, and 3.2 ms after the shear waves were generated by body force

Figs. 3(a) and 3(c) present the shear modulus map of cross-section A calculated from the displacement field of shear waves using the 2D AHI equation. Figs. 3(b) and 3(d) illustrate calculations obtained by the 3D AHI equation.

Fig. 3.A shear modulus map of cross-section A calculated from the displacement field of shear waves using the 2D algebraic Helmholtz inversion (AHI) equation from the early stage of cancer; (a) and advanced cancer (c); A shear modulus map of cross-section A obtained by the 3D AHI equation from the early stage of cancer; (b) and advanced cancer (d)

Figs. 3(a) and 3(b) were obtained from the early stage of cancer model, and Fig. 3(c) and 3(d) were computed from the advanced cancer model. The inclusion and surrounding tissue were marked by black dashed circles. The higher shear modulus estimated within the inclusion indicates that the embedded inclusion is stiffer than the surrounding tissue. The shear wave excitation region was marked by a black dashed rectangle as there was no shear wave propagation in this area. The results of these figures show that the shear modulus of the surrounding tissue obtained by 2D AHI equation was similar to that obtained by the 3D AHI equation. On the other hand, the shear modulus of inclusion reconstructed by 3D AHI equation exhibited higher values than that obtained by 2D AHI equation. Each shear modulus for the inclusion and the surrounding tissue was spatially averaged within the area of the inclusion size as depicted in Fig. 3.

Figs. 4(a) and 4(b) shows the shear modulus of inclusion and surrounding tissue computed by the 2D AHI and 3D AHI equations. Fig. 4(a) presents the shear modulus values obtained from the early stage of cancer model. Fig. 4(b) displays results computed from the advanced cancer model. The average shear modulus of the inclusion and surrounding tissue are marked by an open red circle and an open blue square, respectively. The average and standard deviation of shear modulus of inclusion and surrounding tissue are summarized in Table 1.

Fig. 4.The shear modulus of inclusion and surrounding tissue computed by 2D algebraic Helmholtz inversion (AHI) and 3D AHI equations for early stage of cancer (a) and advanced cancer (b)

Table 1.Shear modulus of inclusion and surrounding tissue computed by 2D AHI and 3D AHI equations

For the early stage of cancer model, input shear modulus values, which are marked by black dashed lines, were set at 150 kPa for the inclusion and 60 kPa for the surrounding tissue. The average shear modulus of the inclusion and surrounding tissue obtained by the 2D AHI equation was 130 kPa and 59 kPa, respectively. The average shear modulus of the inclusion and surrounding tissue obtained by the 3D AHI equation was 150 kPa and 60 kPa, respectively For the advanced cancer model, input shear modulus values were set at 600 kPa for the inclusion and 60 kPa for the surrounding tissue. The average shear modulus of the inclusion and surrounding tissue obtained by the 2D AHI equation was 576 kPa and 59 kPa, respectively. The average shear modulus of the inclusion and surrounding tissue obtained by the 3D AHI equation was 600 kPa and 60 kPa, respectively.

 

4. Discussion

The 3D AHI equation exhibited more accurate shear modulus results with less variation for the 3D FE tissue model than the 2D AHI equation (Figs. 4(a) and 4(b)). For the early stage of cancer model, the average shear modulus of inclusion reconstructed by the 3D AHI equation was the same as the input value of 150 kPa. On the other hand, the average shear modulus of inclusion computed by the 2D AHI equation was 130 kPa and it was lower than the input value. In addition, the standard deviation of the shear modulus in the inclusion obtained by the 2D AHI equation was about 2 times higher (19 kPa compared to 10 kPa) than the shear modulus reconstructed by the 3D AHI equation (Fig. 4(a)). For the advanced cancer model, the average shear modulus of inclusion reconstructed by the 3D AHI equation was the same as the input value of 600 kPa. On the other hand, the average shear modulus of inclusion computed by the 2D AHI equation was 576 kPa and it was lower than the input value (Fig. 4(b)). The underestimation and high variation of shear modulus computed by the 2D AHI equation may have been caused by the attenuation of shear wave in the x direction (lateral direction). The shear wave attenuation occurred about two times more within the inclusion than the surrounding tissue as shown in Fig. 2(b). This attenuation reduced SNR and critically affected the shear modulus reconstruction [2]. Although shear waves in the x direction were attenuated, some shear waves propagated to the y direction (elevational direction) in the inclusion. The 3D AHI equation reconstructed shear modulus using shear wave propagations in both x and y directions, and could compute the shear modulus in the inclusion more accurately than the 2D AHI equation.

A distorted inclusion shape was observed in the reconstructed shear modulus map obtained by the 2D AHI equation for both early stage of cancer model and advanced cancer model (Figs. 3(a) and 3(c)). The reconstructed inclusion shape was shorter in the x direction than the input inclusion shape. On the other hand, the inclusion shape in the shear modulus map reconstructed by the 3D AHI equation was close to the input inclusion shape as shown in Figs. 3(b) and 3(d). The distortion of the inclusion shape may have been caused by a curved shear wave front in the x direction [12]. The 2D AHI equation computed the shear modulus by assuming a planar shear wave front; however, after the shear wave advanced into the inclusion, the wave form was no longer a planar wave, as seen in Fig. 2(b). The 3D AHI equation uses shear wave propagations in both x and y directions, and computations were less affected by the shear wave front shape than in 2D AHI computations.

In this study, the shear modulus of inclusion was selected 150 kPa for the early stage of cancer model and 600 kPa for the advanced cancer model, but the same shear modulus of surrounding tissue was selected 60 kPa for both models. We only increased the shear modulus of inclusion since pathological changes such as progress of malignant cancer result in substantial increase of cancer tissue stiffness in contrast to the surrounding healthy tissues [13-15].

The error of shear modulus of inclusion was remarkably reduced by using 3D AHI equation for the early stage of cancer model as shown in Fig. 4(a). On the other hand, for the advanced cancer model, the error of shear modulus of inclusion obtained by the 2D AHI equation was not noticeable as shown in Fig. 4(b). For the early stage of cancer model, the shear wave speed was close between inside inclusion and surrounding tissue, and the shear modulus computation may have been susceptible to the shear wave attenuation and shear wave diffraction. For the advanced cancer model, the shear wave speed was much faster inside inclusion than surrounding tissue, and shear modulus computation within the inclusion may have been less affected by the shear wave attenuation and shear wave diffraction.

In this study, speckle noise was not considered and z displacements (axial displacements) field of shear wave were obtained by FE simulation. Practical US images contain speckle noise [16], and the axial displacements of shear wave are usually estimated using 1-D autocorrelation function [17]. The speckle noise affects the axial displacements estimation, but it would not impact on shear wave propagation. This speckle noise effect on shear modulus computation may need to be further investigated for practical US shear wave imaging.

Tissue viscosity effects for shear wave speed were not considered within the scope of the study. The 3D AHI equation assumes that the effect of viscosity for shear wave propagation is minute and the shear modulus was calculated by ignoring viscosity effects for shear wave propagation. However, the proportional relationship between shear wave speeds and wave frequencies for viscous tissues has been previously reported [18], and the viscosity effect for shear modulus reconstruction may need further investigation.

The 3D AHI equation substantially improved accuracy of shear modulus computation for early stage of cancer model. The shear wave imaging is not sensitive enough to detect small stiffness changes, especially in early stage of pathological changes [19]. The 3D AHI equation would be able to improve the accuracy of US shear wave imaging for detecting early stage of liver fibrosis [1], liver cirrhosis [7, 8], and breast cancer [4].

 

5. Conclusion

This study evaluated the accuracy of the 3D AHI equation in assessing the shear modulus in heterogeneous tissues using in silico with a FE simulation. The 3D AHI equation produced more accurate shear modulus results with less variation for the 3D FE heterogeneous tissue model than the original 2D AHI equation. Therefore, we expect that the 3D AHI equation would be able to improve the accuracy of diagnosis for 3D US shear wave imaging detecting early stage of liver fibrosis and breast cancer. The accuracy of 3D AHI equation should be further validated in future ex vivo and in vivo experiments for clinical application.