Double pulsed field gradient diffusion MRI to assess skeletal muscle microstructure

Preliminary study to determine whether double pulsed field gradient (PFG) diffusion MRI is sensitive to key features of muscle microstructure related to function.


INTRODUCTION
Skeletal muscle function is related to muscle microstructure. For example, the maximum isometric force-generating capacity of muscle fibers is directly related to muscle fiber cross-sectional area. 1 Changes in muscle microstructure in response to, for example, injury, training, pathology, and age, can alter muscle function and performance. Muscle microstructure is typically assessed with histology, which is highly invasive, destructive to the tissue, and only samples a small volume of a muscle. This has driven interest in developing noninvasive methods to accurately assess changes in muscle microstructure associated with muscle function. Diffusion MRI (dMRI) is a noninvasive approach that quantifies the restricted diffusion of water in tissues with anisotropic microstructure such as skeletal muscle. 2 The sarcolemma is considered to be the primary barrier to diffusion; thus, water is more restricted radially across a muscle fiber (cell) than along its longitudinal axis. The most widely used dMRI pulse sequences used to assess muscle microstructure are single pulsed field gradient (sPFG) sequences, in which dephasing of the diffusion signal along a defined diffusion direction in a defined time period is measured and fit to a third-order tensor model called DTI. 3,4 Of these sequences, spin-echo DTI is the most common, in which diffusion-encoding gradients are applied before and after a 180 • refocusing pulse. However, spin-echo DTI pulse sequences suffer from lower diffusion sensitivity in muscle because of the relatively short diffusion times (Δ) that need to be used in the presence of the short T 2 of muscle. 5,6 Although expected differences in the diffusion tensor (DT) are often observed (i.e., more restriction with muscle atrophy [7][8][9] ), the sensitivity of these measurements to detect small changes in muscle microstructure is generally not known and not validated. Stimulated echo DTI pulse sequences are less common but allow for prolonged Δ to be sampled. However, stimulated echo DTI sequences result in longer scan times and decreased SNR, thus resulting in less accurate diffusion measurements than spin-echo DTI, given the same time constraints.

F I G U R E 1
Pulse sequence diagram of a double pulsed field gradient sequence. , gradient duration; , gradient spacing; G, gradient strength; t m , mixing time.
There are several challenges associated with dPFG pulse sequences that have limited its application and use: (1) There is no standard analytical model-analogous to the DT in sPFG experiments-that can be used to quantify features of microstructure; (2) the pulse sequence is not widely available; and (3) they have inherently longer TEs potentially limiting the SNR. In fact, until recently it has primarily been studied with numerical simulations. However, we have recently developed an analytical approach to quantify restricted diffusion measured with dPFG dMRI pulse sequences called diffusion tensor subspace imaging (DiTSI) and demonstrated how this technique can enhance sensitivity to measuring microstructural anisotropy in neural tissues. 26 The dPFG pulse sequences analyzed with DiTSI can be used to assess the distribution of axon fiber sizes 27 as well as mitigating the problem of gyral bias. 26 Furthermore, this pulse sequence has been assessed with both numerical simulation and validated on both 3T and 7T MRI scanners. 16,[28][29][30] Although dPFG pulse sequences may have enhanced sensitivity to tissue microstructure compared with sPFG sequences, this application has not been previously explored in skeletal muscle. The goal of this study is to evaluate whether dPFG dMRI pulse sequences are predictive of microstructural features of skeletal muscle that are related to muscle function via numerical simulation. Additionally, to validate that the findings from numerical simulation can be replicated in a real experiment, a rat model of muscle hypertrophy was imaged using these tools on a 7T MRI. This study presents preliminary experiments designed to determine whether dPFG provides useful new information on skeletal muscle microstructure and whether the technique warrants further investigation. We demonstrate that dPFG dMRI with DiTSI analysis can potentially be used to assess key features of muscle microstructure related to muscle function (i.e., fiber size, fiber area).

