Paper The following article is Open access

3D elastic tensor imaging in weakly transversely isotropic soft tissues

, , , , and

Published 25 July 2018 © 2018 Institute of Physics and Engineering in Medicine
, , Citation M Correia et al 2018 Phys. Med. Biol. 63 155005 DOI 10.1088/1361-6560/aacfaf

0031-9155/63/15/155005

Abstract

We present herein 3D elastic tensor imaging (3D ETI), an ultrasound-based volumetric imaging technique to provide quantitative volumetric mapping of tissue elastic properties in weakly elastic anistropic media. The technique relies on (1) 4D ultrafast shear wave elastography (SWE) at very high volume rate (e.g.  >  8000 Hz, depending only on the imaging depth), (2) a volumetric estimation of shear wave velocity using the eikonal equation and (3) a generalized 3D elastic tensor-based approach. 3D ETI was first evaluated using numerical simulations in homogeneous isotropic and transverse isotropic media. Results showed that 3D ETI can accurately assess tissue stiffness and tissue anisotropy in weakly transversely isotropic media (elastic fractional anisotropy coefficient  <  0.34). Experimental feasibility was shown in vitro in a transverse isotropic phantom. Quantification of the elastic properties by 3D ETI was in good agreement with 2D SWE results performed at different orientations using a clinical ultrafast ultrasound scanner. 3D ETI has the potential to provide a volumetric quantitative map of tissue elastic properties in weakly transversely isotropic soft tissues within less than 20 ms of acquisition for the entire imaged volume.

Export citation and abstract BibTeX RIS

Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

1. Introduction

Elastography techniques either based on magnetic resonance imaging (MRI) (Fowlkes et al 1995, Muthupillai et al 1995, Manduca et al 2001, Rump et al 2007) or on ultrasound (Sarvazyan et al 1998, Catheline et al 1999, Fatemi and Greenleaf 1999, Nightingale et al 2002, Sandrin et al 2003, Bercoff et al 2004, Song et al 2012) have been proposed to evaluate non-invasively the stiffness of soft tissues, in order to provide clinically relevant information (Glaser et al 2012, Gennisson et al 2013). Yet, these techniques are based on the same general principle, i.e. the generation of mechanical shear waves within the tissue by a variety of methods and the tracking of their propagation to assess tissue stiffness. Ultrasound based techniques present several advantages compared to MRI techniques, including portability and real-time capability at a lower cost. However, MRI techniques allow the possibility to acquire the shear wave propagation in 3D.

Shear wave elastography (SWE) has been extensively developed over the past decade since it can provide a quantitative mapping of soft tissues in real-time, and it has been implemented on a clinical ultrasound scanner (Bercoff et al 2004). This technique is based on a push-beam creation within the tissue by acoustic radiation force, and on ultrafast ultrasound imaging to track the induced shear wave propagation. The shear wave velocity is directly related to the elastic properties of the tissue, such as, the shear and Young's modulus. Although clinical SWE is commonly applied to soft tissues with isotropic elastic properties, such as the liver or the breast, the feasibility of SWE to non-invasively characterize tissue anisotropy was also shown in vivo in skeletal muscles (Gennisson et al 2010), arteries (Shcherbakova et al 2017) and in the myocardium (Couade et al 2011, Lee et al 2012). SWE can indeed quantify the shear wave velocities parallel (i.e. higher shear wave velocity) and perpendicular (i.e. lower shear wave velocity) to the tissue fibers, by setting the ultrasound probe parallel or perpendicular to the fibers direction, respectively.

Additionally, elastic tensor imaging (ETI) was proposed (Lee et al 2012) to generalize a tensor-based approach of SWE imaging. By using several SWE acquisitions along different directions, a least-squares inversion can be performed to estimate tissue fibers orientation and stiffness. This technique was initially based on the rotation of a 2D ultrasound probe to create shear waves propagating along different directions. This approach, however, remains challenging to perform in a clinical setting as it requires acquisition repetition during the probe rotation.

Recently these technical challenges were overcome by the development of 4D Ultrafast imaging. Provost et al (2014) introduced a 4D ultrafast ultrasound system that can perform volumetric imaging at 5000 Hz. Gennisson et al (2015) showed that 4D ultrafast shear wave elastography can be performed using a 2D matrix-array probe using milliseconds acquisitions. Therefore, shear wave propagation can be induced and imaged in three-dimensions, and tissue stiffness can be estimated in a full volume.

In this study, we introduce 3D elastic tensor imaging (3D-ETI) to quantify non-invasively the elastic properties of weakly transversely isotropic soft tissues in three-dimensions. This technique is based on 4D ultrafast shear wave elastography imaging to access the elastic tensor and therefore provide a quantitative mapping of the anisotropic elastic properties. Using the Eikonal equation, we developed a method to quantify the elastic properties under the assumption of weak anisotropic and weak heterogeneous medium. The objectives of this study were therefore: (1) to evaluate and validate the method in isotropic and transverse isotropic media through numerical simulations, and (2) to demonstrate the experimental feasibility in a transverse isotropic model.

2. Material and methods

2.1. 3D Elastic Tensor in isotropic and transverse isotropic media

