Abstract

This paper describes a new method for estimating anisotropic mechanical properties of fibrous soft tissue by imaging shear waves induced by focused ultrasound (FUS) and analyzing their direction-dependent speeds. Fibrous materials with a single, dominant fiber direction may exhibit anisotropy in both shear and tensile moduli, reflecting differences in the response of the material when loads are applied in different directions. The speeds of shear waves in such materials depend on the propagation and polarization directions of the waves relative to the dominant fiber direction. In this study, shear waves were induced in muscle tissue (chicken breast) ex vivo by harmonically oscillating the amplitude of an ultrasound beam focused in a cylindrical tissue sample. The orientation of the fiber direction relative to the excitation direction was varied by rotating the sample. Magnetic resonance elastography (MRE) was used to visualize and measure the full 3D displacement field due to the ultrasound-induced shear waves. The phase gradient (PG) of radially propagating “slow” and “fast” shear waves provided local estimates of their respective wave speeds and directions. The equations for the speeds of these waves in an incompressible, transversely isotropic (TI), linear elastic material were fitted to measurements to estimate the shear and tensile moduli of the material. The combination of focused ultrasound and MR imaging allows noninvasive, but comprehensive, characterization of anisotropic soft tissue.

1 Introduction

1.1 Motivation and Background: Anisotropy in Soft Tissues and Biomaterials.

Muscles, tendons, and white matter of the brain are important examples of fibrous soft tissue. These biological tissues are structurally and mechanically anisotropic, meaning that their structure and response to loading depends on direction. Because anisotropy may influence injury mechanisms or reflect tissue health, an understanding of the anisotropic behavior of these materials is important. Despite the prevalence of anisotropy in soft tissues like brain and muscle, and the potential importance of anisotropy in engineered scaffolds and other biomaterials, anisotropic mechanical properties of soft materials are still an active area of research. This is largely due to the technical challenges of anisotropic property estimation.

Shear wave behavior can be used to probe the mechanical properties of elastic anisotropic materials. This approach has been explored in theoretical studies (e.g., Ref. [1]) and experimental studies using ultrasound elastography [25], and magnetic resonance elastography (MRE) [615]. Fibrous materials are often modeled as transversely isotropic (TI) materials (materials with a single axis of symmetry, usually associated with a dominant fiber direction). TI materials can exhibit shear and tensile anisotropy, both of which are important in determining shear-wave behavior [16]. In general, TI materials require five independent parameters to specify their behavior; e.g., two shear moduli, two tensile moduli, and a bulk modulus. If the TI material is incompressible, the number of independent parameters is reduced to three because incompressibility constrains the relationships between the moduli [17].

The use of imaging (magnetic resonance (MR) or ultrasound) to measure shear waves and estimate parameters is known as elastography. Elastography was first used to estimate two elastic parameters of a TI material model: the shear moduli governing shear in planes parallel and perpendicular to the fiber direction. Gennisson et al. [3], used ultrasound elastography to study transversely isotropic phantoms and measured two shear moduli parallel and perpendicular to the fibers. Related experimental ultrasound studies [4,5] uncovered two different shear-wave speeds in transversely isotropic phantoms, but the intrinsic material parameters were not obtained. MRE studies have been performed to estimate these two shear moduli in breast tissue [6], muscle tissue [711], anisotropic phantoms [10], and aligned fibrin gels [12]. However, these two-parameter models are incomplete, as TI materials can also have different tensile moduli [17]. Other theoretical studies include finite element simulations of wave behavior in TI material models [1], analytical studies [2], and models that incorporate anisotropy for applications such as injury prediction [18].

Magnetic resonance elastography can also be used to estimate three parameters (e.g., a baseline shear modulus, a shear anisotropy ratio, and a tensile anisotropy ratio) to fully parameterize incompressible (ITI) or nearly incompressible (NITI) TI material models [9,13,14]. It is also possible to try to estimate five parameters for general TI material models [15], but in soft biological materials, the ratio between bulk modulus and other moduli is extremely large so that estimating the additional parameters of the five-parameter model is both more difficult and less important for predicting deformations.

For any anisotropic material model, a sufficiently rich dataset is needed to obtain accurate parameter estimates. Tweten et al. [17], showed, using simulations, that two types of shear waves (“pure transverse” or slow, and “quasi-transverse” or fast) must be measured, with propagation of both waves in different directions, to accurately estimate all three material parameters. Illustrating this idea, Anderson et al. [19] used multiple excitation methods and showed that parameter estimates of an isotropic material model depended on the direction and location of actuation. Schmidt et al. [16] actuated tissue in two different directions to induce and measure the two shear wave types described by Tweten et al. [17].

1.2 Basic Theory: Waves in ITI or Nearly Incompressible Transversely Isotropic Materials.

Methods to estimate intrinsic properties of anisotropic material from slow and fast shear wave speeds rely on concepts from elasticity and acoustics. Specifically, an ITI material can be described by three independent parameters: shear modulus (μ), shear anisotropy (ϕ), and tensile anisotropy (ζ). These three parameters can be estimated from shear-wave speeds, with known propagation direction, polarization direction, and fiber direction. The slow and fast shear-wave polarization directions (ms and mf, respectively) are obtained from the eigenvectors of the acoustic tensor [17]. They are
ms=n×a|n×a|
(1)
mf=n×ms
(2)

where n is the propagation direction and a is the fiber direction. These polarization direction unit vectors can be used to isolate slow and fast shear waves in the displacement field [16,17].

Shear wave speeds are found from the eigenvalues of the acoustic tensor of the ITI material. Slow shear waves do not induce stretch along the fiber direction in the ITI material. Therefore, the slow shear wave speed (cs) depends only on μ, ϕ, density (ρ), and the angle between the fiber direction and the propagation direction (θ)
ρcs2=μ(1+ϕcos2θ)
(3)
Fast shear waves induce stretch in the fibers of the ITI material, so the fast shear wave speed (cf) also depends on tensile anisotropy (ζ)
ρcf2=μ(1+ϕcos22θ+ζsin22θ)
(4a)
or
ρcf2=μ(1+ϕ)+μ(ζϕ)sin22θ
(4b)
Multiple shear wave speeds (cs,cf) can be measured for a range of propagation angles (θ), which can be obtained, in principle, by varying actuation direction or location. To estimate these three unknown material parameters, measurements of both shear-wave speeds (cs,cf) at multiple different angles of propagation (θ) provide an overdetermined set of equations for the material parameters
[1 cos2θ101cos2θn01cos22θn+1sin22θn+11cos22θn+msin22θn+m][μμϕμζ]=ρ[cs,12cs,n2cf,n+12cf,n+m2]
(5)

where n is the number of slow shear wave measurements and m is the number of fast shear-wave measurements. The system is solved in the least-squares sense to estimate the material properties μ,ϕ,andζ.

1.3 Practical Challenges in Anisotropic Elastography.

Comprehensive characterization of anisotropic tissue by wave behavior requires both slow and fast shear waves with significant amplitudes in multiple propagation and polarization directions [17]. Meeting these requirements is challenging because 3D displacement fields are necessary to observe waves with multiple propagation and polarization directions. (Note that 3D displacement fields are also needed to characterize anisotropic stress–strain relationships by quasi-static methods, as in the recent study by Dargar et al. [20]) Previous experimental studies that investigated anisotropy using slow and fast shear waves employed multiple samples in different experimental setups [16]. To perform anisotropic MRE using only one sample, multiple shear wave directions can be induced by varying frequency, actuator placement, or actuation direction [19]. External (boundary) actuation is noninvasive, but produces shear waves that are largely uncontrolled in direction and vulnerable to attenuation. Direct (invasive) actuation, for example, with an embedded needle or rod (stinger), can produce higher amplitude waves deep in tissue, but it is destructive, so actuating in multiple directions is not possible due to cumulative damage to the sample [16]. In principle, pure ultrasound elastography (ultrasound imaging of ultrasound-induced displacement) is noninvasive and incorporates the ability to actuate in multiple directions, but it does not provide the 3D displacement fields necessary to fully characterize material behavior.

It can also be beneficial to perform anisotropic parameter estimation locally (within a small regional volume) to address tissue heterogeneity and attenuation. In some biological tissues, like white matter in the brain, the length scales of the tissue heterogeneities are relatively small. For anisotropic MRE, accurate parameter estimation requires a small wavelength (high frequency). In MRE, external actuation at high frequencies is susceptible to attenuation, so waves may not reach deep tissue far from the boundary. Ultrasound elastography is also not ideal for small sample volumes because of limits in spatial resolution.

Challenges with existing actuation methods for anisotropic parameter estimation can be summarized thus: (i) boundary actuation is noninvasive, but it involves uncontrolled propagation and polarization directions and high attenuation, especially at high frequencies. (ii) Direct internal actuation is invasive and does not allow for multiple propagation directions per sample. (iii) Ultrasound elastography has low resolution and does not provide a complete 3D displacement field. These limitations motivate the pursuit of an MRE approach with a noninvasive, local, and nondestructive actuator to study anisotropy.

1.4 Potential Solutions Using Focused Ultrasound.

Focused ultrasound (FUS) can be used to address some of these issues. FUS concentrates ultrasound energy into a region of 0.5–3 mm3 [21,22] deep inside tissue (Fig. 1). The ultrasound wave propagation generates a force called “acoustic radiation force,” which, in simplified terms, can be is described by the equation
F(t)=2αI(t)c
(6)