METHODS
The simulation-based approach used in this study is similar to the approach used in previous investigations of relationships between the diffusion tensor and muscle microstructure assessed with multi-echo spin-echo dMRI and time-dependent stimulated-echo dMRI. 31,32 The "Overview of DifSim" and "Histology Informed Model Generation" sections have been previously reported. 31,32 All studies involving animals were approved by the University of California, San Diego, Institutional Animal Care and Use Committee.

Overview of DifSim
DifSim is a numerical simulation software capable of simulating entire dMRI experiments with user-defined sequence parameters in tissues with arbitrarily complex geometry. The location of molecules is tracked using MCell, a validated Monte Carlo simulator for simulating diffusion in cellular microphysiology. As molecules diffuse during a simulated experiment, particle location, magnetization amplitude, and phase are monitored, and the resulting signal over all molecules is calculated. A detailed description of DifSim can be found in Balls et al. 33 and Baxter et al. 34

Histology-informed model generation
To relate dMRI measurements to physiologically accurate models of muscle, we created geometric models based on histology from previously performed animal experiments. A detailed explanation of how models were generated can be found in Berry et al. 31 The models used in this study include healthy control, Cardiotoxin (muscle regeneration, botulinum toxin Botox, chemical denervation/atrophy), surgical denervation (atrophy), and surgical tenotomy (acute hypertrophy, chronic atrophy). A single model was generated from histology from each model at 1, 3, 7, 14, and 30 days after injury. These injury models represent a wide spectrum of common microstructural changes that are found in skeletal muscle associated with injury or pathology. From each model, feret diameter, fiber area, and surface area-to-volume ratio (SA/V) was calculated. Extracellular water volume fractions (corresponding to the magnitude of edema) were approximated from histologic and MRI studies of these tissues and assigned to each model (see Berry et al. 31 ). No model was made with geometry from a Day 1 cardiotoxin model, as there are no clearly histologically defined muscle fibers at this timepoint.

Simulation details
No noise was added to the simulated dMRI signal to measure the exact relationship between muscle microstructure and quantitative measures of restricted diffusion under ideal conditions. Intracellular and extracellular particles were assigned different diffusion coefficients and magnetic relation rates (T 1 , T 2 ) based on literature values of these tissue at 7 T: intracellular: T 1 /T 2 : 1740/25 ms, 35-37 D: 1.8 × 10 −3 mm 2 /s [38][39][40] ; extracellular: T 1 /T 2 : 2500/95 ms, 35,37 D: 2.2 × 10 −3 mm 2 /s. 41,42 Particles were defined as impermeable to the sarcolemma. A minimum of 200 000 particles were simulated to accurately converge on an analytical solution based on the diffusion coefficients and b-values chosen for this experiment. 33 Myofilaments within a muscle fiber, or extracellular matrix proteins outside of muscle fibers, were not physically defined in this model, although it is at least partially reflected in the assigned diffusion coefficients, taken from previous studies of diffusion of small molecules in these tissues. All simulations were run on a Linux cluster with an Intel Xeon E-2697 CPU (2.60 GHz), with one node with 56 cores. The amount of time to run each simulation was less than 60 s.

Experimental data collection
To demonstrate that the dPFG technique can be reproduced to assess microstructure in real skeletal muscle, a well-described rat model of muscle hypertrophy (synergist ablation) was generated and scanned ex vivo on a 7T system at the University of Arizona. Additionally, to compare the inherent dynamic range of diffusion measurements in the same tissue given a standard sPFG and dPFG protocol, an sPFG dMRI sequence was also acquired. Briefly, under anesthesia, an incision was made on the medial calf of a rat, and the gastrocnemius and soleus muscle were surgically removed from the Achilles tendon, leaving the plantaris as the only functional plantar flexor. After 14 days, the animal was sacrificed, the legs harvested, and then fixed in formalin for 3 days. Then, the legs were placed in Galden and  44 From an animal that had undergone the same procedure but was not scanned and remained at the University of California, San Diego, histologic sections were obtained and stained with wheat germ agglutinin to stain the basement membrane, from which the distribution of fiber minimum Feret diameter, fiber area, and SA/V were calculated using the MATLAB command regionprops.

