Fast Open-Source Toolkit for Water T2 Mapping in the Presence of Fat From Multi-Echo Spin-Echo Acquisitions for Muscle MRI

Imaging has become a valuable tool in the assessment of neuromuscular diseases, and, specifically, quantitative MR imaging provides robust biomarkers for the monitoring of disease progression. Quantitative evaluation of fat infiltration and quantification of the T2 values of the muscular tissue's water component (wT2) are two of the most essential indicators currently used. As each voxel of the image can contain both water and fat, a two-component model for the estimation of wT2 must be used. In this work, we present a fast method for reconstructing wT2 maps obtained from conventional multi-echo spin-echo (MESE) acquisitions and released as Free Open Source Software. The proposed software is capable of fast reconstruction thanks to extended phase graphs (EPG) simulations and dictionary matching implemented on a general-purpose graphic processing unit. The program can also perform more conventional biexponential least-squares fitting of the data and incorporate information from an external water-fat acquisition to increase the accuracy of the results. The method was applied to the scans of four healthy volunteers and five subjects suffering from facioscapulohumeral muscular dystrophy (FSHD). Conventional multi-slice MESE acquisitions were performed with 17 echoes, and additionally, a 6-echo multi-echo gradient-echo (MEGE) sequence was used for an independent fat fraction calculation. The proposed reconstruction software was applied on the full datasets, and additionally to reduced number of echoes, respectively, to 8, 5, and 3, using EPG and biexponential least-squares fitting, with and without incorporating information from the MEGE acquisition. The incorporation of external fat fraction maps increased the robustness of the fitting with a reduced number of echoes per datasets, whereas with unconstrained fitting, the total of 17 echoes was necessary to retain an independence of wT2 from the level of fat infiltration. In conclusion, the proposed software can successfully be used to calculate wT2 maps from conventional MESE acquisition, allowing the usage of an optimized protocol with similar precision and accuracy as a 17-echo acquisition. As it is freely released to the community, it can be used as a reference for more extensive cohort studies.

Introduction

Neuromuscular disorders encompass genetic and acquired diseases of lower motor neurons, peripheral nerves, neuromuscular junction, or skeletal muscle, generally causing different degrees of motor impairment in the affected patients. In particular, muscular dystrophies are hereditary degenerative disorders of skeletal muscles, causing progressive replacement of muscle tissue by fat. This happens through disparate pathological processes and molecular mechanisms that are, at least to some extent, disease-specific and related to the peculiar genetic defect that characterizes each of them (1).

However, some of these broad pathological processes are shared by most muscular dystrophies. They can be followed on muscle imaging as an early phase of muscle damage and intramuscular edema, often corresponding to inflammatory/necrotic changes (2, 3), and subsequent stages characterized by progressive deposition of fat and connective tissue (4–6). Different pathological processes are generally present simultaneously in the same patient or even in the same muscle group.

The good soft tissue characterization capabilities of MRI can be exploited to quantify the status of these ongoing processes in the muscle; for example, global T2 contrast can be used to highlight edema, and fat/water separation methods can show fat infiltration.

When moving toward quantitative imaging, global T2, extracted by a monoexponential fitting is a sensitive disease indicator when fat infiltration is not present. However, in neuromuscular diseases where adipocytes significantly replaced muscular tissue, it highly correlates with the fat content of the muscle. In this case, it is therefore rather an indicator of long-term changes in the musculature, albeit an indirect one with respect to fat fraction quantification; conversely, the T2 relaxation of the water component of the tissue (wT2) well correlates with acute “disease activity” (7).

Various acquisition methods have been proposed to quantify wT2 independently of fat infiltration (8–15). In current clinical practice, multi-echo spin-echo (MESE) sequences are typically used, and various types of exponential fitting (biexponential, triexponential) are used for T2 calculation (16, 17). However, these methods are sensitive to multiple confounding factors, such as B1 inhomogeneities.

Marty et al. (14) presented a method based on an extended phase graph (EPG) fitting that addresses many of these issues. Besides, in contrast to other methods such as the ones proposed by Klupp et al. (12), Sousa et al. (11), or Koolstra et al. (13), the EPG-based approach has the advantage of using a conventional spin-echo sequence for the quantification, which is broadly available and therefore does not require any particular sequence modification or hardware. In recent years, this method has been extensively used for wT2 quantification in several studies (18–22). Although there are small differences in the exact protocol used depending on the MR scanner vendor, the differences in the implementation of the fitting process can be more substantial.

As a result of the EPG fitting, a relatively precise estimation of the fat fraction is also obtained. The resulting fat fractions approximate the results of a three-point Dixon acquisition, which has been proven to differentiate well between patients of various neuromuscular diseases and healthy subjects (16, 23). However, fat fraction estimation is not the primary quantitative target of this method, and it can be an issue when accurate fat fraction maps are needed, as weighting factors [such as magnetization transfer effects in multi-slice MESE acquisitions (24–26)] can bias the estimation. An alternative approach to obtain this information is to use a dedicated sequence (usually based on small flip-angle multi-echo gradient-echo imaging) for a purely proton-density-weighted fat fraction estimation. This approach also allows using more accurate fat models that incorporate the complex chemical composition of fat (27). In practice, a typical muscle MRI examination usually comprises a dedicated volume acquisition for single muscle identification (segmentation) and for accurate fat fraction acquisition.