where F(t) is the volumetric radiation force (kg/(s2m2)), α is the tissue absorption coefficient (m1), I(t)is the average acoustic intensity (W/m2), and c is the speed of sound in the material (m/s).

Fig. 1
(a) Focusing of ultrasound waves creates a focal region of increased acoustic radiation force (F). (b) Periodic modulation of the ultrasound waveform creates an oscillatory acoustic radiation force and generates shear waves that propagate radially outward from the focus.
Fig. 1
(a) Focusing of ultrasound waves creates a focal region of increased acoustic radiation force (F). (b) Periodic modulation of the ultrasound waveform creates an oscillatory acoustic radiation force and generates shear waves that propagate radially outward from the focus.

Acoustic radiation force can be used to remotely interrogate tissue mechanical properties by generating impulsive radiation force or harmonic radiation force. Shear wave displacement and velocity from acoustic radiation force can be imaged by ultrasound [23] or MRI [2426]. Prior work in this area includes acoustic radiation force imaging (MR-ARFI) [26,27], transient MRE (t-MRE) [22], harmonic motion imaging (HMI) [21,28], and multiple point ARFI (mpARFI) [29]. In MR-ARFI [26,27] and mpARFI [29], MRI is used to map the displacement field induced by focused ultrasound. In HMI [21,28], ultrasound is used for both harmonic actuation and imaging of displacement. In general, MRI provides greater higher resolution and sensitivity to motion. Isotropic elastic parameters have been previously been estimated from MRI measurement of ultrasound-induced harmonic [2426] or transient [22] shear waves.

The approach used in this study, which we describe as MR imaging of harmonic ultrasound-induced motion (MR-HUM), can overcome several challenges of anisotropic MRE. Harmonic acoustic radiation forces induce oscillations at the focal point and the resulting shear waves propagate radially outward from the focus with a temporal frequency equal to the frequency of the amplitude (square wave, or “on–off”) modulation. By noninvasively producing and imaging shear waves with different propagation and polarization directions, in a region of interest, MR-HUM provides rich data sets that can be inverted to obtain local estimates of anisotropic material parameters. In this study, shear wave speeds of slow and fast shear waves were estimated from their respective phase gradients (PGs). The combined approach (experiment and analysis) is demonstrated by estimating parameters from experimental data from muscle tissue (chicken breast) ex vivo, and using data from simulations of the experiment.

Magnetic resonance imaging of harmonic ultrasound-induced motion can overcome several challenges of anisotropic parameter estimation by noninvasively producing and imaging shear waves with multiple propagation and polarization directions, with small enough wavelengths to produce local estimates of material parameters in a single sample. This study is the first to use MR-HUM to comprehensively and quantitatively characterize anisotropic material properties of a fibrous, soft biological tissue. This is also the first study to model MR-HUM to simulate anisotropic wave propagation and to compare these results to experiments.

2 Methods

Magnetic resonance imaging of harmonic ultrasound-induced motion was performed on muscle tissue (chicken breast) ex vivo. Cylindrical samples were punched from the tissue. The sample could be rotated with respect to the ultrasound transducer. These ensured waves with multiple propagation and polarization directions were obtained in the same sample. Scans were performed at room temperature (∼21 °C) on an Agilent/Varian DirectDrive 4.7 T small-bore animal MRI scanner using a custom-built FUS system (Image Guided Therapy, Pessac, France). Experiments were performed in eleven samples to obtain valid statistical estimates of the expected values of parameters and their variation. Three samples were eliminated from analysis due to problems in scanning, including air bubbles in the sample (which attenuate ultrasound and degrade MR image quality), scanner malfunction, or lack of motion perpendicular to fibers (needed to characterize slow shear waves).

2.1 Sample Preparation.

Chicken breast purchased from a local grocery store was frozen within one day of purchase. The chicken breast was removed from the freezer to thaw at room temperature for ∼1 h. Once the tissue was partially thawed, a 1″ circular punch (McMaster Carr, part 3427A24) was used to cut cylindrical samples (Fig. 2(a)). Each sample was embedded in a gelatin-glycerol gel [30] inside a 50 mL conical tube with an inside diameter of 27 mm and lubricated with canola oil to prevent adhesion to the container (Fig. 2(b)). The sample was then refrigerated until testing. Prior to testing, the chicken/gel sample was removed from the 50 mL tube and inserted into another 50 mL conical tube that had been modified by removing a 25 mm axial section (110 deg of the circumference), creating a window that allowed ultrasound penetration without attenuation. This tube was then placed in a 30-mm diameter coil for scanning (Fig. 2(c)). The ultrasound transducer was placed above the surface of the sample (Fig. 2(d)). A water-filled bladder attached to the transducer provided an acoustic coupling between the transducer and the sample. The focus of the ultrasound transducer was electronically positioned to excite an ellipsoidal region of tissue within the chicken breast. Figure 2(e) shows a schematic diagram of the experimental setup.

Fig. 2
Sample preparation and schematic for MR-HUM. (a) 1″ diameter cylindrical sample of muscle tissue (chicken breast). (b) Sample embedded in gelatin/glycerol mixture for testing. (c) Sample transferred to a 50 mL tube with a cutout window. The tube is placed in the 30-mm diameter RF coil with the cutout facing upward. (d) The ultrasound transducer is placed above the sample. A water bladder covering the transducer provides an acoustic coupling between the transducer and the sample. The sample can be rotated while maintaining contact between the transducer and the sample. (e and f) Schematic of the setup (two cross-sectional views). Black scale bar in (a–c and f) is 10 mm.
Fig. 2
Sample preparation and schematic for MR-HUM. (a) 1″ diameter cylindrical sample of muscle tissue (chicken breast). (b) Sample embedded in gelatin/glycerol mixture for testing. (c) Sample transferred to a 50 mL tube with a cutout window. The tube is placed in the 30-mm diameter RF coil with the cutout facing upward. (d) The ultrasound transducer is placed above the sample. A water bladder covering the transducer provides an acoustic coupling between the transducer and the sample. The sample can be rotated while maintaining contact between the transducer and the sample. (e and f) Schematic of the setup (two cross-sectional views). Black scale bar in (a–c and f) is 10 mm.

During testing, the sample could be rotated in the coil while the ultrasound transducer remained stationary. The water-filled bladder coupling the sample to the transducer was required to remain in the tube cutout area (Fig. 3). The sample rotation controlled the angle between the muscle fibers and the direction of the acoustic radiation force. For this study, two MR-HUM scans were performed on each chicken sample, with actuation at different angles relative to the fiber direction (Fig. 3). The two angles were chosen to be (i) roughly perpendicular to the fiber axis and (ii) at an angle near 45 deg; the method is robust to variations in the precise angle, which is only coarsely controlled by sample rotation. In practice the differences between angles ranged from 23 deg to 53 deg (37 deg ± 11 deg). A total of 11 samples of chicken were imaged in this study.

Fig. 3
The sample was rotated within the coil to change the angle between the fibers and direction of actuation (β). The transducer and focal region of the ultrasound beam remained stationary. This schematic shows samples undergoing actuation at angles of approximately β = 90 deg and β = 45 deg to the muscle fibers. Blue curves indicate the tube with the cutout window.
Fig. 3
The sample was rotated within the coil to change the angle between the fibers and direction of actuation (β). The transducer and focal region of the ultrasound beam remained stationary. This schematic shows samples undergoing actuation at angles of approximately β = 90 deg and β = 45 deg to the muscle fibers. Blue curves indicate the tube with the cutout window.

2.2 Imaging Procedures

2.2.1 MR-HUM.

The tissue was harmonically oscillated by modulating the acoustic radiation force of the focused ultrasound beam. The ultrasound transducer generated a 1500 kHz focused beam modulated by a square wave at 400 Hz (Fig. 1), which produced shear waves in the sample at the frequency of the modulation with a peak positive pressure of 1 MPa. The frequency, 400 Hz, was chosen to induce waves with multiple voxels per wavelength and multiple wavelengths in the image domain. MRE data were acquired with a modified 2D multislice spin-echo sequence [31] with 1 mm isotropic voxels, TR = 1000 ms, and TE = 33 − 34 ms covering a volume of either 32 × 32 × 12 mm3 or 24 × 24 × 12 mm3. Sinusoidal motion-encoding gradients (four cycles) of amplitude 20 G/cm were synchronized with the harmonic FUS excitation to induce phase contrast proportional to displacement.

Magnetic resonance elastography data (phase-contrast images) were phase-unwrapped and converted to displacement, and the effects of rigid body motion were subtracted. During analysis, imaging data were masked at 10-mm radius sphere from the center of actuation (data outside this region were excluded) because MR-HUM shear waves at this frequency dissipate quickly with distance from the focus.

2.2.2 Diffusion Tensor Imaging.

Diffusion tensor imaging (DTI) was performed for all samples at all orientations tested. Diffusion tensors were estimated from diffusion-weighted images obtained using 31 diffusion-weighting gradient directions and two averages [32]. The scan used 2-mm isotropic voxel resolution over an imaging volume of 48 × 48 × 15 mm3. Fractional anisotropy (FA) was estimated from diffusion tensor eigenvalues, and fiber direction (a) was estimated from the first principal eigenvector [33].

2.3 Modeling and Simulation.

A finite element model (comsolmultiphysics; v. 5.3a, Stockholm, Sweden) of a NITI cylinder was used to simulate MR-HUM in anisotropic tissue similar to the experimental sample. The data from these ideal situations were used to evaluate the anisotropic parameter estimation methods (described below) based on the PG of the shear wave field.

