# Modelling the impulse diffraction field of shear waves in transverse isotropic viscoelastic medium

Simon Chatelin, Jean-Luc Gennisson, Miguel Bernal, Mickael Tanter, Mathieu Pernot

► To cite this version:

Simon Chatelin, Jean-Luc Gennisson, Miguel Bernal, Mickael Tanter, Mathieu Pernot. Modelling the impulse diffraction field of shear waves in transverse isotropic viscoelastic medium. Physics in Medicine and Biology, 2015, 60 (9), pp.3639-3654. 10.1088/0031-9155/60/9/3639. hal-02324714

HAL Id: hal-02324714

https://hal.science/hal-02324714

Submitted on 25 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Modelling the impulse diffraction field of shear waves in transverse isotropic viscoelastic medium

To cite this article: Simon Chatelin et al 2015 *Phys. Med. Biol.* **60** 3639

View the [article online](#) for updates and enhancements.

You may also like

- *In vivo* evaluation of the elastic anisotropy of the human Achilles tendon using shear wave dispersion analysis
J Brum, M Bernal, J L Gennisson et al.

- *An in silico* framework to analyze the anisotropic shear wave mechanics in cardiac shear wave elastography
Annette Caenen, Mathieu Pernot, Mathias Peirlinck et al.

- *A finite element model to study the effect of tissue anisotropy on ex vivo arterial shear wave elastography measurements*
D A Shcherbakova, N Debusschere, A Caenen et al.
PAPER • OPEN ACCESS

# Modelling the impulse diffraction field of shear waves in transverse isotropic viscoelastic medium

To cite this article: Simon Chatelin et al 2015 *Phys. Med. Biol.* **60** 3639

View the [article online](#) for updates and enhancements.

You may also like

- *In vivo* evaluation of the elastic anisotropy of the human Achilles tendon using shear wave dispersion analysis
J Brum, M Bernal, J L Gennisson et al.

- *An in silico* framework to analyze the anisotropic shear wave mechanics in cardiac shear wave elastography
Annette Caenen, Mathieu Pernot, Mathias Peirlinck et al.

- *A finite element model to study the effect of tissue anisotropy on ex vivo arterial shear wave elastography measurements*
D A Shcherbakova, N Debusschere, A Caenen et al.
# Modelling the impulse diffraction field of shear waves in transverse isotropic viscoelastic medium

Simon Chatelin, Jean-Luc Gennisson, Miguel Bernal, Mickael Tanter and Mathieu Pernot

Institut Langevin, ESPCI ParisTech, PSL Research University, UMR 7587 CNRS, U979 INSERM, Paris, France

E-mail: simon.chatelin@espci.fr

Received 24 November 2014, revised 3 March 2015
Accepted for publication 3 March 2015
Published 16 April 2015

<figure><img src="image_1.png" /></figure>

CrossMark

## Abstract

The generation of shear waves from an ultrasound focused beam has been developed as a major concept for remote palpation using shear wave elastography (SWE). For muscular diagnostic applications, characteristics of the shear wave profile will strongly depend on characteristics of the transducer as well as the orientation of muscular fibers and the tissue viscoelastic properties. The numerical simulation of shear waves generated from a specific probe in an anisotropic viscoelastic medium is a key issue for further developments of SWE in fibrous soft tissues. In this study we propose a complete numerical tool allowing 3D simulation of a shear wave front in anisotropic viscoelastic media. From the description of an ultrasonic transducer, the shear wave source is simulated by using Field's II software and shear wave propagation described by using the Green's formalism. Finally, the comparison between simulations and experiments are successively performed for both shear wave velocity and dispersion profile in a transverse isotropic hydrogel phantom, *in vivo* forearm muscle and *in vivo* biceps *brachii*.

Keywords: shear wave imaging, viscoelasticity, anisotropy, Green's functions, transverse isotropy, fibrous tissue, elastography

(Some figures may appear in colour only in the online journal)

<figure><img src="image_2.png" /></figure>

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.
# Modelling the impulse diffraction field of shear waves in transverse isotropic viscoelastic medium

Simon Chatelin, Jean-Luc Gennisson, Miguel Bernal, Mickael Tanter and Mathieu Pernot

Institut Langevin, ESPCI ParisTech, PSL Research University, UMR 7587 CNRS, U979 INSERM, Paris, France

E-mail: simon.chatelin@espci.fr

Received 24 November 2014, revised 3 March 2015
Accepted for publication 3 March 2015
Published 16 April 2015

<figure><img src="image_2.png" /></figure>

CrossMark

## Abstract

The generation of shear waves from an ultrasound focused beam has been developed as a major concept for remote palpation using shear wave elastography (SWE). For muscular diagnostic applications, characteristics of the shear wave profile will strongly depend on characteristics of the transducer as well as the orientation of muscular fibers and the tissue viscoelastic properties. The numerical simulation of shear waves generated from a specific probe in an anisotropic viscoelastic medium is a key issue for further developments of SWE in fibrous soft tissues. In this study we propose a complete numerical tool allowing 3D simulation of a shear wave front in anisotropic viscoelastic media. From the description of an ultrasonic transducer, the shear wave source is simulated by using Field's II software and shear wave propagation described by using the Green's formalism. Finally, the comparison between simulations and experiments are successively performed for both shear wave velocity and dispersion profile in a transverse isotropic hydrogel phantom, *in vivo* forearm muscle and *in vivo* biceps *brachii*.

Keywords: shear wave imaging, viscoelasticity, anisotropy, Green's functions, transverse isotropy, fibrous tissue, elastography

(Some figures may appear in colour only in the online journal)

<figure><img src="image_3.png" /></figure>

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

Understanding and simulating shear wave propagation in the biological soft tissues are key features for the development of shear wave elastography (SWE) (Fatemi and Greenleaf 1998, Sarvazyan et al 1998). In SWE, the mechanical properties of the viscoelastic soft tissue are deduced from properties of shear wave propagation, as group or phase velocity, attenuation or dispersion. Each property can be studied differently in order to retrieve tissue mechanical properties. As an example, shear wave dispersion is linked to elasticity and viscosity of tissues and can be characterized using appropriate rheological models (Chen et al 2004, Deffieux et al 2009). Shear wave group velocity would give access to tissue stiffness (Fatemi and Greenleaf 1998, Sarvazyan et al 1998) and shear wave attenuation to viscosity (Catheline et al 2004). These parameters were studied in different organs often considered as isotropic (Palmeri et al 2008, Tanter et al 2008, Sandrin et al 2003, Bavu et al 2011, Cosgrove et al 2011, Deffieux et al 2011, Berg et al 2012). In fibrous media, such as muscle, biological tissues can be considered as transversely isotropic. Many studies with different techniques were performed regarding shear wave velocities as a function of the main tissue axis (Hoffmeister et al 1996, Gennisson et al 2010, Aubry et al 2013).

By combining acoustic radiation force with ultrafast imaging (Tanter et al 2014), supersonic shear imaging (SSI) was proposed a decade ago as a remote generation of low-frequency shear waves in biological soft tissues (Bercoff et al 2004a). The anisotropic shear viscoelasticity was previously investigated experimentally by using SSI technique in *in vivo biceps brachii* by rotating the ultrasonic probe with respect to the main fiber's direction. (Gennisson et al 2010). SSI was also proposed on other strong anisotropic tissues as *ex vivo myocardium* or *in vivo Achille's* tendon (Lee et al 2012, Brum et al 2013) or arteries (Shcherbakova et al 2014). However, understanding wave anisotropy of shear viscosity and elasticity is still difficult. Recently, a theory of shear wave propagation in a soft tissue with anisotropy of both elastic and dissipative properties has been proposed with no fewer than four independent components for these two characteristics (Rudenko and Sarvazyan 2014). Due to the complex dependence of the shear wave propagation in such tissues, there is a need to have models and numerical simulations of shear wave propagation in incompressible anisotropic viscoelastic tissues. Recently, a theoretical model was proposed and compared to finite element simulations to estimate the elastic tensor of an incompressible transverse isotropic material (Rouze et al 2013).

For over half a century, the waves generated by point sources have been described mainly for seismology in homogeneous isotropic and anisotropic elastic media, such as piezoelectric, thermoelastic and poroelastic solids (Buchwald 1959, Aki and Richards 1980, Burridge 1967, Willis 1980, Norris 1994, Wang and Achenbach 1994). The theoretical Green's functions in an infinite transverse isotropic media have been calculated by Vavryčuk (2001). This description was extended some years later to viscoelastic media by including a frequency dependence in the Green's functions (Vavryčuk 2007, 2008). From these studies, velocity and attenuation of waves propagating from a point source in homogeneous media of arbitrary anisotropy and attenuation strength can be calculated using the Green's formalism. At the same time, the use of the Green's functions approach for low-frequency elastic waves through semi-infinite isotropic and homogeneous soft media was proposed for the development of transient elastography in the biological tissues (Catheline et al 2004, Sandrin et al 2004). The extension to transverse isotropic elastic media has been developed to support the evidence of the polarization of low frequency shear for transient elastography in muscular soft tissues (Gennisson et al 2003). The application of the Green's formalism in an Infinite viscous and elastic medium was then developed for analyzing elastic waves induced by acoustic radiation force imaging (Bercoff et al 2004b).
# 1. Introduction

Understanding and simulating shear wave propagation in the biological soft tissues are key features for the development of shear wave elastography (SWE) (Fatemi and Greenleaf 1998, Sarvazyan et al 1998). In SWE, the mechanical properties of the viscoelastic soft tissue are deduced from properties of shear wave propagation, as group or phase velocity, attenuation or dispersion. Each property can be studied differently in order to retrieve tissue mechanical properties. As an example, shear wave dispersion is linked to elasticity and viscosity of tissues and can be characterized using appropriate rheological models (Chen et al 2004, Deffieux et al 2009). Shear wave group velocity would give access to tissue stiffness (Fatemi and Greenleaf 1998, Sarvazyan et al 1998) and shear wave attenuation to viscosity (Catheline et al 2004). These parameters were studied in different organs often considered as isotropic (Palmeri et al 2008, Tanter et al 2008, Sandrin et al 2003, Bavu et al 2011, Cosgrove et al 2011, Deffieux et al 2011, Berg et al 2012). In fibrous media, such as muscle, biological tissues can be considered as transversely isotropic. Many studies with different techniques were performed regarding shear wave velocities as a function of the main tissue axis (Hoffmeister et al 1996, Gennisson et al 2010, Aubry et al 2013).

By combining acoustic radiation force with ultrafast imaging (Tanter et al 2014), supersonic shear imaging (SSI) was proposed a decade ago as a remote generation of low-frequency shear waves in biological soft tissues (Bercoff et al 2004a). The anisotropic shear viscoelasticity was previously investigated experimentally by using SSI technique in *in vivo biceps brachii* by rotating the ultrasonic probe with respect to the main fiber's direction. (Gennisson et al 2010). SSI was also proposed on other strong anisotropic tissues as *ex vivo myocardium* or *in vivo Achille's* tendon (Lee et al 2012, Brum et al 2013) or arteries (Shcherbakova et al 2014). However, understanding wave anisotropy of shear viscosity and elasticity is still difficult. Recently, a theory of shear wave propagation in a soft tissue with anisotropy of both elastic and dissipative properties has been proposed with no fewer than four independent components for these two characteristics (Rudenko and Sarvazyan 2014). Due to the complex dependence of the shear wave propagation in such tissues, there is a need to have models and numerical simulations of shear wave propagation in incompressible anisotropic viscoelastic tissues. Recently, a theoretical model was proposed and compared to finite element simulations to estimate the elastic tensor of an incompressible transverse isotropic material (Rouze et al 2013).

For over half a century, the waves generated by point sources have been described mainly for seismology in homogeneous isotropic and anisotropic elastic media, such as piezoelectric, thermoelastic and poroelastic solids (Buchwald 1959, Aki and Richards 1980, Burridge 1967, Willis 1980, Norris 1994, Wang and Achenbach 1994). The theoretical Green's functions in an infinite transverse isotropic media have been calculated by Vavryčuk (2001). This description was extended some years later to viscoelastic media by including a frequency dependence in the Green's functions (Vavryčuk 2007, 2008). From these studies, velocity and attenuation of waves propagating from a point source in homogeneous media of arbitrary anisotropy and attenuation strength can be calculated using the Green's formalism. At the same time, the use of the Green's functions approach for low-frequency elastic waves through semi-infinite isotropic and homogeneous soft media was proposed for the development of transient elastography in the biological tissues (Catheline et al 2004, Sandrin et al 2004). The extension to transverse isotropic elastic media has been developed to support the evidence of the polarization of low frequency shear for transient elastography in muscular soft tissues (Gennisson et al 2003). The application of the Green's formalism in an Infinite viscous and elastic medium was then developed for analyzing elastic waves induced by acoustic radiation force imaging (Bercoff et al 2004b).
In the present study we propose another approach using analytical description based on the transverse isotropic, also called hexagonal, viscoelastic Green's formalism. This approach is coupled with a description of ultrasound transducer field simulated by Field's II software (Jensen and Svendsen 1992, Jensen 1996). This numerical approach is successively compared to SSI experiments on a transverse isotropic tissue-mimicking phantom, in *in vivo* forearm muscles and in *in vivo biceps brachii*. The aim of this study is to demonstrate that this simple approach leading to a Green's function formalism can be used to reproduce shear wave propagation in transverse isotropic viscoelastic biological tissue, subject to the prior knowledge of the mechanical properties of the medium.

## 2. Materials and methods

The objective of this study is to propose a new method combining the numerical simulations of the ultrasonic fields and the analytical description of the shear wave propagation. Its application and its comparison to SSI experimental measurements are then successively performed in an anisotropic phantom (non-biological tissue but homogeneous, perfectly transverse isotropic and completely characterized), a forearm muscle (*in vivo* biological tissue supposed homogeneous and transverse isotropic for the analysis of the shear wave field and anisotropic stiffness) and five biceps muscles (*in vivo* biological tissue supposed homogeneous and transverse isotropic for the analysis of the dispersion profile).

### 2.1. Analytical simulation of shear wave propagation in a transverse isotropic viscoelastic medium

The numerical shear wave propagation is developed in two successive parts: after modelling the spatiotemporal distribution of the radiation force field, the displacement and axial velocity fields are calculated using Green's formalism.

#### 2.1.1. Simulation of the radiation force field.
First the pressure focused spot, generating the radiation force field and acting as a shear wave source, is calculated by using Field II ultrasound simulation software (Jensen and Svendsen 1992, Jensen 1996). The linear transducer array used experimentally is designed numerically with the following parameters: geometry (64 elements, 9.6 mm length), focus (40 mm), elements in emission (64 elements), pitch (150 µm), central frequency (5 MHz) and elevation focus (40 mm). The radiation force field $F(\mathbf{x})$ generated in the dissipative medium is mainly oriented along the focus axis (named Z axis) and is deduced from the pressure-field pattern $P(\mathbf{x})$, as defined in equation (1) (Torr 1984).

$$F(\mathbf{x}) = \frac{2\alpha P^2(\mathbf{x})}{\zeta c} \qquad (1)$$

where $\alpha$ is the ultrasound attenuation coefficient (0.6 dB MHz$^{-1}$ cm$^{-1}$ (Culjat et al 2010)), $c$ the sound celerity in the medium and $x = (x, y, z)$ the Cartesian coordinates of the observation point. $\zeta$ is the acoustic impedance. Then each point of the resulting focal spot is used as a shear wave source point in the following part of the simulation. For each point, amplitude of the shear wave is modulated by the force field amplitude.

#### 2.1.2. Calculation of the displacement field.
Secondly, the shear wave propagation is based on Green's function formalism developed by Vavryčuk for seismology (Vavryčuk 2007). The
In the present study we propose another approach using analytical description based on the transverse isotropic, also called hexagonal, viscoelastic Green's formalism. This approach is coupled with a description of ultrasound transducer field simulated by Field's II software (Jensen and Svendsen 1992, Jensen 1996). This numerical approach is successively compared to SSI experiments on a transverse isotropic tissue-mimicking phantom, in *in vivo* forearm muscles and in *in vivo biceps brachii*. The aim of this study is to demonstrate that this simple approach leading to a Green's function formalism can be used to reproduce shear wave propagation in transverse isotropic viscoelastic biological tissue, subject to the prior knowledge of the mechanical properties of the medium.

## 2. Materials and methods

The objective of this study is to propose a new method combining the numerical simulations of the ultrasonic fields and the analytical description of the shear wave propagation. Its application and its comparison to SSI experimental measurements are then successively performed in an anisotropic phantom (non-biological tissue but homogeneous, perfectly transverse isotropic and completely characterized), a forearm muscle (*in vivo* biological tissue supposed homogeneous and transverse isotropic for the analysis of the shear wave field and anisotropic stiffness) and five biceps muscles (*in vivo* biological tissue supposed homogeneous and transverse isotropic for the analysis of the dispersion profile).

### 2.1. Analytical simulation of shear wave propagation in a transverse isotropic viscoelastic medium

The numerical shear wave propagation is developed in two successive parts: after modelling the spatiotemporal distribution of the radiation force field, the displacement and axial velocity fields are calculated using Green's formalism.

#### 2.1.1. Simulation of the radiation force field.
First the pressure focused spot, generating the radiation force field and acting as a shear wave source, is calculated by using Field II ultrasound simulation software (Jensen and Svendsen 1992, Jensen 1996). The linear transducer array used experimentally is designed numerically with the following parameters: geometry (64 elements, 9.6 mm length), focus (40 mm), elements in emission (64 elements), pitch (150 µm), central frequency (5 MHz) and elevation focus (40 mm). The radiation force field $F(x)$ generated in the dissipative medium is mainly oriented along the focus axis (named Z axis) and is deduced from the pressure-field pattern $P(x)$, as defined in equation (1) (Torr 1984).

$$F(x) = \frac{2\alpha P^2(x)}{\zeta c} \qquad (1)$$

where $\alpha$ is the ultrasound attenuation coefficient (0.6 dB MHz$^{-1}$ cm$^{-1}$ (Culjat et al 2010)), $c$ the sound celerity in the medium and $x = (x, y, z)$ the Cartesian coordinates of the observation point. $\zeta$ is the acoustic impedance. Then each point of the resulting focal spot is used as a shear wave source point in the following part of the simulation. For each point, amplitude of the shear wave is modulated by the force field amplitude.

#### 2.1.2. Calculation of the displacement field.
Secondly, the shear wave propagation is based on Green's function formalism developed by Vavryčuk for seismology (Vavryčuk 2007). The
medium is assumed to be semi-infinite, transversely isotropic, viscoelastic and homogeneous. As written in equation (2), the generalized Hooke's law expresses the constitutive relation of a heterogeneous, anisotropic and viscoelastic solid.

$$ \sigma_{ij} = C_{ijkl}\epsilon_{ij}, \quad \text{for } (i, j, k, l) = (1, 2, 3) \qquad (2) $$

where $\sigma_{ij}$ and $\epsilon_{ij}$ are the Cartesian components of the stress and strain tensors, respectively, and $C_{ijkl}$ is the complex component of the fourth-order viscoelastic tensor.

The linear equation of motion for wave propagation derives from the momentum conservation combined with the constitutive relation and can be written as:

$$ \rho \ddot{u}(\mathbf{x}, t) = \nabla \cdot \Sigma(\mathbf{x}, t) + F(\mathbf{x}, t) \qquad (3) $$

where $u(\mathbf{x},t)$ is the displacement field, $\Sigma(\mathbf{x},t)$ is the stress tensor, $F(\mathbf{x},t)$ represents the impulse body force (i.e. the radiation force field defined by the equation (1)), $\rho$ is the density and $\nabla$ is the divergence operator. $\mathbf{x}$ and $t$ correspond to the position vector and time, respectively.

The exact Green's function $U$ gives a time-harmonic plane wave description of $u$ by following the linear equation of motion (equation (3)). The exact elastodynamic Green's function in a semi-infinite, homogeneous, anisotropic, viscoelastic medium can be expressed from a point source in the frequency domain as the sum of three regular and two singular terms. This sum corresponds to a superposition of three homogeneous plane waves propagating in an arbitrary direction with an arbitrary phase velocity and attenuation: one compressional wave (P-wave), one fast shear wave (SV-wave) and one slow shear wave (SH-wave).

The calculation of the displacement vector from one point source to N point sources is developed in the frequency domain in equation (4). As illustrated in figure 1 for temporal representation, the terms $G_{kl}^{(1)}$, $G_{kl}^{(2)}$, $G_{kl}^{(3)}$, $G_{kl}^{(4)}$ and $G_{kl}^{(5)}$ correspond to the contribution of the P-wave, SV-wave, SH-wave, coupling between P-wave and SV-wave and coupling between P-wave and SH-wave, respectively (for more details on this coefficient, the reader can refer to Vavryčuk (2007)).

$$ U_k(\mathbf{x}, \omega) = \sum_{l}^{N} \frac{F_l}{4\pi\rho} \left\{ G_{kl}^{(1)}(\mathbf{x}, \omega) + G_{kl}^{(2)}(\mathbf{x}, \omega) + G_{kl}^{(3)}(\mathbf{x}, \omega) + G_{kl}^{(4)}(\mathbf{x}, \omega) + G_{kl}^{(5)}(\mathbf{x}, \omega) \right\} \qquad (4) $$

where $\omega$ is the frequency, $(k, l)$ the direction index numbers of displacement and source, respectively, $F_l$ the radiation force amplitude and $\rho$ the density. The other terms $G_{kl}$ are defined as:

$$ \begin{aligned} G_{kl}^{(1)}(\mathbf{x}, \omega) &= \frac{g_k^{(1)} g_l^{(1)}}{\tau^{(1)} \sqrt{a_{11}^3}} e^{i\omega\tau^{(1)}}, & G_{kl}^{(2)}(\mathbf{x}, \omega) &= \frac{g_k^{(2)} g_l^{(2)}}{\tau^{(2)} \sqrt{a_{44}^3}} e^{i\omega\tau^{(2)}}, & G_{kl}^{(3)}(\mathbf{x}, \omega) &= \frac{g_k^{(3)} g_l^{(3)}}{\tau^{(3)} a_{66} \sqrt{a_{44}}} e^{i\omega\tau^{(3)}} \\ G_{kl}^{(4)}(\mathbf{x}, \omega) &= \frac{1}{\sqrt{a_{44}}} \frac{g_k^{(3)\perp} g_l^{(3)\perp} - g_k^{(3)} g_l^{(3)}}{R^2} \left(\frac{-1}{i\omega}\right) \left[e^{i\omega\tau^{(3)}} - e^{-i\omega\tau^{(2)}}\right] \\ G_{kl}^{(5)}(\mathbf{x}, \omega) &= \frac{3g_k^{(1)} g_l^{(1)} - \delta_{kl}}{r} \left\{ \frac{1}{a_{11}} e^{i\omega\tau^{(1)}} \left[ \left( \frac{-1}{i\omega\tau^{(1)}} \right) + \left( \frac{-1}{i\omega\tau^{(1)}} \right)^2 \right] \right. \\ & \qquad \left. - \frac{1}{a_{44}} e^{i\omega\tau^{(2)}} \left[ \left( \frac{-1}{i\omega\tau^{(2)}} \right) + \left( \frac{-1}{i\omega\tau^{(2)}} \right)^2 \right] \right\}. \end{aligned} \qquad (5) $$

The travel times $\tau^{(1)}$, $\tau^{(2)}$ and $\tau^{(3)}$ of the compressional wave P, the shear wave parallel to the fibers SV and shear wave perpendicular to the fibers SH respectively, are defined by:
medium is assumed to be semi-infinite, transversely isotropic, viscoelastic and homogeneous. As written in equation (2), the generalized Hooke's law expresses the constitutive relation of a heterogeneous, anisotropic and viscoelastic solid.

$$ \sigma_{ij} = C_{ijkl}\epsilon_{ij}, \quad \text{for } (i, j, k, l) = (1, 2, 3) \qquad (2) $$

where $\sigma_{ij}$ and $\epsilon_{ij}$ are the Cartesian components of the stress and strain tensors, respectively, and $C_{ijkl}$ is the complex component of the fourth-order viscoelastic tensor.

The linear equation of motion for wave propagation derives from the momentum conservation combined with the constitutive relation and can be written as:

$$ \rho \ddot{u}(\mathbf{x}, t) = \nabla \cdot \Sigma(\mathbf{x}, t) + F(\mathbf{x}, t) \qquad (3) $$

where $u(\mathbf{x}, t)$ is the displacement field, $\Sigma(\mathbf{x}, t)$ is the stress tensor, $F(\mathbf{x}, t)$ represents the impulse body force (i.e. the radiation force field defined by the equation (1)), $\rho$ is the density and $\nabla$ is the divergence operator. $\mathbf{x}$ and $t$ correspond to the position vector and time, respectively.

The exact Green's function $U$ gives a time-harmonic plane wave description of $u$ by following the linear equation of motion (equation (3)). The exact elastodynamic Green's function in a semi-infinite, homogeneous, anisotropic, viscoelastic medium can be expressed from a point source in the frequency domain as the sum of three regular and two singular terms. This sum corresponds to a superposition of three homogeneous plane waves propagating in an arbitrary direction with an arbitrary phase velocity and attenuation: one compressional wave (P-wave), one fast shear wave (SV-wave) and one slow shear wave (SH-wave).

The calculation of the displacement vector from one point source to N point sources is developed in the frequency domain in equation (4). As illustrated in figure 1 for temporal representation, the terms $G_{kl}^{(1)}$, $G_{kl}^{(2)}$, $G_{kl}^{(3)}$, $G_{kl}^{(4)}$ and $G_{kl}^{(5)}$ correspond to the contribution of the P-wave, SV-wave, SH-wave, coupling between P-wave and SV-wave and coupling between P-wave and SH-wave, respectively (for more details on this coefficient, the reader can refer to Vavryčuk (2007)).

$$ U_k(\mathbf{x}, \omega) = \sum_l^N \frac{F_l}{4\pi\rho} \left\{ G_{kl}^{(1)}(\mathbf{x}, \omega) + G_{kl}^{(2)}(\mathbf{x}, \omega) + G_{kl}^{(3)}(\mathbf{x}, \omega) + G_{kl}^{(4)}(\mathbf{x}, \omega) + G_{kl}^{(5)}(\mathbf{x}, \omega) \right\} \qquad (4) $$

where $\omega$ is the frequency, $(k, l)$ the direction index numbers of displacement and source, respectively, $F_l$ the radiation force amplitude and $\rho$ the density. The other terms $G_{kl}$ are defined as:

$$ \begin{aligned} G_{kl}^{(1)}(\mathbf{x}, \omega) &= \frac{g_k^{(1)} g_l^{(1)}}{\tau^{(1)} \sqrt{a_{11}^3}} e^{i\omega\tau^{(1)}}, \quad G_{kl}^{(2)}(\mathbf{x}, \omega) = \frac{g_k^{(2)} g_l^{(2)}}{\tau^{(2)} \sqrt{a_{44}^3}} e^{i\omega\tau^{(2)}}, \quad G_{kl}^{(3)}(\mathbf{x}, \omega) = \frac{g_k^{(3)} g_l^{(3)}}{\tau^{(3)} a_{66} \sqrt{a_{44}}} e^{i\omega\tau^{(3)}} \\ G_{kl}^{(4)}(\mathbf{x}, \omega) &= \frac{1}{\sqrt{a_{44}}} \frac{g_k^{(3)\perp} g_l^{(3)\perp} - g_k^{(3)} g_l^{(3)}}{R^2} \left(\frac{-1}{i\omega}\right) \left[e^{i\omega\tau^{(3)}} - e^{-i\omega\tau^{(2)}}\right] \\ G_{kl}^{(5)}(\mathbf{x}, \omega) &= \frac{3g_k^{(1)} g_l^{(1)} - \delta_{kl}}{r} \left\{ \frac{1}{a_{11}} e^{i\omega\tau^{(1)}} \left[ \left( \frac{-1}{i\omega\tau^{(1)}} \right) + \left( \frac{-1}{i\omega\tau^{(1)}} \right)^2 \right] \right. \\ & \qquad \left. - \frac{1}{a_{44}} e^{i\omega\tau^{(2)}} \left[ \left( \frac{-1}{i\omega\tau^{(2)}} \right) + \left( \frac{-1}{i\omega\tau^{(2)}} \right)^2 \right] \right\}. \end{aligned} \qquad (5) $$

The travel times $\tau^{(1)}, \tau^{(2)}$ and $\tau^{(3)}$ of the compressional wave P, the shear wave parallel to the fibers SV and shear wave perpendicular to the fibers SH respectively, are defined by:
<figure><img src="image_3.png" /><figcaption>Figure 1. Representation of the axial displacement field at one time and calculated from a source point in one plane using the Green’s functions. The wave front is calculated as the summation of a P-wave, a SV-wave, a SH-wave, a coupling term between P-wave and SV-wave and a coupling term between P-wave and SH-wave.</figcaption></figure>

$$ \tau^{(1)} = \frac{r}{\sqrt{a_{11}}}, \quad \tau^{(2)} = \frac{r}{\sqrt{a_{44}}}, \quad \tau^{(3)} = \frac{r}{\sqrt{a_{66}}} \sqrt{N_1^2 + N_2^2 + \frac{a_{66}}{a_{44}} N_3^2}. \quad (6) $$

The polarization vectors are defined by:

$$ g^{(1)} = \begin{bmatrix} N_1 \\ N_2 \\ N_3 \end{bmatrix}, \quad g^{(2)} = \frac{1}{\sqrt{N_1^2 + N_2^2}} \begin{bmatrix} -N_1 N_3 \\ -N_2 N_3 \\ N_1^2 + N_2^2 \end{bmatrix}, \quad g^{(3)} = \frac{1}{\sqrt{N_1^2 + N_2^2}} \begin{bmatrix} N_2 \\ -N_1 \\ 0 \end{bmatrix}, \\ g^{(3)\perp} = \frac{1}{\sqrt{N_1^2 + N_2^2}} \begin{bmatrix} N_1 \\ N_2 \\ 0 \end{bmatrix}. \quad (7) $$

Then, $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.

At last, assuming a Voigt’s model, $a_{11} = 1/\rho(c_p^2 + i\omega\eta_p)$, $a_{44} = 1/\rho(c_{SV}^2 + i\omega\eta_{SV})$ and $a_{66} = 1/\rho(c_{SH}^2 + i\omega\eta_{SH})$ are the density-normalized viscoelastic coefficients of the complex Christoffel tensor, where $c_P$, $c_{SV}$, $c_{SH}$ and $\eta_P$, $\eta_{SV}$, $\eta_{SH}$ correspond to the ‘zero-viscosity’ celerity and to the dynamic viscosity of the P-wave, SV-wave and SH-wave, respectively. It is important to note that the $c_P$, $c_{SV}$ and $c_{SH}$ parameters correspond to the group velocities in only one case of pure elasticity, i.e. with negligible viscosity. Under viscoelasticity assumption,
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Representation of the axial displacement field at one time and calculated from a source point in one plane using the Green's functions. The wave front is calculated as the summation of a P-wave, a SV-wave, a SH-wave, a coupling term between P-wave and SV-wave and a coupling term between P-wave and SH-wave.</figcaption></figure>

$$ \tau^{(1)} = \frac{r}{\sqrt{a_{11}}}, \quad \tau^{(2)} = \frac{r}{\sqrt{a_{44}}}, \quad \tau^{(3)} = \frac{r}{\sqrt{a_{66}}} \sqrt{N_1^2 + N_2^2 + \frac{a_{66}}{a_{44}} N_3^2}. \quad (6) $$

The polarization vectors are defined by:

$$ g^{(1)} = \begin{bmatrix} N_1 \\ N_2 \\ N_3 \end{bmatrix}, \quad g^{(2)} = \frac{1}{\sqrt{N_1^2 + N_2^2}} \begin{bmatrix} -N_1 N_3 \\ -N_2 N_3 \\ N_1^2 + N_2^2 \end{bmatrix}, \quad g^{(3)} = \frac{1}{\sqrt{N_1^2 + N_2^2}} \begin{bmatrix} N_2 \\ -N_1 \\ 0 \end{bmatrix}, \\ g^{(3)\perp} = \frac{1}{\sqrt{N_1^2 + N_2^2}} \begin{bmatrix} N_1 \\ N_2 \\ 0 \end{bmatrix}. \quad (7) $$

Then, $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.

At last, assuming a Voigt's model, $a_{11} = 1/\rho(c_p^2 + i\omega\eta_p)$, $a_{44} = 1/\rho(c_{SV}^2 + i\omega\eta_{SV})$ and $a_{66} = 1/\rho(c_{SH}^2 + i\omega\eta_{SH})$ are the density-normalized viscoelastic coefficients of the complex Christoffel tensor, where $c_p, c_{SV}, c_{SH}$ and $\eta_p, \eta_{SV}, \eta_{SH}$ correspond to the 'zero-viscosity' celerity and to the dynamic viscosity of the P-wave, SV-wave and SH-wave, respectively. It is important to note that the $c_p, c_{SV}$ and $c_{SH}$ parameters correspond to the group velocities in only one case of pure elasticity, i.e. with negligible viscosity. Under viscoelasticity assumption,
$[\eta_{1,SV} = c_{SV}^2; \eta_{2,SV} = \eta_{SV}]$ and $[\eta_{1,SH} = c_{SH}^2; \eta_{2,SH} = \eta_{SH}]$ correspond to the elasticity and damping parameters of the Voigt's model for the longitudinal and transverse direction, respectively. Assuming homogeneous density of the tissue, these $a_{11}$, $a_{44}$ and $a_{66}$ coefficients of the density-normalized transverse isotropic stiffness tensor are defined by a link to the complex viscoelastic Christoffel tensor, as shown by equation (8).

$$A_{ijkl} = \frac{1}{\rho} C_{ijkl}, \qquad (8)$$

The harmonic $U_k(\mathbf{x}, \omega)$ displacement field is then multiplied by the Fourier transform of a temporal rectangular function with $t = 300 \mu s$ duration (corresponding to the duration for the application of the radiation force field $F(\mathbf{x})$). Finally, the temporal 3D displacement field $U_k(\mathbf{x}, t)$ is deduced using the inverse Fourier transform, as well as the axial velocity field of the shear wave velocity $\dot{U}_z(\mathbf{x}, t)$ after temporal differentiation of $U_z(\mathbf{x}, t)$.

## 2.2. Experimental validation in a transversely isotropic phantom

In this part the probe was a linear phased array dedicated to pediatric cardiac diagnosis (linear phased array, 5 MHz central frequency, 150 µm pitch, 64 elements, Vermon, Tours, France) connected to an ultrafast Aixplorer scanner (8280 frames/s, SuperSonic Imagine, Aix-en-Provence, France).

The experimental measurement of the maximal radiation pressure, as well as the ultrasonic beam pressure pattern, was first performed for this probe by using a laser interferometry technique. A high resolution optical system has been used for imaging pulsed ultrasonic fields transmitted in water by the ultrasonic probe. As described by Barrière and Royer (Barrière and Royer 2001b), the particle displacement was measured by a compact heterodyne interferometer whose probe beam was reflected by a thin membrane immersed in front of the transducer. After modelization of the spatio-temporal response of the membrane to a specific pressure (Royer et al 1992, Barrière and Royer 2001b), the beam pressure pattern was obtained. From these measurements, the mechanical Index ($MI$) and spatial-peak time average ($I_{SPTA}$) were calculated and both $MI_{0.3}$ and $I_{SPTA_{0.3}}$ were inferior to 1.90 and 720 mW cm$^{-2}$, respectively, as imposed by the Food and Drug Administration (FDA) for the muscular soft tissue. The aim of these experimental measurements is both to control the numerical pressure-field pattern obtained from the FieldII software and to insure the innocuousness of the in vivo measurements.

Then the Green's function formalism was validated by performing SSI experiments in a transverse isotropic phantom. The phantom was prepared using a recently described protocol aimed at inducing mechanical transverse isotropy in a cryogel PVA phantom by the alignment of its polymeric fibers along a preferential direction (Millon et al 2006, Wan et al 2009, Chatelin et al 2014). The sample was composed of a 5w% PVA solution (molecular weight 89000–98000, 99+% hydrolyzed, Sigma-Aldrich, St Louis, US) and 1w% cellulose (20 µm in diameter, S3504 Sigmacell, Sigma-Aldrich, St Louis, US) diluted in degassed distilled water. After transfer into a parallelepiped PVC mold (100 × 60 × 50 mm$^3$), the solution underwent two freeze/thaw cycles (i.e. 20 °C during 12 h and −18 °C during the 12 next h) followed by three freeze/thaw cycles, which include stretching to 80% of the initial length using a customized static tensile-test set up. The cross-link was induced in the polymer along a preferential orientation that corresponds to the stretching direction.

SSI experiments were performed using the following protocol: an ultrasonic focused beam was used to remotely generate radiation acoustic forced resulting in low-frequency propagation
$[\eta_{1,SV} = c_{SV}^2; \eta_{2,SV} = \eta_{SV}]$ and $[\eta_{1,SH} = c_{SH}^2; \eta_{2,SH} = \eta_{SH}]$ correspond to the elasticity and damping parameters of the Voigt's model for the longitudinal and transverse direction, respectively. Assuming homogeneous density of the tissue, these $a_{11}$, $a_{44}$ and $a_{66}$ coefficients of the density-normalized transverse isotropic stiffness tensor are defined by a link to the complex viscoelastic Christoffel tensor, as shown by equation (8).

$$A_{ijkl} = \frac{1}{\rho} C_{ijkl}, \qquad (8)$$

The harmonic $U_k(\mathbf{x}, \omega)$ displacement field is then multiplied by the Fourier transform of a temporal rectangular function with $t = 300 \mu s$ duration (corresponding to the duration for the application of the radiation force field $F(\mathbf{x})$). Finally, the temporal 3D displacement field $U_k(\mathbf{x}, t)$ is deduced using the inverse Fourier transform, as well as the axial velocity field of the shear wave velocity $\dot{U}_z(\mathbf{x}, t)$ after temporal differentiation of $U_z(\mathbf{x}, t)$.

## 2.2. Experimental validation in a transversely isotropic phantom

In this part the probe was a linear phased array dedicated to pediatric cardiac diagnosis (linear phased array, 5 MHz central frequency, 150 µm pitch, 64 elements, Vermon, Tours, France) connected to an ultrafast Aixplorer scanner (8280 frames/s, SuperSonic Imagine, Aix-en-Provence, France).

The experimental measurement of the maximal radiation pressure, as well as the ultrasonic beam pressure pattern, was first performed for this probe by using a laser interferometry technique. A high resolution optical system has been used for imaging pulsed ultrasonic fields transmitted in water by the ultrasonic probe. As described by Barrière and Royer (Barrière and Royer 2001b), the particle displacement was measured by a compact heterodyne interferometer whose probe beam was reflected by a thin membrane immersed in front of the transducer. After modelization of the spatio-temporal response of the membrane to a specific pressure (Royer et al 1992, Barrière and Royer 2001b), the beam pressure pattern was obtained. From these measurements, the mechanical Index ($MI$) and spatial-peak time average ($I_{SPTA}$) were calculated and both $MI_{0.3}$ and $I_{SPTA_{0.3}}$ were inferior to 1.90 and 720 mW cm$^{-2}$, respectively, as imposed by the Food and Drug Administration (FDA) for the muscular soft tissue. The aim of these experimental measurements is both to control the numerical pressure-field pattern obtained from the FieldII software and to insure the innocuousness of the in vivo measurements.

Then the Green's function formalism was validated by performing SSI experiments in a transverse isotropic phantom. The phantom was prepared using a recently described protocol aimed at inducing mechanical transverse isotropy in a cryogel PVA phantom by the alignment of its polymeric fibers along a preferential direction (Millon et al 2006, Wan et al 2009, Chatelin et al 2014). The sample was composed of a 5w% PVA solution (molecular weight 89000–98000, 99+% hydrolyzed, Sigma-Aldrich, St Louis, US) and 1w% cellulose (20 µm in diameter, S3504 Sigmacell, Sigma-Aldrich, St Louis, US) diluted in degassed distilled water. After transfer into a parallelepiped PVC mold (100 × 60 × 50 mm$^3$), the solution underwent two freeze/thaw cycles (i.e. 20 °C during 12 h and −18 °C during the 12 next h) followed by three freeze/thaw cycles, which include stretching to 80% of the initial length using a customized static tensile-test set up. The cross-link was induced in the polymer along a preferential orientation that corresponds to the stretching direction.

SSI experiments were performed using the following protocol: an ultrasonic focused beam was used to remotely generate radiation acoustic forced resulting in low-frequency propagation
of shear waves (figure 2). Then by using ultrafast imaging mode, a movie of the propagating shear wave was recorded, giving access to shear wave speed (Bercoff et al 2004a). The shear wave group velocity in the longitudinal ($\theta = 0^{\circ}$) and transverse ($\theta = 90^{\circ}$) directions were used as input for the numerical simulations. Because of the results from the first step of the experiments, the input parameters used for the simulations are related in the result subsection for the tests on both phantom, forearm and biceps muscle. The viscosity value was measured in a previous study (0.1 Pa s) by different rheological technique in the transverse isotropic PVA phantom (Chatelin et al 2014). It has to be noticed that the use of viscosity inferior to 0.4 Pa s produces negligible effects on the shear wave profile for the field of view and stiffness considered in this study. The profiles of the shear wave speed in function of the probe angle $\theta$ (from -90° to 90° by 10° steps) were then compared between experiments and simulations.

## 2.3. Experimental validation in vivo on human muscle

To present a distribution of fibers oriented along a preferential direction, muscle is the best candidate for a transverse isotropic elastic medium. In this step comparison of the simulations with the measurements of anisotropic shear wave group velocity and shear wave dispersion was performed in vivo on the contracted forearm and in the biceps brachii respectively.

**On the forearm muscle,** the shear wave speed was estimated using a dedicated flight-time algorithm that measures the shear wave group velocity for various angles $\theta$ (from -90° to 90° by 10° steps) between the probe and the fibers. The maximum speed was reached by aligning the probe with the muscle ($\theta = 0^{\circ}$), and the minimum when the probe was perpendicular to the muscle fibers ($\theta = \pm 90^{\circ}$). These values were used as inputs for the Green's function numerical simulations as the longitudinal and transverse shear wave speed, respectively. After simulation, the profiles of the shear wave speed were then compared between experiments and simulations over the probe angle $\theta$ (from -90° to 90° by 10° steps).

**On the biceps brachii,** numerical simulations were compared to experimental data already published in Gennisson et al (2010). The reader can refer to this article for more details about the experimental material and methods. In this study, the phase velocity profiles were estimated as a function of the frequency using the algorithm described by Deffieux et al (2009). The shear wave phase for each frequency was computed by calculating the Fourier transform of experimental velocity field acquired as a function of time. The phase velocity $c_{\varphi}$ of shear wave is given by equation (9) and the real part of the complex wave number $Re[k]$ can be estimated from a linear fit of the phase function along the axis of propagation of the wave, where $\omega$ is the angular frequency.

$$c_{\varphi} = \omega/\text{Re}[k] \qquad (9)$$

Assuming a Voigt's model, the estimation of the shear wave speed at multiple frequencies leads to the global estimation of both shear modulus and viscosity. The longitudinal $\eta_{SV}$ and transverse $\eta_{SH}$ viscosities were identified from the numerical simulations by fitting the phase velocity as a function of the frequency (equations (10) and (11)) in accordance with the Voigt's model. This model is essentially applicable to the highest frequencies and for the far-field from the push area. The shear wave profile has a bandwidth from 120 to 800 with a central frequency close to 450 Hz.

$$c_{\varphi SV} = \sqrt{\frac{2(\mu_{SV}^2 + \omega^2\eta_{SV}^2)}{\rho(\mu_{SV} + \sqrt{\mu_{SV}^2 + \omega^2\eta_{SV}^2})}} \qquad (10)$$
of shear waves (figure 2). Then by using ultrafast imaging mode, a movie of the propagating shear wave was recorded, giving access to shear wave speed (Bercoff et al 2004a). The shear wave group velocity in the longitudinal ($\theta = 0^{\circ}$) and transverse ($\theta = 90^{\circ}$) directions were used as input for the numerical simulations. Because of the results from the first step of the experiments, the input parameters used for the simulations are related in the result subsection for the tests on both phantom, forearm and biceps muscle. The viscosity value was measured in a previous study (0.1 Pa s) by different rheological technique in the transverse isotropic PVA phantom (Chatelin et al 2014). It has to be noticed that the use of viscosity inferior to 0.4 Pa s produces negligible effects on the shear wave profile for the field of view and stiffness considered in this study. The profiles of the shear wave speed in function of the probe angle $\theta$ (from -90° to 90° by 10° steps) were then compared between experiments and simulations.

## 2.3. Experimental validation in vivo on human muscle

To present a distribution of fibers oriented along a preferential direction, muscle is the best candidate for a transverse isotropic elastic medium. In this step comparison of the simulations with the measurements of anisotropic shear wave group velocity and shear wave dispersion was performed in vivo on the contracted forearm and in the biceps brachii respectively.

**On the forearm muscle,** the shear wave speed was estimated using a dedicated flight-time algorithm that measures the shear wave group velocity for various angles $\theta$ (from -90° to 90° by 10° steps) between the probe and the fibers. The maximum speed was reached by aligning the probe with the muscle ($\theta = 0^{\circ}$), and the minimum when the probe was perpendicular to the muscle fibers ($\theta = \pm 90^{\circ}$). These values were used as inputs for the Green's function numerical simulations as the longitudinal and transverse shear wave speed, respectively. After simulation, the profiles of the shear wave speed were then compared between experiments and simulations over the probe angle $\theta$ (from -90° to 90° by 10° steps).

**On the biceps brachii,** numerical simulations were compared to experimental data already published in Gennisson et al (2010). The reader can refer to this article for more details about the experimental material and methods. In this study, the phase velocity profiles were estimated as a function of the frequency using the algorithm described by Deffieux et al (2009). The shear wave phase for each frequency was computed by calculating the Fourier transform of experimental velocity field acquired as a function of time. The phase velocity $c_{\varphi}$ of shear wave is given by equation (9) and the real part of the complex wave number $Re[k]$ can be estimated from a linear fit of the phase function along the axis of propagation of the wave, where $\omega$ is the angular frequency.

$$c_{\varphi} = \omega/\text{Re}[k] \qquad (9)$$

Assuming a Voigt's model, the estimation of the shear wave speed at multiple frequencies leads to the global estimation of both shear modulus and viscosity. The longitudinal $\eta_{SV}$ and transverse $\eta_{SH}$ viscosities were identified from the numerical simulations by fitting the phase velocity as a function of the frequency (equations (10) and (11)) in accordance with the Voigt's model. This model is essentially applicable to the highest frequencies and for the far-field from the push area. The shear wave profile has a bandwidth from 120 to 800 with a central frequency close to 450 Hz.

$$c_{\varphi SV} = \sqrt{\frac{2(\mu_{SV}^2 + \omega^2\eta_{SV}^2)}{\rho(\mu_{SV} + \sqrt{\mu_{SV}^2 + \omega^2\eta_{SV}^2})}} \qquad (10)$$
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The basic principle of SSI (a): an acoustic beam (grey area) was generated and focused (blue spot) from a probe to induce shear waves in the medium (red lines and arrows). These waves were caught by using ultrafast imaging mode. The setup for the experimental measurements in a transverse isotropic PVA hydrogel phantom (b) and in the forearm muscle (c) are based on this principle. The in-plane stiffness constants for different orientations are deduced from the SSI shear wave speed measurements by rotation of the probe (θ angle).</figcaption></figure>

$$c_{\phi SH} = \frac{\sqrt{\frac{2(\mu_{SH}^2 + \omega^2\eta_{SH}^2)}{\rho(\mu_{SH} + \sqrt{\mu_{SH}^2 + \omega^2\eta_{SH}^2})}}}{\sqrt{\rho(\mu_{SH} + \sqrt{\mu_{SH}^2 + \omega^2\eta_{SH}^2})}} \quad (11)$$

$\mu_{SV}$ and $\mu_{SH}$ are the longitudinal and transverse shear elasticity, respectively, and correspond to the elastic parameters of the Voigt's model. The damping part of the longitudinal and transverse Voigt's models is given by the longitudinal $\eta_{SV}$ and transverse $\eta_{SH}$ dynamic viscosities.

The $\mu_{SV}$, $\mu_{SH}$, $\eta_{SV}$ and $\eta_{SH}$ parameters have been identified by fitting equations (10) and (11) to the longitudinal and transverse shear wave dispersion profiles in Gennisson et al (2010). These parameters have been used in the present study as input parameters for the anisotropic viscoelastic Green's simulations. The numerical dispersion profiles were compared to the experimental ones presented by Gennisson et al in 2010.

Because the parameters used as input for the simulations came mostly from the first part of the experiments, we decided to detail these parameters in tables in the section 3.

## 3. Results

### 3.1. Comparison between simulation and SSI experiments in phantom

The pressure focal spot size at full width at −6dB attenuation at the focus was 0.8 mm in the azimuthal direction, 0.8 mm in the elevation direction and 10 mm in the axial beam (Z) direction, with a maximal value of 1.205 MPa.

Figure 3 represents the experimental shear wave propagation in the transverse isotropic phantom in both longitudinal (a and c) and transverse (b and d) direction. The resulting axial (Z) velocity field of the shear wave velocity over time in the longitudinal (a) and transverse
<figure><img src="image_5.png" /><figcaption>**Figure 2.** The basic principle of SSI (a): an acoustic beam (grey area) was generated and focused (blue spot) from a probe to induce shear waves in the medium (red lines and arrows). These waves were caught by using ultrafast imaging mode. The setup for the experimental measurements in a transverse isotropic PVA hydrogel phantom (b) and in the forearm muscle (c) are based on this principle. The in-plane stiffness constants for different orientations are deduced from the SSI shear wave speed measurements by rotation of the probe (θ angle).</figcaption></figure>

$$c_{\phi SH} = \sqrt{\frac{2(\mu_{SH}^2 + \omega^2\eta_{SH}^2)}{\rho(\mu_{SH} + \sqrt{\mu_{SH}^2 + \omega^2\eta_{SH}^2})}} \quad (11)$$

$\mu_{SV}$ and $\mu_{SH}$ are the longitudinal and transverse shear elasticity, respectively, and correspond to the elastic parameters of the Voigt's model. The damping part of the longitudinal and transverse Voigt's models is given by the longitudinal $\eta_{SV}$ and transverse $\eta_{SH}$ dynamic viscosities.

The $\mu_{SV}$, $\mu_{SH}$, $\eta_{SV}$ and $\eta_{SH}$ parameters have been identified by fitting equations (10) and (11) to the longitudinal and transverse shear wave dispersion profiles in Gennisson et al (2010). These parameters have been used in the present study as input parameters for the anisotropic viscoelastic Green's simulations. The numerical dispersion profiles were compared to the experimental ones presented by Gennisson et al in 2010.

Because the parameters used as input for the simulations came mostly from the first part of the experiments, we decided to detail these parameters in tables in the section 3.

## 3. Results

### 3.1. Comparison between simulation and SSI experiments in phantom

The pressure focal spot size at full width at −6dB attenuation at the focus was 0.8 mm in the azimuthal direction, 0.8 mm in the elevation direction and 10 mm in the axial beam (Z) direction, with a maximal value of 1.205 MPa.

Figure 3 represents the experimental shear wave propagation in the transverse isotropic phantom in both longitudinal (a and c) and transverse (b and d) direction. The resulting axial (Z) velocity field of the shear wave velocity over time in the longitudinal (a) and transverse
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The experimental axial (Z) velocity field of the shear wave velocity over time in the longitudinal (a) and transverse (b) directions (superposed to B-mode images) are represented at the same sampling times (1 ms after the push generation) for SSI in the anisotropic hydrogel phantom. The spatio-temporal experimental imaging of the shear wave in the longitudinal (c) and transverse (d) directions in a focal area (rectangle in red color) are compared to the numerical ones for the same orientations (e and f, respectively).</figcaption></figure>

(b) directions is represented at the same sampling times (1 ms after the push generation). By fitting the slope, the spatio-temporal imaging (c and d) of the shear wave gives an evaluation of the longitudinal and transverse velocity. $C_{SV}$ and $c_{SH}$ were found to be 3.24 ± 0.853 m s⁻¹ and 1.95 ± 0.326 m s⁻¹, respectively. Based on these experimental values and the pressure field calculated using the FieldII software, the parameters summarized in table 1 are implemented in the numerical simulations of the shear wave propagation in the anisotropic phantom.

The resulting numerical longitudinal (e) and transverse (f) spatio-temporal imaging is illustrated in figure 4 for comparison with the experimental results (c and d). Despite some differences in the shear wave amplitudes, similar shear wave speeds are observed between experimental and numerical data, validating the simulations for these two particular directions.

Figure 4 shows the shear wave speed variations for both experimental and numerical results in the fibrous medium by variation of the in-plane angle $\theta$ between the fibers of the phantom and the shear wave propagation direction (a). The shear wave profile propagating over different time steps is also represented (b).

## 3.2. Simulation versus in vivo SSI experiments

Similarly to the SSI data for the anisotropic phantom, the experimental shear wave propagation in the forearm muscle is illustrated in figure 5 for both longitudinal (a and c) and transverse (b and d) directions, i.e. along ($\theta$ = 0°) or perpendicularly ($\theta$ = 90°) to the muscular fibers. The resulting axial (Z) velocity fields in both longitudinal (a) and transverse (b) directions are represented at the same sampling times (1 ms after the push generation). The results are
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The experimental axial (Z) velocity field of the shear wave velocity over time in the longitudinal (a) and transverse (b) directions (superposed to B-mode images) are represented at the same sampling times (1 ms after the push generation) for SSI in the anisotropic hydrogel phantom. The spatio-temporal experimental imaging of the shear wave in the longitudinal (c) and transverse (d) directions in a focal area (rectangle in red color) are compared to the numerical ones for the same orientations (e and f, respectively).</figcaption></figure>

(b) directions is represented at the same sampling times (1 ms after the push generation). By fitting the slope, the spatio-temporal imaging (c and d) of the shear wave gives an evaluation of the longitudinal and transverse velocity. $C_{SV}$ and $c_{SH}$ were found to be 3.24 ± 0.853 m s⁻¹ and 1.95 ± 0.326 m s⁻¹, respectively. Based on these experimental values and the pressure field calculated using the FieldII software, the parameters summarized in table 1 are implemented in the numerical simulations of the shear wave propagation in the anisotropic phantom.

The resulting numerical longitudinal (e) and transverse (f) spatio-temporal imaging is illustrated in figure 4 for comparison with the experimental results (c and d). Despite some differences in the shear wave amplitudes, similar shear wave speeds are observed between experimental and numerical data, validating the simulations for these two particular directions.

Figure 4 shows the shear wave speed variations for both experimental and numerical results in the fibrous medium by variation of the in-plane angle $\theta$ between the fibers of the phantom and the shear wave propagation direction (a). The shear wave profile propagating over different time steps is also represented (b).

## 3.2. Simulation versus in vivo SSI experiments

Similarly to the SSI data for the anisotropic phantom, the experimental shear wave propagation in the forearm muscle is illustrated in figure 5 for both longitudinal (a and c) and transverse (b and d) directions, i.e. along ($\theta$ = 0°) or perpendicularly ($\theta$ = 90°) to the muscular fibers. The resulting axial (Z) velocity fields in both longitudinal (a) and transverse (b) directions are represented at the same sampling times (1 ms after the push generation). The results are
<table><caption>**Table 1.** Experimental SSI data used to determine the numerical parameters as input for the shear wave simulations in the transverse isotropic phantom.</caption><thead><tr><th></th><th>c<sub>SV</sub><br>[m s<sup>-1</sup>]</th><th>c<sub>SH</sub><br>[m s<sup>-1</sup>]</th><th>c<sub>P</sub> [m s<sup>-1</sup>]</th><th>η<sub>P</sub><br>[Pas]</th><th>η<sub>SV</sub><br>[Pas]</th><th>η<sub>SH</sub><br>[Pas]</th><th>ρ<br>[g cm<sup>-3</sup>]</th></tr></thead><tbody><tr><td>Experimental measurements</td><td>3.24 ± 0.85</td><td>1.95 ± 0.33</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Numerical parameters</td><td>3.24</td><td>1.95</td><td>1500</td><td>0.1</td><td>0.1</td><td>0.1</td><td>1000</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** The shear wave speed variations (b) are extracted from the simulation (blue dots and line) as a function of the probe/fibers angle $\theta$ and compared to the experimental data (orange dots, including standard deviation) in the anisotropic phantom. The longitudinal Z-displacements induced in the focal plane by the acoustic radiation force field are represented at different sampling times from the Green's function formalism simulation in a transversely isotropic viscoelastic medium (a).</figcaption></figure>

analyzed in a Region Of Interest (ROI) defined manually to obtain data in a homogeneous part of the muscle, as far as possible from muscular layer interfaces. By fitting the slope, the spatio-temporal imaging (c and d) of the shear wave gives experimental measurements for the longitudinal and transverse group velocity of $6.69 \pm 0.234 \text{ m s}^{-1}$ and $2.16 \pm 0.187$, respectively. For this simulation, the viscosity has been implemented in accordance with values from the literature, supposing a similar viscous behavior for both forearm and biceps muscle (Gennisson et al 2010). We consider these group velocities as equal to the ‘zero-viscosity’ velocities, i.e. giving a direct estimation of the elastic part of the Voigt’s model. This assumption will be further discussed in the section 4. Consequently, these parameters are used as $c_{SV}$ and $c_{SH}$ input parameters for the numerical Green’s simulations. Based on these experimental values and the pressure field exposed in the 3.1 subsection, the parameters summarized in table 2 are implemented in the numerical simulations of the shear wave propagation in the forearm muscle.

Figure 5(e) shows the shear wave speed variations for both experimental and numerical results in the fibrous medium by variation of the in-plane angle $\theta$ between the muscular fibers
<table><caption>**Table 1.** Experimental SSI data used to determine the numerical parameters as input for the shear wave simulations in the transverse isotropic phantom.</caption><thead><tr><th></th><th>c<sub>SV</sub><br>[m s<sup>-1</sup>]</th><th>c<sub>SH</sub><br>[m s<sup>-1</sup>]</th><th>c<sub>P</sub> [<br>m s<sup>-1</sup>]</th><th>η<sub>P</sub><br>[Pas]</th><th>η<sub>SV</sub><br>[Pas]</th><th>η<sub>SH</sub><br>[Pas]</th><th>ρ<br>[g cm<sup>-3</sup>]</th></tr></thead><tbody><tr><td>Experimental<br>measurements</td><td>3.24 ± 0.85</td><td>1.95 ± 0.33</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Numerical<br>parameters</td><td>3.24</td><td>1.95</td><td>1500</td><td>0.1</td><td>0.1</td><td>0.1</td><td>1000</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** The shear wave speed variations (b) are extracted from the simulation (blue dots and line) as a function of the probe/fibers angle $\theta$ and compared to the experimental data (orange dots, including standard deviation) in the anisotropic phantom. The longitudinal Z-displacements induced in the focal plane by the acoustic radiation force field are represented at different sampling times from the Green's function formalism simulation in a transversely isotropic viscoelastic medium (a).</figcaption></figure>

analyzed in a Region Of Interest (ROI) defined manually to obtain data in a homogeneous part of the muscle, as far as possible from muscular layer interfaces. By fitting the slope, the spatio-temporal imaging (c and d) of the shear wave gives experimental measurements for the longitudinal and transverse group velocity of 6.69 ± 0.234 m s⁻¹ and 2.16 ± 0.187, respectively. For this simulation, the viscosity has been implemented in accordance with values from the literature, supposing a similar viscous behavior for both forearm and biceps muscle (Gennisson et al 2010). We consider these group velocities as equal to the ‘zero-viscosity’ velocities, i.e. giving a direct estimation of the elastic part of the Voigt’s model. This assumption will be further discussed in the section 4. Consequently, these parameters are used as $c_{SV}$ and $c_{SH}$ input parameters for the numerical Green’s simulations. Based on these experimental values and the pressure field exposed in the 3.1 subsection, the parameters summarized in table 2 are implemented in the numerical simulations of the shear wave propagation in the forearm muscle.

Figure 5(e) shows the shear wave speed variations for both experimental and numerical results in the fibrous medium by variation of the in-plane angle $\theta$ between the muscular fibers
<figure><img src="image_7.png" /><figcaption>**Figure 5.** The experimental axial (Z) velocity fields in both longitudinal (a) and transverse (b) directions (superposed to B-mode images) are represented at the same sampling times for SWE in the *in vivo* forearm muscle. The spatio-temporal experimental imaging of the shear wave in the longitudinal (c) and transverse (d) directions confirms the anisotropic profile by a significant variation of the slope. The shear wave speed variations (e) are extracted from these experimental measurements (green dots, including standard deviation) as a function of the probe / fibers angle θ and compared to the anisotropic data (blue dots and line) from the numerical simulations.</figcaption></figure>

<table><caption>**Table 2.** Experimental SWE data used to determine the numerical parameters as input for the shear wave simulations in the *in vivo* forearm muscle.</caption><thead><tr><th></th><th>c<sub>SV</sub><br>[m s<sup>-1</sup>]</th><th>c<sub>SH</sub><br>[m s<sup>-1</sup>]</th><th>c<sub>P</sub><br>[m s<sup>-1</sup>]</th><th>η<sub>SV</sub><br>[Pa s]</th><th>η<sub>SH</sub><br>[Pa s]</th><th>η<sub>P</sub><br>[Pa s]</th><th>ρ<br>[g cm<sup>-3</sup>]</th></tr></thead><tbody><tr><td>Experimental measurements</td><td>6.69 ± 0.23</td><td>2.16 ± 0.19</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Numerical parameters</td><td>6.50</td><td>2.16</td><td>1500</td><td>0.65 ± 0.10<sup>a</sup></td><td>0.92 ± 0.06<sup>a</sup></td><td>0.1</td><td>1000</td></tr></tbody></table>

and the shear wave propagation direction. The highest speed is well correlated with the fiber orientation for both approaches. Whatever the fiber orientation is, the experimental results and the viscoelastic Green's function simulations are in very good agreement. This helps to validate the elastic behavior of the shear wave propagation in the anisotropic simulations.

The capability of the numerical method proposed in this study to simulate the anisotropic viscoelastic behavior of fibrous soft tissues is here compared to experimental data on the *in vivo* relaxed biceps from Gennisson *et al* (2010). The experimental values from this study, as well as the numerical input parameters for the simulations in the biceps muscle, are summarized in table 3.

For θ = 0° and θ = 90°, the numerical longitudinal and transverse shear wave velocity fields are extracted and compared to the experimental data, as illustrated in figure 6. As expected, while the longitudinal propagation is non dispersive, the transverse shear wave propagation is dispersive. Both curves are in good agreement with the experimental data. This helps to show that the previously described transverse isotropic viscoelastic Green's

<footnote><sup>a</sup> *In vivo* data on five healthy volunteers from Gennisson *et al* (2010).</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** The experimental axial (Z) velocity fields in both longitudinal (a) and transverse (b) directions (superposed to B-mode images) are represented at the same sampling times for SWE in the *in vivo* forearm muscle. The spatio-temporal experimental imaging of the shear wave in the longitudinal (c) and transverse (d) directions confirms the anisotropic profile by a significant variation of the slope. The shear wave speed variations (e) are extracted from these experimental measurements (green dots, including standard deviation) as a function of the probe / fibers angle θ and compared to the anisotropic data (blue dots and line) from the numerical simulations.</figcaption></figure>

<table><caption>**Table 2.** Experimental SWE data used to determine the numerical parameters as input for the shear wave simulations in the *in vivo* forearm muscle.</caption><thead><tr><th></th><th>c<sub>SV</sub><br>[m s<sup>-1</sup>]</th><th>c<sub>SH</sub><br>[m s<sup>-1</sup>]</th><th>c<sub>P</sub><br>[m s<sup>-1</sup>]</th><th>η<sub>SV</sub><br>[Pa s]</th><th>η<sub>SH</sub><br>[Pa s]</th><th>η<sub>P</sub><br>[Pa s]</th><th>ρ<br>[g cm<sup>-3</sup>]</th></tr></thead><tbody><tr><td>Experimental measurements</td><td>6.69 ± 0.23</td><td>2.16 ± 0.19</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Numerical parameters</td><td>6.50</td><td>2.16</td><td>1500</td><td>0.65 ± 0.10<sup>a</sup></td><td>0.92 ± 0.06<sup>a</sup></td><td>0.1</td><td>1000</td></tr></tbody></table>

and the shear wave propagation direction. The highest speed is well correlated with the fiber orientation for both approaches. Whatever the fiber orientation is, the experimental results and the viscoelastic Green's function simulations are in very good agreement. This helps to validate the elastic behavior of the shear wave propagation in the anisotropic simulations.

The capability of the numerical method proposed in this study to simulate the anisotropic viscoelastic behavior of fibrous soft tissues is here compared to experimental data on the *in vivo* relaxed biceps from Gennisson *et al* (2010). The experimental values from this study, as well as the numerical input parameters for the simulations in the biceps muscle, are summarized in table 3.

For θ = 0° and θ = 90°, the numerical longitudinal and transverse shear wave velocity fields are extracted and compared to the experimental data, as illustrated in figure 6. As expected, while the longitudinal propagation is non dispersive, the transverse shear wave propagation is dispersive. Both curves are in good agreement with the experimental data. This helps to show that the previously described transverse isotropic viscoelastic Green's

<footnote><sup>a</sup> *In vivo* data on five healthy volunteers from Gennisson *et al* (2010).</footnote>
<table><caption>**Table 3.** Experimental SWE data used to determine the numerical parameters as input for the shear wave simulations in the viscoelastic biceps.</caption><thead><tr><th></th><th>μ<sub>SV</sub><br>[kPa]</th><th>μ<sub>SH</sub><br>[kPa]</th><th>c<sub>P</sub><br>[m s<sup>-1</sup>]</th><th>η<sub>SV</sub><br>[Pas]</th><th>η<sub>SH</sub><br>[Pas]</th><th>η<sub>P</sub><br>[Pas]</th><th>ρ<br>[g cm<sup>-3</sup>]</th></tr></thead><tbody><tr><td>Experimental<br>measurements</td><td>5.86 ± 0.20<sup>a</sup></td><td>1.58 ± 0.15<sup>a</sup></td><td>—</td><td>0.65 ± 0.10<sup>a</sup></td><td>0.92 ± 0.06<sup>a</sup></td><td>—</td><td>—</td></tr><tr><td>Numerical<br>parameters</td><td>5.86</td><td>1.58</td><td>1500</td><td>0.65</td><td>0.92</td><td>0.1</td><td>1000</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Velocity field of the shear wave velocity along the fibers and perpendicularly to the fibers from numerical simulations and experimental SSI acquisitions* in five *in vivo* viscoelastic biceps, using the values presented in table 3. The minimal, mean and maximal numerical curves are calculated using the mean as well as mean ± standard deviation values from Gennisson *et al* (2010) as input parameters.</figcaption></figure>

function can be used to simulate not only the elastic but also the viscous anisotropic behavior of *in vivo* fibrous tissue.

Finally, while the experiments in the phantom contribute towards validating the approach in a well-known, perfectly transverse isotropic, homogeneous medium, the tests in the forearm and biceps muscle illustrate the potential of combining numerical simulations and analytical formalism to determine the transverse isotropic shear wave profile and the anisotropic dispersion profile, respectively.

## 4. Discussion and conclusions

Understanding shear wave propagation in fibrous tissue is one of the most challenging issues for elastography. Through the comparison with phantom-based as well as *in vivo* experimental data on muscle, this paper proposes a simple numerical tool for the simulation of shear wave elastography in fibrous viscoelastic medium, with numerous perspectives in terms of mechanical anisotropy and viscosity measurements on *in vivo* fibrous soft tissues.

As illustrated in figure 4(a), the shear wave speed profile is similar in simulations and experiments and is significantly faster in the longitudinal direction than in the transverse direction. Whatever the fiber orientation is, the experimental results and the viscoelastic Green's

<footnote>ª *In vivo* data on five healthy volunteers from Gennisson *et al* 2010).</footnote>
<table><caption>**Table 3.** Experimental SWE data used to determine the numerical parameters as input for the shear wave simulations in the viscoelastic biceps.</caption><thead><tr><th></th><th>μ<sub>SV</sub><br>[kPa]</th><th>μ<sub>SH</sub><br>[kPa]</th><th>c<sub>P</sub><br>[m s<sup>-1</sup>]</th><th>η<sub>SV</sub><br>[Pas]</th><th>η<sub>SH</sub><br>[Pas]</th><th>η<sub>P</sub><br>[Pas]</th><th>ρ<br>[g cm<sup>3</sup>]</th></tr></thead><tbody><tr><td>Experimental<br>measurements</td><td>5.86 ± 0.20<sup>a</sup></td><td>1.58 ± 0.15<sup>a</sup></td><td>—</td><td>0.65 ± 0.10<sup>a</sup></td><td>0.92 ± 0.06<sup>a</sup></td><td>—</td><td>—</td></tr><tr><td>Numerical<br>parameters</td><td>5.86</td><td>1.58</td><td>1500</td><td>0.65</td><td>0.92</td><td>0.1</td><td>1000</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Velocity field of the shear wave velocity along the fibers and perpendicularly to the fibers from numerical simulations and experimental SSI acquisitions* in five *in vivo* viscoelastic biceps, using the values presented in table 3. The minimal, mean and maximal numerical curves are calculated using the mean as well as mean ± standard deviation values from Gennisson *et al* (2010) as input parameters.</figcaption></figure>

function can be used to simulate not only the elastic but also the viscous anisotropic behavior of *in vivo* fibrous tissue.

Finally, while the experiments in the phantom contribute towards validating the approach in a well-known, perfectly transverse isotropic, homogeneous medium, the tests in the forearm and biceps muscle illustrate the potential of combining numerical simulations and analytical formalism to determine the transverse isotropic shear wave profile and the anisotropic dispersion profile, respectively.

## 4. Discussion and conclusions

Understanding shear wave propagation in fibrous tissue is one of the most challenging issues for elastography. Through the comparison with phantom-based as well as *in vivo* experimental data on muscle, this paper proposes a simple numerical tool for the simulation of shear wave elastography in fibrous viscoelastic medium, with numerous perspectives in terms of mechanical anisotropy and viscosity measurements on *in vivo* fibrous soft tissues.

As illustrated in figure 4(a), the shear wave speed profile is similar in simulations and experiments and is significantly faster in the longitudinal direction than in the transverse direction. Whatever the fiber orientation is, the experimental results and the viscoelastic Green's

<footnote>ª *In vivo* data on five healthy volunteers from Gennisson *et al* 2010).</footnote>
function simulations are in very good agreement. For both experimental and numerical results, the shear wave velocity is gradually dependent on the $\theta$ angle between the fibers and the shear wave propagation direction, which is in good agreement with the assumption of mechanical transverse isotropy. These variations of the shear wave speed can be easily understood by observing the shear wave profile propagating over different time steps (figure 4(b)). The shear waves are propagating from the center of the radiation force field and significantly faster along than orthogonally to the fiber direction. It can also be observed that shear wave amplitude is significantly higher in the fiber main direction than in the orthogonal one. However, due to technological limitations, the actual probes give the possibility of visualizing simultaneously the here-described shear wave profile in only one bi-dimensional plane. The numerical simulations are still one of the best ways to illustrate this phenomenon. Recent development of 3D ultrafast ultrasound imaging *in vivo* would give a promising perspective for the complete anisotropic shear wave field experimental visualization (Provost et al 2014).

