Real‐time MR elastography for viscoelasticity quantification in skeletal muscle during dynamic exercises

To develop and test real‐time MR elastography for viscoelastic parameter quantification in skeletal muscle during dynamic exercises.


| INTRODUCTION
Muscle activity involves the generation of mechanical forces, which deform or strain the muscle. This in turn results in mechanical stress, or force per area, acting on tendons, blood vessels, or surrounding tissue. According to the basic laws of elasticity, a measured strain is related to stress by the underlying tissue stiffness: the higher the elastic modulus at a given strain, the higher the stress that is generated. 1 For that reason, elastography, which depicts the elastic modulus of soft tissues in medical imaging, allows quantification of muscle activity in clinical settings. 2,3 For example, elastography in rehabilitation medicine can help in the design of assistive technologies because it noninvasively provides information on muscle status that is directly linked to muscle structure and function. 4 Ultrasound-based elastography (USE) techniques have been used extensively to investigate skeletal muscles and to quantify changes in muscle stiffness induced by activity. [5][6][7][8][9][10] USE is an excellent tool for studying muscle action because it allows repeated evaluation with high frame rates from 25 to 40 Hz. 7,11 However, USE is limited in the measurement of different components of mechanical strain, as required to account for anisotropic elastic tensor components such as transverse isotropy in muscle tissue. 12 Furthermore, most USE techniques have a limited FOV and cannot visualize several muscles at the same time. 13 By contrast, MR elastography (MRE) can quantify viscoelastic tissue parameters in several organs across extended FOVs. 14 MRE can also encode any strain components, which is especially favorable for studying anisotropic elastic coefficients in muscle tissue. [15][16][17][18][19][20] However, MRE is limited by long acquisition times of the order of several minutes, during which organs move or muscle load changes. 21 Therefore, so far MRE has only been used to study muscle mechanical parameters either in a relaxed state, 19,[22][23][24] under sustained loading with the assumption that submaximal contraction was maintained at a constant level, 2,15,[25][26][27] or-as in the case of myocardial tissue-by assuming periodic changes in elastic parameters to which MRE data acquisition can be repetitively synchronized. [28][29][30] Direct measurement of dynamic viscoelasticity changes without the need for synchronization would enhance the usability of MRE as a tool for quantifying muscle function. However, in classical MRE, externally induced harmonic motion is repeatedly encoded for acquiring phases, components, or frequencies of the wave fields. Many MRE sequences require precise resynchronization to the external harmonic motion by adjusting the vibration phase or the start of motion encoding. Furthermore, 2-or 3D k-spaces are normally repeatedly acquired either as single lines or as full image slices, again imposing the need for resynchronization of sequence and motion. 21 We here introduce an MRE method that permits the monitoring of dynamic viscoelasticity parameters in 2 dimensions without resynchronization of sequence and motion. To this end, the method captures 2 components of harmonic motion fields in an interleaved fashion by utilizing the principle of stroboscopic undersampling of harmonic vibrations (SUV). SUV was recently proposed for MRE of the brain as a method of sampling relatively fast tissue oscillations by a sequence of relatively slow encoding steps. 31 Sampling of harmonic vibrations with a frequency larger than the frame rate of the sequence causes aliasing. However, in SUV the aliased frequency is well defined because both sequence timing and the harmonic vibration period are precisely known. In our earlier work on SUV in the brain, we performed k-space acquisition across several shots in conjunction with cardiac triggering, which required additional resynchronization of the encoding sequence and external motion. 31 In the present work, we extend the concept of SUV to ungated single-shot MRE, allowing us to capture nonperiodic viscoelasticity changes during dynamic voluntary activation of skeletal muscles. Our hypothesis is that we can measure the changes in viscoelastic parameters induced in different groups of skeletal muscles of the lower extremity upon plantar flexion and dorsiflexion. The gastrocnemius (GAST), soleus (SOL), peroneus (PER), and tibialis posterior (TIBP) belong to the posterior compartment of the leg, which acts as plantar flexor. In contrast, the tibialis anterior (TIBA) and extensor digitorum longus belong to the anterior compartment of the lower leg, which acts as dorsiflexor. However, activation in the lower leg normally invokes a more complex pattern of muscle activity for the following reasons: 1) the muscles can have more than 1 function (e.g., extensor digitorum longus is a dorsiflexor that also stretches the toes 32 ); 2) the function of the muscles depends on leg posture (e.g., GAST function is limited when the knee is in a bent position 33 ); and 3) muscle activation always reflects the vector sum of counteracting forces, which varies with extension of the muscle in time and space. Consequently, muscle activity patterns vary greatly from subject to subject, over time, and with force and anatomy. The concert of lower leg muscle activation including GAST, SOL, PER, TIBA, and TIBP has never been addressed in a synchronous manner by elastography in real time.
The method described here provides approximately 6 frames per second of shear-wave speed (SWS) maps and phase-angle maps of the complex shear modulus (φ). This high frame rate of viscoelasticity imaging is considered to represent real-time relative to the duration of the loading cycle of skeletal muscles, which is of the order of seconds. Hereinafter, we refer to our method of single-shot stroboscopic sampling MRE as real-time MRE (rt-MRE). The feasibility of the method is demonstrated in lower-extremity muscles of healthy participants performing voluntary exercises of isometric plantar flexion and dorsiflexion.

