Assessment of fitting methods and variability of IVIM parameters in muscles of the lumbar spine at rest

Intravoxel incoherent motion (IVIM) MRI provides insight into tissue diffusion and perfusion. Here, estimates of perfusion fraction (f), pseudo-diffusion coefficient (D*), and diffusion coefficient (D) obtained via different fitting methods are compared to ascertain (1) the optimal analysis strategy for muscles of the lumbar spine and (2) repeatability of IVIM parameters in skeletal muscle at rest. Diffusion-weighted images were acquired in the lumbar spine at rest in 15 healthy participants. Data were fit to the bi-exponential IVIM model to estimate f,D* and D using three variably segmented approaches based on non-linear least squares fitting, and a Bayesian fitting method. Assuming that perfusion and diffusion are temporally stable in skeletal muscle at rest, and spatially uniform within a spinal segment, the optimal analysis strategy was determined as the approach with the lowest temporal or spatial variation and smallest residual between measured and fit data. Inter-session repeatability of IVIM parameters was evaluated in a subset of 11 people. Finally, simulated IVIM signal at varying signal to noise ratio were evaluated to understand precision and bias. Experimental results showed that IVIM parameter values differed depending on the fitting method. A three-step non-linear least squares fitting approach, where D,f, and D* were estimated sequentially, generally yielded the lowest spatial and temporal variation. Solving all parameters simultaneously yielded the lowest residual between measured and fit data, however there was substantial spatial and temporal variability. Results obtained by Bayesian fitting had high spatial and temporal variability in addition to a large residual between measured and fit data. Simulations showed that all fitting methods can fit the IVIM data at signal to noise ratios >35, and that D* was the most challenging to accurately obtain. Overall, this study motivates use of a three-step non-linear least squares fitting strategy to quantify IVIM parameters in skeletal muscle.


Introduction
Measurement of blood flow to skeletal muscle provides insight into muscle metabolism, and activation in response to exercise.Non-invasive magnetic resonance imaging (MRI) methods to measure microvascular blood flow in skeletal muscle such as arterial spin labeling (ASL) and intra-voxel incoherent motion (IVIM) have been employed in a variety of muscles, including those in the leg (1-4), shoulder (5,6), and lumbar spine (7)(8)(9).While ASL can quantify perfusion in physiologically relevant units, its spatial extent is somewhat limited and the perfusion response to stimuli such as exercise or induced ischemia dissipates quickly (10).Moreover, it remains difficult to quantify the low perfusion state of muscle at rest due to the limited signal-to-noise ratio (SNR) of ASL data.
Alternatively, IVIM is based on the incoherent translational motion of intravascular and intracellular water, and therefore is sensitive to the entire capillary blood volume within a voxel.Building off the theory of molecular diffusion, the motion of water in randomly oriented capillaries will also result in incoherent phase accrual, leading to signal decay in the presence of diffusion-encoding gradients (11,12).Due to the difference in scale of motion in the intracellular (e.g., molecular diffusion, D) vs. intravascular (so called pseudo-diffusion, D*) spaces, signal in intravascular space decays much more quickly than the intracellular component (13).The molecular and pseudo-diffusion coefficients, along with the volume fraction of the intravascular space (i.e., perfusion fraction, f) can be estimated by fitting data acquired at multiple diffusion encoding strengths to a bi-exponential model.
To estimate the IVIM parameters f, D, and D* from diffusion-weighted images, a variety of fitting strategies have been employed (14)(15)(16).However, in skeletal muscle, where resting perfusion is quite low, the rationale behind the choice of a particular analysis approach has not been clearly and methodically evaluated.Currently, there is not a single agreed upon analysis strategy, resulting in a lack of comparability of results across studies.Additionally, given that studies using IVIM to evaluate skeletal muscle perfusion may interrogate the response to exercise either acutely (same day) or over the course of some therapeutic intervention, it is also important to understand the repeatability of IVIM parameters in skeletal muscle at rest.Thus, the purpose of this work was to compare the estimates of IVIM parameters obtained via different fitting methods to ascertain the optimal strategy and normative values for muscles of the lumbar spine and to determine the variability of IVIM parameters across sessions.The findings of this study are essential for enhancing comparability across studies and establishing a consensus analysis strategy for future research.