In a similar way to finite differences or finite elements simulations, the analytical approach using Green's formalism gives the opportunity of understanding shear wave propagation in a complex medium. As observed in this study, the mechanical properties measured by bi-dimensional SSI in a fibrous tissue is greatly dependent on the observation direction of the shear wave propagation, i.e. on the orientation of the ultrasound probe, as shown in figures 3 and 5. However, it has been shown in this study that the numerical simulations are powerful tools to predict or interpret experimental SSI measurements. It has to be noticed that the analytical formalism is subject to an *a priori* knowledge of some the mechanical of the soft tissue, mainly the longitudinal and transverse stiffness and viscosity coefficients.

Both experimental and numerical profiles of the shear wave speed in the function of the angle illustrated in figures 5(a) and 6(b) are similar to the analytical expressions for an homogeneous transverse isotropic medium (Royer and Dieulesaint 2000, Wang et al 2013). It is also important to show that the velocity of the longitudinal wave as well as the profile of the propagation velocity for a transverse (shear) wave in the function of the angle of incidence $\theta$ (figures 4(a) and 5(e)) have a similar shape to the analytical expression recently described by Rudenko and Sarvazyan (Rudenko and Sarvazyan). As described by the authors, the velocity of the transverse wave polarized orthogonally to the plane of incidence as a function of angle $\theta$ between fiber axis and wave vector can be determined by two independent moduli $a_{44}$ and $a_{66}$, similarly to the here-developed approach for the Green's formalism. Like Rudenko and Sarvazyan using the fourth rank viscosity tensor, both descriptions for the dissipative properties are determined by two independent components. The present approach provides, in addition, more detailed aspects on wave dispersion and implies specific parameters of the ultrasonic probe.