In this work, we build upon this concept to present a fast and open software for wT2 fitting, which can be used as a reference implementation for reproducibility studies. The postprocessing performance is optimized by the usage of GPU processing and the creation of a dictionary (of adjustable size) of simulated signals incorporating slice profile information (14, 28). Additionally, this work further extends the initial concept as it implements multiple fitting methods for the wT2 estimation, and it can incorporate the information coming from a separate fat/water acquisition (of arbitrary resolution and field of view) to constrain the fitting. With this constraint, the possibility of reducing the number of acquired echoes for the fitting (thus reducing the scan time) is also analyzed.

Materials and Methods

Software Implementation

A stand-alone, command-line application was developed in Python, using standard mathematical extension libraries (NumPy, SciPy) and, for hardware acceleration, the pycuda extension for interfacing with the general-purpose graphical processing unit.

The core of the application consists of a fast, GPU-accelerated implementation of signal simulation based on the extended phase graph concept (29, 30). As the method is intended for multi-echo spin-echo acquisitions, a Carr-Purcell-Meiboom-Gill (CPMG) (31, 32) simulation is performed at every run of the program, adapting the timing and the number of echoes to the actual sequence parameters. The slice profile of the radiofrequency pulses is taken into account in the simulation, assuming hanning-windowed sinc pulses. The slice profile is calculated through the application of a Shinnar-LeRoux transform (33) and the slice width of the refocusing pulse is assumed to be a factor of 1.2 larger than the excitation pulse, as per characteristics of the pulse sequence of the used MR scanner. For different vendors and acquisition protocols, the user can provide custom slice profiles and refocusing width factors as parameters to the program. The simulation logic is directly written in GPU-specific C++ language for maximum efficiency. Thanks to the high parallelization of the GPU tasks, signals corresponding to multiple combinations of wT2, fat fraction, and B1 inhomogeneity factors can be simulated concurrently and placed in a dictionary, whose size can be chosen by the user. The fat model is assumed to have a single spectral peak with a fixed T2 value, which can either be given a priori or estimated from the input data. T1 values for both water and fat are assumed to be constant (1,400 and 365 ms, respectively), as in (14).

Subsequently, the time course of each voxel in the input series is compared to each entry in the dictionary by using the correlation metric. The parallelization of the GPU is again leveraged in this step, as the correlation operation is equivalent to the following matrix multiplication:

$\begin{array}{l}\text{C}=\text{S}\times \text{D}\end{array}$

where S is the signal matrix, having a separate voxel on each row and time in the column direction; D is the dictionary matrix, having time in the row direction and dictionary entries in the row directions; C is the result matrix holding the signal correlations with the dictionary. After multiplication, the index of the maximum value of each column identifies the parameter combination that best fits the measured signal.

The above procedure describes the unconstrained fitting. When an externally derived fat fraction map is also provided, it is aligned (based on the slice orientation and position) with each slice of the original spin-echo data. The information is then taken into account by restricting the maximum search to the given fat fraction for each voxel.

In addition to EPG fitting, the software can also perform double-exponential fitting. This fitting can neither correct for B1 inhomogeneity nor slice profile and is meant to be used when the EPG fitting is unstable or the true slice profile is not known.

The code is released under a free software license (GNU General Public License v3) at the website: https://www.github.com/fsantini/MyoQMRI.

Acquisition Protocol

A conventional 2D multi-slice multi-echo spin-echo (MESE) acquisition protocol was prepared for a commercial whole-body 3T MR scanner (MAGNETOM Skyra, Siemens Healthcare, Erlangen, Germany) equipped with an 18-channel body array coil and integrated spine coil with the following acquisition parameters: number of echoes 17, number of slices 7, TR 4,100 ms, first TE and echo spacing 10.9 ms, bandwidth 250 Hz/px, matrix size 192 × 384, resolution 1.2 × 1.2 mm², slice thickness 10 mm, gap between slices 30 mm.

For fat fraction quantification, a 3D multi-echo gradient-echo (MEGE) acquisition using a custom sequence was prepared with the following parameters: number of echoes 6, TR 35 ms, first TE/echo spacing 1.7/1.5 ms, flip angle 7°, bandwidth 1,050 Hz/px, matrix size 396 × 432 × 52, resolution 1.0 × 1.0 × 5.0 mm3. The sequence had a monopolar readout with interleaved echo spacing (even and odd echoes acquired in subsequent repetitions). The acquired images were postprocessed with the publicly-available algorithm FattyRiot (27) to obtain the fat fraction maps.

Human Experiments

