Association of Cervical and Lumbar Paraspinal Muscle Composition Using Texture Analysis of MR-Based Proton Density Fat Fraction Maps

In this study, the associations of cervical and lumbar paraspinal musculature based on a texture analysis of proton density fat fraction (PDFF) maps were investigated to identify gender- and anatomical location-specific structural patterns. Seventy-nine volunteers (25 men, 54 women) participated in the present study (mean age ± standard deviation: men: 43.7 ± 24.6 years; women: 37.1 ± 14.0 years). Using manual segmentations of the PDFF maps, texture analysis was performed and texture features were extracted. A significant difference in the mean PDFF between men and women was observed in the erector spinae muscle (p < 0.0001), whereas the mean PDFF did not significantly differ in the cervical musculature and the psoas muscle (p > 0.05 each). Among others, Variance(global) and Kurtosis(global) showed significantly higher values in men than in women in all included muscle groups (p < 0.001). Not only the mean PDFF values (p < 0.001) but also Variance(global) (p < 0.001), Energy (p < 0.001), Entropy (p = 0.01), Homogeneity (p < 0.001), and Correlation (p = 0.037) differed significantly between the three muscle compartments. The cervical and lumbar paraspinal musculature composition seems to be gender-specific and has anatomical location-specific structural patterns.


Introduction
Muscle structure and composition change over the human life span. Not only increasing body mass index (BMI) and age are associated with structural alterations, including fatty infiltration and atrophy, but also katabolic conditions, such as sarcopenia and cachexia [1][2][3][4][5].
There is evidence for connections between the overall muscle quality and muscular performance [6]. Further, the loss of functional muscle tissue due to fatty replacement leads to increasing disability and mortality [7]. Recently, it has also been shown that preoperative sarcopenia is associated with poorer overall survival in pancreatic cancer patients following pancreaticoduodenectomy [8]. Thus, non-invasive muscle status assessments beyond quantification of mere volumetric changes are of increasing interest with regard to personalized medicine and an individual, multi-parametric patient evaluation.
In clinical practice, established imaging methods such as dual-energy x-ray absorptiometry (DXA) and computed tomography (CT) are used to assess muscle status [9]. In certain settings, more advanced imaging techniques such as single-voxel proton magnetic resonance spectroscopy (MRS) and chemical shift encoding-based water-fat magnetic resonance imaging (CSE-MRI) are important diagnostic tools to extract biomarkers such as the proton density fat fraction (PDFF) or analyze the chemical structure of fatty acids within different tissue types [10]. There are plenty of data to support the role of water-fat MRI in the context of muscle composition analysis in healthy adults [11][12][13][14][15][16][17]. Recently, studies using texture analysis (TA) of PDFF maps derived from CSE-MRI showed promising results with regard to structural changes of vertebral bone marrow and paravertebral musculature with good functional correlation [18,19]. As demonstrated in these studies, TA can generate novel information which captures structural changes in osseous compartments and musculature beyond the mean PDFF [18,19].
However, there is still a lack of normative studies investigating the structural composition and morphological association of large muscle groups. The paravertebral musculature exhibits age-and gender-specific degeneration patterns, which are relevant in the context of musculoskeletal disorders [2,20,21]. The comparison of structural muscular changes with regard to gender and anatomical location has the potential to contribute valuable information on muscle (patho) physiology.
Therefore, this study investigated the associations of different cervical and lumbar paravertebral muscle groups based on TA of PDFF maps to identify gender-and regionspecific structural patterns.

Subjects
Seventy-nine volunteers with self-reported absence of any musculoskeletal or metabolic diseases (male = 25, female = 54) were recruited from a large cohort recruited for evaluation of determinants of basal metabolic rate [22]. The study protocol and procedures were approved by the local ethics committee of the Faculty of Medicine. Inclusion criteria were: age of 18 years or older, no history of severe diseases or surgery within the last three months, no acute physical impairment. Exclusion criteria were standard contraindications for MRI (implanted pacemaker or MRI-incompatible implants or devices). All subjects provided written informed consent prior to study inclusion.