Study design
This study was approved by the Institutional Review Board of the University of California, San Diego.Sixteen healthy people were recruited for the study after providing written informed consent.Participants were excluded if they had contraindications to MRI or had a history of low back pain.Individuals were instructed to avoid any physical activity on the day of the testing session.MR imaging was performed at 3T (MR750, GE Healthcare, Milwaukee, Wisconsin, USA) using a CTL phased-array spine coil for signal reception.The imaging protocol consisted of localizers, T1-weighted anatomical reference images, followed by collection of IVIM data while the subject rested in the supine position.Twelve participants returned for a second, identical visit, allowing for the evaluation of inter-session repeatability of the IVIM parameters.

Image analysis
Image acquisition and analysis procedures are summarized in Figure 1.Specifically, the single IVIM acquisition described above was used for all sub-analyses, which are described in detail subsequently.Briefly, the four repeats were analyzed independently to evaluate intra-session temporal stability across the entire paraspinal muscle ROI.To assess intrasession spatial consistency the repeats were averaged across the entire paraspinal muscle ROI, and the mean diffusion-weighted signal in each slice was used for analysis.To assess Inter-session repeatability, repeated exams on separate days were compared in the subset of participants in which data were available.

2.3.1
Preprocessing-Phase and distortion correction of the diffusion-weighted images was performed via TOPUP based on blip-up and blip-down calibration scans acquired at the end of the IVIM acquisition (17,18).Denoising was conducted using a local principal component analysis filter (19).Rigid body motion correction over the dynamic image series was performed in AFNI (20).

2.3.2
Image processing-Preprocessed data were imported into MATLAB (R2017a, Natick, MA), and subsequent analysis was performed using in-house written code, which can be found at https://github.com/ekenglund/Muscle_IVIM_analysis.First, a 2D 3 × 3 median filter was applied.IVIM data were then averaged over the three diffusion encoding directions for each b-value, yielding four series of low-to-high b-value images (one series per repeat).The four repeats were initially analyzed separately to assess temporal stability, then were averaged together prior to the calculation of slice-wise and voxelwise IVIM parameter maps (described below).SNR was calculated from a b = 0 s/mm 2 image for each subject by dividing the mean of the signal of a region of interest (ROI) placed within the muscle by the standard deviation of the signal of an equally sized ROI placed outside of the body (21,22).

IVIM parameter fitting-
For each analysis, IVIM parameters f, D*, and D were obtained assuming bi-exponential signal decay as a function of b-value (13): where S b /S 0 is the measured diffusion-weighted data normalized by the non-diffusion weighted acquisition, f is the perfusion fraction, representing the relative fraction of signal from the intravascular vs. extravascular (i.e., intracellular) spaces and serving to weight the bi-exponential model by the two coefficients, D* the pseudo-diffusion coefficient, and D the diffusion coefficient.The product of f and D* fD* was also calculated.Four separate approaches were evaluated: • 1-step non-linear least squares (NLLS) fitting; all three IVIM parameters were computed simultaneously using an NLLS fitting approach.= A e −bD where A is an offset term.Once D was obtained, the value was fixed and all data were used to determine f and D* from the full bi-exponential equation via NLLS fitting.
• 3-step NLLS fitting; first the natural log of the high b-value data were fit to a first order polynomial to determine D. This line was extrapolated to b = 0 s/mm 2 to obtain f, and finally the NLLS fit of the full bi-exponential model was used to calculate D*.
Bayesian fitting; a Bayesian-probability based approach was used to estimate D, f, and D* based on the joint posterior probability as described in (23,24).See supplemental methods for specific details of the Bayesian fitting analysis.
In addition to the obtained parameter values, the norm of the residual between the measured and fit data was recorded for each fitting method.
ROIs encompassing in the left and right erector spinae and multifidus (together referred to as paraspinal) muscles were manually defined on all slices of the T1-weighted images in Horos (v3.3.6, www.horosproject.org) and transferred to the IVIM data as previously described (25).In the final mask, voxels were excluded if noise or motion caused the signal intensity of the first diffusion-weighted acquisition to exceed the non-diffusion-weighted acquisition.Additionally, voxels were excluded if the solution for any of the voxelwise fits of the IVIM parameters was equal to the set fit constraint to prevent bias or overfitting of the signal to the IVIM equation.