The acquisitions described above were performed on the thighs of 5 (three male, median age 46 y, range 29–61 y) subjects with diagnosed facioscapulohumeral muscular dystrophy (FSHD), which is a particular form of muscular dystrophy characterized by progressive muscle wasting and fatty replacement often preceded by muscle edema on MRI (34), and of four subjects without a history of neuromuscular diseases (two male, median age 54.5 y, range 26–72 y). The acquisitions were performed according to the local ethics regulations and informed consent was obtained from the participants.

Data Analysis

The computer used for the postprocessing is a current mid-range personal computer (Ryzen 2600, Advanced Micro Devices, Santa Clara, CA, equipped with 32GB of RAM and a GeForce 1060 GPU, NVIDIA corporation, Santa Clara, CA).

The proposed algorithm was applied to the MESE acquisition after generation of a dictionary containing 60 linearly spaced values for wT2 (range 20–80 ms), 20 values for B1 factor (range 40–140%), and 101 values for the fat fraction (range 0–100%), for a total of 121,200 parameter combinations. The fat T2 was estimated from a subsample of the subjects (through a single-component EPG fitting in regions of subcutaneous fat) and subsequently assumed constant at 151 ms. Maps derived from EPG and double exponential fitting, with and without the constraint of the external proton-density-weighted fat fraction as calculated from the MEGE acquisition, were produced (an overview of the assumed and fitted parameters for each method is given in Table 1). The fitting was repeated by discarding later spin echoes (and only retaining 8, 5, or 3) to evaluate the robustness of the algorithm with respect to the number of echoes.

TABLE 1

Table 1. Overview of the presented fitting methods.

Regions of interest (ROIs) were manually drawn by one reader bilaterally on every MESE slice over the cross section of the following muscles: Vastus Lateralis (VM), Vastus Medialis (VM), Rectus Femoris (RF).

The average and standard deviation for wT2 values and fat fraction were extracted, and the following indicators were calculated:

• Average error between fat fraction calculated from MESE and the one deriving from MEGE (considered the “gold standard”) - only for nonconstrained reconstructions, calculated as:

$\begin{array}{l}ER{R}_{ff}=\underset{pop,subj,roi}{\text{mean}}(F{F}_{MESE})-\underset{pop,subj,roi}{\text{mean}}(F{F}_{MEGE} ),\end{array}$

where $\underset{pop,subj,roi}{\text{mean}}$ represents the averaging operation over the whole population, over all ROIs of each subject, and over the voxels of each ROI.

• Pooled standard deviation across the ROIs—this indicator is similar to an “average of the standard deviations,” and it indicates the variability of the fitted wT2 values over a small spatial region; it is related to noise or tissue inhomogeneities due to disease activity and is calculated as follows:

$\begin{array}{l}S{D}_{pooled}=\sqrt{\underset{pop,subj}{\text{mean}}\left(\underset{roi}{{\text{sd}}^{\text{2}}}(T_2) \right)},\end{array}$

where $\underset{roi}{{\text{sd}}^{\text{2}}}$ represents the squared standard deviation over each ROI.

• Standard deviation of the averages over the ROIs (hereby termed “intrasubject standard deviation”) - this indicator relates to the variability over larger areas, for example arising from field inhomogeneity, and is calculated as follows:

$\begin{array}{l}S{D}_{intrasubject}=\sqrt{\underset{pop}{\text{mean}}\left(\underset{subj}{{\text{sd}}^{\text{2}}}\left(\underset{roi}{\text{mean}}(T_2)\right)\right)} \end{array}$

This indicator was only calculated in the volunteer cohort, because patients might have variability in the wT2 values due to their pathology.

• Global average wT2 for patients and control groups, calculated as follows:

$\begin{array}{l}{\text{mean}}_{T2,(patients|controls)}=\underset{patients|controls}{\text{mean}}\left(\underset{subj,roi}{\text{mean}}(T_2)\right)\\ {\text{sd}}_{T2,(patients|controls)}=\underset{patients|controls}{\text{sd}}\left(\underset{subj,roi}{\text{mean}}(T_2)\right)\end{array}$

Statistical analysis was performed with R (35). To quantify the influence of the fat fraction on the fitting of wT2, the correlation between the average fat fraction and the average wT2 for each ROI was calculated. For the calculation of the correlation coefficients the log-transformed fat fractions were used to compensate for the skewness of the fat fraction distribution. However, the visualization was done on the original fat fraction axis (0–100%) for ease of interpretation.

In addition, correlation of repeated measurements (rmcorr function) instead of the standard Pearson coefficient was used to compensate for dependence of measurements of the same subject. In this case the subject is introduced as a variable and both wT2 and fat fraction are considered as measures.

Validation

In order to evaluate the absolute accuracy of the fitting method, the same data was analyzed using the EPG wT2 fitting procedure of the QMRTools software package (36). The slice profile was adapted to match the one used in the current acquisitions and the fitting was applied to the full 17-echo dataset.

The same ROIs as in the previous analysis were considered, and the average wT2 in each ROI was calculated and compared to the corresponding ROI of our method. Average and standard deviation of the errors were calculated, and the agreement was visualized in a Bland-Altman plot.