The simulation domain was a linear, elastic, nearly incompressible, homogenous cylinder of 27 mm diameter and 50 mm length (Figs. 4(a) and 4(b); dimensions were chosen to match experimental samples). The domain was discretized into 100,505 quadratic Lagrange elements, corresponding to 432,883 degrees-of-freedom. The boundaries of the cylinder were rigid. A harmonic body force was applied to a spherical (1 mm radius) region at a single location; the force was oriented at different angles with respect to fiber direction in each of five simulations by rotating the fiber direction (Figs. 4(a)4(c)). For each simulation, the cylinder material had one fiber direction with an angle of β=0deg,30deg,45deg,60deg,or90deg relative to the actuation direction (Fig. 4(c)). The harmonic body load was applied at a single frequency in the z-direction. The solution for the steady-state frequency response was found using comsol's frequency domain solver. Displacement data from the simulation were exported into matlab and interpolated onto a 3D grid with 1 mm3 voxel resolution for analysis using the LiveLink feature of comsol (“mphinterp” command). The harmonic body load produced shear waves propagating with approximately spherical wave fronts outward from the center of the cylinder.

Fig. 4
Simulation of MR-HUM. (a and b) A body load is applied to the small spherical region in the center of the cylinder (a,x-y view; b, x-z view). The β = 90 deg case is shown. (c) Five models for simulation of MR-HUM, showing the fiber direction at β = 90, 60, 45, 30, and 0 deg to the actuation (z-) direction.
Fig. 4
Simulation of MR-HUM. (a and b) A body load is applied to the small spherical region in the center of the cylinder (a,x-y view; b, x-z view). The β = 90 deg case is shown. (c) Five models for simulation of MR-HUM, showing the fiber direction at β = 90, 60, 45, 30, and 0 deg to the actuation (z-) direction.

Analysis of data from the simulations was performed using only data from the spherical region within a 10-mm radius of the center of the cylinder (location of the harmonic body load) to eliminate effects of wave dissipation and reflections from boundaries. Figure 4 shows the spherical region enclosed in dotted lines (a and b) and the five simulation orientations (c).

Two sets of material properties were modeled in separate sets of simulations. One simulation set incorporated approximately brain-like shear modulus, with material parameters μ=2000Pa,ϕ=1,andζ=2, and a body force density at the focal region of 50 kN/m3 at 300 Hz. A second simulation set incorporated stiffer, approximately muscle-like, shear modulus, with parameters of μ=7500Pa,ϕ=1,andζ=1, and a force density of 150 kN/m3 at 400 Hz. The second simulation set approximates the MR-HUM chicken breast experiment. The actuation body forces were chosen to produce micron-level displacements in the simulation, similar to those seen in the experiment. Data from these simulations, which are discretized but noise-free, representing an idealized scenario, were used to evaluate the phase-gradient method for anisotropic property estimation.

2.4 Anisotropic Parameter Estimation

2.4.1 Estimation of Wave Speeds and Apparent Shear Moduli From PG.

The wave propagation direction (n) of the data was assumed to be purely radial, emanating from the ultrasound focus. Slow and fast polarization directions (ms and mf,respectively)and propagation-fiber angle, θ, were calculated using the assumed propagation direction and the fiber direction (a) of the sample/simulation using Eqs. (1) and (2). The curl of the displacement was then sorted into slow and fast components (note that curl polarizations are perpendicular to corresponding displacement polarizations)
Γs=Γ·mf
(7)
Γf=Γ·ms
(8)
Next, the phase angles of the slow and fast (ψsandψf) waves were calculated
ψs=Γs
(9)
ψf=Γf
(10)
The wave numbers (ks and kf) were estimated from the gradients of phase
ks=ψs
(11)
kf=ψf
(12)
Wavelength for slow and fast waves (λs and λf) were calculated from the wave numbers
λs=2π|ks|
(13)
λf=2π|kf|
(14)
Apparent shear modulus was calculated from wavelength, using the density value of ρ=1000 (kg/m3) and the frequency of the actuation (f)
μapp=μs=ρ(λsf)2
(15)
μapp=μf=ρ(λff)2
(16)

2.4.2 Voxel Inclusion Criteria and Classification of Voxels as Slow or Fast.

For a voxel to be included in the analysis, multiple conditions must be met to ensure that approximations and assumptions are accurate in that voxel: (i) The voxel must experience a shear wave amplitude above a specified threshold; (ii) the voxel must be within a specified radius of the center of actuation; (iii) the propagation direction within the voxel must be close to that of radially propagating waves, and (iv) the voxel must have a fractional (diffusion) anisotropy above a specified threshold. Table 1 summarizes the inclusion criteria for the analysis.

Table 1

Inclusion criteria for analysis of anisotropic parameter estimation for both simulations and experiments. Parameters were chosen to be consistent with experimental studies, which had lower wave amplitude and generally low FA.

Inclusion criteriaEquationParameter
Amplitude|U|>A|U|medianA=0.1
Propagation directionn·er>propthreshpropthresh=0.75
Radial distanceRmin<r<RmaxRmin=2mmRmax=10mm
Fraction anisotropyFA>FAthreshFAthresh=0.01
Inclusion criteriaEquationParameter
Amplitude|U|>A|U|medianA=0.1
Propagation directionn·er>propthreshpropthresh=0.75
Radial distanceRmin<r<RmaxRmin=2mmRmax=10mm
Fraction anisotropyFA>FAthreshFAthresh=0.01
After voxels are selected based on these inclusion criteria, they must also meet classification criteria to be identified as either a slow or fast voxel, as determined by the dominant shear wave polarization in that voxel. A voxel was classified as a fast or slow shear wave voxel if the normalized curl component in the fast or slow polarization direction exceeded a minimum “polarization threshold” (polthresh). The normalized fast and slow shear wave curl components are
Γ̂f=Γ·ms|Γ|
(17)
Γ̂s=Γ·mf|Γ|
(18)

(directions of curl polarization are perpendicular to displacement polarization directions). Thus, a voxel would be classified as fast if Γ̂f>polthresh and as slow if Γ̂s>polthresh. Voxels that did not meet either of these criteria were excluded from the analysis. Table 2 outlines the classification criteria used for the PG method.

Table 2

Classification criteria for PG analysis

Classification criteria for PGEquationParameter
Polarization direction—slowΓ̂s>polthreshpolthresh=0.75
Polarization direction—fastΓ̂f>polthreshpolthresh=0.75
Classification criteria for PGEquationParameter
Polarization direction—slowΓ̂s>polthreshpolthresh=0.75
Polarization direction—fastΓ̂f>polthreshpolthresh=0.75

2.4.3 Parameter Estimation by Multiple Linear Regression.

Phase gradient was used to estimate the material properties of the chicken breast samples. This approach was used, as described above, to separate the waves by polarization direction (slow and fast) and to approximate the wave speeds (cs and cf) and apparent shear modulus (μapp=ρc2) for each type. The three unknown parameters of an NITI material were then estimated from the equations for slow and fast shear waves using a multiple linear regression model of the form
y=β0+β1x1+β2x2
(19)
The unknown parameters are β0=μ, β1=μϕ and β2=μζ. The dependent variable is the apparent shear modulus: y=μapp (μapp=μs for slow voxels and μapp=μf for fast voxels). The independent variables in the multiple regression are defined in terms of the angle θ as follows:
x1=cos2θ(slowvoxels)cos22θ(fastvoxels)andx2=0(slowvoxels)sin22θ(fastvoxels)
(20)

Voxels classified as either slow or fast were included in the fit based on the criteria in Table 2. These values and the corresponding values of the independent variables for slow and fast voxels were used in the data sets fitted by the linear regression equation (Eq. (19)) to estimate values of μ, μϕ, and μζ. Figure 5 outlines all the steps of this PG-based estimation method.

Fig. 5
Flow chart outlining the steps of shear wave separation and anisotropic parameter estimation using PG
Fig. 5
Flow chart outlining the steps of shear wave separation and anisotropic parameter estimation using PG

2.5 Direct Mechanical Testing.

For comparison with parameter estimates from MR-HUM, estimates of shear moduli were obtained by dynamic shear testing (DST) in a limited number of samples, and estimates of direction-dependent tensile moduli were obtained (in different samples) by biaxial testing. These techniques are well-established; our implementation is described briefly in Secs. 2.5.1 and 2.5.2 [13].

2.5.1 Dynamic Shear Testing.

Seven thin, round samples (15.6 mm dia, 2.9–3.9 mm thick) were tested by DST using a custom-built instrument [13]. The sample was placed in a holder on a flexure-mounted platform, with the dominant fiber axis (determined visually) aligned either parallel or perpendicular to the axis of motion. The flexure was oscillated by a voice coil actuator at frequencies from 20–100 Hz (chirp). An upper platen instrumented with two force transducers (209C11, PCB Piezotronics, Depew, NY) was lowered until contact was established, and then lowered further to achieve 10% compression. Shear stress was calculated from the measured horizontal traction force on the sample, divided by sample area, and shear strain from the horizontal displacement. The shear modulus was computed from the average ratio of shear stress to shear strain at 45 Hz. The sample was then rotated 90 deg and the shear modulus in the other direction was measured similarly.

2.5.2 Biaxial Testing.