2.3.4
Defining the optimal fitting strategy-Two assumptions were made in order to define the optimal fitting strategy for IVIM data collected in skeletal muscle: (1) IVIM parameters are temporally constant in skeletal muscle at rest ( 26); (2) the variability of IVIM parameters between adjacent slices within a spinal segment in healthy paraspinal muscle is likely negligible due to its segmental innervation and function as an intervertebral stabilizer (27,28).Based on these assumptions, we considered the optimal method as one that minimized the residual between the measured and fit data and had the least variability over time and between adjacent slices (e.g., intra-session repeatability).Variability was assessed by the within-subject coefficient of variation (WS-CV) (29).

Intra-session temporal stability-Average signal intensity was computed
across the entire paraspinal ROI for each b-value and repeat (Figure 1B), leading to four temporally-resolved signal decay curves per participant.Data were normalized by the nondiffusion weighted signal intensity, and IVIM parameters were computed via the four fitting methods for each repeat of the IVIM acquisition.The mean and standard deviation (SD) of f, D*, D, fD*, and the squared residual norm were computed for each fitting method, and variability of IVIM parameters was assessed by the WS-CV.

Intra-session spatial consistency-Diffusion
weighted data were averaged across the four repeats, and within the ROI for each of the 22 slices, yielding one signal decay curve per slice.IVIM parameters were then computed with each fitting method for each slice (Figure 1B).Four slices spanning the superior and inferior margins of the vertebral body were assigned to each lumbar spine segment (L1-L5), and the variability within a spinal segment was evaluated.Mean and SD of IVIM parameters and mean squared residual norm for each fitting method were determined.Assuming that physiologically, perfusion and diffusion are relatively constant within a spinal segment (e.g., across all four slices located within L1), the WS-CV was computed for each of the five spinal segments as a metric of spatial consistency.

2.3.7
Inter-session repeatability-For inter-session repeatability analysis, voxelwise parameter maps were computed within the defined paraspinal muscle ROI for the baseline and repeat scans using each of the four previously defined fitting methods (Figure 1C).Following voxelwise quantification of f, D*, D, and fD*, IVIM parameter means and SD were computed over the entire paraspinal ROI as well as within the erector spinae and multifidus muscles separately, and results were compared between fitting methods.The SD of IVIM parameters for all voxels within the ROI was computed for both the baseline and repeat data to assess dispersion.The average percent of voxels included in the final ROI was also determined and compared between analysis methods.WS-CV was computed to assess the repeatability of ROI-averaged IVIM parameters across scan sessions.

Simulation
To understand the bias and precision of these fitting methods under relevant signal to noise ratio (SNR) conditions, a Monte Carlo based simulation was performed.Step NLLS, and Bayesian fitting approaches.An interrater correlation coefficient (2, 1) was calculated for each fitting method for each SNR level to assess the bias between fit parameters of f, D, and D* and the original values.fD* was not explicitly evaluated in the simulation model, as it the product of two o coefficients that were directly in the model.Therefore, its accuracy depends on the individual accuracy of these parameters.

Statistical analyses
Repeated measures one-way Analysis of Variance (ANOVA) tests were used to compare the IVIM parameter values between fitting methods for each analysis.

Results
Image artifacts precluded analysis in one subject, leaving 15 participants (mean ± SD, age = 44.4± 13.0 years, 8 male) for intra-session and 11 (42.5 ± 14.6 years, 6 male) for inter-session repeatability assessments.The average time between scan sessions was 76 ± 37 days.SNR of the images was 46.6 ± 10.6.

Intra-session temporal stability
Averaged over the entire paraspinal muscle ROI, it was observed that overall, the choice of analysis method resulted in significantly different IVIM parameters (p < 0.0001 for all parameters, Figure 2A), and three-step fitting had the least temporal variability (Table 1).
The average squared norm of the residual was highest for Bayesian fitting (0.0016) and was improved for all NLLS fitting methods.WS-CV was minimized for all IVIM parameters with the three-step fitting approach.