Results

The software fitted the MESE datasets with EPG dictionary matching in an average time of 2 s/slice, plus 12 s per dataset for the generation of the dictionary (operation performed at every run of the program, which was faster than caching the results to disk). The double exponential fitting was not GPU-optimized and took 390 s (6 m 30 s) per slice.

Exemplary outputs from an unconstrained 17-echo reconstruction for a patient and a volunteer are shown in Figure 1. No noticeable artifacts are visible in the wT2 and fat fraction maps, whereas the B1 maps show some inconsistencies, mostly located in the fascia, thus remaining mostly masked in the other maps.

FIGURE 1

Figure 1. Parameter maps derived from an unconstrained 17-echo EPG dictionary matching for an FSHD patient (top row) and a healthy volunteer (bottom row). The rightmost panel is the registered fat fraction map deriving from a multi echo gradient echo acquisition (MEGE) for reference.

WT2 maps had a homogeneous appearance with EPG matching both for 17 and 8 echoes (Figure 2). However, the unconstrained reconstruction showed a noticeable difference in quantitative values in areas with heavy fat infiltration, whereas the fat-fraction-constrained reconstructions appeared similar between 17, 8, and 5 considered echoes. The visual quality of the map was insufficient in any combination when only three echoes were used for the reconstruction (Figure 2). Quantitatively, it can be observed that a reduced number of echoes in the unconstrained fitting has an effect on wT2 fitting when fat infiltration is present, resulting in a highly significant positive correlation (p < 0.001) for EPG fitting with 5 and 8 echoes (correlation coefficient r = 0.73 and 0.61, respectively, Figure 3). Conversely, the correlation was significantly negative (r = −0.71, p < 0.001) for the unconstrained fitting with 17 echoes, and for the constrained fitting of 5 and 8, but not 17 echoes (r = −0.60, −0.41, −0.19; p < 0.001, p = 0.005, 0.2, respectively).

FIGURE 2

Figure 2. wT2 maps for an FSHD patient using a varying number of acquired echoes by EPG simulation with and without an external fat fraction map as a reconstruction constraint.

FIGURE 3

Figure 3. wT2 values over each ROI vs. fat fraction in the same ROI, for all ROIs drawn in patients and healthy volunteers. The solid lines represent the linear regression with shaded gray areas indicating the 95% confidence intervals.

The double exponential fitting showed decreased wT2 with any number of echoes when unconstrained fitting was used, in addition to artifacts arising from B1 inhomogeneities (see Figure 4 and Table 2). Correlation between wT2 and the fat fraction was always negative and significant (p = 0.02 for the constrained 8 and 17 echoes, p < 0.001 otherwise) in all cases except, notably, the unconstrained fitting of five echoes (p = 0.29), with negative correlation coefficients ranging from r = −0.34 for the constrained 8-echo reconstruction to r = −0.62 for the unconstrained 8-echo reconstruction (Figure 3).

FIGURE 4

Figure 4. wT2 maps for an FSHD patient using a varying number of acquired echoes by double exponential fitting with and without an external fat fraction map as a reconstruction constraint.

TABLE 2

Table 2. Summary of results for the different fitting methods, with and without an external proton-density-weighted fat fraction (FF) constraint.

According to the quantitative quality metrics (Table 2), considering a higher number of echoes improved both the noise (pooled standard deviation) and the homogeneity across the ROIs for each subject (intrasubject variability). For a reduced number of echoes, the constrained reconstruction retained good homogeneity down to five considered echoes for the EPG reconstruction (3.1 ms for the constrained case vs. 6.7 for the unconstrained). The double exponential fitting showed good homogeneity (intrasubject standard deviation ranging from 2.7 to 4.6 milliseconds), but more noise (pooled standard deviation ranging from 5.9 to 10.5 ms).

The constrained reconstruction resulted in a generally higher pooled standard deviation due to the artifacts introduced by the misregistration of the images.

The accuracy of the wT2 values was close to the QMRTools implementation, showing an average difference of 0.9 ms (5.2 ms standard deviation, Figure 5).

FIGURE 5

Figure 5. Bland Altman plot of the agreement of the wT2 values obtained by the QMRTools software package and the presented software, for an unconstrained 17-echo EPG reconstruction. The bias and the 95% confidence intervals are depicted in the plot.

Concerning the fat fraction accuracy, although it is not the primary focus of this method, the EPG matching resulted in consistent underestimation compared to the FattyRiot algorithm, with a larger number of echoes providing results which were closer in average to the gold standard (from −4.2 percentage points (p.p.) for 17 echoes to -7.8 p.p. for eight echoes). The double exponential fitting provided an overestimation ranging from +3.2 p.p. for 17 echoes to +8.6 p.p. for five echoes.

Discussion