Four approximately square samples (25–31 mm square, 10–20 mm thick) were tested in a planar biaxial tensile test machine (TestResources, Shakopee, MN) with custom grips [34]. Samples were loaded in strip biaxial configuration with 5% strain applied either parallel or perpendicular to the fiber axis (estimated visually) while the other dimension (perpendicular or parallel, respectively) was held fixed. An equibiaxial test, in which both dimensions were subjected to 5% strain, was also performed on all samples. A preload of 0.5 N was applied before testing. Cyclic loading for ten cycles at 0.5 Hz to the maximum strain of 5% was applied, and force was recorded from each of the four load cells on the instrument. Estimates of the ratio of the two tensile moduli were obtained from the results, using the equations of linear, elastic, transversely isotropic materials under planar biaxial loading.

3 Results

3.1 Fiber Direction and Displacement Fields

3.1.1 Experiment.

Fiber directions, estimated using DTI, are shown for one sample in Fig. 6. In this figure and subsequent figures, data fields are masked to show only the region within 10 mm radius of the focus. Figure 6 shows the sample area (dotted line) and fiber directions from one sample with fibers at 51 deg and 87 deg to the actuation direction. This sample contains fibers with a consistent, clear orientation.

Fig. 6
Fiber directions estimated by DTI from one sample at two different angles relative to the actuation direction. (a) Schematic diagram. The region of the sample outlined by dotted lines (top) is the partial sphere of 10 mm radius centered about the focal region used in the analysis. The sample was rotated 36 deg between the two experiments. (b–e) DTI estimates of fiber direction are displayed for multiple views for the two orientations: (b and c) β = 51 deg and (d and e) β = 87 deg.
Fig. 6
Fiber directions estimated by DTI from one sample at two different angles relative to the actuation direction. (a) Schematic diagram. The region of the sample outlined by dotted lines (top) is the partial sphere of 10 mm radius centered about the focal region used in the analysis. The sample was rotated 36 deg between the two experiments. (b–e) DTI estimates of fiber direction are displayed for multiple views for the two orientations: (b and c) β = 51 deg and (d and e) β = 87 deg.

All three components (U,V,W) of the shear-wave displacement field are shown in panels (ac) of Figs. 7 and 8, for the same sample depicted in Fig. 6. Displacement fields were consistent with expected wave behavior in an ITI or NITI material. Noncircular waves were observed for all samples; typically, the wave fronts were elliptical with the major semi-axis aligned with the fiber direction from DTI. Propagation direction was estimated from the wave fields using an array of directional filters [31,35]. Slow and fast polarization directions were calculated from propagation direction and fiber direction. Figures 7 and 8 (panels dg) show the results from this directional filtering analysis for the chicken breast sample of Fig. 6, on a slice near the center of actuation. As noted above, actuation was 51 deg and 87 deg, respectively, to the fiber direction at 400 Hz. The fiber direction (a), propagation direction (n), slow polarization direction (ms), and fast polarization direction (mf) are shown. Voxels within the specified distance from the focus are displayed.

Fig. 7
MR-HUM data for chicken breast sample with actuation direction 87 deg to the fiber direction. (a–c) Shear wave displacement components (U, V, W) for a slice near the center of actuation. (d–g) Fiber direction (a), propagation direction (n), slow polarization direction (ms), and fast polarization direction (mf) are shown for the sample. Components are encoded by color: red = x; green = y; blue = z. Data within a 10 mm radius from the focus are shown. Scale bar in (c) is 5 mm.
Fig. 7
MR-HUM data for chicken breast sample with actuation direction 87 deg to the fiber direction. (a–c) Shear wave displacement components (U, V, W) for a slice near the center of actuation. (d–g) Fiber direction (a), propagation direction (n), slow polarization direction (ms), and fast polarization direction (mf) are shown for the sample. Components are encoded by color: red = x; green = y; blue = z. Data within a 10 mm radius from the focus are shown. Scale bar in (c) is 5 mm.
Fig. 8
MR-HUM data for chicken breast sample with actuation direction 51 deg to the fiber direction for directional filtering analysis. (a–c) Shear wave displacement components (U, V, W) for a slice near the center of actuation. (d–g) Fiber direction (a), propagation direction (n), slow polarization direction (ms), and fast polarization direction (mf) are shown for the sample. Components are encoded by color: red = x; green = y; blue = z. Data within a 10 mm radius from the focus are shown. Scale bar in (c) is 5 mm.
Fig. 8
MR-HUM data for chicken breast sample with actuation direction 51 deg to the fiber direction for directional filtering analysis. (a–c) Shear wave displacement components (U, V, W) for a slice near the center of actuation. (d–g) Fiber direction (a), propagation direction (n), slow polarization direction (ms), and fast polarization direction (mf) are shown for the sample. Components are encoded by color: red = x; green = y; blue = z. Data within a 10 mm radius from the focus are shown. Scale bar in (c) is 5 mm.

3.1.2 Fiber Direction and Displacement Fields in Simulation.

Figures 9 and 10 show displacement fields, fiber directions, and propagation and polarization directions of shear waves in simulations of muscle-like tissue (analogous to experimental Figs. 7 and 8). As in the experiment, noncircular wave fronts propagate radially outward from the location of the oscillatory body force. Both the displacement fields and the propagation and polarization direction fields closely resemble those observed in the experiment.

Fig. 9
MR-HUM data for simulation of muscle-like tissue with actuation direction 90 deg to the fiber direction. (a–c) Shear wave displacement components (U, V, W) for a slice at the center of actuation. (d–g) Fiber direction (a), propagation direction (n), slow polarization direction (ms), and fast polarization direction (mf) are shown for the simulation. Components are encoded by color: red = x; green = y; blue = z. Data within a 10-mm radius from the focus are shown. Scale bar in (c) is 5 mm.
Fig. 9
MR-HUM data for simulation of muscle-like tissue with actuation direction 90 deg to the fiber direction. (a–c) Shear wave displacement components (U, V, W) for a slice at the center of actuation. (d–g) Fiber direction (a), propagation direction (n), slow polarization direction (ms), and fast polarization direction (mf) are shown for the simulation. Components are encoded by color: red = x; green = y; blue = z. Data within a 10-mm radius from the focus are shown. Scale bar in (c) is 5 mm.
Fig. 10
MR-HUM data for simulation of muscle-like tissue with actuation direction 45 deg to the fiber direction. (a–c) Shear wave displacement components (U, V, W) for a slice at the center of actuation. (d–g) Fiber direction (a), propagation direction (n), slow polarization direction (ms), and fast polarization direction (mf) are shown for the simulation. Components are encoded by color: red = x; green = y; blue = z. Data within a 10-mm radius from the focus are shown. Scale bar in (c) is .5 mm.
Fig. 10
MR-HUM data for simulation of muscle-like tissue with actuation direction 45 deg to the fiber direction. (a–c) Shear wave displacement components (U, V, W) for a slice at the center of actuation. (d–g) Fiber direction (a), propagation direction (n), slow polarization direction (ms), and fast polarization direction (mf) are shown for the simulation. Components are encoded by color: red = x; green = y; blue = z. Data within a 10-mm radius from the focus are shown. Scale bar in (c) is .5 mm.

3.2 Parameter Estimates From Phase Gradient of Slow and Fast Shear Waves.

3.2.1 Parameter Estimates From Simulation.

Parameter estimates obtained from simulation data are shown first to demonstrate the performance of the method under close-to-ideal conditions. Voxels were first classified into slow and fast categories based on polarization direction (c). Figure 11 shows the results of initial voxel classification for displacement (U: row 1, panels a and b, k and l) and curl fields (Γ: row 2, panels c and d, m and n) with amplitude thresholding for the two cases β=45deg and β=90deg at 400 Hz. The phase angle (ψ) of each shear wave component was calculated from the curl component. Figure 11 row 3, panels (e and f) and (o and p), shows the phase fields for the two cases, with arrows representing the propagation direction for the voxels that meet all the criteria for slow or fast waves. (Note: there are no arrows on voxels categorized as fast for the β=90deg case because those voxels did not meet criteria for inclusion). The angle (θ) between propagation direction and fiber direction, and the apparent shear modulus (μapp) are shown for slow and fast voxels that met inclusion criteria, in Fig. 11 Rows 4 and 5, panels (gj) and (qt). The majority of the voxels included in the β=90deg case (all those in the slice depicted here) were classified as slow.

