Shear Wave Propagation and Estimation of Material Parameters in a Nonlinear, Fibrous Material

Elastographic techniques, including both ultrasound elastography and magnetic resonance elastography (MRE) have great potential for noninvasive evaluation of the mechanics of soft tissues. Harmonic MRE is based on MR (magnetic resonance) imaging of shear waves induced by external vibration of the tissue, followed by inversion of the displacement fields to estimate material parameters. MRE has been used to quantify noninvasively the material properties of many biological tissues, such as skeletal muscle [1], liver [2,3], and brain [4]. Most MRE studies use linear elastic or viscoelastic material models, and typically the material is assumed to be isotropic. Recently, MRE has been extended to use anisotropic material models, such as a three-parameter model [5–7] for nearly incompressible, transversely isotropic (TI), fibrous materials. However, these models still rely on the assumptions of linear elasticity, while many biological tissues exhibit nonlinear strain–stress relationships [8].

Nonlinear hyperelastic models have been successfully applied to describe the mechanics of soft biological materials [9,10]. The Holzapfel–Gasser–Ogden (HGO) model in particular is straightforward and widely used to model fibrous soft tissues [11,12]; it contains separate terms to describe the contributions of fiber deformation to the strain energy, and can model an isotropic nonlinear material (with $κ$ = 1/3), or a strongly anisotropic material with single or multiple “families” of fibers. Estimating the parameters of hyperelastic models from experimental data remains an important challenge. Here, we demonstrate that wave speed data, such as those available from MRE studies, can be used for this purpose.

Shear waves in MRE consist of infinitesimal dynamic deformations, which may be superimposed on larger, quasi-static “predeformations.” Shear wave speeds in a nonlinear material are determined by both its mechanical properties and its deformation state. In this study, closed-form expressions for shear wave speeds in the HGO model are obtained in terms of the model parameters and imposed predeformations. Analytical expressions for wave speeds were confirmed by performing finite element (FE) simulations of shear waves in a predeformed cube of HGO material with a single fiber family. Local frequency estimation (LFE) was used to estimate speeds of shear waves with various polarization and propagation directions from simulated displacement fields. Finally, we demonstrate, using simulated data, the feasibility of estimating the material parameters of the HGO model from shear wave speeds.

where $A$ is a fourth-order elasticity tensor which relates the incremental strain tensor, $ε̃$⁠, and incremental stress tensor, $σ̃$⁠. In Cartesian coordinates, this relationship can be expressed in indicial notation, $σ̃pi=Apiqjε̃qj$⁠. For nonlinear models, such as the HGO model, the components of the elasticity tensor can be obtained from the relationship $Apiqj=FpαFqβ((∂2W)/(∂Fiα∂Fjβ))$⁠, where $F$ is the deformation gradient tensor which accounts for the effects of predeformation [13,14], and $W$ is the strain energy function. Thus, in principle, shear wave speeds can be used to estimate material parameters.

Since the acoustic tensor, $Q$⁠, depends on the propagation direction, $n$⁠, in general wave speeds depend on $n$⁠. Also, $Q$⁠, may have up to three distinct eigenvalues (wave speeds) and three corresponding eigenvectors (polarization directions) so that there may be three plane waves that propagate in the same direction. However, material symmetries and constraints reduce the number of possible wave speeds. In an isotropic linear elastic material with shear modulus, $μ$ and bulk modulus, $K$⁠, the acoustic tensor is the same for all propagation directions, and only two wave speeds exist: one longitudinal and one transverse (shear). Longitudinal waves in isotropic materials have $c2=(K+4μ/3)/ρ$ [15], and polarization parallel to the propagation direction (⁠$m=n$); corresponding shear waves have $c2=μ/ρ$ and polarization $m⊥n$⁠. In an isotropic, incompressible linear elastic material, the longitudinal wave speed becomes infinite, and only one parameter, $μ$⁠, remains to estimate.

where K and μ are the bulk modulus and the isotropic shear modulus, respectively, $I¯1$ is the modified first invariant defined by $I¯1=J−2/3I1, (J=detF)$⁠, where $I1$ is the first variant of Cauchy-Green strain tensor $C$⁠. Many soft materials have shear moduli roughly between 10² and 10⁵ Pa, spanning a cellular collagen and fibrin gels [18–21], brain tissue [22,23], and muscle [24]. Anisotropic terms in the strain energy density function can have different forms depending on fiber arrangement