In this work, we presented a software application that can quickly and reliably calculate wT2 maps in the presence of spectrally inhomogeneous voxels containing both non-fatty tissue and lipids, while at the same time estimating fat fraction from conventional multi-echo spin-echo acquisitions. This work follows the concept introduced by Marty et al. (14). Still, it additionally provides the possibilities of performing double exponential (and, optionally, single exponential) fitting and, more interestingly, of incorporating external fat fraction information to improve the accuracy of the fitting. There are no specific hardware requirements for this program; however, for better performance, a reasonable Cuda-compatible graphic card should be used. Consumer-grade GPUs deliver good performance, but more extensive dictionaries than the one tested in this current setup require increased memory on the device. The code is released with a free open source license and, in contrast to existing available implementations (36), this software package is exclusively based on free software and is platform-independent. This implementation can thus be considered ready to be widely used in a clinical research context and as a reference by other researchers.

In the comparison with an existing implementation, the absolute wT2 values obtained by this method are close, but not identical, to the ones obtained from QMRTools, although both methods are based on the same conceptual framework. One explanation could be the usage of different metrics while performing the dictionary matching. This finding highlights the necessity of having consistent acquisition and data processing pipelines, and the necessity of characterizing a quantitative method primarily in terms of reproducibility and precision. It is generally true also in other fields of quantitative MRI, that the absolute values obtained are method-dependent and thus comparisons need to be carefully considered (37, 38).

During the optimization of the EPG matching algorithm, it appeared clear that an accurate slice profile was very important to obtain absolute values of wT2 close to the literature. Small changes to the profile, or selecting a too large refocusing slice width factor, could introduce a bias of a few milliseconds in the estimated values. The intersubject and intrasubject variabilities, however, minimally changed, so the user should be advised to obtain the exact sequence characteristics and to use a coherent parameter set when reconstructing multiple datasets.

The introduction of the constrained reconstruction appears beneficial for the fitting of wT2, providing reliable and homogeneous results, not correlated with the amount of infiltration of the muscle, with considering as few as five acquired echoes. This is a valuable result in light of providing better spatial coverage by the MESE sequence: reducing the number of acquired echoes allows exciting multiple interleaved slices in a simple repetition time and thus lessen the interslice gaps. However, the alignment of the fat fraction map can introduce some artifacts if patient motion occurs between the two acquisitions. In the current implementation, no image registration is performed. The alignment of the MESE and MEGE datasets is currently only performed based on the orientation information provided in the image headers. Image registration is, however, implementable using free python libraries and could be added if needed.

The EPG matching, as expected, could account for the B1 inhomogeneities and therefore lead to lower artifacts in the areas where flip angles deviate from the nominal value in comparison to double-exponential fitting.

In general, unconstrained double exponential fitting with reduced echoes resulted in underestimating wT2 values relative to other methods with unconstrained reconstruction. On the other hand, thanks to its fewer degrees of freedom in the fitting, it showed remarkable homogeneity across the various ROIs even with few echoes, suggesting that it might still be considered a robust and straightforward approach when no other methods are available. The constrained double exponential reconstruction, conversely, produced results farthest from the expected values (12, 14). An analysis of the data shows that the unconstrained double exponential fitting consistently overestimates the fat fraction. The possible explanation is that the longer exponential decay due to the stimulated echoes gets assigned to the fat component during the fitting process in the unconstrained case, whereas fixing the fat fraction in the constrained reconstruction produces a longer apparent exponential decay constant in what is practically a monoexponential fitting of the residuals.

Generally, most of the obtained wT2 values negatively correlate with fat fraction, which is in line with the previous findings (15, 39), with the exception of the EPG matching of the reduced echoes, showing a positive correlation, suggesting a failure to separate the wT2 from fat.

As an overall comparison, both fitting methods (double exponential fitting and EPG matching) have advantages and disadvantages. EPG matching appears accurate and precise even when a lower number of echoes is used, especially when paired with an external fat fraction constraint, and it has a higher insensitivity to B1 inhomogeneity, but it requires precise knowledge of the acquisition parameters to obtain unbiased values. Double exponential fitting fails in regions of B1 inhomogeneity, its accuracy is poor when few echoes are used, and the external fat fraction constraints introduces a bias in the obtained values. However, it requires very little knowledge of the acquisition parameters, and it can therefore be chosen when the sequence characteristics are unknown. When the full 17-echo acquisitions are used, the two fitting methods have similar characteristics when averaged across the muscles; however, the sensitivity of the double exponential fitting to B1 inhomogeneity might mask local intramuscular changes in patients with neuromuscular disorders.

Although not explored in the present validation for the results to be more comparable, the program assumes a constant value for fat T2. Still, it gives the possibility of indicating it as a runtime parameter or estimating it from the image itself. Similarly, T1 values for water and fat are fixed, although this can be assumed to have a small effect on the final result (14). Other effects like j-coupling are also not explicitly introduced in the model but rather incorporated in the assumed value of fat T2, which is dependent on the characteristics of the MESE sequence (40).