DiTSI analysis
For a full description of DiTSI analysis, please see Frank et al. 26 Briefly, from the second sampling scheme (Figure 2), we can construct a quantity of lower dimensionality but physical significance, the radial and angular SDs of , as follows: where = ( , ) and is the angular average of the spin density function (r, The angular variance 1 , 2 characterizes diffusion properties at different angular scales and contains a great deal of information about the microstructural features of the tissue of interest. Although it is a multidimensional tensor and a practicality, it is useful to derive average quantities of lower dimensionality, particularly scalar quantities that can be overlayed to easily discern spatial variations in diffusion properties. Averaging the radial variances matrix over all N × N radii provides the scalar quantity spherical anisotropy (SA) as follows:

Experimental data analysis
The dPFG data were analyzed using the DiTSI method described previously. From the spin-echo sPFG data, the diffusion tensor was calculated using the AFNI command 3dDWItoDT. 45  Regions of interest were drawn for the plantaris muscle using the AFNI draw data set tool. Regions of interest were reformatted to the dPFG, DTI, and T 2 map sequences using the AFNI command 3dresample.

Statistics
To determine if there is a relationship between muscle microstructure and SA, separate linear regressions were performed for fiber area, fiber diameter, and SA/V. To understand the distribution of SA, diameter, area, and SA/V measurements within a tissue of interest, histograms were calculated. All statistics were performed in Prism 9 (GraphPad, La Jolla, CA, USA).

Simulation: Radial dMRI
From the radial sampling scheme, the dPFG signal plotted over a range of relative orientations (G 2 ) for each of the gradient vectors (G 1 ) demonstrates an isotropic diffusion signal for normal healthy muscle ( Figure 3). There is no apparent angular dependence ( × ) on the diffusion signal for any of the diffusion sensitivities investigated (q-values). However, for the surgically denervated (atrophic) model, an anisotropic diffusion signal as a function of and was observed ( Figure 3). The disparity between diffusion signals for different θ increased with increased diffusion sensitivity (q-value). With denervation

F I G U R E 3
Simulated experiments. Simulated double pulsed field gradient experiment using a radial sampling scheme demonstrates differences in the diffusion signal as a function of . Control (top) and denervated (bottom) models of skeletal muscle were investigated for = 0 • (pink), 30 • (red), 60 • (green), and 90 • (blue). Increased asymmetry of the diffusion signal is observed in the denervated skeletal muscle compared with control, a result of the highly angular nature of denervated skeletal muscle fibers. injuries, skeletal muscle fibers adopt a more angular shape, which is uniquely detected using a radial dPFG sampling scheme.

Simulation: DiTSI analysis
The models with histology-informed muscle microstructure covered a range of fiber diameters from 28.9 to 65.9 μm, areas from 986 to 5097 μm 2 , and SA/V from 0.13 to 0.28/μm. From the second sampling scheme, SA for each model was successfully calculated. Overall, SA was found to be an excellent predictor of muscle fiber diameter ( Figure 4A; r 2 = 0.83; p < 0.0001), fiber area ( Figure 4B; r 2 = 0.71; p < 0.0001), and SA/V (Figure 4c; r 2 = 0.97; p < 0.0001).

Ex vivo high-field MRI: DiTSI analysis
To demonstrate that this technique is translatable to real experiments, an animal model of muscle hypertrophy was scanned ( Figure 5A,B). In the hypertrophied muscle, a decrease in SA was observed compared with the contralateral control muscle ( Figure 5C). Although the mean difference in SA was not large (hypertrophy = 0.201 ± 0.082 vs. control = 0.224 ± 0.077), the voxel

F I G U R E 4
Simulated experiments. Correlations between spherical anisotropy (SA) and fiber diameter (A), fiber area (B), and surface area to volume ratio (C) from simulated models of muscle microstructure. Significant correlations with good agreement were found for each relationship.  distribution of SA for the hypertrophied muscle was right-skewed in comparison with a normal distribution for the contralateral control ( Figure 5D). Histologic examination of hypertrophy versus control muscle revealed increased fiber diameter (42.1 ± 9.7 μm vs. 31.0 ± 8.0 μm), increased fiber area (1865 ± 792 μm 2 vs. 1044 ± 458 μm 2 ), and decreased SA/V (0.224 ± 0.049/μm vs. 0.321 ± 0.095/μm) in hypertrophied muscle versus control, respectively ( Figure 6). Overall, the distribution of these microstructural values was broad and demonstrated that there is a wide variance in the microstructural features observed, similar to the SA distributions. To compare these findings with current DTI techniques, a spin-echo sPFG pulse sequence was also performed. Gross visualization of the fractional anisotropy (FA) versus SA maps for the same tissue indicates that SA has much greater contrast and dynamic range than FA ( Figure 5C,E). In addition, the distribution of FA for voxels within the hypertrophic and control region of interest is much smaller for FA than SA. To verify that changes in the diffusion signal were not due to inflammation from the surgical process, T 2 relaxation of the tissue was calculated ( Figure 5G). The mean T 2 relaxation of the control and hypertrophic muscle was 27.8 ms and 30.8 ms, respectively, and the distribution of T 2 for voxels within the regions of interest was similar ( Figure 5H), suggesting that the differences in SA were driven by microstructure and not increased water content due to edema.

DISCUSSION AND CONCLUSIONS
The goal of this study was to perform a preliminary investigation into the use of dPFG pulse sequences to study muscle microstructure. In addition, the sensitivity of a novel analysis technique to quantify restricted diffusion from dPFG dMRI experiments to key features of muscle microstructure was investigated. The key findings from this study are as follows: (1) Radial dPFG shows sensitivity to a phenotype of muscle injury that is difficult to observe without histology (i.e., angular fibers); (2) SA-a scalar value from DiTSI analysis-is highly sensitive to muscle microstructural features predictive of function (i.e., fiber size, fiber area); and (3) these techniques and analysis tools can be translated to real experiments in skeletal muscle. The dynamic range of SA derived from the dPFG data in its relation to fiber size is greater than what we can measure using sPFG pulse sequences with a standard DT model. Increased dynamic range of this measurement allows for increased sensitivity to detecting changes in underlying tissue microstructure, which allows for more refined applications of this technique to support noninvasive assessment of muscle physiology.
The dMRI pulse sequences have the potential to be a highly valuable clinical tool if they can be developed to assess subtle changes in muscle microstructure, which is required to evaluate muscle health, monitor disease progression, diagnose musculoskeletal disease, or monitor efficacy of a treatment. Currently, there is no consensus on the optimal parameters for dPFG sequences in clinical settings, but generally they will vary with application. Our objective in the current study is to investigate the relationship of the dPFG parameters to muscle microstructure-a necessary preliminary step to any determination of parameters for clinical optimality. As a first step, careful simulations have been used to elucidate sensitivity to various pulse sequences, pulse sequence parameters, SNRs, and states of tissue health. However, translating these findings to real-world experiments presents a challenge, as there are limitations associated with real-world experiments such as image artifacts, scan time, movement, and SNR. Furthermore, there is a need to carefully validate these measurements with phantoms or histology in preclinical animal models to support potential clinical adoption of dMRI in the clinic. Currently, dPFG pulse sequences are not standardized on preclinical or clinical MRI scanners. The pulse sequences used in this study using a preclinical 7T animal MRI, and in a previous study using a clinical 3T MRI, 26 had to be programmed by an experienced MRI pulse sequence programmer. If dPFG sequences can be reliably shown to be a useful research and clinical tool, these sequences may become stock sequences from MRI vendors.
Healthy muscle fibers have a polygonal shape and normal distribution of fiber sizes. With physical or chemical denervation, fibers tend to adopt a more angular shape in addition to undergoing atrophy. Denervation has been widely studied using dMRI. Generally, increased FA and decreased radial diffusivity are observed in chronically denervated muscle, which is consistent with smaller muscle fibers. No difference in λ1 is typically observed, suggesting that the microstructural changes are primarily related to the size of muscle fibers. 7 Acutely after denervation injury, increased mean diffusivity and decreased FA have been observed, which suggests enlargement of the extracellular fluid space. 46 Although the ratio between λ2 and λ3 has been proposed to be a metric of fiber ellipticity, 47 this theory has never been histologically validated, although its popularity has propagated through the literature since its proposal. However, using standard spin-echo DTI, it is unclear whether atrophy due to denervation can be separated from age-related atrophy (e.g., sarcopenia), disuse atrophy (e.g., immobilization), or non-nerve-related pathology (e.g., muscular dystrophy). Although it is currently not clear whether each of these categories of atrophy are biologically distinct and have differential effects on muscle microstructure that can be detected histologically or with dMRI. The radial dPFG analysis demonstrated the ability to detect angular changes in muscle fiber shape associated with denervation. Specifically, increasing anisotropy was observed with increasing q-value. As the maximum q-value that can be obtained experimentally is based on maximum gradient strength of an MRI scanner, translating the experiments to a clinical scanner will be hardware-dependent and thus may vary from vendor to vendor. Although these results are preliminary, they suggest that the radial sampling scheme used has high angular sensitivity, which is sensitive to the presence of angulated fibers. Future work will further investigate the sensitivity of this method to detecting acute versus chronic denervation associated atrophy.
A major limitation to the use of dPFG pulse sequences is the lack of standard analysis tools to quantify restricted diffusion. The dMRI methods based on sPFG pulse sequences are the current standard on clinical scanners, which typically include rather standard DTI analysis methods on the MRI console itself (e.g., fitting diffusion data to the tensor model). Unfortunately, the shortcoming of the de facto standard sPFG sequences is their insensitivity to anisotropy variations, whose average over a voxel is zero or low enough to obscure microscopic effects. This is particularly problematic in a tissue with mixture of elongated structures of different widths and orientations, and cells with a range of sizes and shapes. These all induce local anisotropy that is poorly detected (if at all) by sPFG, 26 and whose sources are multivariate, making estimation of microscopic tissue properties an exceedingly poorly posed inverse problem. Despite numerous compartment modeling approaches, 35,[48][49][50] these sPFG methods are ultimately limited in both their sensitivity and specificity. This problem is generally insurmountable without a higher-order diffusion encoding such as dPFG. DiTSI is an analysis technique that decomposes dPFG collected with a Cartesian product of n gradient directions into an n 4 -dimensional hyper-volume Q. The analysis of high-dimensionality data is difficult, and the difficulty can be reduced by averaging over appropriate radial or angular subspaces to produce lower dimensional anisotropy metrics. The two useful scalars are (1) radial variances of Q averaged over all angles (i.e., the radial anisotropy) and (2) angular variances averaged over all (or range of) radii (i.e., the SA). The lowest-order component of radial anisotropy is approximately equal to the SA; therefore, for simplicity, we can limit our study using SA only. This simple measure produces detailed local anisotropy maps caused by microstructural variations with sensitivity and specificity not attainable with sPFG analysis methods.
Assessing changes in fiber size is one of the most pressing clinical applications for dMRI. Currently, there are two main models that attempt to quantify muscle fiber size from dMRI data: (1) the random permeable barrier model 51 and (2) the apparent muscle fiber diameter (AFD) model. 8 Random permeable barrier model uses stimulated-echo sPFG dMRI to measure restricted diffusion at a number of diffusion times ( = 20 −> 1000 ms) that is fit to a model which considers water diffusion hindered by a network of randomly oriented semipermeable membranes. This approach has been used to quantify changes in microstructure associated with postnatal growth, 52 postexercise fiber dilation in healthy patients but not in patients with chronic extracellular compartment syndrome, 53 rotator cuff repair, 54 and disuse atrophy and recovery. 55 Overall, this technique is highly accurate at measuring the SA/V of a muscle fiber and can detect changes in fiber permeability, which may be an early marker of muscle injury. Both RPBM and the DiTSI analyses presented here have demonstrated high sensitivity to SA/V. However, SA/V is not a standard muscle microstructural measurement, and its relationship to force generating capacity is the focus of ongoing work. The AFD model is relatively new and relies on fitting acquired data to a forward simulated model that relates restricted diffusion measurements from a spin-echo, multishell DTI pulse sequence to a cylinder model of muscle fibers. 8 This dictionary-based approach has not yet been validated against histology, and the predicted fiber sizes are larger than what has been previously reported in the literature. 56 Although the AFD model has predicted fiber atrophy associated with denervation versus control (60.7 vs. 89.7 μm), this same model has also predicted fiber atrophy in non-denervated muscles (71.6 μm). While this model can assess expected changes in AFD associated with injury, it still must be validated against real tissue microstructures in order to have clinical utility. However, in their use of sPFG acquisitions, both RPBM and AFD are intrinsically limited to investigations of macroscopic anisotropy, unlike the DiTSI method, which is based on dPFG acquisitions.
Given the relationship between SA and fiber diameter derived in Figure 4A, the predicted mean fiber size of the hypertrophy and control muscles was 46.0 and 42.6 μm. These predicted measurements were larger than the measured fiber sizes for hypertrophy (42.1 μm) and control (31.0 μm) muscle in a different sample. Histology is inherently very localized, as the scale at which fiber size is measured can be orders of magnitude smaller than achieved voxel dimensions, depending on the imaging protocol. Rectifying mean fiber sizes with mean diffusion measurements averaged over the entire volume of a muscle is difficult, but it is currently the best way to validate dMRI sensitivity to fiber size. Spatially matching histologic measurements to dMRI measurements using histologic images registered to MRI data would allow for more specific regional variations in fiber sizes to be assessed and is the subject of ongoing work in the lab.
This study has a few limitations. Intracellular structures were not directly modeled, as they cannot be clearly defined from histology. However, the presence of these structures is partially accounted for by the assigned diffusion coefficients for molecules in the intracellular and extracellular regions, which was derived from careful measurements in these structures. Additionally, a systematic evaluation of the effect of additional muscle microstructure components on SA and additional quantitative measurements from DiTSI analysis was not performed. An extended phase graph algorithm was not used to calculate T 2 relaxation, as the effect of B 1 field inhomogeneity was assumed to be negligible, and the analysis technique that was used (Levenberg-Marquardt fitting) is commonly used to assess T 2 changes associated with inflammation in skeletal muscle. 57,58 Moreover, a computational model of muscle hypertrophy was not simulated in this study, because the model of muscle hypertrophy that was scanned has a muscle fiber size (42 μm) that was in the dynamic range of the existing muscle models (28-65 μm). This is because the hypertrophy model was performed on 3-month-old rats that have not yet reached skeletal maturity, whereas the histology-informed models that were simulated were generated from skeletally mature 6-month-old and 12-month-old rats. Furthermore, there are several pulse-sequence parameters that can be evaluated, including b-value(s), primary and secondary diffusion directions, diffusion time, voxel size, TE, and SNR, which can all potentially affect the resulting restricted diffusion measurement. This study was designed to be a preliminary analysis and feasibility study to determine whether dPFG pulse sequences should potentially be applied to study skeletal muscle microstructure. Future in silico experiments will be performed to (1) determine the sensitivity of dPFG to muscle microstructure in idealized models in which key microstructural features of skeletal muscle (fiber size, fibrosis, edema, and permeability) will be systematically simulated over physiologically relevant dimensions individually and in combination; and (2) determine the optimal pulse-sequence parameters required to measure muscle fiber size in physiologically relevant situations (e.g., inflammation, SNR). After these relationships have been established using numerical simulation, these findings can be validated using 3D-printed anisotropic diffusion phantoms with user-defined microstructure. 9