1. Single fiber-family model (TI): Terms in the strain energy density function reflecting the effects of fibers with a distribution of orientations centered on the unit vector $a0$⁠, which is the fiber direction before predeformation:

   $Waniso HGO1=k12k2expk2κI¯1+1−3κI¯4−12−1, for I¯4>1$

   (7)

where $I¯4$ is the modified pseudo-variant defined by $I¯4=J−2/3I4$⁠, $I4=a0⋅C⋅a0$ is the squared stretch in the fiber direction.

1. Multiple fiber-family model (orthotropic): Additional fiber families (with a principal direction of $a0i$⁠, and the same properties k₁, k₂, and κ) can be modeled by adding contributions from $I4i=a0i⋅C⋅a0i$ to the strain energy, as

   $Waniso HGON=∑ik12k2expk2E¯i2−1; E¯i=κI¯1+1−3κI¯4i−1, for I¯4i>1$

   (8)

The effects of $k1$ and $k2$ on stress–strain behavior in simple shear are shown in Fig. 2. For example, for simple shear $γYZ$ in a plane containing fibers at an angle of $π/4$ (⁠$a=(j+k)/2$⁠), $k1$ describes the initial slope of the curve (Fig. 2(a)), $k2$ describes the nonlinearity of the curve (Fig. 2(b)).

Fiber distributions corresponding to different values of $κ$ are shown in Fig. 3, and $κ$ captures the distribution of the fiber orientations, ranging from alignment in a single direction (κ = 0) to no preferred direction (κ = 1/3). When κ = 0, all fibers are assumed to be perfectly aligned, and when κ = 1/3, the material is isotropic. We note that, formally, fibers in the HGO model do not contribute to the stress or to the strain energy when they are in compression $(I4<1)$⁠. We did not model the bilinearity between fiber tension and compression for wave propagation in the HGO model in the undeformed case (assuming that the fibers can resist an infinitesimal compressive strain in wave propagation, or equivalently, an infinitesimal tensile prestrain).

To verify the analytical results, FE simulations of shear wave propagation were performed using finite element software (comsolmultiphysics v5.3, Burlington, MA). A static predeformation step (either (i) simple shear or (ii) tension) and a frequency-domain perturbation step were performed in a cubic domain (50 $×$ 50 $×$ 50 mm³, Figs. 4 and 5). The HGO model was implemented in comsol to model the elastic behavior; an isotropic loss factor of 0.1 was used to provide a small amount of viscoelastic damping. We set the frequency of excitation equal to 200 Hz, in order to provide multiple wavelengths in the model domain. The domain was discretized into 5000 hexahedral elements. To demonstrate convergence, the results were confirmed at higher resolution. In order to compare the undeformed case to the cases with finite predeformation, we assume the fibers can resist infinitesimal compressive loads during wave motion. A periodic boundary condition was applied on the $XZ$-plane and $YZ$-plane, for fast and slow shear waves, respectively. The periodic boundary conditions eliminate boundary effects on the vertical sides of the cube, allowing for a closer comparison of the analytical and numerical results. The assigned default parameters are as follows: predeformation $γXZ=γYZ=0.2$⁠; initial isotropic shear modulus, $μ0=1$ kPa; density $ρ$ = 1000 kg/m³; initial anisotropy ratio, $k1/μ0=2$⁠; nonlinearity parameter, $k2=5$⁠; fiber dispersion parameter, $κ=1/12;$ and ratio of bulk modulus to initial shear modulus, $K/μ0=104$⁠. To obtain either slow or fast shear waves, a harmonic displacement (simple shear case) or harmonic force (lengthening/shortening cases) was imposed to the top surface in the corresponding polarization direction, defined by the $ms$ or $mf$ unit vector, respectively. The LFE method [25] was applied to estimate the wavelength (and thus wave speed) from simulated data. LFE provides an estimate of wave speed at each “voxel” in a discretized version of the 3D wave field. The mean values and standard deviations of wave speeds from all voxels in a central region of interest are used to generate the symbols and error bars in plots [25,26]. The LFE parameters used in this study are $ρ0=1$ for the center frequency and $L0=11$ for the number of filters [26].