Fig. 11
PG analysis of data from simulation of NITI cylinder with focal actuation oriented at 45 deg (a–j) and 90 deg (k–t) to fiber direction at 400 Hz. Voxels were included based on criteria from Table 1. Images are from the center slice normal to the z-axis. Row 1: Displacement fields. (a, k) Displacement-field component (Us) due to shear waves with slow polarization. (b, l) Displacement-field component (Uf) due to shear waves with fast polarization. Row 2: Curl fields. (c, m) Curl-field component (Γs) due to shear waves with slow polarization. (d, n) Curl-field component (Γf) due to shear waves with fast polarization. Row 3: Phase fields. (e, o) Phase angle (ψ) of slow shear wave curl field, Γs. Black arrows represent the propagation direction. Arrows only appear over voxels that meet the criteria for inclusion in the analysis (see Table 2). (f, p) Phase angle (ψ) of fast shear-wave curl field, Γf. In panel (p), there are no arrows showing the propagation direction because no fast voxels for this case meet the criteria for inclusion (see Table 2). Row 4: Fiber-propagation angle. (g, q) Angle between propagation direction and fiber direction (θ) for slow voxels that met inclusion criteria. (h, r) Angle between propagation direction and fiber direction (θ) for fast voxels that met inclusion criteria. Row 5: Apparent shear modulus. (i, s) Apparent shear modulus (μapp) in slow voxels that met inclusion criteria. (j, t) Apparent shear modulus (μapp) categorized by fast polarization. In panels (r) and (t), no fast voxels met inclusion criteria. Scale bar in (l) is 5 mm.
Fig. 11
PG analysis of data from simulation of NITI cylinder with focal actuation oriented at 45 deg (a–j) and 90 deg (k–t) to fiber direction at 400 Hz. Voxels were included based on criteria from Table 1. Images are from the center slice normal to the z-axis. Row 1: Displacement fields. (a, k) Displacement-field component (Us) due to shear waves with slow polarization. (b, l) Displacement-field component (Uf) due to shear waves with fast polarization. Row 2: Curl fields. (c, m) Curl-field component (Γs) due to shear waves with slow polarization. (d, n) Curl-field component (Γf) due to shear waves with fast polarization. Row 3: Phase fields. (e, o) Phase angle (ψ) of slow shear wave curl field, Γs. Black arrows represent the propagation direction. Arrows only appear over voxels that meet the criteria for inclusion in the analysis (see Table 2). (f, p) Phase angle (ψ) of fast shear-wave curl field, Γf. In panel (p), there are no arrows showing the propagation direction because no fast voxels for this case meet the criteria for inclusion (see Table 2). Row 4: Fiber-propagation angle. (g, q) Angle between propagation direction and fiber direction (θ) for slow voxels that met inclusion criteria. (h, r) Angle between propagation direction and fiber direction (θ) for fast voxels that met inclusion criteria. Row 5: Apparent shear modulus. (i, s) Apparent shear modulus (μapp) in slow voxels that met inclusion criteria. (j, t) Apparent shear modulus (μapp) categorized by fast polarization. In panels (r) and (t), no fast voxels met inclusion criteria. Scale bar in (l) is 5 mm.

After classification into slow or fast voxels, the apparent shear modulus (μapp) and propagation-fiber angle (θ) were used to estimate the anisotropic material parameters (μ,μϕ,μζ) using the multiple linear regression model (Eq. (19)). The multiple linear regression analysis was performed using the linear regression function (“fitlm”) in the matlab statistics and machine learning toolbox. Figure 12 shows plots of apparent shear modulus versus angle in slow and fast voxels, for all cases of the simulation, for the case of brain-like stiffness at 300 Hz (A-slow voxels and B-fast voxels) and for the case of the muscle-like stiffness at 400 Hz (C-slow voxels and D-fast voxels).

Fig. 12
Apparent shear modulus from all NITI cylinder simulations for all cases, estimated by the PG method. Each dot represents one voxel that met criteria for slow (panels a and c) or fast (panels b and d) voxels in the PG analysis. The black solid line represents the linear relationship expected for the input parameters for brain-like tissue (a and b): μ = 2 kPa, ϕ = 1, ζ = 2 and muscle like tissue (c and d): μ = 7.5 kPa, ϕ = 1, ζ = 1. The black dashed line represents the linear regression model for the estimated material parameters found using PG and multiple linear regression. (a) Apparent shear modulus in slow voxels for all simulation cases for brain-like tissue. (b) Apparent shear modulus in fast voxels for all simulation cases for brain-like tissue. (c) Apparent shear modulus in slow voxels for all simulation cases for muscle-like tissue. (d) Apparent shear modulus in fast voxels for all simulation cases for muscle-like tissue.
Fig. 12
Apparent shear modulus from all NITI cylinder simulations for all cases, estimated by the PG method. Each dot represents one voxel that met criteria for slow (panels a and c) or fast (panels b and d) voxels in the PG analysis. The black solid line represents the linear relationship expected for the input parameters for brain-like tissue (a and b): μ = 2 kPa, ϕ = 1, ζ = 2 and muscle like tissue (c and d): μ = 7.5 kPa, ϕ = 1, ζ = 1. The black dashed line represents the linear regression model for the estimated material parameters found using PG and multiple linear regression. (a) Apparent shear modulus in slow voxels for all simulation cases for brain-like tissue. (b) Apparent shear modulus in fast voxels for all simulation cases for brain-like tissue. (c) Apparent shear modulus in slow voxels for all simulation cases for muscle-like tissue. (d) Apparent shear modulus in fast voxels for all simulation cases for muscle-like tissue.

Table 3 shows anisotropic parameter estimates from PG analysis for both simulation cases: brain-like stiffness and muscle-like stiffness. The inputs, estimated values, and error (with respect to the exact value) are shown.

Table 3

Comparison between parameter values estimated by PG and the exact parameter (input value) for simulations of MR-HUM in brain-like tissue and muscle-like tissue. The “input” column shows the material parameters used for the simulation. The “estimated” column shows the parameter estimates from PG analysis. For the brain-like stiffness simulation, 7867 voxels were used in the linear model fit (R2 = 0.897). For the muscle-like stiffness simulation, 9594 voxels were used in the linear model fit (R2 = 0.816). The parameters μ,μϕ,andμζ are in units of kPa; ϕ and ζ are unitless.

InputEstimatedError (%)
Brain-like stiffnessμ (kPa)22.09+4.7
μϕ (kPa)22.50+25.0
μζ (kPa)43.14−21.4
ϕ11.19+19.4
ζ21.50−25.0
Muscle-like stiffnessμ (kPa)7.58.19+9.2
μϕ (kPa)7.59.44+25.9
μζ (kPa)7.510.7+43.1
ϕ11.15+15.3
ζ11.31+31.1
InputEstimatedError (%)
Brain-like stiffnessμ (kPa)22.09+4.7
μϕ (kPa)22.50+25.0
μζ (kPa)43.14−21.4
ϕ11.19+19.4
ζ21.50−25.0
Muscle-like stiffnessμ (kPa)7.58.19+9.2
μϕ (kPa)7.59.44+25.9
μζ (kPa)7.510.7+43.1
ϕ11.15+15.3
ζ11.31+31.1

3.2.2 Parameter Estimates From Experimental Data.

A total of eight samples were analyzed. Three samples were excluded because of problems with scans, including air bubbles in the sample, scanner malfunction, or if the scan did not include at least one test in which |β90deg|< 30 deg (motion must be close to perpendicular to fibers to excite slow shear waves). Figure 13 shows the results of PG analysis for the chicken breast sample of Fig. 6 (actuation directions of 51 deg and 87 deg to the fiber direction, 400 Hz). Images are from a slice near the center of actuation. Slow and fast components are shown for displacement (U: row 1, panels a and b, k and l) and curl fields (Γ: row 2, panels c and d, m and n), masked by inclusion criteria (Table 1). Phase angle (ψ) was calculated from the slow (Γs) and fast (Γf) curl components. Figure 13 row 3, panels (e and f) and (o and p) show the phase fields of the shear-wave components for the two fiber direction cases; arrows represent the propagation direction for the voxels that meet all inclusion criteria. The angle between propagation and fiber directions (θ), and the apparent shear modulus (μapp) were found for all the included slow or fast voxels (per Tables 1 and 2). Figure 13 Rows 4 and 5, panels (gj) and (qt) show the fiber-propagation angles and apparent shear moduli at one slice for the two cases. The majority of the voxels were classified as fast for the β=51deg case and most voxels were classified as slow for the β=87deg case.

Fig. 13
Results from PG analysis of MR-HUM data from a chicken breast sample with actuation oriented at β = 51 deg (a–j) and β = 87 deg (k–t) to fiber direction at 400 Hz. Voxels were included based on criteria from Table 1. Images are from a central slice normal to the z-axis. Row 1: Displacement fields. (a, k) Displacement-field component (Us) due to shear waves with slow polarization. (b, l) Displacement-field component (Uf) due to shear waves with fast polarization. Row 2: Curl fields. (c, m) Curl-field component (Γs) due to shear waves with slow polarization. (d, n) Curl-field component (Γf) due to shear waves with fast polarization. Row 3: Phase fields. (e, o) Phase angle (ψ) of slow shear wave curl field, Γs. Black arrows represent the propagation direction. Arrows only appear over voxels that meet the criteria for inclusion in the analysis (see Table 2). (f, p) Phase angle (ψ) of fast shear wave curl field, Γf. Row 4: Fiber-propagation angle. (g, q) Angle between propagation direction and fiber direction (θ) for slow voxels that met inclusion criteria. (h, r) Angle between propagation direction and fiber direction (θ) for fast voxels that met inclusion criteria. Row 5: Apparent shear modulus. (i, s) Apparent shear modulus (μapp) in slow voxels that met inclusion criteria. (j, t) Apparent shear modulus (μapp) categorized by fast polarization. Scale bar in (l) is 5 mm.
Fig. 13
Results from PG analysis of MR-HUM data from a chicken breast sample with actuation oriented at β = 51 deg (a–j) and β = 87 deg (k–t) to fiber direction at 400 Hz. Voxels were included based on criteria from Table 1. Images are from a central slice normal to the z-axis. Row 1: Displacement fields. (a, k) Displacement-field component (Us) due to shear waves with slow polarization. (b, l) Displacement-field component (Uf) due to shear waves with fast polarization. Row 2: Curl fields. (c, m) Curl-field component (Γs) due to shear waves with slow polarization. (d, n) Curl-field component (Γf) due to shear waves with fast polarization. Row 3: Phase fields. (e, o) Phase angle (ψ) of slow shear wave curl field, Γs. Black arrows represent the propagation direction. Arrows only appear over voxels that meet the criteria for inclusion in the analysis (see Table 2). (f, p) Phase angle (ψ) of fast shear wave curl field, Γf. Row 4: Fiber-propagation angle. (g, q) Angle between propagation direction and fiber direction (θ) for slow voxels that met inclusion criteria. (h, r) Angle between propagation direction and fiber direction (θ) for fast voxels that met inclusion criteria. Row 5: Apparent shear modulus. (i, s) Apparent shear modulus (μapp) in slow voxels that met inclusion criteria. (j, t) Apparent shear modulus (μapp) categorized by fast polarization. Scale bar in (l) is 5 mm.