The implementation has some further limitations regarding the accuracy of the results. Specifically, Keene et al. demonstrated that this method can be improved by introducing a correction for the chemical shift in the slice profile and a better estimation for fat T2 (15). These corrections are relatively newly introduced and not routinely used in the current studies, and thus not currently implemented. However, the chemical shift correction requires knowledge of the actual implementation of the pulse sequence that might be difficult to obtain in a clinical setting. It would therefore result in potential loss of generalizability. Similarly, a multi-peak spectral model for the fat could be incorporated into the algorithm, which could improve the results' accuracy. As this method is based on multi spin echoes, the multiple spectral peaks do not result in different chemical shifts as in gradient-echo images. The effect would only be seen in the different T2 values of the fat components; therefore, this functionality is not currently implemented, in line with existing spin-echo-based methods for muscle imaging (14, 15, 17, 23).

One limitation of this work is the relatively small number of datasets. However, both healthy and different diseased subjects were included. While a larger subject cohort would be of interest, one of the goals of this study was to offer the tools for such studies and not to draw conclusions on wT2 measurements of dystrophic muscles. For this goal, clinical studies that focus on more homogenous patient cohorts (grouping, for example, the different stages of neuromuscular disease, where similar pathological changes are expected) are required.

The validation of this work was performed by comparing the proposed fitting with an existing implementation of the same concept on the same data, thus lacking an external reference standard for the values obtained. While such an external reference could be desirable, the scope of this work was not to assess the validity of MESE acquisitions for the estimation of wT2, but only to evaluate the efficacy of the proposed implementation. The availability of an independent fitting procedure allowed a direct comparison of these characteristics without other sources of error deriving from the physics of different acquisition methods.

In conclusion, in this work, we have presented a fast and open implementation of an algorithm for the T2 mapping of the water component of muscle tissues in the presence of fat, based on conventional multi-echo spin-echo sequences, capable of incorporating prior knowledge of the fat fraction. Thanks to the additional information obtained from the multi-echo gradient-echo images, a reduced number of echoes can be used for the spin-echo acquisition with while retaining similar inter- and intrasubject variability and similar absolute values, when and EPG model is used. The EPG matching method is in general to be preferred to a double exponential fitting, provided that the characteristics of the MR sequence (especially in terms of RF pulses) are sufficiently known.

Data Availability Statement

The data analyzed in this study is subject to the following licenses/restrictions: Data included in this work can be shared upon research agreement between the Institutions. Requests to access these datasets should be directed to Anna Pichiecchio, YW5uYS5waWNoaWVjY2hpbyYjeDAwMDQwO21vbmRpbm8uaXQ=.

Ethics Statement

The studies involving human participants were reviewed and approved by Comitato Etico Area Referente Pavia Fondazione IRCCS Policlinico San Matteo and Policlinico Gemelli's Ethics Committe. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

FSa study conception, method and software development, data analysis, and manuscript drafting. XD data analysis, statistics, and manuscript drafting. MP, FSo, and AP data acquisition and patient management. MW and PS method development. MM, ER, and GT patient recruitment and clinical evaluations. OB, AP, and NB study conception and supervision. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the Swiss National Science Foundation (SNSF) grant no. 172876 and by the Italian Ministry of Health (RC 2017-2019 and RF-2016-02362914).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

1. Mercuri E, Muntoni F. Muscular dystrophies. Lancet Lond Engl. (2013) 381:845–60. doi: 10.1016/S0140-6736(12)61897-2

CrossRef Full Text | Google Scholar

2. Tasca G, Pescatori M, Monforte M, Mirabella M, Iannaccone E, Frusciante R, et al. Different molecular signatures in magnetic resonance imaging-staged facioscapulohumeral muscular dystrophy muscles. PLoS ONE. (2012) 7:e38779. doi: 10.1371/journal.pone.0038779

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Lassche S, Küsters B, Heerschap A, Schyns MVP, Ottenheijm CAC, Voermans NC, et al. Correlation between quantitative MRI and muscle histopathology in muscle biopsies from healthy controls and patients with IBM, FSHD and OPMD. J Neuromuscul Dis. (2020) 7:495–504. doi: 10.3233/JND-200543

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Monforte M, Laschena F, Ottaviani P, Bagnato MR, Pichiecchio A, Tasca G, et al. Tracking muscle wasting and disease activity in facioscapulohumeral muscular dystrophy by qualitative longitudinal imaging. J Cachexia Sarcopenia Muscle. (2019) 10:1258–65. doi: 10.1002/jcsm.12473

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Klingler W, Jurkat-Rott K, Lehmann-Horn F, Schleip R. The role of fibrosis in Duchenne muscular dystrophy. Acta Myol Myopathies Cardiomyopathies. (2012) 31:184–95.

PubMed Abstract | Google Scholar