Figure 6 shows the simulation results for slow and fast shear waves propagating in the negative Z-direction in the undeformed configuration and with shear predeformation in the $YZ$- or $XZ$-plane. Shear predeformation in either the $YZ−$ or $XZ$-plane (Figs. 6(e) and 6(f)) increases the wavelength of the fast shear wave compared to the undeformed configuration (Fig. 6(d)), corresponding to an increase in the fast wave speed.

Shear wave speeds are compared in analytical predictions and simulations estimations for the three configurations of Figs. 4 and 5, shown in Figs. 7–10 below. In each configuration, one parameter was varied while holding the remaining parameters at the default parameter values given above. The vertical axes of the panels in the top row of each figure display $cf/c0$⁠, the normalized ratio between the fast shear wave speed and the initial wave speed $c0=μ0/ρ$⁠, where $μ0=1000 Pa$ is the initial isotropic shear modulus, $ρ$ is the density of material. Similarly, the vertical axes of the panels on the bottom row of each figure depict the ratio $cs/c0$ between the slow shear wave speed and the initial wave speed. Results are shown for ranges of the HGO parameters isotropic shear modulus $μ$⁠, HGO model parameters $k1 and k2,$ dispersion parameter $κ,$ and the imposed shear, $γ$⁠. In all three figures, solid lines without error bars (orange in online version) depict the analytical predictions, and solid lines with error bars (blue in online version) display corresponding wave speeds estimated from FE simulations.

Figure 7 shows the relationships between shear wave speeds and the parameters in the HGO model in the undeformed configuration. The horizontal axis represents three normalized parameters (⁠$μ/μ0, k1/μ0, κ)$ in HGO model. There is no effect of changing $k2$ or $γ$ because no predeformation is applied. The fast wave speed increases with increasing $k1$ and $μ$ and decreases with increasing $κ$⁠. In contrast, the slow wave speed is affected only by $μ$⁠.

Figure 8 shows the dependence of wave speed on HGO parameters when simple shear predeformation is imposed in the $YZ$-plane, i.e., in the direction that induces the fiber stretch. The horizontal axis of each panel displays values of one of the four parameters (⁠$μ/μ0, k1/μ0, κ,k2$) in the HGO model or the magnitude of shear $γYZ$⁠. The slow and fast wave speeds all increase with the increasing $k1,k2$⁠, $μ$⁠, and $γYZ$⁠, and decrease with the fiber dispersion, $κ$⁠. The fast wave speed is larger than the slow wave speed due to the stiffening effect of the fibers.

The effects of the HGO parameters on shear wave speeds are illustrated in Fig. 9 for the configuration in which predeformation is applied perpendicular to the original fiber axis. The vertical axis of each panel on the top row shows the (normalized) fast wave speed, and on the bottom row depicts slow wave speed. The horizontal axis of each panel shows the value of the HGO model parameter or the magnitude of shear. Wave speed is influenced by all five parameters $μ/μ0, k1/μ0, κ,k2,$ and $γYZ$⁠. Similar to predeformation in the $YZ$-plane, the fast wave speed increases with the increasing $μ/μ0, k1/μ0,k2,γXZ$ and decreases with the increasing $κ$⁠. Slow wave speeds follow the same trend as the fast wave speeds, but to a lesser extent.

The effect of stretch ratio on wave speed is shown in Fig. 10 for the case of imposed extension. Both fast and slow wave speeds increase with stretch ratio and the simulation results agree well with the analytical predictions.

In Sec. 2.3, we demonstrated that shear wave speeds can be calculated analytically from parameter values of the HGO model, for specific propagation direction and polarization directions. We also confirmed that the analytical solutions agree with simulated wave speeds in a finite cube-shaped domain. Conversely, the parameters of the HGO model can, in principle, be estimated from measured shear wave speeds, for given propagation and polarization directions, in predeformed specimens. In Sec. 4.1, we demonstrate the feasibility of this approach to parameter estimation.

In the proposed experiment, the density $ρ$ of the material is a known value, and the simple shear ratio $γYZ$(or stretch ratio $λ$⁠) can be controlled and measured. For a single value of the predeformation, the four independent linear equations can be solved simultaneously to determine the four independent parameters. If more data are available, the over-determined system can be solved in the least-squares sense. The matlab optimization tool lsqnonlin for solving nonlinear least square problems was used to find parameters that minimized the difference between predicted and measured values of wave speed.