Intra-session spatial consistency
IVIM parameters calculated from slice-wise mean signal intensity, averaged within a spinal segment similarly differed based on the choice of analysis method (p < 0.0001), and additionally all parameters differed as a function of spinal segment (p < 0.0001) (Figure 2B).The mean residual was smallest when using one-step fitting (0.0007), however this led to substantial variation in the quantified parameters.The average WS-CV for one-step fitting was 24.2% for f, 45.2% for D*, 5.9% for D, and 40.2% for fD*.Three-step fitting had the least variability within a spinal segment.The mean WS-CV for three-step fitting was 15.4% for f, 24.5% for D*, 4.2% for D, and 29.2% for fD* (Table 1).Regardless of fitting strategy, it was generally observed that across spinal segments, greater variability in IVIM parameters was observed in L1 and L5 compared to L2-L4 (data are included as Supporting Information, Supplementary Table S1 and Figure S1).

Inter-session repeatability
Figure 3 shows quantified IVIM parameter maps for a representative subject using each of the four fitting strategies.Comparing the two visits, when voxelwise IVIM parameters were averaged across the entire paraspinal muscle ROI, the WS-CV values were in general lower than those observed for the intra-session repeatability assessment (Table 1).Table 2 summarizes the results for the IVIM parameters at the first MRI scan session.In general, two-and three-step fitting methods yielded very similar results.One-step fitting led to higher f and D* values and lower D, and D* was significantly lower for Bayesian fitting.The spread of IVIM parameters (assessed by the mean SD across voxels) was generally lower for two-step and three-step fitting.The final mask comprised 62%, 72%, 71%, and 76% of voxels for one-, two-, three-step NLLS, and Bayesian fitting methods, respectively.

Simulation results and bias evaluation
110,245 different combinations of f, D, and D* were evaluated in this analysis.Overall, excellent agreement was found f and D for SNR levels over 35 for all fitting approaches (ICC > 0.8, Figure 4).The ability to fit D* was best for 1-Step NLLS fitting, and fair for 2-Step NLLS, 3-Step NLLS and Bayesian fitting.