6. Bonati U, Hafner P, Schädelin S, Schmid M, Naduvilekoot Devasia A, Schroeder J, et al. Quantitative muscle MRI: a powerful surrogate outcome measure in Duchenne muscular dystrophy. Neuromuscul Disord. (2015) 25:679–85. doi: 10.1016/j.nmd.2015.05.006

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Carlier PG. Global T2 versus water T2 in NMR imaging of fatty infiltrated muscles: different methodology, different information and different implications. Neuromuscul Disord. (2014) 24:390–2. doi: 10.1016/j.nmd.2014.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Janiczek RL, Gambarota G, Sinclair CDJ, Yousry TA, Thornton JS, Golay X, et al. Simultaneous T2 and lipid quantitation using IDEAL-CPMG. Magn Reson Med. (2011) 66:1293–302. doi: 10.1002/mrm.22916

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Mankodi A, Bishop CA, Auh S, Newbould RD, Fischbeck KH, Janiczek RL. Quantifying disease activity in fatty-infiltrated skeletal muscle by IDEAL-CPMG in Duchenne muscular dystrophy. Neuromuscul Disord. (2016) 26:650–8. doi: 10.1016/j.nmd.2016.07.013

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Schlaeger S, Weidlich D, Klupp E, Montagnese F, Deschauer M, Schoser B, et al. Water T2 mapping in fatty infiltrated thigh muscles of patients with neuromuscular diseases using a T2-prepared 3D turbo spin echo with SPAIR. J Magn Reson Imaging. (2020) 51:1727–36. doi: 10.1002/jmri.27032

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Sousa PL de, Vignaud A, Araújo EC de A, Carlier PG. Factors controlling T2 mapping from partially spoiled SSFP sequence: optimization for skeletal muscle characterization. Magn Reson Med. (2012) 67:1379–90. doi: 10.1002/mrm.23131

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Klupp E, Weidlich D, Schlaeger S, Baum T, Cervantes B, Deschauer M, et al. B1-insensitive T2 mapping of healthy thigh muscles using a T2-prepared 3D TSE sequence. PLoS ONE. (2017) 12:e0171337. doi: 10.1371/journal.pone.0171337

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Koolstra K, Webb AG, Veeger TTJ, Kan HE, Koken P, Börnert P. Water–fat separation in spiral magnetic resonance fingerprinting for high temporal resolution tissue relaxation time quantification in muscle. Magn Reson Med. (2020) 84:646–62. doi: 10.1002/mrm.28143

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Marty B, Baudin P-Y, Reyngoudt H, Azzabou N, Araujo ECA, Carlier PG, et al. Simultaneous muscle water T2 and fat fraction mapping using transverse relaxometry with stimulated echo compensation. NMR Biomed. (2016) 29:431–43. doi: 10.1002/nbm.3459

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Keene KR, Beenakker J-WM, Hooijmans MT, Naarding KJ, Niks EH, Otto LAM, et al. T2 relaxation-time mapping in healthy and diseased skeletal muscle using extended phase graph algorithms. Magn Reson Med. (2020) 84:2656–70. doi: 10.1002/mrm.28290

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Strijkers GJ, Araujo ECA, Azzabou N, Bendahan D, Blamire A, Burakiewicz J, et al. Exploration of new contrasts, targets, and MR imaging and spectroscopy techniques for neuromuscular disease – A Workshop Report of Working Group 3 of the Biomedicine and Molecular Biosciences COST Action BM1304 MYO-MRI. J Neuromuscul Dis. (2019) 6:1–30. doi: 10.3233/JND-180333

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Kan HE, Scheenen TWJ, Wohlgemuth M, Klomp DWJ, van Loosbroek-Wagenmans I, Padberg GW, et al. Quantitative MR imaging of individual muscle involvement in facioscapulohumeral muscular dystrophy. Neuromuscul Disord. (2009) 19:357–62. doi: 10.1016/j.nmd.2009.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Schlaffke L, Rehmann R, Rohm M, Otto LAM, de Luca A, Burakiewicz J, et al. Multi-center evaluation of stability and reproducibility of quantitative MRI measures in healthy calf muscles. NMR Biomed. (2019) 32:e4119. doi: 10.1002/nbm.4119

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Hooijmans MT, Baligand C, Froeling M, Verschuuren JJGM, Webb AG, Niks EH, et al. Multi-parametric MR shows increased T2 heterogeneity in fat infiltrated muscles in Becker Muscular Dystrophy. In: Proc Intl Soc Mag Reson Med 26. Paris (2018). Available from: http://indexsmart.mirasmart.com/ISMRM2018/PDFfiles/0816.html (accessed January 15, 2021).

Google Scholar