An important limitation of our method concerns the assumptions of the mechanical properties of the tissue, which is supposed to be semi-infinite, homogeneous and perfectly transverse isotropic. While the modeled shear wave speeds are in very good agreement for $\theta = 0°$ and $\theta = 90°$, some divergences can be observed *in vivo* in-between angles (such as $\theta = 45°$), as shown in figure 5(e). This remark illustrates not only the limit in considering the muscles as perfectly transverse isotropic, but also in assuming the tissue as homogeneous and semi-infinite in the Green's formalism, while muscles are both multi-layered and multi-directional fibrous medium. One of the major perspectives for this study would consist of the complexification in both the anisotropic modelization of the muscular soft tissues and its multi-layer and multi-phase finite medium assumption.

It has to be noticed that use of the Green's formalism for viscoelastic simulations requires the description of the medium using Voigt's model. This implies to use a 'zero-viscosity' velocities (elastic part of the Voigt models), which can be considered as the group velocities in
function simulations are in very good agreement. For both experimental and numerical results, the shear wave velocity is gradually dependent on the $\theta$ angle between the fibers and the shear wave propagation direction, which is in good agreement with the assumption of mechanical transverse isotropy. These variations of the shear wave speed can be easily understood by observing the shear wave profile propagating over different time steps (figure 4(b)). The shear waves are propagating from the center of the radiation force field and significantly faster along than orthogonally to the fiber direction. It can also be observed that shear wave amplitude is significantly higher in the fiber main direction than in the orthogonal one. However, due to technological limitations, the actual probes give the possibility of visualizing simultaneously the here-described shear wave profile in only one bi-dimensional plane. The numerical simulations are still one of the best ways to illustrate this phenomenon. Recent development of 3D ultrafast ultrasound imaging *in vivo* would give a promising perspective for the complete anisotropic shear wave field experimental visualization (Provost et al 2014).