| Subjects
The study was approved by the ethics committee of the Charité-Universitätsmedizin Berlin in accordance with the World Medical Association Declaration of Helsinki (Ethical Principles for Medical Research Involving Human Subjects), and each volunteer gave written informed consent to participate. Right lower-leg muscles of 15 healthy participants without orthopedic or neuromuscular abnormalities were investigated (2 females, 13 males; mean age 35 ± 10 years, age range 24 to 58 years). Because muscle structure and function can vary significantly between males and females, our study focused on male volunteers. However, we also included 2 female volunteers to explore the feasibility of the method in women. Because no apparent differences in effect size or in measured values were observed between the sexes, we analyzed pooled data.

| MRI scanner and mechanical actuator setup
All measurements were performed on a 1.5 Tesla clinical MRI scanner (Siemens Magnetom Sonata, Erlangen, Germany). Continuous harmonic vibrations were induced by a pressurized air driver comprising a compressed air hose and a 3D-printed flat, bottle-shaped passive actuator placed next to a standard extremity coil (single-element birdcage coil), directly underneath the subject's Achilles tendon ( Figure 1). A forerun of vibrations before the start of the imaging sequence of approximately 3 s was performed in order to establish a mechanical steady state. A vibration frequency of 40 Hz was applied continuously throughout each scan of approximately 30 s.

| Paradigms of muscle activation
To quantify activation of lower-leg muscles by rt-MRE, all participants were asked to perform 2 exercises voluntarily during consecutive scans: 1) isometric plantar flexion and 2) isometric dorsiflexion. Therefore, at the start of MRE data acquisition, the lower-leg muscles were rested for 10 s (baseline phase), followed by a phase of muscle activation for 10 s (exercise) and a second phase of rest for the remaining scan time (recovery). Each participant was given commands over headphones for starting and ending voluntary muscle activation. To allow for isometric muscle activation, a custom-made footplate was mounted to the RF coil and adjusted such that a foot position perpendicular to the lower extremity (90°) was maintained during rest and exercise, as illustrated in Figure 1. The knee was slightly flexed with an angle of approximately 145°, as indicated in Figure 1, to maintain a fully rested position of lower-extremity muscles without any prestretch. The foot and ankle were fixed to the loading plate by hook-and-loop tape to minimize ankle rotations in the plantar and dorsal direction while ensuring sustained isometric conditions and avoiding major displacement of the extremity. Effectively, calf-muscle activation was initiated by the voluntary force imposed onto the loading plate through the ball of the foot without changing position. The loading device is shown in Figure 1.

