Estimation of Anisotropic Material Properties of Soft Tissue by MRI of Ultrasound-Induced Shear Waves

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 [2–5], and magnetic resonance elastography (MRE) [6–15]. 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 [7–11], 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].

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

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 ζ$⁠.

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.

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

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 [24–26]. 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 [24–26] 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.

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

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.

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.

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 mm³ or 24 × 24 × 12 mm³. 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.

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 mm³. Fractional anisotropy (FA) was estimated from diffusion tensor eigenvalues, and fiber direction (a) was estimated from the first principal eigenvector [33].

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 $β=0 deg, 30 deg, 45 deg, 60 deg, or 90 deg$ 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 mm³ 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.

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 $μ=2000 Pa, ϕ=1, and ζ = 2,$ and a body force density at the focal region of 50 kN/m³ at 300 Hz. A second simulation set incorporated stiffer, approximately muscle-like, shear modulus, with parameters of $μ=7500 Pa, ϕ=1, and ζ = 1,$ and a force density of 150 kN/m³ 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.

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.

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

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.

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

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.

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.

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.

All three components (⁠$U, V, W$⁠) of the shear-wave displacement field are shown in panels (a–c) 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 d–g) 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.

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.

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 $β=45 deg$ and $β=90 deg$ 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 $β=90 deg$ 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 (g–j) and (q–t). The majority of the voxels included in the $β=90 deg$ case (all those in the slice depicted here) were classified as slow.

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

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.

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 $|β−90 deg|<$ 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 (g–j) and (q–t) 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 $β= 51 deg$ case and most voxels were classified as slow for the $β=87 deg$ case.

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

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.

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.71 kPa$ (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$⁠.

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 $μ ∼ 33 kPa, ϕ ∼ 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 $μ ∼ 4 kPa and ϕ ∼ 0.6$ at 20–40 Hz [16]. DST of chicken breast samples (n = 7) in this study provided estimates of $μ= 6.19±1.71 kPa and ϕ= 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$⁠.

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 $μ ∼ 12 kPa$ at 200 Hz to $μ ∼ 35 kPa$ 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.

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.

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.

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

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

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