In a similar way to finite differences or finite elements simulations, the analytical approach using Green's formalism gives the opportunity of understanding shear wave propagation in a complex medium. As observed in this study, the mechanical properties measured by bi-dimensional SSI in a fibrous tissue is greatly dependent on the observation direction of the shear wave propagation, i.e. on the orientation of the ultrasound probe, as shown in figures 3 and 5. However, it has been shown in this study that the numerical simulations are powerful tools to predict or interpret experimental SSI measurements. It has to be noticed that the analytical formalism is subject to an *a priori* knowledge of some the mechanical of the soft tissue, mainly the longitudinal and transverse stiffness and viscosity coefficients.

Both experimental and numerical profiles of the shear wave speed in the function of the angle illustrated in figures 5(a) and 6(b) are similar to the analytical expressions for an homogeneous transverse isotropic medium (Royer and Dieulesaint 2000, Wang et al 2013). It is also important to show that the velocity of the longitudinal wave as well as the profile of the propagation velocity for a transverse (shear) wave in the function of the angle of incidence $\theta$ (figures 4(a) and 5(e)) have a similar shape to the analytical expression recently described by Rudenko and Sarvazyan (Rudenko and Sarvazyan). As described by the authors, the velocity of the transverse wave polarized orthogonally to the plane of incidence as a function of angle $\theta$ between fiber axis and wave vector can be determined by two independent moduli $a_{44}$ and $a_{66}$, similarly to the here-developed approach for the Green's formalism. Like Rudenko and Sarvazyan using the fourth rank viscosity tensor, both descriptions for the dissipative properties are determined by two independent components. The present approach provides, in addition, more detailed aspects on wave dispersion and implies specific parameters of the ultrasonic probe.