| rt-MRE pulse sequence
A 2D single-shot, steady-state gradient echo MRE pulse sequence with spiral readout (spiral-out trajectories starting from the center of k-space) was developed. TR was 85 ms, including RF excitation (25° flip angle), motion encoding within 16 ms TE, and signal spoiling. According to the principle of fractional encoding, the bipolar, flow-compensated motion-encoding gradient (0th and 1st moment nulling, F I G U R E 1 Experimental setup with a custom-built footplate mounted on an extremity coil to apply voluntary isometric tension to the lower extremity muscles. A single pressurized-air driver was placed underneath the Achilles tendon (white arrow) 12 ms duration, 35 mT/m amplitude) was shorter than the vibration period of external mechanical stimulation. 34 Two in-plane motion components were continuously encoded in an interleaved fashion, resulting in an overall MRE frame rate of 2 × TR = 170 ms, or 5.9 Hz. Steady-state continuous vibrations were consecutively acquired over 360 TRs, resulting in a total acquisition time of 30.6 s. The first of the 180 phases of both in-plane motion components was acquired without a motion-encoding gradient and treated as a static offset phase (θ 0 ) in order to perform voxel-wise phase corrections in all subsequent images by rotating their MRI signal phases by θ 0 . This phase correction eliminates any static phase offsets and thus reduces the number of phase wraps in the image, which helps to unfold phase wraps resulting from elastic deformation. Data were acquired in a transverse view through the center of the lower leg. Slice thickness, FOV, and matrix size were 5 mm, 160 × 160 mm 2 , and 128 × 128, respectively. A sequence timing diagram including the timing of dynamic muscle activation is provided in Figure 2.

| Parameter reconstruction
In a first step, image registration was performed to account for minor lower-leg displacements during exercise. This was done using the open-source elastix toolbox 35 to register the 360 signal magnitudes of the rt-MRE images acquired over time to the first image of the series. Rigid in-plane transformations were then applied to the real parts and imaginary parts of all 360 rt-MRE images. The phases of these motioncorrected rt-MRE data were further processed for parameter recovery.
To this end, all phase images were Hilbert-transformed, yielding 180 complex-valued wave images for each of the 2 wave-field components. Hilbert transformation was performed using a Gaussian bandpass of σ = 0.3 Hz centered at fundamental vibration frequency, as proposed in Ref. 31. Parameter maps were reconstructed 180 times and were based on 2 consecutive complex-valued wave images depicting the 2 Hilbert-transformed in-plane wave field components. We thus obtained 180 parameter maps generated with a 5.9-Hz frame rate. Representative shear-wave images are given in Figure 3.
Two parameters were reconstructed for each of the 180 time points: 1) SWS (in m/s) as a surrogate parameter for stiffness (SWS in m/s); and 2) phase angle of the complex shear modulus (φ), which measures a material's fluidity between the 2 limits of 0 and π/2 (pure solids have φ = 0; pure fluids have φ = π/2). 36 Hereinafter, φ is referred to as tissue fluidity.
A dual reconstruction pipeline was employed: SWS was reconstructed by wave number-based (k-) multicomponent dual elasto-visco (MDEV) inversion, 37 whereas φ was retrieved from MDEV as detailed previously. 38 Both algorithms are publicly available on https ://bioqic-apps.chari te.de

| Parameter analysis and statistical tests
Time-dependent SWS and φ parameters were further analyzed for the 6 muscle groups investigated in the lower leg: GAST, SOL, PER, extensor digitorum longus, TIBA, and TIBP. Regions of interest were selected based on an additionally acquired higher-resolution T 1 -weighted image, which shared the slice geometry and orientation with rt-MRE. A representative T 1 -weighted image with overlaid regions of interest is shown in Figure 4. SWS and φ were averaged within the respective regions of interest and analyzed for temporal variations and group statistics. For calculation of summarized group statistics, time-course data were temporally averaged across each of the 3 phases (baseline, exercise, recovery), omitting the first and last 2.5 s of each phase to minimize transient effects.
To test for significant changes in SWS and φ during exercise, a linear mixed-effect model was applied with SWS and φ as dependent variables and the 3 phases as independent variables. Participants were treated as a random effect, accounting for non-independency due to repeated measurements in the same participants and muscle conditions across the 3 phases. P values were calculated from multiple comparisons using Tukey method. Statistical significance was assumed for P values below 0.05. All statistical analysis was done in R (version 3.5.1). SNR analysis was performed for all rt-MRE magnitude images and for all wave-displacement images over the entire acquisition time, as detailed in Ref. 39. Unless otherwise specified, errors are given as group SD.