We assume the propagation of small amplitude shear wave that can be described by the linear elasticity theory. The wave propagation in elastic media is governed by the wave equation, that can be solved leading to the Christoffel's equation (Royer and Dieulesaint 2000):

Equation (1)

where ${{G}_{ik}}$ is the Christoffel's matrix given by ${{G}_{ik}}={{C}_{ijkl}}{{n}_{j}}{{n}_{i}}$ (with ${{n}_{i}}$ the direction vector of the wave propagation), ${{C}_{ijkl}}$ , the elastic tensor, $V$ the phase wave velocity, ${{\delta }_{ik}},$ the identity matrix, $\rho $ the density of the medium and ${{U}_{k}}$ the particle displacement vector.

Considering that in SWE, shear waves remotely generated by the acoustic radiation force are polarized along the ultrasound probe axial direction (i.e. z-direction in a Cartesian coordinate system), the particle displacement vector is given by ${{U}_{k}}=[0~\,0\,~1]$ , with $k=1,2,3$ . Also, considering an orthotropic or transverse isotropic medium, the representative matrix of the Christoffel stiffness tensor, using the Voigt notation, can be given by:

Equation (2)

In shear wave propagation, considering SH-propagation, the stiffness tensor can be further simplified to:

Equation (3)

By developing equation (1), considering (3) to the SWE case in transverse isotropic media, and deriving the group velocity from the phase wave velocity and angle, the following relationship can be obtained:

Equation (4)

Where ${{V}_{s}}$ is the group velocity (usually easier to measure in SWE) and $\beta $ the group velocity angle. Further information on linear elasticity theory and the applied mathematical developments, interested reader may refer to appendix.

The equation (4) is valid for any $\beta $ and the stiffness tensor, G, can be solved using the least-squares method for $N$ shear wave propagation measurements. Therefore, equation (4) can be extended to:

Equation (5)

The least-square solution of $ \newcommand{\e}{{\rm e}} {{\left[ \begin{array}{@{}cccccccccccccccccccc@{}} \frac{1}{{{\mu }_{11}}} & \frac{1}{{{\mu }_{22}}} \\\end{array}\!\! \right]}^{T}}$ is obtained by a pseudo-inversion operation. In an isotropic and transverse isotropic medium, the shear wave velocity along and across the tissue fibers, can be further retrieved by:

Equation (6)