An important limitation of our method concerns the assumptions of the mechanical properties of the tissue, which is supposed to be semi-infinite, homogeneous and perfectly transverse isotropic. While the modeled shear wave speeds are in very good agreement for $\theta = 0°$ and $\theta = 90°$, some divergences can be observed *in vivo* in-between angles (such as $\theta = 45°$), as shown in figure 5(e). This remark illustrates not only the limit in considering the muscles as perfectly transverse isotropic, but also in assuming the tissue as homogeneous and semi-infinite in the Green's formalism, while muscles are both multi-layered and multi-directional fibrous medium. One of the major perspectives for this study would consist of the complexification in both the anisotropic modelization of the muscular soft tissues and its multi-layer and multi-phase finite medium assumption.

It has to be noticed that use of the Green's formalism for viscoelastic simulations requires the description of the medium using Voigt's model. This implies to use a 'zero-viscosity' velocities (elastic part of the Voigt models), which can be considered as the group velocities in
only one case of negligible viscosity. This point is important because using the group velocities values as input for the Green's function will be completely correct for the low-viscosity media. It corresponds to the configurations and assumptions used in the tests in the PVA phantom as well as in the forearm muscle. However, this assumption is useless when considering at the same time a viscoelastic medium and the dispersion curves, as for our tests in the biceps muscle. It comes from the use of the Voigt's model in both experimental and analytical approaches. The use of the group velocity is no more useful and the description of the material using the parameters of the Voigt's model will give the best results. Again, this aspect illustrates the necessity of a good *a priori* knowledge of the mechanical properties of the medium for the use of the formalism proposed in this study.