Here, τ is a random value in the standard normal distribution (mean = 0, standard deviation = 1), and $ψ$ is defined as a noise factor to control the range of noise variance. In this paper, the noise is defined on three levels. Values of $ψ$ = 0.01, 0.02, and 0.03 indicate wave speed variance ranges of $±3.3%,±6.6%, and ±10%$ from the expected values, respectively.

For wave speed data without noise, the material parameters can be determined by four equations corresponding to two configurations, the undeformed configuration and one value of predeformation. However, if wave speed data are noisy, more data are needed. In a Monte Carlo approach, ten additional simulated experiments with different predeformations were added to the original two simulated experiments, and these simulated experiments were repeated 1000 times with different random values. For various noise levels, the mean values (⁠$±$ standard deviation) of all four parameters were calculated (Tables 1 and 2). To improve the accuracy, outliers (greater than three standard deviations from the mean) were excluded from wave speed data.

As expected, the standard deviation of each parameter estimate increases with noise level. The mean value of some parameters also deviates from the expected value as noise increases. The parameters $k2$ and $μ0$ are relatively insensitive to the noise level; estimates of $k1$ and $κ$ deviate more when noise increases.

In materials that can be modeled as nonlinear, anisotropic, and nearly incompressible, slow and fast wave speeds can be measured from MRE and used to estimate parameters of the material model. In the examples above, theoretical predictions of shear wave speed values in different configurations (undeformed configuration, simple shear in the $XZ$-plane or $YZ$-plane) agreed well with simulation results.

