## Abstract

This paper describes the propagation of shear waves in a Holzapfel–Gasser–Ogden (HGO) material and investigates the potential of magnetic resonance elastography (MRE) for estimating parameters of the HGO material model from experimental data. In most MRE studies the behavior of the material is assumed to be governed by linear, isotropic elasticity or viscoelasticity. In contrast, biological tissue is often nonlinear and anisotropic with a fibrous structure. In such materials, application of a quasi-static deformation (predeformation) plays an important role in shear wave propagation. Closed form expressions for shear wave speeds in an HGO material with a single family of fibers were found in a reference (undeformed) configuration and after imposed predeformations. These analytical expressions show that shear wave speeds are affected by the parameters ($μ0, k1, k2, κ$) of the HGO model and by the direction and amplitude of the predeformations. Simulations of corresponding finite element (FE) models confirm the predicted influence of HGO model parameters on speeds of shear waves with specific polarization and propagation directions. Importantly, the dependence of wave speeds on the parameters of the HGO model and imposed deformations could ultimately allow the noninvasive estimation of material parameters in vivo from experimental shear wave image data.

## 1 Introduction

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 [57] 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.

## 2 Theoretical Methods

### 2.1 Wave Speeds in Transversely Isotropic Elastic Materials.

Shear wave speeds in an elastic material are calculated from the eigenvalues of the acoustic tensor [13,14], as in the following equation:
$ρc2m=Qn·m$
(1)
where $ρc2$ is the eigenvalue of the acoustic tensor $Q$, $ρ$ is the density of material, $c$ is the shear wave speed, $n$ is the propagation direction of wave, and $m$ is the polarization direction vector of the shear wave. The acoustic tensor $Q$ for a specific propagation direction, $n$, is obtained from Eq. (2) [13,14] below:
$Q=n·A·n$
(2)

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.

Fig. 1
Fig. 1
Close modal
A linear elastic, TI material requires five parameters to describe its constitutive behavior. These can be, for example, two tensile moduli ($E1$, $E2$), two shear moduli ($μ1$, $μ2$) and a bulk modulus ($K$). If the material is incompressible, three parameters are sufficient, for example, baseline shear modulus ($μ2)$ and two ratios: shear anisotropy $ϕTI=μ1/μ2−1$ and tensile anisotropy $ζTI=E1/E2−1$. In a linear elastic, nearly incompressible, TI material, in which anisotropy is due to a single family of aligned fibers, the shear wave speeds depends in a relatively simple fashion on the material properties and the wave propagation direction relative to the fibers [15]. Shear waves can be separated into slow and fast shear waves in different polarization directions (Fig. 1). The slow and fast shear wave polarization directions under no predeformation are defined by the following relationship:
$ms=n×a|n×a| mf=n×ms$
(3)
where $n$ is the propagation direction of the shear wave,$a$ is the fiber orientation after deformation, and $ms$ and $mf$ are the slow and fast polarization directions, respectively. The corresponding wave speeds in this linear elastic material, which depend on the angle, $θ$, between fiber and propagation directions are
$cs2=μ2ρ1+ ϕTI cos2θ; cf2=μ2ρ1+ ϕTI cos22θ+ ζTIsin22θ$
(4)

### 2.2 The Holzapfel–Gasser–Ogden Model.

The Holzapfel–Gasser–Ogden model, which is an influential recent model for fibrous soft tissues, was proposed by Holzapfel [17]. The strain energy density function is a sum of isotropic and anisotropic terms
$W=Wiso+Waniso$
(5)
The isotropic part of its strain energy density function contains both volumetric and isochoric terms
$Wiso=Wvol+WisochoricWvol=K2J−12, Wisochoric=μ2I¯1−3$
(6)

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 102 and 105 Pa, spanning a cellular collagen and fibrin gels [1821], 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 k1, k2, 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)).

Fig. 2
Fig. 2
Close modal

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).

Fig. 3
Fig. 3
Close modal

### 2.3 Closed-Form Expressions for the Relationships Between Model Parameters and Wave Speeds

#### 2.3.1 Closed-Form Expressions for Speeds of Waves Superimposed on Simple Shear.