We have shown that *in vivo* application is possible, for both anisotropic shear wave velocity and dispersion profile simulation. It plays a major role in the design and the optimization of new ultrasonic transducers dedicated to remote elastography in isotropic or fibrous tissue, such as muscle or myocardium. One of the main advantages of the numerical approach is to have a relatively low computational requirement in comparison with finite element and differences approaches. By coupling the Field II and the shear wave Green's formalism simulations, the total pressure field generated by a specific ultrasonic probe is considered as source of wave propagation. It gives the possibility of investigating the shear wave propagation dependence not only with the mechanical properties of the targeted tissue (stiffness, viscosity, transverse isotropy, compressibility, and density) but also with the probe used for elastography. This remark has great perspective for the design of new probes dedicated for SSI.

In this study, the description of wave propagation has been focused on transverse isotropy with fiber direction systematically parallel to the ultrasound beam focal plane. This configuration is certainly the closest to the most common one used for cardiac or muscular SSI. However, other configurations and local fibers-crossings would not be modeled in a satisfactory way by in-plane fibers. In this way, the extension of the Green's formalism to orthotropic medium would be one of the first perspectives of prolongation for this study. Elsewhere, the feasibility of simulating shear waves in 3D anisotropic space has been shown. However, the actual SSI technique does not allow us to visualize a shear wave front out from the 2D imaging plane (with shear waves polarized in the axial direction). This limitation would especially concern SSI with either out-of focal plane fibers or fibers crossing areas. The main suggestion to solve it would consist of the development of a 3D probe for SSI.