| RESULTS
SNR of rt-MRE signal magnitudes averaged across all muscle tissue visible in the image slice changed only slightly (approximately 3%; baseline 12.5 ± 2.0 dB, exercise 12.9 ± 2.0 dB, recovery 12.6 ± 2.0 dB; P < .001). SNR of wave displacement did not change during dorsiflexion, whereas a slight change (of approximately 2%) was observed upon plantar flexion (baseline 54.3 ± 2.5 dB, exercise 55.5 ± 2.5 dB, recovery 54.4 ± 4 dB; P < .01), probably due to minor changes in Achilles tendon position relative to the actuator. Average shear-wave amplitudes were 43.2 ± 14.7 µm. Figure 4 presents representative SWS and φ maps from all 3 phases of muscle status. During plantar flexion, the SWS maps demonstrate a major increase in SWS values in the SOL muscle. In contrast, during dorsiflexion, a major increase is seen in the area of the TIBA muscle. Whereas muscle activity is well perceivable from SWS maps, corresponding φ maps do not show clearly visible activation patterns. However, group mean data reveal a significant increase in φ during exercise in several muscles. Figure 5 shows group mean time courses for SWS and φ. SWS increased in the SOL by 20.0% ± 3.6% (P < .0001) and GAST by 9.0% ± 6.0% (P = .03) due to plantar flexion; whereas φ increased predominantly in the GAST by 13.5% ± 6.1% (P < .01), PER by 14.3% ± 10.4% (P < .01), SOL by 41.3% ± 12.0% (P < .0001), and TIBP by 12.3% ± 4.5% (P = .02). During dorsiflexion, we observe an increase in SWS in the extensor of 14.2% ± 5.2% (P < .01), GAST of 15.4% ± 10.2% (P < .0001), SOL of 6.3% ± 2.4% (P < .01), and TIBA of 41.5% ± 10.2% (P < .0001); whereas φ increased primarily in the extensor by 27.7% ± 0.7% (P < .01), GAST by 15.3% ± 5.6% (P < .0001), SOL by 17.2% ± 6.4% (P < .0001), TIBA by 27.9% ± 2.8% (P < .0001), and TIBP by 21.8% ± 2.8% (P < .0001). No significant difference was found between baseline and recovery for SWS and φ in any of the investigated muscles. Figure 6 plots group mean values of SWS versus φ at baseline and during exercise, demonstrating the viscoelastic property changes of lower extremity muscles due to plantar flexion and dorsiflexion. The black arrow indicates the direction of viscoelastic change from baseline to exercise. Of the 6 muscle groups investigated, 2 were significantly activated by plantar flexion and 4 by dorsiflexion based on SWS, whereas φ changed significantly in 5 muscles during the 2 exercises.
All group mean SWS and φ values and effect sizes are summarized in Table 1.