Magnetic Resonance Imaging
All subjects underwent MRI at 3 Tesla (Ingenia, Philips Healthcare, Best, Netherlands) using the built-in 12-channel posterior coil and a 16-channel anterior coil (dStream Torso coil, Philips Healthcare, Best, Netherlands). Subjects were positioned head-first in supine position.
The imaging protocol comprised an axially prescribed, six-echo three-dimensional (3D) spoiled gradient echo sequence in two stacks for chemical shift encoding-based water-fat separation at the cervical and lumbar spine, respectively.

Muscle Fat Quantification
The gradient echo imaging data were processed online using the fat quantification routine of the MRI vendor (Philips Healthcare, Best, The Netherlands). After phase error correction, complex-based water-fat decomposition was performed using a precalibrated seven-peak fat spectrum model with a single T2 * [23,24]. The PDFF maps were generated voxel by voxel by computing the ratio of the fat signal and the sum of fat and water signals. Mean PDFF values were computed by averaging the PDFF maps over the ROIs of the segmented muscle groups. Segmentations were performed by a radiologist using the free open-source software Medical Imaging Interaction Toolkit (MITK, developed by the Division of Medical and Biological Informatics, German Cancer Research Center, Heidelberg, Germany; www.mitk.org, accessed on 18 October 2021).
The cervical musculature (CE), the erector spinae muscle (ES) and the psoas muscle (PS) were manually segmented bilaterally in the PDFF maps at the level of C5 and L4, respectively ( Figure 1). As a standard, 10 axial slices were segmented for each muscle group. Reproducibility of the segmentation method was described in previous publications. The intraclass correlation coefficient (ICC) for intra-reader and inter-reader reliability was reported to be 0.966 and 0.942 at C5 level, respectively [2]. Further, an excellent reproducibility was reported for PDFF measurements at L5 level, with a root mean square absolute precision error of 0.48% [21].