As shown in figure 6, the numerical approach is thought to reproduce the *in vivo* dispersion in both longitudinal and transverse directions and, by extension, the complete shear wave dispersion in a viscoelastic fibrous tissue. In isotropic as well as anisotropic mediums, shear wave attenuation is not only due to viscosity but also to diffraction, being consequently a source of bias. The use of numerical simulations, such as the analytical approach developed in this study for anisotropic SSI, allows us to separate these two contributions and opens the perspective to experimentally extract viscous losses in anisotropic mediums. In this way, the comparison between numerical and experimental shear wave attenuation in transverse isotropic viscoelastic medium is one of the continuations of this study.

The simulations have been designed and experimentally validated to model the shear wave propagation in a viscoelastic transverse isotropic medium. The results show the feasibility of numerically combining both focal spot simulation and tridimensional displacement field calculation. The comparison with *in vivo* experimental data in a viscoelastic transverse isotropic medium shows the efficiency of the method to model shear wave elasticity imaging a replicated fibrous biological soft tissue using a specific ultrasonic probe.
only one case of negligible viscosity. This point is important because using the group velocities values as input for the Green's function will be completely correct for the low-viscosity media. It corresponds to the configurations and assumptions used in the tests in the PVA phantom as well as in the forearm muscle. However, this assumption is useless when considering at the same time a viscoelastic medium and the dispersion curves, as for our tests in the biceps muscle. It comes from the use of the Voigt's model in both experimental and analytical approaches. The use of the group velocity is no more useful and the description of the material using the parameters of the Voigt's model will give the best results. Again, this aspect illustrates the necessity of a good *a priori* knowledge of the mechanical properties of the medium for the use of the formalism proposed in this study.

We have shown that *in vivo* application is possible, for both anisotropic shear wave velocity and dispersion profile simulation. It plays a major role in the design and the optimization of new ultrasonic transducers dedicated to remote elastography in isotropic or fibrous tissue, such as muscle or myocardium. One of the main advantages of the numerical approach is to have a relatively low computational requirement in comparison with finite element and differences approaches. By coupling the Field II and the shear wave Green's formalism simulations, the total pressure field generated by a specific ultrasonic probe is considered as source of wave propagation. It gives the possibility of investigating the shear wave propagation dependence not only with the mechanical properties of the targeted tissue (stiffness, viscosity, transverse isotropy, compressibility, and density) but also with the probe used for elastography. This remark has great perspective for the design of new probes dedicated for SSI.

In this study, the description of wave propagation has been focused on transverse isotropy with fiber direction systematically parallel to the ultrasound beam focal plane. This configuration is certainly the closest to the most common one used for cardiac or muscular SSI. However, other configurations and local fibers-crossings would not be modeled in a satisfactory way by in-plane fibers. In this way, the extension of the Green's formalism to orthotropic medium would be one of the first perspectives of prolongation for this study. Elsewhere, the feasibility of simulating shear waves in 3D anisotropic space has been shown. However, the actual SSI technique does not allow us to visualize a shear wave front out from the 2D imaging plane (with shear waves polarized in the axial direction). This limitation would especially concern SSI with either out-of focal plane fibers or fibers crossing areas. The main suggestion to solve it would consist of the development of a 3D probe for SSI.

As shown in figure 6, the numerical approach is thought to reproduce the *in vivo* dispersion in both longitudinal and transverse directions and, by extension, the complete shear wave dispersion in a viscoelastic fibrous tissue. In isotropic as well as anisotropic mediums, shear wave attenuation is not only due to viscosity but also to diffraction, being consequently a source of bias. The use of numerical simulations, such as the analytical approach developed in this study for anisotropic SSI, allows us to separate these two contributions and opens the perspective to experimentally extract viscous losses in anisotropic mediums. In this way, the comparison between numerical and experimental shear wave attenuation in transverse isotropic viscoelastic medium is one of the continuations of this study.

The simulations have been designed and experimentally validated to model the shear wave propagation in a viscoelastic transverse isotropic medium. The results show the feasibility of numerically combining both focal spot simulation and tridimensional displacement field calculation. The comparison with *in vivo* experimental data in a viscoelastic transverse isotropic medium shows the efficiency of the method to model shear wave elasticity imaging a replicated fibrous biological soft tissue using a specific ultrasonic probe.
## Acknowledgments

The research leading to these results received funding from the European Research Council under the European Union's Seventh Framework Program (FP/2007–2013) / ERC Grant Agreement n°311025.

## References

Aki L and Richards P G 1980 *Quantitative Seismology, Theory and Methods* vol 1 (New York: Freeman) chapter 4

Aubry S, Risson J R, Kastler A, Barbier-Brion B, Siliman G, Runge M and Kastler B 2013 Biomechanical properties of the calcaneal tendon *in vivo* assessed by transient shear wave elastography *Skeletal Radiol* **42** 1143–50

Barrière C and Royer D 2001a Measurement of large ultrasonic displacements with an optical probe *IEEE Ultrasonics Symp. Proc.* vol 2 pp 1105–8

Barrière C and Royer D 2001b Optical measurement of large transient mechanical displacements *Appl. Phys. Lett.* **79** 878–80

Bavu E *et al* 2011 Noninvasive *in vivo* liver fibrosis evaluation using supersonic shear imaging: a clinical study on 113 hepatitis C virus patients *Ultrason. Med. Biol.* **37** 1361–73

Bercoff J, Tanter M and Fink M 2004a Supersonic shear imaging: a new technique for soft tissue elasticity mapping *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **51** 396–409

Bercoff J, Muller M, Tanter M and Fink M 2004b The role of viscosity in the impulse diffraction field of elastic waves induced by the acoustic radiation force *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **51** 1523–36

Berg W *et al* 2012 Shear-wave elastography improves the specificity of breast US: the BE1 multinational study of 939 masses *Radiology* **262** 435–49

Bernal M, Nenadic I, Urban M W and Greenleaf J F 2011 Material property estimation for tubes and arteries using ultrasound radiation force and analysis of propagating modes *Acoust. Soc. Am.* **129** 1344–54

Brum J, Bernal M, Gennisson J-L and Tanter M 2013 *In vivo* evaluation of the elastic anisotropy of human achilles tendon using shear wave dispersion analysis *Phys. Med. Biol.* **55** 505–23

Buchwald V T 1959 Elastic waves in anisotropic media *Proc. R. Soc. A* **253** 563–80

Burridge R 1967 The singularity on the plane lids of the wave surface of elastic media with cubic symmetry *Q. J. Mech. Appl. Math.* **20** 41–56

Catheline S, Gennisson J-L, Delon G, Fink M, Sinkus R, Abouelkaram S and Culioli J 2004 Measurement of viscoelastic properties of homogeneous soft solid using transient elastography: an inverse problem approach *J. Acoust. Soc. Am.* **116** 3734–41

Chatelin S, Bernal M, Deffieux T, Papadacci C, Flaud P, Nahas A, Boccara C, Gennisson J-L, Tanter M and Pernot M 2014 Anisotropic polyvinyl alcohol hydrogel phantom for shear wave elastography in fibrous biological soft tissue: a multimodality characterization *Phys. Med. Biol.* **59** 6923–40

Chen S, Fatemi M and Greenleaf J F 2004 Quantifying elasticity and viscosity from measurement of shear wave speed dispersion *Acoust. Soc. Am.* **115** 2781–5

Cosgrove D O, Berg W A, Doré C J, Skyba D M, Henry J-P, Gay J, Cohen-Bacrie C and BE1 Study Group 2011 Shear wave elastography for breast masses is highly reproducible *Eur. Radiol.* **22** 1023–32

Culjat M O, Goldenberg D, Tewari P and Singh R S 2010 A review of tissue substitutes for ultrasound imaging *Ultrasound Med. Biol.* **36** 861–73

Deffieux T, Gennisson J-L, Bercoff J and Tanter M 2011 On the effects of reflected waves in transient shear wave elastography *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **58** 2032–5

Deffieux T, Montaldo G, Tanter M and Fink M 2009 Shear wave spectroscopy for *in vivo* quantification of human soft tissues visco-elasticity *IEEE Trans. Med. Imaging* **28** 313–22

Fatemi M and Greenleaf J 1998 Ultrasound-stimulated vibro-acoustic spectrography *Science* **280** 82–5

Gennisson J L, Catheline S, Chaffaï S and Fink M 2003 Transient elastography in anisotropic medium: application to the measurement of slow and fast shear wave speeds in muscles *J. Acoust. Soc. Am.* **114** 536–41
## Acknowledgments

The research leading to these results received funding from the European Research Council under the European Union's Seventh Framework Program (FP/2007–2013) / ERC Grant Agreement n°311025.

## References

Aki L and Richards P G 1980 *Quantitative Seismology, Theory and Methods* vol 1 (New York: Freeman) chapter 4

Aubry S, Risson J R, Kastler A, Barbier-Brion B, Siliman G, Runge M and Kastler B 2013 Biomechanical properties of the calcaneal tendon *in vivo* assessed by transient shear wave elastography *Skeletal Radiol* **42** 1143–50

Barrière C and Royer D 2001a Measurement of large ultrasonic displacements with an optical probe *IEEE Ultrasonics Symp. Proc.* vol 2 pp 1105–8

Barrière C and Royer D 2001b Optical measurement of large transient mechanical displacements *Appl. Phys. Lett.* **79** 878–80

Bavu E *et al* 2011 Noninvasive *in vivo* liver fibrosis evaluation using supersonic shear imaging: a clinical study on 113 hepatitis C virus patients *Ultrason. Med. Biol.* **37** 1361–73

Bercoff J, Tanter M and Fink M 2004a Supersonic shear imaging: a new technique for soft tissue elasticity mapping *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **51** 396–409

Bercoff J, Muller M, Tanter M and Fink M 2004b The role of viscosity in the impulse diffraction field of elastic waves induced by the acoustic radiation force *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **51** 1523–36