where ${{v}_{//}}$ and ${{v}_{\bot }}$ correspond to the shear wave propagation velocities along (i.e. parallel) and across (i.e. perpendicular) the fibers orientation, respectively.

In order to evaluate the degree of elastic anisotropy, we used the elastic fractional anisotropy (FA) defined in a 2D plane (Lee et al 2012). The fractional anisotropy (FA) is a normalized scalar value. A FA value of zero corresponds to an isotropic shear wave propagation (i.e. it is unrestricted in all directions), and a value of one means that shear wave propagation occurs along one specific direction, being completely restricted along all the others:

Equation (7)

To perform 3D elastic tensor imaging, shear wave elastography is used to induce, track and estimate the shear wave velocities parallel and perpendicular to the tissue fibers alignment in three-dimensions.

2.2. 4D ultrafast shear wave elastography

4D ultrafast shear wave elastography relies on the following steps: (1) a spherical delay law was applied in all ultrasound probe transducers to create a push-beam with a cylindrical shape at the focal spot; (2) 4D ultrafast imaging was used to track the shear wave propagation through 3D plane-wave coherent compounding (Provost et al 2014), (figure 1); (3) the Kasai et al autocorrelation method (Kasai et al 1985) was then applied to obtain shear wave propagation cine-loop; (5) a 3D directional filtering approach was proposed to filter the artifacts coming from reflected shear waves; (6) the shear wave velocity was estimated in a full volume through the Eikonal equation application. For more details please refer to Gennisson et al (2015).

Figure 1.

Figure 1. 4D ultrafast shear wave elastography concept using 2D plane-wave transmission. Push-beam definition (A) and Shear wave creation and 4D ultrafast plane-wave imaging (B).

Standard image High-resolution image

The 3D directional filtering approach was proposed, as the extension of the method proposed by Deffieux et al (2012) to the case of 3D cylindrical waves. The method consisted in: (1) converting the shear wave propagation cine-loop from the Cartesian (x, y, z, t) to the cylindrical (r, $\theta ,$ z, t) coordinates system with the center fixed at the push-beam location; (2) applying a Fourier transform along the radius (r) and time (t); (3) keeping, for each $\theta $ direction, the two quadrants of the Fourier space that show ${{k}_{r}}$ and $\omega $ with opposite signs; (4) applying the inverse Fourier Transform; and (5) reconverting the cylindrical to the Cartesian coordinate system.

The Eikonal equation was used to assess the local shear wave velocity by linking the medium phase velocity to the wave phase gradient, and assuming a high-frequency wave propagating in weakly heterogeneous and linear isotropic media. In addition, we assumed that the group velocity and the phase velocity are almost equal under the assumption of low anisotropy and non-dispersive medium. Under these assumptions, the Eikonal equation can be written as (Slawinski 2010):

Equation (8)

Equation (8) allows then for the estimation of the shear wave velocity without prior knowledge of the propagation direction, by calculating the time delays along the Cartesian coordinate system axis. For each direction, the time delays are calculated by the conventional cross-correlation applied on the shear wave profiles for points at a given distance ($dr$ ) and centered on the voxel-of-interest. A Parabolic interpolation and a multiscale approach are also used, by varying the distance between the correlated waveforms within a fixed kernel size, in order to minimize the variance of the time delay estimation.

In order to recover a full volume information of the shear wave propagation pattern in a region-of-interest, N push-beams at a fixed depth were induced in different lateral locations.

To evaluate the performance of 3D ETI we performed the following: (1) numerical simulations of the shear wave propagation using the Green's formalism (Aki and Richards 1980) on isotropic, transverse isotropic (as 1st order approximation), elastic, homogeneous media, and (2) in vitro measurements on transverse isotropic phantom.

To summarize 3D ETI technique, a methodology workflow is presented in figure 2.

Figure 2.

Figure 2. Methods workflow.

Standard image High-resolution image

2.3. Numerical simulation methods

The numerical simulations of the shear wave propagation field, based on temporal Green functions formalism (Aki and Richards 1980), were performed in isotropic and transverse isotropic media (Gennisson et al 2003, Chatelin et al 2015).

2.3.1. Radiation force field simulation

The shear wave propagation field was based on a Green's function formalism (Aki and Richards 1980, 2001). The simulated medium was assumed to be semi-infinite, isotropic or transverse isotropic, elastic and homogeneous. For all the details in Green's function formalism, interested reader may refer to the appendix.

The radiation force field was calculated using the Field II ultrasound software (Jensen and Svendsen 1992, Jensen 1996). The matrix-array transducer used to simulate fived push-beams was designed numerically with the parameters: 32  ×  32 transducer elements, 9.6 mm2 size, 0.3 mm pitch, 5 MHz central-frequency. Push-beams were simulated at five different lateral and elevational locations (i.e. at the x- and y-directions parallel to the surface of the matrix-array probe): at (0, 0) mm, i.e. the center, (3, 3) mm, (−3, 3) mm, (3, −3) mm and (−3, −3) mm, i.e. the four corners of the numerical matrix-array probe. The beam was perpendicular to the probe surface for the central push. The beam was tilted by 9.6° for all the other pushes. For all the five push-beams, the depth focus was fixed at 25 mm (i.e. z-direction and axial direction) and the duration was 125 µs.

Using the described parameters, the radiation force field generated was simulated following the methodology described in Chatelin et al (2015).

2.3.2. Shear wave propagation velocity field

For each push-beam (i.e. five in total, different x- and y-locations), the displacement vector was calculated. A displacement field-of-view of (9.6  ×  9.6  ×  30) mm, a density of 1038 kg m−3 and Poisson ratio 0.49 were considered as parameters. The displacement velocity vector was then calculated by the first-derivative of the described displacement field, considering an imaging volume rate of 8 KHz and step time of 125 µs.

In our study, the shear wave propagation velocity fields were computed for different stiffness and fractional anisotropy coefficients to evaluate the performance of the 3D ETI technique on different media. A summary of the stiffness and fractional anisotropy parameters used is presented in table 1.

Table 1. Parameters used in the numerical simulations.

Mediuma 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
Shear wave velocity parallel to fibers (m s−1) 5.0 5.0 5.0 5.0 5.0 3.5 3.5 3.5 3.5 3.5 2.0 2.0 2.0 2.0 2.0
FA coefficient 0.000 0.074 0.156 0.246 0.343 0.000 0.074 0.156 0.246 0.343 0.000 0.074 0.156 0.246 0.343

aEach Medium is composed of five different simulations, where the push-beams were simulated at different x- and y-locations (i.e. central and four corners of the matrix-array simulated probe).

2.4. Experimental methods

2.4.1. Transverse isotropic phantom experiments

Millon et al (2006) introduced an approach to develop orthotropic hydrogel mimicking soft tissues that were found suitable for elastography in anisotropic medium (Chatelin et al 2014). In this method, mechanical anisotropy is induced by stretching the physical crosslinks of polyvinyl alcohol (PVA) cryogel during freezing/thaw cycles along a preferential direction. Based on this approach, we developed a new sample constituted by 5% PVA solution (molecular weight 89 000–98 000, 99+% hydrolyzed, Sigma-Aldrich, St Louis, US), 1% cellulose (20 µm in diameter, S3504 Sigmacell, Sigma-Aldrich, St Louis, US) and 94% distilled water. First, the sample was transferred into a parallelepiped PVC mold (100  ×  60  ×  50 mm3) and then underwent two isotropic freezing/thaw cycles (i.e. approximately 20 °C during 12 h followed by  −18 °C during 12 h) followed by one anisotropic freezing/thaw cycle (i.e. stretching the sample to 80% of the initial length through a customized static tensile-test apparatus). The cross-link was subsequently induced in the polymer along the stretching direction. The in vitro setup, including the transverse isotropic phantom developed, is presented in figure 3(A).

Figure 3.

Figure 3. The experimental setup with transverse isotropic phantom and 2D matrix-array probe (A) and push-beam locations scheme (B).

Standard image High-resolution image

In order to assess the transverse isotropic phantom stiffness and compare the results with the ones obtained by 3D ETI technique, experiments using a 2D ultrafast ultrasound scanner (Aixplorer scanner, Supersonic Imagine, Aix-en-Provence, France) were performed, using a clinical shear wave elastography mode, with a 1D linear array (SL10-2, 6 MHz central frequency, Supersonic Imagine, Aix-en-Provence, France). In this case, the shear wave velocity was assessed in two orientations, i.e. parallel and perpendicular to the fibers orientations of the transverse isotropic phantom.

2.4.2. 4D ultrafast imaging hardware

A customized, programmable, 1024-channel ultrafast ultrasound system (Provost et al 2014) was used to drive a 32-by-32 matrix-array probe centered at 8 MHz with a 90% bandwidth at  −6 dB, a 0.3 mm pitch and a 0.3 mm element size (Vermon, Tours, France), (figure 3(A)). The system was composed of four Aixplorer systems (Supersonic Imagine, Aix-en-Provence, France), each one with 256-transmit and 128-receive channels units. The four systems were assembled and synchronized, defining a device containing 1024-channels in transmission and 512-channels in receive. Since in receive the channels were multiplexed to 1 of 2 transducers elements, each emission was repeated twice to synthetize 1024-channels. Specifically, data from two consecutive emissions were concatenated to create synthetically one RF data set, corresponding approximately to an acquisition in which all transducers have received simultaneously. Consequently, the effective volume rate is half the pulse-repetition frequency (${\rm PRF}$ ).

2.4.3. 4D ultrafast shear wave elastography imaging sequence

The sequence consisted of a push-beam followed by 4D ultrafast plane-wave imaging. To perform a push-beam by acoustic radiation force, a spherical law was programmed in transmission using 784 transducer elements located in the center of the matrix-array probe at 5 MHz, at 25 mm depth and during 100 µs.

As shown by Bercoff et al (2004) and Gennisson et al (2015), the shear wave propagation information cannot be assessed in the immediate region around the push-beam location. To access a full elastic map of the region-of-interest in the media, we can however rely on another push-beam performed at a different location. In our study, five push-beams were transmitted at different locations (i.e. at the x- and y-directions of the matrix-array probe), (see figure 3(B)). Specifically, the push-beams were arbitrarily performed at (0, 0) mm, i.e. the center, (3, 3) mm, (−3, 3) mm, (3, −3) mm and (−3, −3) mm, i.e. the four corners of the matrix-array probe.

Followed by each push-beam, 4D ultrafast imaging using 2D plane-wave insonification was performed to track the consequent shear wave propagation. Each plane-wave was defined by a pair of angles between the normal vector of the plane wave and the x- and y-directions of the matrix-array probe., In this study, only one 2D plane wave was used, with a (0°, 0°) angles. This choice was made to achieve the largest possible field of view (i.e. 9.6  ×  9.6 mm2), in which the 2D plane wave propagates, and the highest possible volume rate frequency (i.e. ${\rm PRF}$ divided by $N$ ). The effective ${\rm PRF}$ and the volume rate were both equal to 8050 Hz for a 25 mm imaging depth.

Each acquisition consisting in one push-beam followed by the acquisition of 40 ultrafast imaging volumes lasted 5.1 ms. For each push-beam location, four consecutive acquisitions were performed that were temporally averaged, in order to improve the signal-to-noise ratio of the tracked shear wave propagation. The total acquisition time (i.e five push-beams and ultrafast imaging, each repeated four times) was then 102 ms.

2.4.4. Post-processing analysis

A 3D delay-and-sum beamforming using spatial coherent compounding was performed, as described by Provost et al (2014), to reconstruct the ultrasound volumes (i.e. beamformed IQ demodulated data). In in vitro and in vivo results, volumes of (9.6  ×  9.6  ×  25) mm3 and (19.2  ×  19.2  ×  25) mm3 were reconstructed, respectively.

Tissue velocity estimates were obtained through the Kasai autocorrelation method (Kasai et al 1985). Prior to phase estimation and the 3D directional filtering, the data was filtered using a spatiotemporal clutter filter (Demene et al 2015). The final tissue velocity movie was then composed of 38 imaging volumes of (9.6  ×  9.6) mm2 and 25 mm depth. Tensor estimation was then performed using the method described in section 2.2.

Finally, the 3D volume rendering of the shear wave propagation and a tissue fiber tractography representation (i.e. the 3D elastic tensor representation) were performed using the volren and illuminated streamlines functions of the Amira software (Visualization Sciences Group, USA), respectively.

All beamforming and elastic tensor assessment processing were performed offline using a graphical processing unit (Tesla K40, NVidia, USA) within a Matlab (2014b, The MathWorks Inc., USA) interface. The post-processing time was approximately 10 min. This duration included all 4D ultrafast data transfers and calculations, i.e. the beamforming, directional filtering, shear wave velocity vector components estimation, fiber tractography.

3. Results

3.1. Numerical simulations

Figure 4 presents the shear velocities estimated from the numerical simulations data.

Figure 4.

Figure 4. Numerical simulation results of the application of 3D ETI to estimate shear wave velocities parallel and perpendicular to the fibers media. In (A) the estimated parallel shear wave propagation velocities (i.e. parallel to the medium fibers) compared to the theoretical velocities (i.e. 5.0, 3.5 and 2.0 m s−1) and in function of the fractional anisotropic coefficients. In (B) the estimated perpendicular shear wave propagation velocities (i.e. perpendicular to the medium fibers) for the parallel shear wave propagation velocities of 5.0, 3.5 and 2.0 m s−1, and in function of the fractional anisotropic coefficients. In (C) and (D) the correlation and Bland–Altman plots between the theoretical and estimated parallel shear wave velocities. In (E) and (F) the correlation and Bland–Altman plots between the theoretical and estimated perpendicular shear wave velocities.

Standard image High-resolution image

The averaged estimated parallel shear wave velocities (i.e. parallel to the medium fibers, figure 4(A)) were in good agreement with the theoretical velocities (i.e. 5.0, 3.5 and 2.0 m s−1).

Perpendicular shear wave velocities (i.e. perpendicular to the medium fibers, figure 4(B)) were also in a very good agreement with the theoretical values, but only for weak anisotropy (FA  <  0.15). In contrast, for higher fractional anisotropy, the estimation errors increased up to 20% (FA  =  0.34).

Figures 4(C) and (D) present the correlation and Bland-Altman analysis between the theoretical and estimated shear wave velocities parallel to the fibers alignment. The strong correlation corroborated by the Bland-Altman plot, showed that parallel shear velocities were well estimated.

Figures 4(E) and (F) present the correlation and Bland-Altman analysis between the theoretical and estimated shear wave velocities perpendicular to the fibers alignment. It shows, that most of the perpendicular shear velocities were well estimated, but highest errors were found in moderately anisotropic medium.

The fiber orientation was also evaluated and the respective standard-deviation estimation is also presented. These results are presented in table 2 and a 3D ETI fiber tractography case is presented in figure 5.

Table 2. Fiber alignment orientation results.

Medium 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
Shear wave velocity parallel to fibers (m s−1) 5.0 5.0 5.0 5.0 5.0 3.5 3.5 3.5 3.5 3.5 2.0 2.0 2.0 2.0 2.0
Estimated angle  +  stda (°) 177.9  ±  23.9 175.9  ±  26.5 172.7  ±  28.2 169.6  ±  36.3 178.7  ±  11.6 174.7  ±  12.2 173.4  ±  10.9 172.6  ±  19.5 181.8  ±  6.6 182.5  ±  5.6 182.7  ±  6.2 183.5  ±  6.3

aStd stands for the standard-deviation of the spatial-averaged orientation estimation.

Figure 5.

Figure 5. 3D-ETI fiber tractography results on numerical simulations of medium #8.

Standard image High-resolution image

3.2. In vitro experiments: transverse isotropic media

The shear wave propagation results are presented in figure 6 in a top and side views. Five pushes were induced at different spatial locations (i.e. different locations in the x- and y-directions of the matrix array probe). The shear waves showed an elliptical propagation pattern, demonstrating an anisotropic behavior that is directly relate d to the fibers alignment orientation.

Figure 6.

Figure 6. Shear wave propagation results in the transverse isotropic phantom for the five different pushes in a specific time-frame. The first row corresponds to the shear wave propagation view seen from the top (parallel to the matrix-array plane). The second row corresponds to a lateral shear wave propagation view.

Standard image High-resolution image

The 3D elastic tensor was reconstructed in a region of interest of 2.0  ×  2.0  ×  5 mm3 in the center of the imaging volume by using the information of the 4 shear waves propagating from the corners (i.e. pushes located at 2, 3, 4 and 5 in figure 6). The reconstruction included: the 3D directional filtering, the 3D time-of-flight calculation, and the 3D elastic tensor estimation algorithms. The averaged parallel shear wave velocity calculated was $3.35\pm 0.34\,{\rm m}\,{{{\rm s}}^{-1}}$ , the averaged perpendicular shear wave velocity calculated was $1.90\pm 0.17\,{\rm m}\,{{{\rm s}}^{-1}}$ , and the averaged orientation of the fibers estimated was $170.7\pm 10.3{}^\circ $ . The elliptical orientation of shear wave propagation of the central push (i.e. push 1 of figure 6) was also assessed and was found as $163~\pm 5{}^\circ $ . The fractional anisotropy coefficient was calculated and it was $F{{A}_{3{\rm D}\,{\rm ETI}}}=0.38\pm 0.08$ .

The shear wave velocities in two orthogonal 2D imaging planes, i.e. parallel and perpendicular to the fibers alignment, were also assessed by the Aixplorer ultrasound scanner. These results are presented in figure 7. In a region of interest of (9  ×  9) mm2 and we found $3.10\pm 0.10\,{\rm m}\,{{{\rm s}}^{-1}}$ parallel to fibers and $2.20\pm 0.10$ m s−1 perpendicular to the fibers. The fractional anisotropy coefficient was calculated and it was $F{{A}_{2{\rm D}\,{\rm SWE}}}=0.24\pm 0.05$ .

Figure 7.

Figure 7. Aixplorer measurements results in two different planes of the transverse isotropic phantom. In (A) the results in a region-of-interest parallel to the medium fiber direction. In (B) the results in a region-of-interest perpendicular to the medium fiber direction.

Standard image High-resolution image

Finally, the elastic fiber tractography approach was also applied and the results are shown in figure 8. Through this representation approach, it could be noted that the fibers are mainly orientated in the same direction and orientation.

Figure 8.

Figure 8. 3D-ETI fiber tractography results of transverse isotropic phantom.

Standard image High-resolution image

4. Discussion

In this study, we have developed a new approach, 3D elastic tensor imaging, for the assessment of the Elastic Tensor in three-dimensions that is based on 4D ultrafast shear wave elastography imaging. This technique was first evaluated in a numerical study and in an in vitro study, i.e. in a transverse isotropic phantom.

Our method was first evaluated using synthetic data of shear wave propagation in transverse isotropic medium generated by a numerical model. The technique was therefore evaluated in isotropic and transverse isotropic media with different stiffness (i.e. different shear wave propagation velocities) and different anisotropy degrees (i.e. different fractional anisotropy coefficients). These results showed, that 3D ETI technique could accurately estimate the shear wave velocities parallel to the fibers alignment. Likewise, 3D ETI technique could accurately estimate shear wave velocities perpendicular to the fibers alignment, with estimation error increasing for higher FA coefficients (i.e. up to 20% error for 5 m s−1 and up to 13.7% for 3.5 m s−1). Also, the standard-deviation of the estimation shear wave velocities values increased up 1.3 m s−1. These numerical results demonstrated that 3D ETI can accurately estimate anisotropic shear wave velocities in weak transversely isotropic media (FA  <  0.35), but may provide inaccurate measurements at higher FA (i.e. estimation errors higher than 20%). This was however expected since the time-of-flight algorithm based on the Eikonal equation (8), which was used to compute the shear velocity, assumes low anisotropy and low heterogeneity. In what concerns the fibers angle estimations, results have shown that 3D ETI can accurately access the fibers angles direction in weak transversely isotropic media. The angle direction error was found to be  ±  6° (i.e. for FA  =  0.34 and shear wave velocities parallel to the fibers of 2 m s−1), and increased up to 36° (i.e. for FA  =  0.34 and shear wave velocities parallel to the fibers of 5 m s−1).

Based on previous 2D ultrafast shear wave elastography studies (Gennisson et al 2010), 5 push-beams were expected to be enough to sample the entire field-of-view in 3D. Indeed, three pushes were used in 2D to cover an imaging field-of-view of 20 mm, with a distance between those pushes of 6.7 mm. In this study, the distance between the pushes was smaller (i.e. 4.2 mm between the central and the corner pushes, and 6 mm between the corner pushes). Moreover, these settings were used (Gennisson et al 2015) to perform 3D shear wave elastography in vivo in the breast, and 5 pushes were sufficient to reconstruct a similar volume.

3D elastic tensor imaging was also evaluated in vitro in a transverse isotropic phantom. The shear wave velocity vector was assessed in three dimensions and the fiber tractography representation was implemented. The results showed a good agreement between 3D ETI and the control measurements performed using a conventional and clinical SWE ultrafast ultrasound scanner results. Since, the control measurement was within the numerical simulation estimations range limits, we considered that the in vitro 3D ETI allowed to accurately estimate the tissue stiffness and elastic anisotropy (i.e. estimation error approximately 20%). Also, the fiber tractography representation showed that the estimated fibers alignment orientation was as expected, with approximately 11% of estimation error.

The matrix probe offers the possibility to generate complex push beams, however the quantification of the anisotropy requires to generate isotropic push beams. In this study we created isotropic push beams by transmitting the push on a circular window (f/D ratio equal to 3, with a 5 MHz transmission).These parameters were chosen as a compromise between the transmitted energy and frequency, and the isotropy of the spherical push beam.

In terms of imaging parameters, we used one 2D plane-wave transmission, at 0° angle, to achieve the largest field of view and the highest possible volume rate. This choice was yielded to effectively track the shear wave propagation in our field of view (i.e. (9.6  ×  9.6  ×  25) mm3). A higher number and tilt angles could have been used to increase the imaging quality by coherently compounding the receive signals (Provost et al 2014), but the volume rate would have been reduced (i.e. the volume rate is defined by the effective PRF divided by the number of transmitted 2D plane-waves). Moreover, the choice of the angle set may affect the imaging isotropy. Therefore, in order to maintain it at least three independent angles should be used for the coherent compounding. With the current setup this would result in 2500 volumes s−1 which is quite low to perform shear wave imaging. Nevertheless, using electronics with 1024 channel in receive instead of 512 would allow to reach 5000 volumes s−1 even with 3 angles which should be enough to perform SWE.

One limitation of the technique, is the use of 2D plane-wave insonification which limits the field of view. By using 2D plane-waves transmission, the imaging field of view is reduced to a portion of a parallelepiped that is equal to the matrix-array size. Therefore, to track shear wave propagation in stiff tissues in vivo (i.e. shear wave velocities higher than 4 m s−1) a higher volume rate or higher field of view should be achieved. This last restriction could be overcome by the use of 2D diverging-waves or with higher number of electronic channels. Though, in 2D diverging-wave imaging, the transmitted ultrasound energy is considerable lower compared with the plane-wave approach, and in this case the energy decreases with the inverse of the depth square-root (Papadacci et al 2014). Thus, the 3D shear wave elastography implementation using 2D diverging-waves may be more challenging. Concerning the number of electronic channels, in our study we used a custom made prototype of 3D ultrafast imaging scanner with 1024 electronic channels. To the best of our knowledge, only three research systems in the world have currently such a high channel count (Jensen et al 2013, Provost et al 2014, Petrusca et al 2018). We can envision however, that such systems will be increasingly available at the best cost in the near future, with the recent development of high channel count electronics commercially available.

Another limitation of the method, is the assumption of weak heterogeneity in the eikonal equation which implies that the velocity varies spatially at a larger scale than the wavelength. As typical shear wavelength are in the range of a few millimeters to centimeters, our method is not valid for medium with strong heterogeneities such as small nodules. Moreover, we assumed that the anisotropy was weak. We validated our method for fractional anisotropy lower than FA $\sim $ 0.3 which is in the range of typical FA in passive skeletal muscles or in the myocardium. At higher fractional anisotropy, e.g. active muscles, however, a more complex derivation will be required to quantify accurately the elastic properties. Finally, we measured fibers orientation only in the x- and y-direction, parallel to the matrix array probe, when in fact the fibers direction could be also in the z-direction. This limits the application of the technique to imaging settings where the ultrasound beam can be positioned perpendicularly to the fiber direction, e.g. skeletal muscles (Gennisson et al 2010) and myocardium (e.g. in short- and long-axis echographic views), (Villemain et al 2018). In order to overcome this limitation, an approach with steered plane wave could be consider to estimate vector tissue velocities.

3D ETI could indeed be applied in medicine by characterizing non-invasively the elastic properties of skeletal muscles at rest or at low load, such as brachii biceps. Shear wave elastography has been shown to provide important clinical diagnostic information on the musculoskeletal system (Shinohara et al 2010, Koo et al 2014). The precise measurement of the shear modulus in muscles is challenging to perform due to the muscle anisotropy. The probe must be positioned accurately along the fibers direction, which is difficult to ensure during muscle contraction. Reproducibility of the measurement is therefore a real challenge. 3D ETI could potentially offer a solution to improve reproducibility and accuracy of shear wave velocity measurements and provide new diagnostic biomarkers.

Likewise in cardiology, we have shown recently (Villemain et al 2018) that 2D SWE can be performed non-invasively in the passive myocardium (i.e. during diastole) and at two orientations (e.g. in short- and long- axis echographic views). The fractional anisotropy coefficient was then accessed, and it was found to be altered in hypertrophic cardiomyopathy patients. However, the FA quantification was limited by the low number of probe orientations available (i.e. two echographic views). Thus, 3D ETI could potentially be upgraded for cardiology to access simultaneously and accurately the diastolic myocardial stiffness in all directions and to quantify the fractional anisotropy for the diagnostic of cardiomyopathies.

5. Conclusion

3D elastic tensor imaging, a 4D ultrafast shear wave elastography imaging based technique can assess elastic tissue stiffness and anisotropy simultaneously, in three-dimensions, at high volume rate (e.g.  >8000 Hz) and in a few minutes of post-processing. The accuracy of the technique was evaluated numerically in simulated isotropic and transverse isotropic model with different degrees of stiffness and anisotropy. Results showed that 3D ETI can accurately assess tissue stiffness and tissue anisotropy in weakly transversely isotropic media (elastic fractional anisotropy coefficient  <0.34). The in vitro feasibility was also assessed in a transverse isotropic phantom. 3D elastic tensor imaging shows a big potential for in vivo application to assess tissue stiffness and weak elastic anisotropy.

Acknowledgments

This study was supported by the European Research Council under the European Union's Seventh Framework Programme (FP/2007–2013)/ERC Grant Agreement n° 311025 and the ANR-10-IDEX-0001-02 PSL* Research University.

Appendix

A.1. Linear elasticity theory in transverse isotropic media: phase and group velocities

We assume the propagation of small amplitude shear wave that can be described by the linear elasticity theory. The wave propagation in an elastic media is governed by the wave equation:

Equation (A.1)

Here, $\rho $ is the density of the medium, ${{u}_{i}}$ is the displacement vector, $t$ is the time, ${{x}_{i}}$ are the spatial coordinates, and ${{C}_{ijkl}}$ is the stiffness tensor, a fourth-rank tensor. Considering a harmonic plane-wave, equation (A.1) can be solved leading to the Christoffel's equation (Royer and Dieulesaint 2000):

Equation (A.2)

where ${{G}_{ik}}$ is the Christoffel's matrix given by ${{G}_{ik}}={{C}_{ijkl}}{{n}_{j}}{{n}_{i}}$ (with ${{n}_{i}}$ the direction vector of the wave propagation), $V$ the phase wave velocity, ${{\delta }_{ik}},$ the identity matrix. The Christoffel's equation has three independent solutions, corresponding to three different wave-modes. The direction of the particle displacement vector $U$ , or the polarization of each wave-mode, corresponds to an eigenvector of the Christoffel's matrix, while the phase velocity $V$ can be deduced from the associated eigenvalue.

Considering that in SWE, shear waves remotely generated by the acoustic radiation force are polarized along the ultrasound probe axial direction (i.e. z-direction in a Cartesian coordinate system), the particle displacement vector is given by ${{U}_{k}}=[0~\,0\,~1]$ , with $~k=1,2,3$ . Simultaneously, shear waves propagate in the transverse and longitudinal planes (i.e. x- and y-directions in a Cartesian coordinate system) which correspond to a direction vector given by ${{n}_{i}}=[{{n}_{1}}~{{n}_{2}}~0]$ , with $i=1,2,3$ . Thus, the equation (A.2) can be developed into:

Equation (A.3)

with ${{C}_{3111}}=0;{{C}_{3112}}={{C}_{3121}}=0;{{C}_{3211}}=0;{{C}_{3212}}=~{{C}_{3221}}=0;{{C}_{3122}}=0;{{C}_{3222}}=0$ . These six independent components demonstrate that the media can present orthotropic, tetragonal, hexagonal, cubic or isotropic symmetry. Shear wave propagation in equation (A.3) can be rewritten as (Dellinger 1991, Lee et al 2012, Wang et al 2013):

Equation (A.4)

Equation (A.5)

with $\alpha $ the phase velocity angle with respect to the fibers in the propagation plane.

In a nondispersive isotropic medium, the phase and group velocity are equal, but in an anisotropic medium they can differ (figure A1). As a matter of fact, the group velocity can exhibit different values depending on the observation point, while the phase velocity is always a single value dependent on the direction of the transmitted plane-wave (Wang et al 2013). The differences between these two velocities are then dependent of the properties of the medium, such as, viscosity and anisotropy. In SWE, the group velocity (${{V}_{s}}$ ) is usually easier to measure, which can be derived from the phase velocity using equation (A.6), (Dellinger 1991, Wang et al 2013), as:

Equation (A.6)
Figure A1.

Figure A1. Phase (V) and group velocity (Vs) for a plane-wave transmitted in a nondispersive anisotropic medium.

Standard image High-resolution image

The group velocity angle can also be derived from the phase velocity and the phase velocity angle, as given by equation (A.7):

Equation (A.7)

Assuming an orthotropic or transverse isotropic medium, the stiffness tensor is then reduced to:

Equation (A.8)

And, equations (A.4) and (A.5) can be rewritten as:

Equation (A.9)

By developing (A.6), (A.7) and (A.9), the following relationship can be obtained:

Equation (A.10)

A.2. The elastic Green's function in a transversely isotropic medium

The shear wave propagation field was based on a Green's function formalism (Aki and Richards 1980, 2001). The simulated medium was assumed to be semi-infinite, isotropic or transverse isotropic, elastic and homogeneous.

Generally, the stiffness tensor $\mathbf{C}$ links the stress and strain tensors, $\boldsymbol{\sigma }$ and $\boldsymbol{\varepsilon }$ , via the via the generalized Hooke's law, given by:

Equation (A.11)

In the specific case of transverse isotropy, the representative matrix of the Christoffel stiffness tensor was further simplified by using the Voigt notation as followed:

Equation (A.12)

The calculation of the displacement vector from one point source in the temporal domain was expressed using the transversely isotropic elastic Green's formalism of equation (A.11). The terms $~G_{kl}^{(P)}$ , $~G_{kl}^{(SV)}$ , $~G_{kl}^{(SH)}$ , $~G_{kl}^{(SV-SH)}$ and $~G_{kl}^{(P-SV)}$ corresponded to the contribution of the P-wave, SV-wave, SH-wave, near-field term between P-wave and SV-wave and near-field term between P-wave and SH-wave, respectively (Vavryčuk 2001, Gennisson et al 2003, Chatelin et al 2015).

Equation (A.13)

where $t$ is the time, ${\bf x}=\left(x,y,z \right)$ the Cartesian coordinates of specific observation point, ${{F}_{l}}$ the radiation force amplitude along the direction of propagation (as described in section 2.3.1), $\rho $ the density and $\left(k,l \right)$ the direction index numbers of the displacement and the source, respectively. The other parameters of the exact formula for the Green's function are:

Equation (A.14)

The travel times of the compressional P-, shear SV- (i.e. parallel to the fibers) and shear SH- (i.e. orthogonal to the fibers) wave were respectively defined by:

Equation (A.15)

The polarization vectors were defined by:

Equation (A.16)

where $r=\sqrt{x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}$ and $~R=\sqrt{x_{1}^{2}+x_{2}^{2}}$ are the distance from the source point to the observation point and the distance from the focal axis to the observation point, respectively. ${{N}_{m}}={{x}_{m}}/r$ is the unit direction vector from the source point to the observation point. ${{c}_{11}}=\rho v_{P}^{2}$ , ${{c}_{44}}=\rho v_{\parallel }^{2}$ and ${{c}_{66}}=\rho v_{\bot }^{2}$ are the coefficients of the Christoffel stiffness tensor ${\bf C}$ (vP, $v_{\parallel }^{2}$ . and $v_{\bot }^{2}$ correspond to the velocity of the P-wave, SV-wave and SH-wave, respectively).

Please wait… references are loading.