Closed-form expressions that relate wave speeds to model parameters are highly desirable. Such expressions for the HGO model were determined from analytical solution of the eigenvalue problem (Eq. (1)) for shear waves propagating in the negative Z-direction ($n=−k$, Fig. 4) in the undeformed configuration (Fig. 4(a)) and with imposed predeformations in simple shear corresponding to the configurations of Figs. 4(b) and 4(c). Symbolic solutions were obtained using Matlab Symbolic Toolbox (Mathworks, Natick, MA)
$csXZ=μρ+2k1ργXZ2M2expM2k2γXZ4$
(9)
(10)
$csYZ=μρ+2k1ργYZMγYZM+NexpγYZM+N2k2γYZ2$
(11)
$cfYZ =μρ+k1ρN2+6γYZMN+6γYZ2N2+2γYZ22γYZM+N2γYZM+N2expγYZM+N2k2γYZ2$
(12)
where $(csXZ$, $cfXZ)$ and $(csYZ$, $cfYZ )$ are the slow and fast shear wave speeds for predeformations $γXZ$ and $γYZ$, respectively. The terms $M$ and $N$ are combinations of the dispersion parameter $κ$ and the angle $ϕ$ between the fiber and the propagation direction:
(13)
Fig. 4
Fig. 4
Close modal
With no imposed predeformation ($γXZ=γYZ=0$), the speeds reduce to
$cs=c0=μρ; cf=μρ+k1N2ρ$
(14)

#### 2.3.2 Closed-Form Expressions for Speeds of Waves Superimposed on Stretching.

Expressions for the speeds of slow and fast shear waves superimposed on isochoric, static lengthening deformations (Fig. 5) were also obtained. In this situation, the maximum stretch ratio, $λ1=λ$, is used to describe the imposed deformation; other stretches are $λ2=λ3=1/λ$. Wave speeds were obtained with the help of matlab Symbolic Toolbox (Mathworks, Natick, MA)
$CfZ=λ2μρ+k1ρλλ2N2+λ2ML+2N2L2k2expk2L2λ2$
(15)
$CsZ=λ2μρ+λMLk1ρexpk2L2λ2$
(16)
where $(csZ$, $cfZ)$ are the slow and fast shear wave speeds for stretch $λ$ in Z-direction. Definitions of $M$ and $N$ are from Eq. (13). $L$ is defined in terms of $M$ and the stretch ratio $λ$
$L=λ−1M2λ2+λ+1−1$
(17)
Fig. 5
Fig. 5
Close modal
With no imposed predeformation ($λ=1$), the wave speed expressions can be simplified, as above, to
$cs=c0=μρ; cf=μρ+k1N2ρ$
(18)
which are identical to Eq. (14).

### 2.4 Computational Modeling and Simulations.

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 mm3, 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/m3; 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].

## 3 Results: Shear Wave Speeds in Undeformed and Deformed Configurations

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.

Fig. 6
Fig. 6
Close modal

Shear wave speeds are compared in analytical predictions and simulations estimations for the three configurations of Figs. 4 and 5, shown in Figs. 710 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.

Fig. 7
Fig. 7
Close modal
Fig. 8
Fig. 8
Close modal
Fig. 9
Fig. 9
Close modal
Fig. 10
Fig. 10
Close modal

### 3.1 One Family of Fibers in the Undeformed Configuration.

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 $μ$.

### 3.2 One Fiber Family With Predeformation by Simple Shear in the $YZ$-Plane.

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.

### 3.3 One Fiber Family With Predeformation by Simple Shear in the $XZ$-Plane.

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.

### 3.4 One Fiber Family With Predeformation Consisting of Imposed Extension.

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.

## 4 Estimation of Parameters in the Holzapfel–Gasser–Ogden Model

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.

### 4.1 Estimation Method.