Berg W *et al* 2012 Shear-wave elastography improves the specificity of breast US: the BE1 multinational study of 939 masses *Radiology* **262** 435–49

Bernal M, Nenadic I, Urban M W and Greenleaf J F 2011 Material property estimation for tubes and arteries using ultrasound radiation force and analysis of propagating modes *Acoust. Soc. Am.* **129** 1344–54

Brum J, Bernal M, Gennisson J-L and Tanter M 2013 *In vivo* evaluation of the elastic anisotropy of human achilles tendon using shear wave dispersion analysis *Phys. Med. Biol.* **55** 505–23

Buchwald V T 1959 Elastic waves in anisotropic media *Proc. R. Soc. A* **253** 563–80

Burridge R 1967 The singularity on the plane lids of the wave surface of elastic media with cubic symmetry *Q. J. Mech. Appl. Math.* **20** 41–56

Catheline S, Gennisson J-L, Delon G, Fink M, Sinkus R, Abouelkaram S and Culioli J 2004 Measurement of viscoelastic properties of homogeneous soft solid using transient elastography: an inverse problem approach *J. Acoust. Soc. Am.* **116** 3734–41

Chatelin S, Bernal M, Deffieux T, Papadacci C, Flaud P, Nahas A, Boccara C, Gennisson J-L, Tanter M and Pernot M 2014 Anisotropic polyvinyl alcohol hydrogel phantom for shear wave elastography in fibrous biological soft tissue: a multimodality characterization *Phys. Med. Biol.* **59** 6923–40

Chen S, Fatemi M and Greenleaf J F 2004 Quantifying elasticity and viscosity from measurement of shear wave speed dispersion *Acoust. Soc. Am.* **115** 2781–5

Cosgrove D O, Berg W A, Doré C J, Skyba D M, Henry J-P, Gay J, Cohen-Bacrie C and BE1 Study Group 2011 Shear wave elastography for breast masses is highly reproducible *Eur. Radiol.* **22** 1023–32

Culjat M O, Goldenberg D, Tewari P and Singh R S 2010 A review of tissue substitutes for ultrasound imaging *Ultrasound Med. Biol.* **36** 861–73

Deffieux T, Gennisson J-L, Bercoff J and Tanter M 2011 On the effects of reflected waves in transient shear wave elastography *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **58** 2032–5

Deffieux T, Montaldo G, Tanter M and Fink M 2009 Shear wave spectroscopy for *in vivo* quantification of human soft tissues visco-elasticity *IEEE Trans. Med. Imaging* **28** 313–22

Fatemi M and Greenleaf J 1998 Ultrasound-stimulated vibro-acoustic spectrography *Science* **280** 82–5

Gennisson J L, Catheline S, Chaffaï S and Fink M 2003 Transient elastography in anisotropic medium: application to the measurement of slow and fast shear wave speeds in muscles *J. Acoust. Soc. Am.* **114** 536–41
Gennisson J-L, Deffieux T, Macé E, Montaldo G, Fink M and Tanter M 2010 Viscoelastic and anisotropic mechanical properties of *in vivo* muscle tissue assessed by supersonic shear imaging *Ultrasound Med. Biol.* **36** 789–801

Hoffmeister B K, Handley S M, Wickline S A and Miller J G 1996 Ultrasonic determination of the anisotropy of Young's modulus of fixed tendon and fixed myocardium *J. Acoust. Soc. Am.* **100** 3933–40

Jensen J A 1996 Field: a program for simulating ultrasound systems *Med. Biol. Eng. Comput.* **34** 351–3 (10th Nordic-Baltic Conf. on Biomedical Imaging)

Jensen J A and Svendsen N B 1992 Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **39** 262–7

Lee W N, Pernot M, Couade M, Messas E, Bruneval P, Bel A and Tanter M 2012 Mapping myocardial fiber orientation using echocardiography-based shear wave imaging *IEEE Trans. Med. Imaging* **31** 554–62

Millon L E, Mohammadi H and Wan W K 2006 Anisotropic polyvinyl alcohol hydrogel for cardiovascular applications *J. Biomed. Mater. Res. Part B: Appl. Biomater.* **79** 305–11

Norris A N 1994 Dynamic Green functions in anisotropic piezoelectric, thermoelastic and poroelastic solids *Proc. R. Soc. A* **447** 175–88

Palmeri M L, Wang M H, Dahl J J, Frinkley K D and Nightingale K R 2008 Quantifying hepatic shear modulus *in vivo* using acoustic radiation force *Ultrasound Med. Biol.* **34** 546–58

Provost J, Papadacci C, Arango J E, Imbault M, Fink M, Gennisson J-L, Tanter M and Pernot M 2014 3D ultrafast ultrasound imaging *in vivo* *Phys. Med. Biol.* **59** L1

Rouze N C, Wang M H, Palmeri M L and Nightingale K R 2013 Finite element modeling of impulsive excitation and shear wave propagation in an incompressible, transversely isotropic medium *J. Biomech.* **46** 2761–8

Royer D and Dieulesaint E 2000 *Elastic Waves in Solids* vol 1 (Heidelberg: Springer)

Royer D, Dubois N and Benoist P 1992 Optical probing of acoustic fields—application to the ultrasonic testing of steam generator tubes *IEEE Ultrason. Symp. Proc.* vol 2 pp 805–8

Rudenko O V and Sarvazyan A P 2014 Wave anisotropy of shear viscosity and elasticity *Acoust. Phys.* **60** 710–8

Sandrin L *et al* 2003 Transient elastography: a new noninvasive method for assessment of hepatic fibrosis *Ultrasound Med. Biol.* **29** 1705–13

Sandrin L, Cassereau D and Fink M 2004 The role of the coupling term in transient elastography *J. Acoust. Soc. Am.* **115** 73–83

Sarvazyan A P, Rudenko O V, Swanson S D, Fowlkes J B and Emelianov S Y 1998 Shear wave elasticity imaging—a new ultrasonic technology of medical diagnostic *Ultrasound Med. Biol.* **20** 1419–36

Shcherbakova D A *et al* 2014 Supersonic shear wave imaging to assess arterial nonlinear behavior and anisotropy: proof of principle via *ex vivo* testing of the horse aorta *Adv. Mech. Eng.* **6**

Tanter M, Bercoff J, Athanasiou A, Deffieux T, Gennisson J-L, Montaldo G, Muller M, Tardivon A and Fink M 2008 Quantitative assessment of breast lesion viscoelasticity: initial clinical results using supersonic shear imaging *Ultrasound Med. Biol.* **34** 1373–86

Tanter M and Fink M 2014 Ultrafast imaging in biomedical ultrasound *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **61** 102–19

Torr G R 1984 The acoustic radiation force *Am. J. Phys.* **52** 402–8

Vavryčuk V 2001 Exact elastodynamic Green functions for simple types of anisotropy derived from higher-order ray theory *Stud. Geophys. Geod.* **45** 67–84

Vavryčuk V 2007 Asymptotic Green's function in homogeneous anisotropic viscoelastic media *Proc. R. Soc. Am.* **463** 2689–707

Vavryčuk V 2008 Velocity, attenuation, and quality factor in anisotropic viscoelastic media: a perturbation approach *Geophysics* **73** 63–73

Wan W, Millon L E and Mohammadi H 2009 Anisotropic hydrogels *United States Patent Application Publication* US 2009/0214623 A1

Wang C Y and Achenbach J D 1994 Elastodynamic fundamental solutions for anisotropic solids *Geophys. J. Int.* **118** 384–92

Wang M, Byram B, Palmeri M, Rouze N and Nightingale K 2013 Imaging transverse isotropic properties of muscle by monitoring acoustic radiation force induced shear waves using a 2D matrix ultrasound array *IEEE Trans. Med. Imaging* **32** 1671–84

Willis J R 1980 A polarization approach to the scattering of elastic waves. I. Scattering by a single inclusion *J. Mech. Phys. Solids* **28** 287–305
Gennisson J-L, Deffieux T, Macé E, Montaldo G, Fink M and Tanter M 2010 Viscoelastic and anisotropic mechanical properties of *in vivo* muscle tissue assessed by supersonic shear imaging *Ultrasound Med. Biol.* **36** 789–801

Hoffmeister B K, Handley S M, Wickline S A and Miller J G 1996 Ultrasonic determination of the anisotropy of Young's modulus of fixed tendon and fixed myocardium *J. Acoust. Soc. Am.* **100** 3933–40

Jensen J A 1996 Field: a program for simulating ultrasound systems *Med. Biol. Eng. Comput.* **34** 351–3 (10th Nordic-Baltic Conf. on Biomedical Imaging)

Jensen J A and Svendsen N B 1992 Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **39** 262–7

Lee W N, Pernot M, Couade M, Messas E, Bruneval P, Bel A and Tanter M 2012 Mapping myocardial fiber orientation using echocardiography-based shear wave imaging *IEEE Trans. Med. Imaging* **31** 554–62

Millon L E, Mohammadi H and Wan W K 2006 Anisotropic polyvinyl alcohol hydrogel for cardiovascular applications *J. Biomed. Mater. Res. Part B: Appl. Biomater.* **79** 305–11

Norris A N 1994 Dynamic Green functions in anisotropic piezoelectric, thermoelastic and poroelastic solids *Proc. R. Soc. A* **447** 175–88

Palmeri M L, Wang M H, Dahl J J, Frinkley K D and Nightingale K R 2008 Quantifying hepatic shear modulus *in vivo* using acoustic radiation force *Ultrasound Med. Biol.* **34** 546–58

Provost J, Papadacci C, Arango J E, Imbault M, Fink M, Gennisson J-L, Tanter M and Pernot M 2014 3D ultrafast ultrasound imaging *in vivo* *Phys. Med. Biol.* **59** L1

Rouze N C, Wang M H, Palmeri M L and Nightingale K R 2013 Finite element modeling of impulsive excitation and shear wave propagation in an incompressible, transversely isotropic medium *J. Biomech.* **46** 2761–8

Royer D and Dieulesaint E 2000 *Elastic Waves in Solids* vol 1 (Heidelberg: Springer)

Royer D, Dubois N and Benoist P 1992 Optical probing of acoustic fields—application to the ultrasonic testing of steam generator tubes *IEEE Ultrason. Symp. Proc.* vol 2 pp 805–8

Rudenko O V and Sarvazyan A P 2014 Wave anisotropy of shear viscosity and elasticity *Acoust. Phys.* **60** 710–8

Sandrin L *et al* 2003 Transient elastography: a new noninvasive method for assessment of hepatic fibrosis *Ultrasound Med. Biol.* **29** 1705–13

Sandrin L, Cassereau D and Fink M 2004 The role of the coupling term in transient elastography *J. Acoust. Soc. Am.* **115** 73–83

Sarvazyan A P, Rudenko O V, Swanson S D, Fowlkes J B and Emelianov S Y 1998 Shear wave elasticity imaging—a new ultrasonic technology of medical diagnostic *Ultrasound Med. Biol.* **20** 1419–36

Shcherbakova D A *et al* 2014 Supersonic shear wave imaging to assess arterial nonlinear behavior and anisotropy: proof of principle via *ex vivo* testing of the horse aorta *Adv. Mech. Eng.* **6**

Tanter M, Bercoff J, Athanasiou A, Deffieux T, Gennisson J-L, Montaldo G, Muller M, Tardivon A and Fink M 2008 Quantitative assessment of breast lesion viscoelasticity: initial clinical results using supersonic shear imaging *Ultrasound Med. Biol.* **34** 1373–86

Tanter M and Fink M 2014 Ultrafast imaging in biomedical ultrasound *IEEE Trans. Ultrason. Ferroelectr. Freq. Control* **61** 102–19

Torr G R 1984 The acoustic radiation force *Am. J. Phys.* **52** 402–8

Vavryčuk V 2001 Exact elastodynamic Green functions for simple types of anisotropy derived from higher-order ray theory *Stud. Geophys. Geod.* **45** 67–84

Vavryčuk V 2007 Asymptotic Green's function in homogeneous anisotropic viscoelastic media *Proc. R. Soc. Am.* **463** 2689–707

Vavryčuk V 2008 Velocity, attenuation, and quality factor in anisotropic viscoelastic media: a perturbation approach *Geophysics* **73** 63–73

Wan W, Millon L E and Mohammadi H 2009 Anisotropic hydrogels *United States Patent Application Publication* US 2009/0214623 A1

Wang C Y and Achenbach J D 1994 Elastodynamic fundamental solutions for anisotropic solids *Geophys. J. Int.* **118** 384–92

Wang M, Byram B, Palmeri M, Rouze N and Nightingale K 2013 Imaging transverse isotropic properties of muscle by monitoring acoustic radiation force induced shear waves using a 2D matrix ultrasound array *IEEE Trans. Med. Imaging* **32** 1671–84

Willis J R 1980 A polarization approach to the scattering of elastic waves. I. Scattering by a single inclusion *J. Mech. Phys. Solids* **28** 287–305