Discussion
This study evaluated the variability of IVIM parameters in skeletal muscle at rest using a variety of fitting methods.With additional use of IVIM in healthy controls, in response to exercise or other interventions, or in patient populations, it is important to understand the typical IVIM parameter values in skeletal muscle at rest in healthy individuals, the repeatability of those parameters across scan sessions, and how fitting methods may influence those values.This methodical evaluation of the various fitting strategies demonstrated that the selected fitting strategy directly influences the calculation and interpretation of output IVIM parameters.Provided that our assumptions that skeletal muscle at rest has relatively temporally stable (over several minutes) and spatially consistent (within one spinal segment) perfusion and diffusion characteristics are valid, the optimal non-linear least squares fitting approach appears to be the segmented three-step fitting method as intra-session temporal and spatial variability was minimized and few voxels were excluded from the final mask.
The simulations performed herein echo the experimental results, finding that the one-step NLLS fitting method had the lowest bias between known and fit parameters, but that all fitting methods provide reasonable estimates of the IVIM parameters at SNR >35.Given that the average SNR of our experimental data was 46.6, we posit that there is negligible difference of bias or precision between fitting methods.A prior study by Damon defined the threshold for accurate diffusion measurements as 40 (30), in relative agreement with what was observed in our simulation.However, in general, image acquisition should be optimized to achieve high SNR in order to ensure accurate measurements of perfusion and restricted diffusion.
In addition to SNR, another important factor to consider for IVIM analysis in skeletal muscle is fat suppression.Both intra-muscular fat depots as well as subcutaneous fat can lead to potential errors IVIM.In this study, spatial-spectral fat suppression was used.Further, the low-bandwidth phase encoding direction was selected to be anterior to posterior, such that any unsuppressed signal from subcutaneous fat would alias posteriorly, outside of the anatomy of interest.Nevertheless, the seminal study by Cameron et al. (31) demonstrates that effective fat suppression is crucial for accurate measurement of diffusion parameters including IVIM.A conclusion from the Cameron, et al. study is that even with triple-fat suppressed data, there may still be reliability issues when estimating D*.Both the Cameron et al. and this paper agree that D* is routinely the most difficult IVIM-based parameter to accurately estimate given currently available pulse sequence parameters and analysis approaches.
While one-step fitting minimized the residual between measured and fit data in both experimental and simulated data, it resulted in substantial temporal and spatial variability and the final mask contained the fewest voxels due to solutions occurring at the set boundary condition for at least one parameter from experimental data.Several prior studies have evaluated fitting methods for analysis of IVIM parameters (15,16,32,33).In general, those studies have focused on the evaluation of IVIM parameters in cancerous tissues and have found, in agreement with the results herein, that parameter estimates for f, D, and D* differ depending on the fitting algorithm selected.A number of additional algorithms have been evaluated recently, including convolutional neural network based parameter estimation strategies (34).Overall, the results presented herein demonstrate the importance of considering the fitting method when comparing results across prior studies.
In general, it was observed that D had the least intra-session spatial and temporal variability, and the lowest inter-session variability of the quantified IVIM parameters.Perfusion fraction was stable with intra-session WS-CVs of ~15% and inter-session WS-CVs <10%.D* and fD* had higher variability, suggesting that it may be more difficult to detect the effect of an intervention with these parameters.However, prior studies have demonstrated changes in D* on the order of 100% (5-9) and as large as 289% (3) in response to acute bouts of exercise.Nevertheless, the variability and sensitivity of parameters needs to be taken into consideration when powering studies and interpreting results.
Comparing these results to prior studies, there has been a wide range of IVIM parameters reported in muscle (8,9,35).This variability may in part be due to the lack of standardization of fitting approach, as results herein show significant differences in parameter values depending on choice of analysis method.A recent study also evaluated the intra-and inter-session repeatability of IVIM parameters in tumor and muscle, finding improved repeatability with Bayesian fitting, however the mean values of IVIM parameters in muscle were not provided (36).Ranges for IVIM parameters in skeletal muscle have spanned 0.5 (37) to 2.19 (4) × 10 −3 mm 2 /s for D; 0.03 (3) to 0.26 (37) for f; and 8 (38) to 287 (39) × 10 −3 mm 2 /s for D*.
Focusing on only investigations of the paraspinal muscles, a recent study by Federau, et al. examined the effect of lumbar extension exercise on IVIM metrics in controls and patients with adolescent idiopathic scoliosis (7).IVIM parameters were quantified via two-step fitting in the paraspinal muscles at rest in healthy volunteers were in reasonable agreement with those reported herein; f, D*, and D: 0.10, 32 × 10 −3 mm 2 /s, and 1.4 × 10 −3 mm 2 /s, (right), and 0.8, 53 × 10 −3 mm 2 /s, and 1.4 × 10 −3 mm 2 /s (left).In addition to evaluating the IVIM parameters at rest, we have previously assessed the response to high intensity lumbar spine extension exercise in healthy individuals and people with low back pain using a three-step fitting approach (9).The response to varied exercise intensity was also evaluated in healthy individuals (8).Results of these prior studies showed that acute, targeted exercise is associated with an increase in IVIM parameters.Importantly, differences in IVIM parameters have been shown to be clinically relevant, for example, changes in D and D* in response to an acute exercise bout were shown to be highly predictive of reductions in disability and pain in response to an exercise-based rehabilitation program in individuals with low back pain, suggesting that this method may have utility in identifying individuals who may benefit from prescribed treatments (9).It also may help better understand the contribution of underlying microvascular impairments in various musculoskeletal conditions (7,40).
Our results here and in a previous, independent study (41) suggest that diffusion characteristics differ in the paraspinal muscles as a function of spinal segment.The observation of changing mean diffusivity over the span of the lumbar spine may draw into question the assumption of spatial uniformity, however we note that spatial consistency was assumed within one spinal segment, not across the entire spine.Indeed, differences in IVIM parameters across the spinal segments are likely to reflect differences in structure and physiology of the paraspinal muscles from cranial to caudal, including differences in size (including the relative contribution of erector spinae vs. multifidus), tissue composition, and metabolic capacity (41)(42)(43).We also note that the IVIM parameters were generally more variable at L1 and L5, potentially due to respiratory motion at the superior extent and coil sensitivity at the inferior extent.
This study has several limitations.First, there is no gold standard to compare the experimental results.Simulations of various bi-exponential fitting approaches in the context of varied measurement noise have been completed previously (15), and were performed here to evaluate bias.In this study, we sought to primarily make use of the acquired data, based on assumptions regarding the underlying physiology to evaluate the variability.However, the lack of a true gold standard with which to compare the experimental data remains a limitation.In addition, we did not evaluate the choice of cutoff value for the high b-value regime (44), which may be important when analyzing IVIM parameters in response to a physiologic challenge such as exercise, nor did we evaluate the influence of pre-processing steps, such as applying a median filter.Last, we did not evaluate any variation in the cutoff ranges of IVIM parameters, but instead chose to use a wide range of physiologically plausible values for f, D, and D*.Given that the mean D was less than 1.4 × 10 −3 mm 2 /s, the minimum D* was set to 1.5 × 10 −3 mm 2 /s.Overlapping ranges between D and D* could lead to issues with fit convergence, though this was not practically observed.