The example system is shown in Fig. 11. The angle $ϕ$ of the fiber axis is chosen to be $π/4$ radians from the base of the specimen in the undeformed configuration, as in Figs. 4 and 5. The propagation direction $n$ is along negative Z-axis, and fiber direction $a$ is in the $YZ$-plane. Experiments are separated into two steps. In the first step, the fast wave speed and slow wave speed are measured without predeformation, by applying horizontal, harmonic displacements to the top surface in the fast or slow polarization directions. The slow wave speed is a function of only the isotropic shear modulus, $μ$, and density ρ, but the fast wave speed is a function of $μ, k1,$ and $κ$ (Eq. (19)). In the second step, the fast wave speed and slow wave speed are measured after applying a predeformation of simple shear in the $YZ−$ plane). In this configuration, both fast and slow wave speeds are functions of all five parameters ($μ, k1,k2,κ$, and $γYZ$) (Eq. (20))
$cs0=μ/ρ; cf0=f(k1,μ,κ)$
(19)
$cs=fk1,k2,μ,κ,γYZ; cf=fk1,k2,μ,κ,γYZ$
(20)
Fig. 11
Fig. 11
Close modal
For the analogous situation using imposed extension (stretch ratio $λ$), Eq. (20) can be written as
$cs=fk1,k2,μ,κ,λ;cf=fk1,k2,μ,κ,λ$
(21)

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.

### 4.2 Sensitivity to Noise.

Because experimental data inevitably contain noise or measurement errors, it is necessary to quantify the robustness of parameter estimates. For each wave speed estimate, random noise was applied from a normal distribution, as shown in the following equation:
$cf(noise)=cf(1+ψτ)$
(22)
$cs(noise)=cs(1+ψτ)$
(23)

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.

Table 1

Comparison of HGO model parameter estimates for different noise levels (imposed shear)

Noise level$μ$(Pa)$k1$(Pa)$κ$$k2$
Expected100020000.0835
($ψ$ = 0.01)$1000 ± 11$$2034 ± 367$$0.083 ± 0.023$$5 ± 0.0$
($ψ$ = 0.02)$1000 ± 23$$2050 ± 739$$0.077 ± 0.044$$4.8 ± 0.1$
($ψ$ = 0.03)$1001 ± 34$$2161 ± 1107$$0.079 ± 0.057$$4.7 ± 0.2$
Noise level$μ$(Pa)$k1$(Pa)$κ$$k2$
Expected100020000.0835
($ψ$ = 0.01)$1000 ± 11$$2034 ± 367$$0.083 ± 0.023$$5 ± 0.0$
($ψ$ = 0.02)$1000 ± 23$$2050 ± 739$$0.077 ± 0.044$$4.8 ± 0.1$
($ψ$ = 0.03)$1001 ± 34$$2161 ± 1107$$0.079 ± 0.057$$4.7 ± 0.2$
Table 2

Comparison of HGO model parameter estimates for different noise levels (imposed extension)

Noise level$μ$(Pa)$k1$(Pa)$κ$$k2$
Expected100020000.0835
($ψ$ = 0.01)$1000 ± 11$$2007 ± 390$$0.081 ± 0.023$$4.9 ± 0.1$
($ψ$ = 0.02)$1000 ± 20$$2044 ± 752$$0.077 ± 0.043$$4.7 ± 0.3$
($ψ$ = 0.03)$1002 ± 32$$2068 ± 1125$$0.075 ± 0.050$$4.5 ± 0.5$
Noise level$μ$(Pa)$k1$(Pa)$κ$$k2$
Expected100020000.0835
($ψ$ = 0.01)$1000 ± 11$$2007 ± 390$$0.081 ± 0.023$$4.9 ± 0.1$
($ψ$ = 0.02)$1000 ± 20$$2044 ± 752$$0.077 ± 0.043$$4.7 ± 0.3$
($ψ$ = 0.03)$1002 ± 32$$2068 ± 1125$$0.075 ± 0.050$$4.5 ± 0.5$

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.

## 5 Discussion

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]).

## 6 Conclusion

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.

## Acknowledgment

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

## Funding Data

• NSF (Funder ID: 10.13039/100000001).

• NIH (Funder ID: 10.13039/100000002).

## Nomenclature

• $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 =

•
• $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

## References