| DISCUSSION
This paper presents rt-MRE, a method to study dynamic viscoelasticity changes in soft tissues without synchronization to MRE data acquisition. rt-MRE was found to be able to quantify the activation of posterior muscles in the lower extremity during plantar flexion and the activation of anterior muscles during dorsiflexion. Consistent with our hypothesis, SOL and TIBA (the largest flexor muscles) in a bent knee position-demonstrated the strongest effects of antagonistic parameter changes in rt-MRE. rt-MRE captured the in-plane viscoelastic shear modulus using the principle of SUV, which yielded relatively high frame rates of approximately 6 Hz. Higher frame rates are feasible, for example, by encoding only 1 motion component-similar to rapid ultrasound-based elastographywhich is predominantly sensitive to the displacement component along the axis of the ultrasound beam. 12 Capturing 2 orthogonal in-plane wave components in the transverse view relative to the principal axis of the muscle (e.g., x 3 ) allows shear-modulus reconstruction of μ 12 , which is related to the plane of isotropy in transversely isotropic materials. 40 In our previous work, we exploited the directional filter capabilities of the curl operator to retrieve 3 independent elasticity parameters and 3 independent viscosity parameters according to the constitutive laws of incompressible transversely isotropic materials. 19  ). 20 The transversely isotropic viscoelastic coefficients related to the plane of symmetry and x 3 -axis have not been addressed in our present study. Instead, we intentionally chose our image plane and directions of motion encoding such that predominantly slow-transverse (ST) wave modes were captured, whereas contributions of fast-transverse wave modes (FT) were minimized or fully canceled. ST modes relate to waves in which both polarization and propagation are directed transversely to the x 3 -axis, whereas FT modes relate to waves in which either polarization or propagation is directed along x 3 . 26 Thus, waves with x 1 -or x 2 -polarization propagating within the plane of isotropy (x 1 to x 2 ) belong to ST modes and were acquired by our experimental setup.
Unlike in our previous work, 31 we captured shear waves only at a single-drive frequency. Therefore, multifrequency inversion as proposed by Guo et al could not be applied. 19 To enhance the stability of the wave inversion against noise, we employed k-MDEV inversion, which has been shown to minimize the effect of noise on SWS parameters by using noise-robust, single-order gradients instead of second-order gradients as required for anisotropic 3-parameter inversion. 37 k-MDEV comprises directional bandpass filtering in 2 dimensions and thus analyzes only those waves that propagate in-plane, which in our scenario are ST waves. To validate the effectiveness of directional bandpass filtering in k-MDEV against the established curl-based method, we preprocessed our wave images by applying the curl operator to in-plane wave components followed by standard k-MDEV processing. The resulting SWS values and effect sizes were very similar to those obtained by our standard k-MDEV pipeline, suggesting that longitudinal waves were effectively suppressed. Specifically, using curl preprocessing, we found that SOL and TIBA had baseline SWS values of 1.03 ± 0.09 m/s and 1.07 ± 0.1 m/s during plantar flexion and 1.01 ± 0.09 m/s and 1.05 ± 0.1 m/s during dorsiflexion (instead of 0.97 ± 0.09 m/s, 1.06 ± 0.15 m/s, 0.96 ± 0.08 m/s, and 1.03 ± 0.13 m/s, respectively; all P > 0.05). In addition, the effect sizes obtained by curl-based preprocessing were statistically not different from those obtained by standard k-MDEV processing (22.26% ± 1.03%, 6.42% ± 7.52% vs. 17.16% ± 4.82%,   18 respectively. The differences might be related to technical differences such as 3D smoothing before inversion in Ref. 19, possibly mixing FT components into ST waves-or shorter wavelengths due to noise as encountered by Green et al. 18 Other reasons might relate to preconditioning. For example, in Guo et al., the lower extremity was stretched by 180° and lower leg muscles were rested in the RF coil. 19 In our current setup, the muscles investigated rested without any prestrainthat is, without being supported by the table or coil-and in an angled knee position (145°), which may have contributed to softer values.
Few data have been reported on changes in μ * 12 induced by muscle activation. The relative increase in stiffnessof the order of 20% and 40% for the SOL and TIBA, respectively-is in a range similar to that found for μ * 13 by Klatt et al 26 in femoral muscles (approximately 40%) or by Barnhill et al 27 in the quadriceps muscle (approximately 30%). In F I G U R E 7 Hypothesized relationship between polar interactions within muscle microstructures and MRE-measured fluidity. A sketch of the general morphology of a skeletal muscle cell (fiber) with and without contraction is shown on the left. A molecular illustration of muscle contraction is shown on the right. The backbone of a muscle fiber is composed of actin and myosin filaments, where actin interacts with troponin and tropomyosin. Actin-myosin interactions are blocked when the muscle rests due to diminished calcium ion concentration and the fact that ATP is bound to myosin. In this state, polar interactions of proteins with water (depicted as bound water) lubricate the sliding filaments, thus reducing friction. Muscle contraction occurs upon a flux of calcium ions that bind to troponin, causing conformational changes and weakening troponin interactions with tropomyosin so that "cross-bridging" between actin and myosin can take place. This process alters the proteins' hydrophilic nature and their polar interactions with ions and water. The sliding motion of interacting filaments, which shortens fiber length, pushes away bound and free water, thus reducing the lubrication zone and increasing molecular friction. Of note, the lubrication zone illustrates the concept of reduced friction due to inhibition of hydrophobic protein-protein interactions and is not locally restricted to the shown actin molecules. Our present results for φ as well as MRI and MRE data reported in Refs. 26, 41, and 42 suggest that increased mechanical friction in skeletal muscle upon contraction is associated with higher fluidity due to the hypothesized molecular mechanisms just outlined. ATP, adenosine triphosphate contrast, larger changes were observed for μ * 13 in the biceps, 15 with slopes of stiffness versus load of up to 34.0 ± 3.5 kPa/kg. USE in the same muscle demonstrated a strong increase in SWS (2.7 to 5.7 m/s) when measured along the fibers, in contrast to a shallow increase (1.2 to 1.8 m/s) perpendicular to the fibers. 7 Furthermore, brachialis muscle tissue was nondispersive along the fibers as opposed to strongly dispersive perpendicular to the fibers. 7 Strong dispersion of perpendicular SWS agrees with our findings of relatively high φ values because φ specifies the slope of the dispersion curve, which is also known as loss angle. φ of the order of 0.65 and 0.75 for SOL and TIBA muscles at rest indicates solid yet highly dispersive tissue properties. Activation results in φ values larger than π/4 for both muscles, suggesting that fluidity of muscle tissue significantly increases upon activation toward predominantly fluid tissue properties. This increase in fluidity agrees with previous reports of increasing viscoelastic power-law coefficients α (α = 2φ/π, assuming that the springpot model is valid 21 ) in skeletal muscles during activation. 26,41 A hypothetical explanation for the high fluidity φ of activated muscle is provided in Figure 7, which illustrates the release of bound water from polar myosin-actin groups in sliding filaments. This possible interpretation is supported by an increase in proton T 2 relaxation time of contracted muscle, suggesting a release of water molecules during transition from the relaxed to the contracted state of muscle tissue. 42 Our study has some limitations. We only investigated lowerextremity muscles. These are known to display fiber pennation angles that are misaligned to the global transversely isotropic coordinate system. Mean pennation angles of approximately 10° are reported for the SOL and TIBA muscles, resulting in an uncertainty of ST versus FT-SWS values of about 1.5%. 43 Higher deviations can occur, for example, in the GAST, where the mean fiber pennation angle is 32°, giving rise to approximately 15% FT contribution to ST waves. Taking additional scans such as diffusion tensor imaging would allow us to avoid assumptions about the local coordinate system [44][45][46] ; however, diffusion tensor imaging cannot be combined with MRE in real time. A further limitation is that we did not quantify muscle activity by electromyography or by measuring the load applied to the fixation plate. Such measurement is technically challenging because the leg and foot positions need to be fixed to preclude displacement during acquisition. Intersubject variation in physique, habituation, and muscle training results in different motion patterns, which act on the loading plate and lead to differences in the stimulation of muscle groups resulting from voluntary motions of plantar flexion or dorsiflexion. Otherwise, this individual variability is exactly what rt-MRE can account for if established as an imaging marker of dynamic muscle activation, which would make the method a promising tool for clinical applications. 47,48 5 | CONCLUSION rt-MRE was introduced for real-time mapping of dynamic muscle activation based on SWS and φ with high frame rates of approximately 6 Hz. SWS and φ, which reflect tissue stiffness and fluidity related to the perpendicular shear modulus μ *

12
, were shown to increase in muscles of the lower extremity during voluntary plantar flexion and dorsiflexion. Specifically, the SOL muscle showed increases in SWS and φ values of around 20% and 40%, respectively, upon plantar flexion. Similarly, the TIBA muscle was activated by dorsiflexion, as reflected by a 40% increase in SWS and 30% increase in φ. Allowing assessment of function-related viscoelasticity changes in soft tissues, rt-MRE is potentially useful for the study of physiological processes and the diagnosis of diseases.