Fast and slow shear wave speeds provide complementary information. The fast shear wave speed is affected by the stiffness of the fibers, while the slow shear wave speed is not. In transversely isotropic materials, displacements in the direction of slow shear wave polarization do not induce fiber stretch. In addition, in the example of this paper, simple shear in the $YZ$-plane (which contains the principal fiber axis, (⁠$a$⁠) directly stretches the fibers and significantly affects fast shear wave speeds. In contrast, simple shear in the $XZ$-plane involves displacements perpendicular to the fibers, and does not stretch the fibers appreciably. The measured wave speed in this condition deviates little from the wave speed in the no predeformation condition. Extension in the Z-direction also stretches the fibers, which increases shear wave speeds.

Using the closed-form expressions for wave speed, and data from either simulations or experiments, we can estimate the parameters of a nonlinear anisotropic material model. Unlike linear elastic materials, predeformation plays an important role in determining wave speeds. Without predeformation, the slow shear wave speed varies only with the isotropic shear modulus $μ$ of material and the fast shear wave depends on $μ, k1,$ and $κ$⁠. If predeformations that stretch the fibers are imposed, the fast shear wave speed depends on all the HGO parameters, $μ,k1,k2,κ,$ as well as the magnitude of the predeformation.

The accuracy of material parameter estimates is affected by the levels of noise in the measured wave speed data. The isotropic shear modulus $μ$ is the least sensitive to noise because it is directly derived from the slow wave speed with no predeformation. For the other three parameters, the nonlinearity $k2$ is less sensitive to noise than $k1$ and $κ$⁠, because $k2$ has fewer interactions with other factors in the experiment.

Among the limitations of this paper, in computing wave speeds we do not impose the bilinearity that excludes fibers from resisting infinitesimal levels of dynamic compressive strain. In practice, this assumption could be avoided by imposing a minimal predeformation greater than the wave amplitude. We have considered the original version of the HGO model, which is still widely used. A new version of HGO model has recently been proposed [27], which might also be analyzed by this approach. Parameter estimates improve, in terms of both increased accuracy and reduced variance, with more MRE experiments. Balancing the desired precision of the result and the cost of experiments must be considered carefully, as in all experimental studies.

Only one fiber family is considered in this paper, but it is plausible to generalize this approach to estimate wave speeds and material parameters in a material with multiple fiber families. Some special cases can be considered qualitatively. For simplicity, consider a second fiber family with $ϕ=−45 deg$⁠. For the situation in which simple shear is imposed in the YZ-plane, one fiber family will be stretched and the other will be compressed. In the original HGO model, fibers in compression $(I4<1)$ do not contribute to the stress or to the strain energy. Therefore, wave speeds in a material with two fiber families would be the same as with one fiber family. For the same reason, if the material is compressed in the Z-direction (⁠$λ < 1$⁠), wave speeds are equal to those in an isotropic material because all fibers are under compression. For the idealized example of imposed extension, adding a second fiber family would simply double the effects of a single fiber family. For other configurations the addition of a second fiber family creates orthotropic material symmetry, instead of the transverse isotropy considered in this paper, and would (in general) require further analysis.

Experimental studies that exploit this approach to characterize fibrous soft tissues are planned for future work; these studies would involve superimposing small amplitude shear waves on large finite deformations. Instead of the idealized simple shear deformations of this paper, dense measurements of actual predeformations would need to be combined with regional measurements of shear wave speed. While challenging, this approach promises the possibility of comprehensive, noninvasive tissue characterization in vivo. Tissue may already be in a predeformed state (like white matter in the brain) [28,29], or quasi-static loading might be imposed by respiration (liver [30]), ocular pressure (eye [31]), or external force (intervertebral disk, muscle or breast [32,33]).

MR elastography can be used, in principle, to estimate parameters of the HGO material model in soft fibrous materials from the speeds of slow and fast shear waves. To demonstrate the ability to obtain accurate results, closed-form expressions for the wave speeds, as functions of predeformation and material parameters, were derived and confirmed by numerical simulations. These results illustrate the feasibility of a new approach to parameter estimation for nonlinear material models of fibrous soft matter.

Financial support for this study was provided by NSF Grant No. CMMI-1727412, NIH Grant Nos. R01/R56 NS055951 and R01 EB027577.

• NSF (Funder ID: 10.13039/100000001).

• NIH (Funder ID: 10.13039/100000002).

 
• $a$ =

  initial fiber direction vector

 
• $A$ =

  elasticity tensor

 
• $c$ =

  shear wave speed

 
• $C$ =

  Cauchy-Green strain tensor

 
• $cf$ =

  fast wave speed under predeformation

 
• $cf0$ =

  fast wave speed under no predeformation

 
• $cs$ =

  slow wave speed under predeformation

 
• $cs0$ =

  slow wave speed under no predeformation

 
• $cf(noise)$ =

  slow wave speed with noise

 
• $cs(noise)$ =

  slow wave speed with noise

 
• $c0$ =

  initial shear wave speed

 
• $E1,E2$ =

  two tensile moduli in TI material

 
• F =

  deformation gradient tensor

 
• $I1$ =

  first variant of Cauchy-Green strain tensor

 
• $I4$ =

  squared stretch in the fiber direction

 
• $I¯1$ =

  modified first variant

 
• $I¯4$ =

  modified pseudo-variant

 
• $K$ =

  bulk modulus of the material

 
• $k1$ =

  initial slope of strain-stress curve in the HGO model

 
• $k2$ =

  nonlinearity of strain-stress curve in the HGO model

 
• $L$ =

  abbreviation in closed-form expression of HGO model

 
• $L0$ =

  number of filters

 
• $m$ =

  polarization direction vector

 
• $M$ =

  abbreviation in closed-form expression of HGO model

 
• $mf$ =

  fast polarization direction

 
• $ms$ =

  slow polarization direction

 
• $n$ =

  wave propagation direction vector

 
• $N$ =

  abbreviation in closed-form expression of HGO model

 
• $Q$ =

  acoustic tensor

 
• $W$ =

  strain energy function

 
• $γ$ =

  shear predeformation

 
• $γXZ$ =

  simple shear ratio in $XZ$-plane

 
• $γYZ$ =

  simple shear ratio in $YZ$-plane

 
• $ζTI$ =

  parameter in TI material model

 
• $θ$ =

  angle between fiber and propagation directions

 
• $κ$ =

  dispersion parameter of fibers in the HGO model

 
• $λ$ =

  stretch ratio in Z-direction

 
• $μ$ =

  isotropic shear modulus in the HGO model

 
• $μ0$ =

  initial value of isotropic shear modulus in the HGO model

 
• $μ1,μ2$ =

  two shear moduli in TI material

 
• $ρ$ =

  density of the material

 
• $ρ0$ =

  LFE parameter

 
• $τ$ =

  random value from normal distribution

 
• $ϕ$ =

  deviation angle between fiber direction and wave propagation direction

 
• $ϕTI$ =

  parameter in TI material model

 
• $ψ$ =

  noise factor in sensitivity analysis