Texture Analysis
Based on the gray-level distribution in an image, TA can be utilized to characterize structural properties of segmentations via the quantification of different texture features (TFs) [25][26][27]. TFs were calculated for each of the segmented muscle groups, using the PDFF maps. Extracted TFs included 3 global features (Variance(global), Skewness(global), and Kurtosis(global)) and 8 second-order features (Energy, Contrast, Entropy, Homogeneity and Correlation, calculated as described in [28]; Variance and Sum-average [29]; Dissimilarity [30]). Equations for the computed TFs can be found in the Appendix A. TFs were averaged over both sides, weighted by muscle volume, to extract bilateral TF values (e.g., Variance (global) CE , Variance (global) ES , and Variance (global) PS ).
Intensity histogram analysis was applied for the calculation of global features. The optimal histogram bin size and number is a point of discussion, and strongly depending on data characteristics and purpose of the histogram analysis. In our study, we chose to calculate the number of bins by taking the median of the three following methods: Sturges' method, Scott's method, and the Freedman-Diaconis method. Using this approach, we obtained the most reasonable results compared to visual inspection of the histograms and best representation of the relevant data characteristics [31][32][33].
Gray-level co-occurrence matrix (GLCM) analysis was applied for the calculation of second-order TFs [28]. The preprocessing included gray-level quantization of the PDFF maps to prevent sparseness by normalizing the image intensities. For this purpose, we used 200 equally sized bins and the minimum and maximum gray levels, which correspond to values of 0% and 100%, respectively.
The entries of the GLCMs at different angular directions θ = (0 • , 45 • , 90 • , and 135 • ) were generated by computing the joint probability of two adjacent voxel intensities at a given offset d = (dx, dy, dz) and given θ, with dx, dy and dz denoting the displacement along the x-, yand z-axis, respectively.
Three-dimensional GLCM analysis was conducted by computation of the co-occurrence probabilities of voxel intensities from the 26 neighbors, which are aligned in 13 directions. In this process, discretization length differences were taken into account and adjusted for. Averaging the 13 directions ensures rotation invariance. All described preprocessing steps (isotropic resampling and gray-level uniform quantization) as well as the actual TA were performed with MATLAB 2021a (MathWorks Inc., Natick, MA, USA) using a radiomics toolbox (https://github.com/mvallieres/radiomics/, accessed on 18 October 2021) [34][35][36].

Statistical Analysis
For the statistical analyses, SPSS (version 20.0; IBM SPSS Statistics for Windows, Armonk, NY, USA) was used. Statistical significance was considered at p < 0.05 (two-sided) in all conducted tests.
The Kolmogorov-Smirnov test indicated non-Gaussian data distribution. Differences in age, BMI, mean PDFF, and TFs of all muscle compartments between males and females were assessed with Wilcoxon-Mann-Whitney tests. Differences in mean PDFF and TFs between the three muscle compartments were analyzed with Wilcoxon signed-rank tests. Furthermore, partial correlations (r), adjusting for age and BMI, were determined in a pairwise manner for mean PDFF and TFs between different muscle groups. For correlation testing, the whole included cohort as well as males and females were calculated separately.

Gender-Specific Results
Mean PDFF showed significant differences between men and women in ES (p < 0.0001), but not in CE (p = 0.149) and PS (p = 0.611).
There were multiple TFs with significant differences between men and women in all muscle groups. Amongst others, Variance (global), Kurtosis (global), and Dissimilarity showed significant differences between men and women in all included muscle groups (Table 1). For instance, Variance (global) demonstrated significantly higher values in men than in women in all included muscle groups (p < 0.001). In contrast to mean PDFF, genderspecific differences in CE and PS could be detected for Kurtosis (global) and Dissimilarity ( Figure 2).

Muscle-Specific Results
The mean PDFF values (p < 0.001), Energy (p < 0.001), Entropy (p = 0.01), Homogeneity (p < 0.001), Sum-average (p = 0.012) and Correlation (p = 0.037) differed significantly for all three muscle groups. There were also other TA, such as Kurtosis (global), Skewness (global) and Dissimilarity, which only differed in two muscle departments. In general, a high muscle specificity for the included TA with regard to all included muscle groups could be detected.

Correlations of PDFF Measurements and Texture Features between Muscle Groups
A partial correlation analysis, adjusted for age and BMI as potential confounders, was performed. No significant correlations could be detected regarding mean PDFF or TFs between the muscle groups in the whole cohort as well as for males and females separately. The TFs of each muscle group correlated with the extracted mean PDFF of the same muscle group, but not with the TFs of other muscle compartments. For CE, the TF correlation showed the highest correlation with mean PDFF (r = 0.786; p < 0.001). For ES, higher values in Variance (global) were associated with increased mean PDFF values (r = 0.863; p < 0.001). For PS, mean PDFF correlated significantly with Homogeneity (r = 0.334; p = 0.03; Figure 3).

Discussion
In the present study, the structural composition of the cervical and lumbar paravertebral musculature was assessed using the TA of PDFF maps generated with CSE-MRI. TA enabled the detection of gender-specific differences beyond the mean PDFF. In all muscle groups, PDFF heterogeneity as assessed by TA was greater in males than females. Within each muscle group, there were several TFs that correlated with the mean PDFF, but no significant inter-muscular correlations could be revealed. This finding suggests quite distinct anatomical location-specific structural patterns of the paraspinal muscles studied. Further, Variance (global) and Kurtosis (global) showed the potential for differentiating male and female structural patterns for the CE and PS, and thereby offered morphological information which could not be detected using PDFF values alone.
For a non-invasive quantitative multi-parametric muscle status assessment, CSE-MRI, MRS, and CT can be used in a preclinical but also in a clinical setting, providing valid biomarkers such as the mean PDFF and fatty acid characterization [9,11,37]. Recently, the literature on structural characterization of musculature using TA of PDFF maps has become available [19,38]. In these studies, an association of quantitative multiparametric muscle imaging and isometric strength was reported [28,38]. Due to an aging population and a predominantly sedentary lifestyle, alterations in muscle composition, such as increasing fatty infiltration, with associated frailty and impaired functional integrity, are valuable markers of functional muscle integrity and may allow an early diagnosis of impaired function and prediction of adverse health consequences such as frailty [39]. Detecting and quantifying structural changes in the initial stages but also differentiating pathologic alterations from age-related degeneration will become a challenge in the near future.
In all three muscle compartments, the TFs Variance (global), Kurtosis (global), and Dissimilarity differed significantly in all three muscle compartments between males and females. These results are in line with prior studies showing the gender-associated structural variations in paravertebral musculature [14,15]. Our findings point to distinct intra-and inter-individual TF distribution patterns and suggest different underlying pathophysiological mechanisms for different muscle groups. By circumscribing the physiological spectrum of TF distribution in healthy subjects, our study is a first step towards the quantitative visualization of the pathophysiologic processes occurring in muscle tissue due to aging and related to gender. As our findings demonstrate, TA offers in-depth insights into paravertebral muscle composition beyond simply quantifying the extent of fatty infiltration. Thus, the present results bring us closer to understanding the intra-and inter-individual variations with regard to differentiating the start of disease manifestations from normal physiological conditions.
To the best of our knowledge, this is the first study to show gender-specific differences in paravertebral muscle composition using TFs in a considerably large cohort of healthy subjects. The results extracted from our analysis could be used in oncologic patients at risk of sarcopenia or cachexia to improve the early detection and monitoring of these severe complications and to facilitate adequate intervention [4,8,40,41]. The comparably short scan time of about two minutes per acquisition for the lumbar spine could be added to the scanning protocol in the course of oncologic disease monitoring without significant additional time expense. Altered TFs could be indicative of the start of structural changes before volumetric changes or functional impairments are evident. In the light of the potential clinical implementation in primary and secondary preventions, the reported data might serve as a reference standard for early pathology detection at the muscle level.
There are some methodological and conceptional limitations to our study. First, the retrospective nature of the study design makes it prone to selection bias. Second, manual segmentation always holds the possibility for inaccuracies. However, freely available and reliable automated segmentation algorithms are just about to be implemented for analyses of paraspinal muscles. Third, although a wide range of different age groups was recruited, we had a 2:1 female-to-male ratio, which has to be considered when interpreting our results. Fourth, there were no follow-up scans, which would have allowed a longitudi-nal analysis. Further longitudinal and interventional studies could also reveal changing muscle composition in association with exercise or caloric restriction. Lastly, the analysis lacks information about potential confounding factors such as physical activity. Future analyses should include information on such confounders.

Conclusions
The results of the present study demonstrate gender-specific and anatomical locationspecific, distinct differences of fatty paraspinal muscle infiltration in a cohort of healthy subjects. The missing inter-muscular correlations after partial correlation testing suggest increasing BMI and age to be driving forces for muscle tissue changes. Thus, the presented TFs provide additional structural information of different muscle groups with potential for clinical implementation in differentiating physiologic from pathologic alterations. Informed Consent Statement: Informed consent was obtained from all subjects involved in the study. Written informed consent has been obtained from the patients to publish this paper.

Data Availability Statement:
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A.
Appendix A. 1

. Global Texture Features
Let P define the first-order histogram of a volume V(x, y, z) with isotropic voxel size. P(i) represents the number of voxels with gray-level i, and N g represents the number of gray-level bins set for P. The i th entry of the normalized histogram is then defined as: The Global texture features are then defined as: Appendix A.2. Second-Order Texture Features Let P define the gray-level co-occurrence Matrix (GLCM) of a quantized volume V(x, y, z) with isotropic voxel size. P(i, j) represents the number of times the voxels of graylevel i were neighbors with voxels of gray-level j in V, and N g represents the pre-defined number of quantized gray-levels set in V. The entry (i, j) of the normalized GLCM was defined as: The following quantities were also defined: The second-order texture features are then defined as: [i p(i, j) + j p(i, j)] (A16)