Estimation of material parameters from slow and fast shear waves in an incompressible, transversely isotropic material

Abstract

This paper describes a method to estimate mechanical properties of soft, anisotropic materials from measurements of shear waves with specific polarization and propagation directions. This method is applicable to data from magnetic resonance elastography (MRE), which is a method for measuring shear waves in live subjects or in vitro samples. Here, we simulate MRE data using finite element analysis. A nearly incompressible, transversely isotropic (ITI) material model with three parameters (shear modulus, shear anisotropy, and tensile anisotropy) is used, which is appropriate for many fibrous, biological tissues. Both slow and fast shear waves travel concurrently through such a material with speeds that depend on the propagation direction relative to fiber orientation. A three-parameter estimation approach based on directional filtering and isolation of slow and fast shear wave components (directional filter inversion, or DFI) is introduced. Wave speeds of each isolated shear wave component are estimated using local frequency estimation (LFE), and material properties are calculated using weighted least squares. Data from multiple finite element simulations are used to assess the accuracy and reliability of DFI for estimation of anisotropic material parameters.

Keywords: Anisotropy; Inversion algorithms; MR elastography; Shear waves; Transversely isotropic material.

Conflict of interest statement

Conflict of Interest Statement

None of the authors has a conflict of interest that could influence the work described in this manuscript.

Copyright © 2015 Elsevier Ltd. All rights reserved.

Figures

Figure 1

a) Transversely isotropic material with fiber reinforcement. Tensile moduli in directions b) parallel and c) perpendicular to the fibers are given by E1 and E2, respectively. Shear moduli in planes d) parallel and c) perpendicular to the fibers are given by μ1 and μ, respectively. The 13-plane (not shown) has the same shear and tensile properties as the 12-plane. The dashed boxes indicate the undeformed case.

Figure 2

A displacement field with a single propagation direction, n⃗, at an angle θ from the fiber direction, a⃗, can be decomposed into two shear waves, (a) “slow” and (b) “fast” with different polarization directions. This is illustrated for the case in which the fiber direction is aligned with the x-axis. (a) The displacements of the slow shear wave are in the m⃗s polarization direction which lies in the shaded plane. (b) The displacements of the fast shear wave are in the m⃗f polarization direction which lies in the shaded (xz) plane. Note that the wavelength of the fast shear waves is longer than that of the slow shear wave for the same frequency.

Figure 3

The effect of tensile modulus ζ and propagation direction θ on the a) slow cs and b) fast cf shear speeds is shown (μ = ρ = ϕ = 1 and κ → ∞). The tensile modulus increases along a radius from the origin with an angle θ from the θ = 0 axis. An increase in ζ increases cf, but has no effect on cs. The effects of ζ and bulk modulus κ on the c) fast shear speed and d) pressure wave speed cp are shown (μ = ρ = ϕ = 1 and θ = 135°). The fast shear speed approaches a constant value for finite κ.

Figure 4

The process of estimating the shear wave speed for DFI begins with the 3D displacement field. The data displayed in this figure is from the cylindrical simulation shown in Fig. 5. All slices are shown in the xy-plane with the slice location and coordinate system indicated in the upper left hand corner of this figure. The U, V, and W displacement fields are in the x, y, and z directions, respectively. The total displacement field is decomposed into slow and fast shear waves and directionally filtered using the propagation and polarization directions shown, resulting in slow and fast shear wave displacement fields for each direction. Next, wave speeds are estimated from the slow and fast shear wave displacement fields using LFE. Inclusion criteria using an amplitude threshold and certainty threshold result in amplitude and certainty masks, respectively for both the slow and fast shear waves. These amplitude and certainty masks are applied to the speed estimates resulting in the slow and fast shear wave speed estimates shown at the end of the process. Outlier wave speeds (> 1 standard deviation from the mean) are not included in the subsequent parameter fitting step.

Figure 5

Finite element (Comsol™) simulation of Case 1 with displacement in the z-direction shown. The location and direction of excitation is shown by the arrow in a), and the resulting propagation is shown in both the b) xz-plane and c) xy-plane. The lines indicate fiber direction. Note that the wavelength is longer in the direction parallel to planes containing the fibers. For all cases, the boundary conditions (BCs) include a 5 μm excitation at 200 Hz on the inner boundary radius = 1.6 mm; fixed displacement on the outer boundary radius = 23 mm; and free displacement on the top and bottom faces. For all cases, the output data was discretized to simulated images with “field of view” of 48 × 48 × 24 mm3 with a 1 mm3 voxel size. d) Propagation direction vector set used for the local and global inversion approaches.

Figure 6

Analytical propagation speeds (lines) and mean estimated propagation speeds from simulation (symbols) of a) slow and b) fast shear waves. Parameters for Case 1 (dotted line, * symbols), Case 2 (dotted line, □ symbols), Case 3 (solid line, ○ symbols), and Case 4 (dashed line, x symbols) are given in Table 1. Mean wave speed estimates are calculated by averaging voxel estimates for each direction using the process outlined in Fig. 4. Note that Cases 1 and 2 have the same theoretical curve, but Case 2 has a wider range of angles, θ.

Figure 7

Local estimates of parameter values for Case 1 (see Table 1) with added noise (SNR=10) using DFI. a) W-displacement field of slice 12 without noise (SNR= ∞) above and with noise (SNR=10) below. The b) shear modulus (μsim = 1000 Pa), c) shear anisotropy (ϕsim = 1), d) tensile anisotropy (ζsim = 2), and e) R 2 are shown for slices 8 through 17. For the parameters μ, ϕ, and ζ, the range shown is ±50% of the true values (this range contains 98% of all estimated values). The full range is shown for R2.