After classification into slow and fast categories, data from the respective voxels were used in the multiple linear regression model (matlablmfit” command). Figure 14 shows the apparent shear modulus (μs and μf) versus angle (θ) in voxels classified as slow or fast, in one sample.

Fig. 14
Apparent shear modulus for voxels categorized as slow (μs) or fast (μf) estimated by PG analysis from one sample (two MR-HUM experiments, panels a and b). Each dot represents one voxel that met the slow-wave (a) or fast-wave (b) criteria (Table 2) for PG analysis. The black dashed line represents the multiple linear regression model for the estimated material parameters found using PG. (a) Apparent shear modulus in slow voxels for one sample (data in Figs. 7 and 8). (b) Apparent shear modulus in fast voxels for one sample (data in Figs. 7 and 8).
Fig. 14
Apparent shear modulus for voxels categorized as slow (μs) or fast (μf) estimated by PG analysis from one sample (two MR-HUM experiments, panels a and b). Each dot represents one voxel that met the slow-wave (a) or fast-wave (b) criteria (Table 2) for PG analysis. The black dashed line represents the multiple linear regression model for the estimated material parameters found using PG. (a) Apparent shear modulus in slow voxels for one sample (data in Figs. 7 and 8). (b) Apparent shear modulus in fast voxels for one sample (data in Figs. 7 and 8).

Figure 15 shows the parameter estimates from PG analysis for all of the chicken breast samples. Estimates from each sample are plotted with the 95% confidence interval for that estimate. Table 4 shows the estimated parameters (mean ± std. dev.) using data from all eight samples.

Fig. 15
Results of PG anisotropic parameter estimation for all eight chicken breast samples included. (a) Estimates of μ, μϕ, and μζ for each of the eight samples (dots) are plotted with their 95% confidence intervals. Black diamonds show the mean and standard deviation of parameter estimates from all eight samples.
Fig. 15
Results of PG anisotropic parameter estimation for all eight chicken breast samples included. (a) Estimates of μ, μϕ, and μζ for each of the eight samples (dots) are plotted with their 95% confidence intervals. Black diamonds show the mean and standard deviation of parameter estimates from all eight samples.
Table 4

Average estimated parameter values from the PG analysis of the eight chicken samples. Values are shown with the standard deviation. μ,μϕ,andμζ are in units of kPa; ϕ and ζ are unitless.

ParameterEstimated
μ (kPa)8.95 ± 1.23
μϕ (kPa)5.71 ± 2.97
μζ (kPa)11.7 ± 4.78
ϕ0.668 ± 0.410
ζ1.37 ± 0.679
ParameterEstimated
μ (kPa)8.95 ± 1.23
μϕ (kPa)5.71 ± 2.97
μζ (kPa)11.7 ± 4.78
ϕ0.668 ± 0.410
ζ1.37 ± 0.679

3.3 Results From Direct Mechanical Testing.

Seven samples were tested by DST. Parameter estimates were limited to the range from 25 to 45 Hz, due to the bandwidth of the instrument. In this range, the average baseline shear modulus (the apparent shear modulus when shear was applied perpendicular to the fiber axis) was μDST=6.19±1.71kPa (mean ± std. dev.), and the average shear anisotropy ratio was ϕDST=0.84±0.30. Four samples were tested in biaxial stretch. In these tests the average value of the tensile modulus perpendicular to the fiber axis was E2=13.4±7.58 kPa, the average value of the tensile modulus parallel to the fiber axis was E1=22.2±7.01 kPa, and the average value of the tensile anisotropy parameter was ζbiax=0.93±0.65.

4 Discussion

This paper introduces the use of MR imaging of harmonic, ultrasound-induced motion (MR-HUM) and phase gradient analysis for anisotropic parameter estimation in fibrous, soft materials. The approach is demonstrated using experimental data from soft tissue (ex vivo chicken breast) and data from analogous simulations. MR-HUM has several key features that make it attractive for material characterization. First, the excitation of shear waves is noninvasive and nondestructive. Second, waves in MR-HUM propagate from the center of the actuation with approximately spherical wave fronts, and attenuate with distance from the focus. This allows more accurate estimation of properties within local regions of heterogeneous tissue, like white matter tracts, by reducing wave interactions with surrounding tissues. Third, the direction of actuation of the focal region can be varied with respect to the dominant fiber axis. In this study, the sample was rotated to vary the fiber direction, so one experimental sample could be imaged with multiple directions of actuation. This allowed controlled excitation of both slow and fast shear waves, which is necessary for accurate anisotropic parameter estimation [17]. MR-HUM provides more control over the direction of shear wave propagation and polarization than excitation of an external boundary. MR-HUM also avoids disruption of sample integrity, such as that caused by an embedded needle or thin rod, which prevents more than one direction of actuation from being explored. Thus, MR-HUM can provide estimates of all three parameters of an NITI material in a single sample.

Focused ultrasound, at high power or prolonged exposure, can heat the focal region. In MR-HUM, the square wave (on–off) modulation of the focused ultrasound reduces heating by cutting the duty cycle in half. MRE sequences were optimized to run quickly (3 min, 37 s), and a time interval between scans of ∼10 min was maintained for sample cooling. In principle, MRI can be used to monitor temperature changes, but sample heating was not quantified in this study. All studies were performed at power levels that did not cause obvious physical changes (color, stiffness, temperature) in the sample.

Phase gradient analysis was used to estimate parameters from simulated and experimental MR-HUM data. Simulations allowed us to rigorously assess the ability of the PG method to estimate parameters in the absence of noise or other sources of error in experimental data. Estimates for the baseline shear modulus varied moderately between samples. Part of this variation may be due to true sample-to-sample variability. However, nonideal features of experimental data likely contributed to some of the variation in experimental parameter estimates. “Wrapping” (due to 2π ambiguity) in phase estimates may introduce some error due to artificially high values of phase gradient at discontinuities. Improved unwrapping techniques might increase the accuracy of the PG method. In addition, the phase gradient is computed by numerical differentiation of the phase field; numerical differentiation has intrinsic error due to discretization, and amplifies the effects of noise. However, despite these potential limitations, estimates of baseline shear modulus calculated from PG were consistent with estimates of shear modulus from established MRE methods such as “local direct inversion” [35] or local frequency estimation (LFE) [36].

Using MR-HUM, chicken breast was observed to be moderately anisotropic in terms of both shear and tensile moduli (Fig. 15, Table 4). The experimental results are consistent with previous studies on turkey breast and cardiac muscle, as well as our current direct testing of chicken breast, as shown in Fig. 16, which compares MR-HUM with anisotropic material parameter estimates obtained using different testing methods (denoted by marker shape) for different muscle tissues (denoted by color). Schmidt et al. estimated the anisotropic parameters of turkey breast using LFE of MRE data and DST of thin samples. For MRE/LFE, they estimated μ33kPa,ϕ1.3,andζ9.2 using piezo-electric direct and surface actuation at 800 Hz (the ζ estimate is suspected to be unreliable due to challenges in estimating wavelength in the sample). For DST, they estimated μ4kPaandϕ0.6 at 20–40 Hz [16]. DST of chicken breast samples (n = 7) in this study provided estimates of μ=6.19±1.71kPaandϕ=0.84±0.30 at 25–45 Hz. Humphrey et al. performed biaxial testing of resting cardiac muscle. From their data, we were able to estimate the tensile anisotropy from the linear elastic region of the equibiaxial test as ζ=0.61±0.25 [37]. Our biaxial testing of chicken breast (n = 4) provided estimates of ζ=0.93±0.65.

Fig. 16
Comparison of estimates (± std. dev.) of anisotropic parameters μ, ϕ, and ζ from various testing methods (denoted by marker shape) and muscle types (denoted by color). MR-HUM is the only method that provided estimates of all three parameters from the same sample. (a) Estimates of baseline shear modulus, μ, from DST (chicken and turkey), MR-HUM (chicken), and MRE using LFE (MRE-LFE, turkey [16]). Muscle tissue is viscoelastic, which means μ is expected to increase with frequency. (b) Estimates of shear anisotropy, ϕ, from DST (chicken and turkey [16]), MR-HUM (chicken), and MRE-LFE (turkey [16]). (c) Estimates of tensile anisotropy, ζ, from biaxial testing (chicken and cardiac muscle [37]), and MR-HUM (chicken).
Fig. 16
Comparison of estimates (± std. dev.) of anisotropic parameters μ, ϕ, and ζ from various testing methods (denoted by marker shape) and muscle types (denoted by color). MR-HUM is the only method that provided estimates of all three parameters from the same sample. (a) Estimates of baseline shear modulus, μ, from DST (chicken and turkey), MR-HUM (chicken), and MRE using LFE (MRE-LFE, turkey [16]). Muscle tissue is viscoelastic, which means μ is expected to increase with frequency. (b) Estimates of shear anisotropy, ϕ, from DST (chicken and turkey [16]), MR-HUM (chicken), and MRE-LFE (turkey [16]). (c) Estimates of tensile anisotropy, ζ, from biaxial testing (chicken and cardiac muscle [37]), and MR-HUM (chicken).