1.
Bensamoun
,
S. F.
,
Charleux
,
F.
,
Debernard
,
L.
,
Themar-Noel
,
C.
, and
Voit
,
T.
,
2015
, “
Elastic Properties of Skeletal Muscle and Subcutaneous Tissues in Duchenne Muscular Dystrophy by Magnetic Resonance Elastography (MRE): A Feasibility Study
,”
IRBM
,
36
(
1
), pp.
4
9
.10.1016/j.irbm.2014.11.002
2.
Chamarthi
,
S. K.
,
Raterman
,
B.
,
Mazumder
,
R.
,
Michaels
,
A.
,
Oza
,
V. M.
,
Hanje
,
J.
,
Bolster
,
B.
,
Jin
,
N.
,
White
,
R. D.
, and
Kolipaka
,
A.
,
2014
, “
Rapid Acquisition Technique for MR Elastography of the Liver
,”
Magn. Reson. Imaging
,
32
(
6
), pp.
679
683
.10.1016/j.mri.2014.02.013
3.
Yang
,
C.
,
Yin
,
M.
,
Glaser
,
K. J.
,
Zhu
,
X.
,
Xu
,
K.
,
Ehman
,
R. L.
, and
Chen
,
J.
,
2017
, “
Static and Dynamic Liver Stiffness: An Ex Vivo Porcine Liver Study Using MR Elastography
,”
Magn. Reson. Imaging
,
44
, pp.
92
95
.10.1016/j.mri.2017.08.009
4.
Kolipaka
,
A.
,
Wassenaar
,
P. A.
,
Cha
,
S.
,
Marashdeh
,
W. M.
,
Mo
,
X.
,
Kalra
,
P.
,
Gans
,
B.
,
Raterman
,
B.
, and
Bourekas
,
E.
,
2018
, “
Magnetic Resonance Elastography to Estimate Brain Stiffness: Measurement Reproducibility and Its Estimate in Pseudotumor Cerebri Patients
,”
Clin. Imaging
,
51
, pp.
114
122
.10.1016/j.clinimag.2018.02.005
5.
Guo
,
J.
,
Hirsch
,
S.
,
Scheel
,
M.
,
Braun
,
J.
, and
Sack
,
I.
,
2016
, “
Three-Parameter Shear Wave Inversion in MR Elastography of Incompressible Transverse Isotropic Media: Application to In Vivo Lower Leg Muscles
,”
Magn. Reson. Med.
,
75
(
4
), pp.
1537
1545
.10.1002/mrm.25740
6.
Labus
,
K. M.
, and
Puttlitz
,
C. M.
,
2016
, “
An Anisotropic Hyperelastic Constitutive Model of Brain White Matter in Biaxial Tension and Structural-Mechanical Relationships
,”
J. Mech. Behav. Biomed. Mater.
,
62
, pp.
195
208
.10.1016/j.jmbbm.2016.05.003
7.
Schmidt
,
J. L.
,
Tweten
,
D. J.
,
,
A. A.
,
Reiter
,
A. J.
,
Okamoto
,
R. J.
,
Garbow
,
J. R.
, and
Bayly
,
P. V.
,
2018
, “
Measurement of Anisotropic Mechanical Properties in Porcine Brain White Matter Ex Vivo Using Magnetic Resonance Elastography
,”
J. Mech. Behav. Biomed. Mater.
,
79
, pp.
30
37
.10.1016/j.jmbbm.2017.11.045
8.
Darvish
,
K. K.
, and
Crandall
,
J. R.
,
2001
, “
Nonlinear Viscoelastic Behavior of Brain Tissue in Oscillatory Shear Deformation
,”
Methods
,
23
, pp.
633
645
.10.1016/S1350-4533(01)00101-1
9.
Panda
,
S. K.
, and
Buist
,
M. L.
,
2018
, “
A Finite Nonlinear Hyper-Viscoelastic Model for Soft Biological Tissues
,”
J. Biomech.
,
69
, pp.
121
128
.10.1016/j.jbiomech.2018.01.025
10.
Ghazy
,
M.
,
Elgindi
,
M. B.
, and
Wei
,
D.
,
2018
, “
Analytical and Numerical Investigations of the Collapse of Blood Vessels With Nonlinear Wall Material Embedded in Nonlinear Soft Tissues
,”
Alexandria Eng. J.
,
57
(
4
), pp.
931
965
.10.1016/j.aej.2018.03.002
11.
Volokh
,
K. Y.
,
2011
, “
Modeling Failure of Soft Anisotropic Materials With Application to Arteries
,”
J. Mech. Behav. Biomed. Mater.
,
4
(
8
), pp.
1582
1594
.10.1016/j.jmbbm.2011.01.002
12.
Shearer
,
T.
,
2015
, “
A New Strain Energy Function for the Hyperelastic Modelling of Ligaments and Tendons Based on Fascicle Microstructure
,”
J. Biomech.
,
48
(
2
), pp.
290
297
.10.1016/j.jbiomech.2014.11.031
13.
Ogden
,
R. W.
,
1997
, “
Wave and Vibrations
,”
Non-Linear Elastic Deformations
,
E. Horwood
,
Chichester, UK
, p.
473
.
14.
Volokh
,
K.
,
2016
, “
Plane Waves in Incompressible Material
,”
Mechanics of Soft Materials
,
Springer
,
Singapore
, p.
96
.
15.
Birch
,
F.
,
1961
, “
The Velocity of Compressional Waves in Rocks to 10 Kilobars—Part 2
,”
J. Geophys. Res.
,
66
(
7
), pp.
2199
2224
.10.1029/JZ066i007p02199
16.
Tweten
,
D. J.
,
Okamoto
,
R. J.
,
Schmidt
,
J. L.
,
Garbow
,
J. R.
, and
Bayly
,
P. V.
,
2015
, “
Estimation of Material Parameters From Slow and Fast Shear Waves in an Incompressible, Transversely Isotropic Material
,”
J. Biomech.
,
48
(
15
), pp.
4002
4009
.10.1016/j.jbiomech.2015.09.009
17.
Holzapfel
,
G.
,
2000
,
Nonlinear Solid Mechanics: A Continuum Approach for Engineering
,
Wiley
,
Chichester, UK
.
18.
Namani
,
R.
,
Feng
,
Y.
,
Okamoto
,
R. J.
,
Sakiyama Elbert
,
S. E.
,
Genin
,
G. M.
, and
Bayly
,
P. V.
,
2012
, “
Elastic Characterization of Transversely Isotropic Soft Materials by Dynamic Shear and Asymmetric Indentation
,”
ASME J. Biomech. Eng.
,
134
(
6
), p.
061004
.10.1115/1.4006848
19.
Sundararaghavan
,
H. G.
,
Monteiro
,
G. A.
,
Lapin
,
N. A.
,
Chabal
,
Y. J.
,
Miksan
,
J. R.
, and
Shreiber
,
D. I.
,
2008
, “
Genipin-Induced Changes in Collagen Gels: Correlation of Mechanical Properties to Fluorescence
,”
J. Biomed. Mater. Res.—Part A
,
87
(
2
), pp.
308
320
.10.1002/jbm.a.31715
20.
Zhang
,
Y.
,
Xu
,
B.
, and
Chow
,
M. J.
,
2011
, “
Experimental and Modeling Study of Collagen Scaffolds With the Effects of Crosslinking and Fiber Alignment
,”
Int. J. Biomater.
,
2011
, p.
172389
.10.1155/2011/172389
21.
Lai
,
V. K.
,
Lake
,
S. P.
,
Frey
,
C. R.
,
Tranquillo
,
R. T.
, and
Barocas
,
V. H.
,
2012
, “
Mechanical Behavior of Collagen-Fibrin Co-Gels Reflects Transition From Series to Parallel Interactions With Increasing Collagen Content
,”
ASME J. Biomech. Eng.
,
134
(
1
), p.
011004
.10.1115/1.4005544
22.
Goriely
,
A.
,
Geers
,
M. G. D.
,
Holzapfel
,
G. A.
,
Jayamohan
,
J.
,
Jérusalem
,
A.
,
Sivaloganathan
,
S.
,
Squier
,
W.
,
van Dommelen
,
J. A. W.
,
Waters
,
S.
, and
Kuhl
,
E.
,
2015
, “
Mechanics of the Brain: Perspectives, Challenges, and Opportunities
,”
Biomech. Model. Mechanobiol.
,
14
(
5
), pp.
931
965
.10.1007/s10237-015-0662-4
23.
Budday
,
S.
,
Sommer
,
G.
,
Holzapfel
,
G. A.
,
Steinmann
,
P.
, and
Kuhl
,
E.
,
2017
, “
Viscoelastic Parameter Identification of Human Brain Tissue
,”
J. Mech. Behav. Biomed. Mater.
,
74
, pp.
463
476
.10.1016/j.jmbbm.2017.07.014
24.
Okamoto
,
R. J.
,
Moulton
,
M. J.
,
Peterson
,
S. J.
,
Li
,
D.
,
Pasque
,
M. K.
, and
Guccione
,
J. M.
,
2000
, “
Epicardial Suction: A New Approach to Mechanical Testing of the Passive Ventricular Wall
,”
ASME J. Biomech. Eng.
,
122
(
5
), pp.
479
487
.10.1115/1.1289625
25.
Knutsson
,
H.
,
Westin
,
C. F.
, and
Granlund
,
G.
,
1994
, “
Local Multiscale Frequency and Bandwidth Estimation
,”
International Conference on Image Processing
(
ICIP
), Austin, TX, Nov. 13–16, pp.
36
40
.10.1109/ICIP.1994.413270
26.
Okamoto
,
R. J.
,
Johnson
,
C. L.
,
Feng
,
Y.
,
,
J. G.
, and
Bayly
,
P. V.
,
2014
, “
MRE Detection of Heterogeneity Using Quantitative Measures of Residual Error and Uncertainty
,”
Proc. SPIE
,
9038
, p.
90381E
. 10.1117/12.2044633
27.
Li
,
K.
,
Ogden
,
R. W.
, and
Holzapfel
,
G. A.
,
2018
, “
Modeling Fibrous Biological Tissues With a General Invariant That Excludes Compressed Fibers
,”
J. Mech. Phys. Solids
,
110
, p.
011004
.10.1016/j.jmps.2017.09.005
28.
Xu
,
G.
,
Bayly
,
P. V.
, and
Taber
,
L. A.
,
2009
, “
Residual Stress in the Adult Mouse Brain
,”
Biomech. Model. Mechanobiol.
,
8
(
4
), pp.
253
262
.10.1007/s10237-008-0131-4
29.
Xu
,
G.
,
Knutsen
,
A. K.
,
Dikranian
,
K.
,
Kroenke
,
C. D.
,
Bayly
,
P. V.
, and
Taber
,
L. A.
,
2010
, “
Axons Pull on the Brain, but Tension Does Not Drive Cortical Folding
,”
ASME J. Biomech. Eng.
,
132
(
7
), p.
071013
.10.1115/1.4001683
30.
Yun
,
M. H.
,
Seo
,
Y. S.
,
Kang
,
H. S.
,
Lee
,
K. G.
,
Kim
,
J. H.
,
An
,
H.
,
Yim
,
H. J.
,
Keum
,
B.
,
Jeen
,
Y. T.
,
Lee
,
H. S.
,
Chun
,
H. J.
,
Um
,
S. H.
,
Kim
,
C. D.
, and
Ryu
,
H. S.
,
2011
, “
The Effect of the Respiratory Cycle on Liver Stiffness Values as Measured by Transient Elastography
,”
J. Viral Hepat.
,
18
(
9
), pp.
631
636
.10.1111/j.1365-2893.2010.01376.x
31.
Nguyen
,
T.-M.
,
Arnal
,
B.
,
Song
,
S.
,
Huang
,
Z.
,
Wang
,
R. K.
, and
O'Donnell
,
M.
,
2015
, “
Shear Wave Elastography Using Amplitude-Modulated Acoustic Radiation Force and Phase-Sensitive Optical Coherence Tomography
,”
J. Biomed. Opt.
,
20
(
1
), p.
016001
.10.1117/1.JBO.20.1.016001
32.
Chan
,
D. D.
,
Gossett
,
P. C.
,
Butz
,
K. D.
,
Nauman
,
E. A.
, and
Neu
,
C. P.
,
2014
, “
Comparison of Intervertebral Disc Displacements Measured Under Applied Loading With MRI at 3.0 T and 9.4 T
,”
J. Biomech.
,
47
(
11
), pp.
2801
2806
.10.1016/j.jbiomech.2014.05.026
33.
Capilnasiu
,
A.
,
,
M.
,
Fovargue
,
D.
,
Patel
,
D.
,
Holub
,
O.
,
Bilston
,
L.
,
Screen
,
H.
,
Sinkus
,
R.
, and
Nordsletten
,
D.
,
2019
, “
Magnetic Resonance Elastography in Nonlinear Viscoelastic Materials Under Load
,”
Biomech. Model. Mechanobiol.
,
18
(
1
), pp.
111
135
.10.1007/s10237-018-1072-1