20. Otto LAM, Pol W-L van der, Schlaffke L, Wijngaarde CA, Stam M, Wadman RI, et al. Quantitative MRI of skeletal muscle in a cross-sectional cohort of patients with spinal muscular atrophy types 2 and 3. NMR Biomed. (2020) 33:e4357. doi: 10.1002/nbm.4357

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Heskamp L, Okkersen K, van Nimwegen M, Ploegmakers MJ, Bassez G, Deux J-F, et al. Quantitative muscle MRI depicts increased muscle mass after a behavioral change in myotonic dystrophy type 1. Radiology. (2020) 297:132–42. doi: 10.1148/radiol.2020192518

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Keene KR, Vught L, Velde NM, Ciggaar IA, Notting IC, Genders SW, et al. The feasibility of quantitative MRI of extra-ocular muscles in myasthenia gravis and Graves' orbitopathy. NMR Biomed. (2021) 34:e4407. doi: 10.1002/nbm.4407

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Carlier PG, Marty B, Scheidegger O, Loureiro de Sousa P, Baudin P-Y, Snezhko E, et al. Skeletal muscle quantitative nuclear magnetic resonance imaging and spectroscopy as an outcome measure for clinical trials. J Neuromuscul Dis. (2016) 3:1–28. doi: 10.3233/JND-160145

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Radunsky D, Blumenfeld-Katzir T, Volovyk O, Tal A, Barazany D, Tsarfaty G, et al. Analysis of magnetization transfer (MT) influence on quantitative mapping of T2 relaxation time. Magn Reson Med. (2019) 82:145–58. doi: 10.1002/mrm.27704

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Weigel M, Helms G, Hennig J. Investigation and modeling of magnetization transfer effects in two-dimensional multislice turbo spin echo sequences with low constant or variable flip angles at 3 T. Magn Reson Med. (2010) 63:230–4. doi: 10.1002/mrm.22145

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Melki PS, Mulkern RV. Magnetization transfer effects in multislice RARE sequences. Magn Reson Med. (1992) 24:189–95. doi: 10.1002/mrm.1910240122

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Welch EB, Smith DS, Avison MJ, Berglund J, Kullberg J, Ahlström H. Fattyriot - Final Winning Entry Of The 2012 Ismrm Challenge On Water-Fat Reconstruction. Zenodo (2015). Available from: https://zenodo.org/record/16741 (accessed October 9, 2020).

Google Scholar

28. Huang C, Altbach MI, El Fakhri G. Pattern recognition for rapid T2 mapping with stimulated echo compensation. Magn Reson Imaging. (2014) 32:969–74. doi: 10.1016/j.mri.2014.04.014

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Weigel M. Extended phase graphs: dephasing, RF pulses, and echoes - pure and simple. J Magn Reson Imaging. (2015) 41:266–95. doi: 10.1002/jmri.24619

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Hennig J. Echoes—how to generate, recognize, use or avoid them in MR-imaging sequences. Part I: fundamental and not so fundamental properties of spin echoes. Concepts Magn Reson. (1991) 3:125–43. doi: 10.1002/cmr.1820030302

CrossRef Full Text | Google Scholar

31. Carr HY, Purcell EM. Effects of diffusion on free precession in nuclear magnetic resonance experiments. Phys Rev. (1954) 94:630–8. doi: 10.1103/PhysRev.94.630

CrossRef Full Text | Google Scholar

32. Meiboom S, Gill D. Modified spin-echo method for measuring nuclear relaxation times. Rev Sci Instrum. (1958) 29:688–91. doi: 10.1063/1.1716296

CrossRef Full Text | Google Scholar

33. Pauly J, Roux PL, Nishimura D, Macovski A. Parameter relations for the Shinnar-Le Roux selective excitation pulse design algorithm (NMR imaging). IEEE Trans Med Imaging. (1991) 10:53–65. doi: 10.1109/42.75611

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Tasca G, Monforte M, Ottaviani P, Pelliccioni M, Frusciante R, Laschena F, et al. Magnetic resonance imaging in a large cohort of facioscapulohumeral muscular dystrophy patients: pattern refinement and implications for clinical trials. Ann Neurol. (2016) 79:854–64. doi: 10.1002/ana.24640

PubMed Abstract | CrossRef Full Text | Google Scholar

35. R Development Core Team. R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing (2008). Available from: http://www.R-project.org (accessed February 15, 2021).

Google Scholar

36. Froeling M. QMRTools: a mathematica toolbox for quantitative MRI analysis. J Open Source Softw. (2019) 4:1204. doi: 10.21105/joss.01204

CrossRef Full Text | Google Scholar

37. Roujol S, Weingärtner S, Foppa M, Chow K, Kawaji K, Ngo LH, et al. Accuracy, precision, and reproducibility of four T1 mapping sequences: a head-to-head comparison of MOLLI. ShMOLLI, SASHA, and SAPPHIRE. Radiology. (2014) 272:683–9. doi: 10.1148/radiol.14140296

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Baeßler B, Schaarschmidt F, Stehning C, Schnackenburg B, Giolda A, Maintz D, et al. Reproducibility of three different cardiac T2-mapping sequences at 1.5T. J Magn Reson Imaging. (2016) 44:1168–78. doi: 10.1002/jmri.25258

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Schlaeger S, Weidlich D, Klupp E, Montagnese F, Deschauer M, Schoser B, et al. Decreased water T2 in fatty infiltrated skeletal muscles of patients with neuromuscular diseases. NMR Biomed. (2019) 32:e4111. doi: 10.1002/nbm.4111

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Hardy PA, Henkelman RM, Bishop JE, Poon ECS, Plewes DB. Why fat is bright in rare and fast spin-echo imaging. J Magn Reson Imaging. (1992) 2:533–40. doi: 10.1002/jmri.1880020511

PubMed Abstract | CrossRef Full Text | Google Scholar