The anisotropic parameter estimates from MR-HUM generally agree with estimates of DST, MRE-LFE, and biaxial testing (Fig. 16). For a viscoelastic tissue like chicken breast, the baseline shear modulus of the material is expected to increase with increased frequency. Riek et al. noted the increase in estimated isotropic shear modulus of bovine muscle ex vivo from μ12kPa at 200 Hz to μ35kPa at 800 Hz using MRE [38]. The baseline shear modulus estimate from MR-HUM at 400 Hz is consistent with the results from DST and MRE-LFE in chicken and turkey breast, when considering the effects of viscoelasticity. The shear anisotropy estimated by MR-HUM is consistent with the estimates from DST in chicken and turkey breast. The estimate of shear anisotropy (ϕ from MR-HUM is within the large standard deviation of the MRE-LFE estimate. The tensile anisotropy (ζ) estimated by MR-HUM is consistent with the estimates from biaxial testing on chicken and cardiac tissue. The large standard deviations and spread of ϕ and ζ estimated from traditional methods (DST and biaxial testing) illustrate the challenges of anisotropic parameter estimation that make “ground truth” parameter values almost impossible to obtain by direct mechanical testing of soft tissue like muscle or brain. Because of this, both (i) verification of the method using data from simulation and (ii) comparison of baseline shear modulus to estimates from other approaches provide important evidence for viability of the anisotropic parameter estimation method.

5 Conclusion

The combination of focused ultrasound, to generate shear waves and MR elastography to image and analyze these waves, provides a noninvasive and nondestructive approach (MR-HUM) to estimate anisotropic material parameters in soft fibrous tissue. Future work will explore the application of this approach to characterize tissues like muscle and white matter of the brain, in vivo.

Acknowledgment

All MRE and DTI experiments were performed in the Small Animal MR Facility, Mallinckrodt Institute of Radiology, Washington University in St. Louis. We are grateful to Ryan Castile for technical assistance with biaxial testing.

Funding Data

  • ONR (DURIP) (Grant No. N00014-17-1-2301; Funder ID: 10.13039/100000006).

  • NSF (Grant No. CMMI-1727412; Funder ID: 10.13039/100000001).

  • NIH (Grant No. R01 EB027577; Funder ID: 10.13039/100000002).

Nomenclature

Mathematical Symbols

    Mathematical Symbols
     
  • a =

    fiber direction

  •  
  • cs,cf =

    wave speed of slow and fast shear waves (respectively)

  •  
  • f =

    frequency

  •  
  • ks,kf =

    wave numbers of slow and fast shear waves (respectively)

  •  
  • ms,mf =

    shear wave polarization for slow and fast shear waves (respectively)

  •  
  • n =

    propagation direction

  •  
  • polthresh =

    polarization threshold

  •  
  • β =

    angle between actuation direction and fiber direction

  •  
  • Γs,Γf =

    curl of slow and fast shear waves (respectively)

  •  
  • ζ =

    tensile anisotropy

  •  
  • θ =

    angle between propagation direction and fiber direction

  •  
  • λs,λf =

    wavelength of slow and fast shear waves (respectively)

  •  
  • μ =

    baseline shear modulus

  •  
  • μapp =

    apparent shear modulus

  •  
  • ρ =

    density

  •  
  • ϕ =

    shear anisotropy

  •  
  • ψs,ψf =

    phase angles of slow and fast shear waves (respectively)

Acronyms

    Acronyms
     
  • DST =

    dynamic shear testing

  •  
  • DTI =

    diffusion tensor imaging

  •  
  • FUS =

    focused ultrasound

  •  
  • HMI =

    harmonic motion imaging

  •  
  • ITI =

    incompressible transversely isotropic

  •  
  • LFE =

    local frequency estimation

  •  
  • mpARFI =

    multiple point acoustic radiation force imaging

  •  
  • MR-ARFI =

    acoustic radiation force imaging

  •  
  • MRE =

    magnetic resonance elastography

  •  
  • MR-HUM =

    magnetic resonance imaging of harmonic ultrasound-induced motion

  •  
  • MRI =

    magnetic resonance imaging

  •  
  • NITI =

    nearly incompressible transversely isotropic

  •  
  • PG =

    phase gradient

  •  
  • t-MRE =

    transient magnetic resonance elastography

  •  
  • TI =

    transversely isotropic

References

1.
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
(
16
), pp.
2761
2768
.10.1016/j.jbiomech.2013.09.008
2.
Royer
,
D.
,
Gennisson
,
J. L.
,
Deffieux
,
T.
, and
Tanter
,
M.
,
2011
, “
On the Elasticity of Transverse Isotropic Soft Tissues
,”
J. Acoust. Soc. Am.
,
129
(
5
), pp.
2757
2760
.10.1121/1.3559681
3.
Gennisson
,
J. L.
,
Catheline
,
S.
,
Chaffai
,
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
(
1
), pp.
536
541
.10.1121/1.1579008
4.
Aristizabal
,
S.
,
Amador
,
C.
,
Qiang
,
B.
,
Kinnick
,
R. R.
,
Nenadic
,
I. Z.
,
Greenleaf
,
J. F.
, and
Urban
,
M. W.
,
2014
, “
Shear Wave Vibrometry Evaluation in Transverse Isotropic Tissue Mimicking Phantoms and Skeletal Muscle
,”
Phys. Med. Biol.
,
59
(
24
), pp.
7735
7752
.10.1088/0031-9155/59/24/7735
5.
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 2-D Matrix Ultrasound Array
,”
IEEE Trans. Med. Imag.
,
32
(
9
), pp.
1671
1684
.10.1109/TMI.2013.2262948
6.
Sinkus
,
R.
,
Tanter
,
M.
,
Catheline
,
S.
,
Lorenzen
,
J.
,
Kuhl
,
C.
,
Sondermann
,
E.
, and
Fink
,
M.
,
2005
, “
Imaging Anisotropic and Viscous Properties of Breast Tissue by Magnetic Resonance-Elastography
,”
Magn. Reson. Med.
,
53
(
2
), pp.
372
387
.10.1002/mrm.20355
7.
Green
,
M. A.
,
Geng
,
G.
,
Qin
,
E.
,
Sinkus
,
R.
,
Gandevia
,
S. C.
, and
Bilston
,
L. E.
,
2013
, “
Measuring Anisotropic Muscle Stiffness Properties Using Elastography
,”
NMR Biomed.
,
26
(
11
), pp.
1387
1394
.10.1002/nbm.2964
8.
Klatt
,
D.
,
Papazoglou
,
S.
,
Braun
,
J.
, and
Sack
,
I.
,
2010
, “
Viscoelasticity-Based MR Elastography of Skeletal Muscle
,”
Phys. Med. Biol.
,
55
(
21
), pp.
6445
6459
.10.1088/0031-9155/55/21/007
9.
Papazoglou
,
S.
,
Rump
,
J.
,
Braun
,
J.
, and
Sack
,
I.
,
2006
, “
Shear Wave Group Velocity Inversion in MR Elastography of Human Skeletal Muscle
,”
Magn. Reson. Med.,
56
(
3
), pp.
489
497
.10.1002/mrm.20993
10.
Qin
,
E. C.
,
Sinkus
,
R.
,
Geng
,
G.
,
Cheng
,
S.
,
Green
,
M.
,
Rae
,
C. D.
, and
Bilston
,
L. E.
,
2013
, “
Combining MR Elastography and Diffusion Tensor Imaging for the Assessment of Anisotropic Mechanical Properties: A Phantom Study
,”
J. Magn. Reson. Imag.
,
37
(
1
), pp.
217
226
.10.1002/jmri.23797
11.
Qin
,
E. C.
,
Jugé
,
L.
,
Lambert
,
S. A.
,
Paradis
,
V.
,
Sinkus
,
R.
, and
Bilston
,
L. E.
,
2014
, “
In Vivo Anisotropic Mechanical Properties of Dystrophic Skeletal Muscles Measured by Anisotropic MR Elastographic Imaging: The Mdx Mouse Model of Muscular Dystrophy
,”
Radiology
,
273
(
3
), pp.
726
735
.10.1148/radiol.14132661
12.
Namani
,
R.
,
Wood
,
M. D.
,
Sakiyama-Elbert
,
S. E.
, and
Bayly
,
P. V.
,
2009
, “
Anisotropic Mechanical Properties of Magnetically Aligned Fibrin Gels Measured by Magnetic Resonance Elastography
,”
J. Biomech.
,
42
(
13
), pp.
2047
2053
.10.1016/j.jbiomech.2009.06.007
13.
Feng
,
Y.
,
Okamoto
,
R. J.
,
Namani
,
R.
,
Genin
,
G. M.
, and
Bayly
,
P. V.
,
2013
, “
Measurements of Mechanical Anisotropy in Brain Tissue and Implications for Transversely Isotropic Material Models of White Matter
,”
J. Mech. Behav. Biomed. Mater.
,
23
, pp.
117
132
.10.1016/j.jmbbm.2013.04.007
14.
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
15.
Romano
,
A.
,
Scheel
,
M.
,
Hirsch
,
S.
,
Braun
,
J.
, and
Sack
,
I.
,
2012
, “
In Vivo Waveguide Elastography of White Matter Tracts in the Human Brain
,”
Magn. Reson. Med.
,
68
(
5
), pp.
1410
1422
.10.1002/mrm.24141
16.
Schmidt
,
J. L.
,
Tweten
,
D. J.
,
Benegal
,
A. N.
,
Walker
,
C. H.
,
Portnoi
,
T. E.
,
Okamoto
,
R. J.
,
Garbow
,
J. R.
, and
Bayly
,
P. V.
,
2016
, “
Magnetic Resonance Elastography of Slow and Fast Shear Waves Illuminates Differences in Shear and Tensile Moduli in Anisotropic Tissue
,”
J. Biomech.
,
49
(
7
), pp.
1042
1049
.10.1016/j.jbiomech.2016.02.018
17.
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
18.
Chatelin
,
S.
,
Deck
,
C.
,
Renard
,
F.
,
Kremer
,
S.
,
Heinrich
,
C.
,
Armspach
,
J. P.
, and
Willinger
,
R.
,
2011
, “
Computation of Axonal Elongation in Head Trauma Finite Element Simulation
,”
J. Mech. Behav. Biomed. Mater.
,
4
(
8
), pp.
1905
1919
.10.1016/j.jmbbm.2011.06.007
19.
Anderson
,
A. T.
,
Van Houten
,
E. E.
,
McGarry
,
M. D.
,
Paulsen
,
K. D.
,
Holtrop
,
J. L.
,
Sutton
,
B. P.
,
Georgiadis
,
J. G.
, and
Johnson
,
C. L.
,
2016
, “
Observation of Direction-Dependent Mechanical Properties in the Human Brain With Multi-Excitation MR Elastography
,”
J. Mech. Behav. Biomed. Mater.
,
59
, pp.
538
546
.10.1016/j.jmbbm.2016.03.005
20.
Dargar
,
S.
,
Rahul
,
Kruger
,
U.
, and
De
,
S.
,
2019
, “
In Vivo Layer-Specific Mechanical Characterization of Porcine Stomach Tissue Using a Customized Ultrasound Elastography System
,”
ASME J. Biomech. Eng.
,
10
(
10
), p.
101004
.10.1115/1.4043259
21.
Konofagou
,
E. E.
,
Maleke
,
C.
, and
Vappou
,
J.
,
2012
, “
Harmonic Motion Imaging (HMI) for Tumor Imaging and Treatment Monitoring
,”
Curr. Med. Imag. Rev.
,
8
(
1
), pp.
16
26
.10.2174/157340512799220616
22.
Liu
,
Y.
,
Liu
,
J.
,
Fite
,
B. Z.
,
Foiret
,
J.
,
Ilovitsh
,
A.
,
Leach
,
J. K.
,
Dumont
,
E.
,
Caskey
,
C. F.
, and
Ferrara
,
K. W.
,
2017
, “
Supersonic Transient Magnetic Resonance Elastography for Quantitative Assessment of Tissue Elasticity
,”
Phys. Med. Biol.
,
62
(
10
), pp.
4083
4106
.10.1088/1361-6560/aa6674
23.
Doherty
,
J. R.
,
Trahey
,
G. E.
,
Nightingale
,
K. R.
, and
Palmeri
,
M. L.
,
2013
, “
Acoustic Radiation Force Elasticity Imaging in Diagnostic Ultrasound
,”
IEEE Trans. Ultrason., Ferroelectr., Freq. Control
,
60
(
4
), pp.
685
701
.10.1109/TUFFC.2013.2617
24.
Sinkus
,
R.
,
Tanter
,
M.
,
Bercoff
,
J.
,
Siegmann
,
K.
,
Pernot
,
M.
,
Athanasiou
,
A.
, and
Fink
,
M.
,
2008
, “
Potential of MRI and Ultrasound Radiation Force in Elastography: Applications to Diagnosis and Therapy
,”
Proc. IEEE
,
96
(
3
), pp.
490
499
.10.1109/JPROC.2007.913536
25.
Wu
,
T.
,
Felmlee
,
J. P.
,
Greenleaf
,
J. F.
,
Riederer
,
S. J.
, and
Ehman
,
R. L.
,
2000
, “
MR Imaging of Shear Waves Generated by Focused Ultrasound
,”
Magn. Reson. Med.
,
43
(
1
), pp.
111
115
.10.1002/(SICI)1522-2594(200001)43:1<111::AID-MRM13>3.0.CO;2-D
26.
McDannold
,
N.
, and
Maier
,
S. E.
,
2008
, “
Magnetic Resonance Acoustic Radiation Force Imaging
,”
Med. Phys.
,
35
(
8
), pp.
3748
3758
.10.1118/1.2956712
27.
Liu
,
Y.
,
Fite
,
B. Z.
,
Mahakian
,
L. M.
,
Johnson
,
S. M.
,
Larrat
,
B.
,
Dumont
,
E.
, and
Ferrara
,
K. W.
,
2015
, “
Concurrent Visualization of Acoustic Radiation Force Displacement and Shear Wave Propagation With 7T MRI
,”
PLoS One
,
10
(
10
), p.
e0139667
.10.1371/journal.pone.0139667
28.
Chen
,
H.
,
Hou
,
G. Y.
,
Han
,
Y.
,
Payen
,
T.
,
Palermo
,
C. F.
,
Olive
,
K. P.
, and
Konofagou
,
E. E.
,
2015
, “
Harmonic Motion Imaging for Abdominal Tumor Detection and High-Intensity Focused Ultrasound Ablation Monitoring: An In Vivo Feasibility Study in a Transgenic Mouse Model of Pancreatic Cancer
,”
IEEE Trans. Ultrason. Ferroelectr. Freq. Control
,
62
(
9
), pp.
1662
1673
.10.1109/TUFFC.2015.007113
29.
Odeen
,
H.
,
de Bever
,
J.
,
Hofstetter
,
L. W.
, and
Parker
,
D. L.
,
2019
, “
Multiple-Point Magnetic Resonance Acoustic Radiation Force Imaging
,”
Magn. Reson. Med.
,
81
(
2
), pp.
1104
1117
.10.1002/mrm.27477
30.
Okamoto
,
R. J.
,
Clayton
,
E. H.
, and
Bayly
,
P. V.
,
2011
, “
Viscoelastic Properties of Soft Gels: Comparison of Magnetic Resonance Elastography and Dynamic Shear Testing in the Shear Wave Regime
,”
Phys. Med. Biol.
,
56
(
19
), pp.
6379
6400
.10.1088/0031-9155/56/19/014
31.
Clayton
,
E. H.
, Garbow, J. R., and
Bayly
,
P. V.
,
2011
, “
Frequency-Dependent Viscoelastic Parameters of Mouse Brain Tissue Estimated by MR Elastography
,”
Phys. Med. Biol.,
56
(
8
), pp. 2391–2406.10.1088/0031-9155/56/8/005
32.
Shimony
,
J. S.
,
McKinstry
,
R. C.
,
Akbudak
,
E.
,
Aronovitz
,
J. A.
,
Snyder
,
A. Z.
,
Lori
,
N. F.
,
Cull
,
T. S.
, and
Conturo
,
T. E.
,
1999
, “
Quantitative Diffusion-Tensor Anisotropy Brain MR Imaging: Normative Human Data and Anatomic Analysis
,”
Radiology
,
212
(
3
), pp.
770
784
.10.1148/radiology.212.3.r99au51770
33.
Mori
,
S.
, and
van Zijl
,
P. C.
,
2002
, “
Fiber Tracking: Principles and Strategies—A Technical Review
,”
NMR Biomed.
,
15
(
7–8
), pp.
468
480
.10.1002/nbm.781
34.
Deeken
,
C. R.
,
Thompson
,
D. M.
,
Castile
,
R. M.
, and
Lake
,
S. P.
,
2014
, “
Biaxial Analysis of Synthetic Scaffolds for Hernia Repair Demonstrates Variability in Mechanical Anisotropy, Non-Linearity and Hysteresis
,”
J. Mech. Behav. Biomed. Mater.
,
38
, pp.
6
16
.10.1016/j.jmbbm.2014.06.001
35.
Okamoto
,
R. J.
,
Romano
,
A. J.
,
Johnson
,
C. L.
, and
Bayly
,
P. V.
,
2019
, “
Insights Into Traumatic Brain Injury From MRI of Harmonic Brain Motion
,”
J. Exp. Neurosci.
,
13
.10.1177/1179069519840444
36.
Manduca
,
A.
,
Muthupillai
,
R.
,
Rossman
,
P. J.
,
Greenleaf
,
J. F.
, and
Ehman
,
R. L.
,
1996
, “
Image Processing for Magnetic-Resonance Elastography
,”
SPIE
Paper No. 2710, pp. 616–623.
10.1117/12.237965
37.
Humphrey
,
J. D.
,
Strumpf
,
R. K.
, and
Yin
,
F. C. P.
,
1990
, “
Determination of a Constitutive Relation for Passive Myocardium—II: Parameter Estimation
,”
ASME J. Biomech. Eng.
,
112
(
3
), pp.
340
346
.10.1115/1.2891194
38.
Riek
,
K.
,
Klatt
,
D.
,
Nuzha
,
H.
,
Mueller
,
S.
,
Neumann
,
U.
,
Sack
,
I.
, and
Braun
,
J.
,
2011
, “
Wide-Range Dynamic Magnetic Resonance Elastography
,”
J. Biomech.
,
44
(
7
), pp.
1380
1386
.10.1016/j.jbiomech.2010.12.031