Conclusion
Overall, this study compares repeatability of IVIM parameters in skeletal muscle at rest using a variety of bi-exponential fitting strategies and helps to define normal values of IVIM parameters in muscles of the lumbar spine at rest in healthy adult volunteers.The three-step NLLS fitting approach yielded results with the least temporal and spatial variability, low inter-session variability, and a large inclusion rate of voxels in the final mask.Future studies evaluating muscles of the lumbar spine in response to acute bouts of exercise, pathology, or long-term rehabilitation therapy may provide additional insight into the physiologic underpinnings of the IVIM parameters in skeletal muscle.Representative IVIM parameter maps for each fitting strategy.In general, D was very spatially uniform across the paraspinal muscles.Relative anatomic contrast is observed in f, D*, and fD* maps for the 2-and 3-step NLLS fitting approaches.Results of simulation experiment.Ideal signal decay curves were generated from varied f, D, and D* values, then Gaussian noise was added to correspond to different levels of image SNR.Each fitting method was used to estimate IVIM parameters, and the interrater correlation coefficient (ICC) was computed between the simulated and fit parameter values.In general, 1 step non-linear least squares fitting had the least bias between simulated and fit parameters.The other fitting methods were in general agreement between simulated and fit data at SNR levels above 35.

FIGURE 1 .
FIGURE 1. Acquisition and analysis overview.(A) Shows acquisition scheme, and the right paraspinal ROI is indicated in yellow on the T1-weighted anatomical reference image.(B) Intra-session analyses were performed on the time-resolved (yellow volumes, whole muscle ROI), spatially averaged (to assess temporal stability) and slice-averaged, time averaged data (to evaluate spatial consistency) separately.(C) Inter-session repeatability was evaluated by first computing the voxelwise IVIM parameter maps for f, D*, D, and fD*, then averaging across the paraspinal ROI and comparing between scan sessions.Representative signal decay curve from a single voxel [blue dot in left panel of (C)] is shown in the middle panel of (C).In all cases, IVIM parameters were obtained by four different approaches including variably segmented, non-linear least squares fitting approaches and a Bayesian fitting method.

FIGURE 2 .
FIGURE 2. Summary of intra-session IVIM parameters showing mean and standard deviation across N = 15 healthy individuals.IVIM parameters are shown for the temporally resolved (A)and slice-wise resolved (B) analyses for each fitting method.Slice-resolved analyses in (B) were evaluated on a slice-by-slice basis then averaged across slices and summarized here.In general, there was a significant difference between IVIM parameter values depending on the choice of fitting method (p < 0.0001 for both time-resolved and slice-resolved analyses on RM-ANOVA).The slice-wise analysis revealed a significant difference of IVIM parameters

TABLE 1
Summary of intra-session and inter-session variability.WS-CV is within-subject coefficient of variation.

TABLE 2
Summary of mean (standard deviation) IVIM parameter values across the paraspinal ROI at visit 1.The dispersion of parameter values in the voxelwise map is presented on bottom.

mean (SD) of the dispersion across full paraspinal muscle ROI
Front Musculoskelet Disord.Author manuscript; available in PMC 2024 August 12.