T2* and quantitative susceptibility mapping in an equine model of post-traumatic osteoarthritis: assessment of mechanical and structural properties of articular cartilage

Objective: To investigate the potential of quantitative susceptibility mapping (QSM) and T2* relaxation time mapping to determine mechanical and structural properties of articular cartilage via univariate and multivariate analysis. Methods: Samples were obtained from a cartilage repair study, in which surgically induced full-thickness chondral defects in the sti ﬂ e joints of seven Shetland ponies caused post-traumatic osteoarthritis (14 samples). Control samples were collected from non-operated joints of three animals (6 samples). Magnetic resonance imaging (MRI) was performed at 9.4 T, using a 3-D multi-echo gradient echo sequence. Biomechanical testing, digital densitometry (DD) and polarized light microscopy (PLM) were utilized as reference methods. To compare MRI parameters with reference parameters (equilibrium and dynamic moduli, proteoglycan content, collagen ﬁ ber angle and -anisotropy), depth-wise pro ﬁ les of MRI parameters were acquired at the biomechanical testing locations. Partial least squares regression (PLSR) and Spearman's rank correlation were utilized in data analysis. Results: PLSR indicated a moderate-to-strong correlation ( r ¼ 0.49 e 0.66) and a moderate correlation ( r ¼ 0.41 e 0.55) between the reference values and T2* relaxation time and QSM pro ﬁ les, respectively (excluding super ﬁ cial-only results). PLSR correlations were noticeably higher than direct correlations between bulk MRI and reference parameters. 3-D parametric surface maps revealed spatial variations in the MRI parameters between experimental and control groups. Conclusion: Quantitative parameters from 3-D multi-echo gradient echo MRI can be utilized to predict the properties of articular cartilage. With PLSR, especially the T2* relaxation time pro ﬁ le appeared to correlate with the properties of cartilage. Furthermore, the results suggest that degeneration affects the QSM-contrast in the cartilage. However, this change in contrast is not easy to quantify.


Introduction
Osteoarthritis (OA) is a progressive disease that leads to restricted mobility and severe joint pain1,2. The onset of OA may be due to joint trauma, such as a ligamental tear or focal cartilage loss. 2 Magnetic resonance imaging (MRI) is among the best non-invasive tools available for diagnosis of OA. Especially, quantitative MRI methods have been proposed for diagnosis of OA. 3 However, quantitative parameters are usually time-consuming to measure in full 3-D imaging due to the requirement of multiple datasets to allow for data quantification (e.g., measurement of T1 relaxation time using gold standard inversion recovery sequence); furthermore, they usually require specific sequence design. 4e7 To overcome the limitations set by desired imaging time and resolution, we propose the use of quantitative susceptibility and T2* relaxation time mapping (QSM and T2*, respectively) that are measurable using standard multi-echo gradient echo sequences. QSM maps the magnetic susceptibility distribution inside the imaging target, which can be altered for example by accumulation of iron or calcifications inside the target region. 8e11 In OA, mineralization of cartilage has been reported 12,13 ; hence QSM might be sensitive in diagnosis of OA. So far, mostly preliminary studies with limited data sets on QSM of articular cartilage have been reported. 14e16 Wei et al. showed that susceptibility of articular cartilage displays orientational anisotropy and that the susceptibility has depth-wise varying contrast in articular cartilage. 14 In another study, Wei et al. utilized the anisotropy of susceptibility to perform collagen fiber orientation tracking using susceptibility tensor imaging. 15,17 Nyk€ anen et al. showed that proteoglycan loss did not affect the susceptibility in ex vivo cartilage. 16 The same study also revealed that anisotropy of QSM was seemingly different from the anisotropy of T2* relaxation, which has been linked to the structural anisotropy of the collagen network in cartilage. 16,18 Furthermore, in a recent in vivo study, Wei et al. have shown that the decrease of depth-wise variation in the susceptibility is linked to osteoarthritis in vivo. 19 In conclusion, previous studies have pointed out that the collagen fiber orientation could be the source of the QSM contrast in the cartilage. Besides the few studies on articular cartilage, some QSM studies have targeted the venous structures of the epiphysealarticular cartilage complex. 20e23 T2* and related T2 relaxation time measurements have been studied more extensively in cartilage imaging and have been related to properties of the collagen network. 24 Our aim was to validate the use of QSM and T2* for prediction of cartilage biomechanical and structural properties using samples from an equine post-traumatic OA study. 25 We hypothesized (1) that cartilage properties can be predicted from depth-wise profiles of quantitative MRI-parameters using partial least squares regression (PLSR) and (2) that this is a more sophisticated way of utilizing MRI data than direct Spearman's rank correlation analysis between bulk values of quantitative MRI and reference parameters.

Samples
Osteochondral samples were obtained from a cartilage repair study involving seven Shetland ponies (N ¼ 7, 6 females and 1 male, Age ¼ 8.8 ± 3.5 years). 25 For the repair study, surgical lesions were induced in the medial trochlear cartilage in both hinds and repaired with a combination of chondrons and mesenchymal stem cells in different carrier hydrogels. The ponies were sacrificed after 1 year. Wedge-shaped samples ( Fig. 1) were collected from the medial trochleae post-mortem and included part of the lesions as well as surrounding cartilage. Stifle joints from three healthy ponies (with matching age range) were acquired from a local abattoir (van de Veen, Nijkerk, The Netherlands). From these stifles six control samples (one from each stifle) were taken, leading to a total of 20 samples (14 experimental, 6 control). The number of animals was based on the power calculations performed for the original cartilage repair study. 26 The ethical permission for animal study was given by the Ethics Committee of Utrecht University for Animal Experiments in compliance with the Institutional Guidelines on the Use of Laboratory Animals. The animal study was carried out in a surgical theatre at the Department of Equine Sciences, Utrecht University, The Netherlands (Permission DEC 2014.III.11.098).
Biomechanical testing, digital densitometry, and polarized light microscopy Prior to MRI, the samples underwent biomechanical indentation testing. The indentation was performed using a system that had a 250 g load cell (accuracy ± 0.25%, Model 31, Honeywell Sensotec Sensors, Columbus, OH, USA) and an actuator (displacement resolution 0.1 mm, PM500-1 A, Newport, Irvine, CA, USA). Several locations (proximal, central and distal with respect to the anatomical location in the joint) were investigated, because cartilage thickness and properties vary throughout the joint. Thus, biomechanical measurements were conducted on a pre-defined grid of 12 testing locations on each sample (Fig. 1), leading to a total number of 235 biomechanical testing locations (two samples had less than 12 locations). 25 Adjacent points could be measured after each other in quick succession since the utilized indender had a small diameter (~0.5 mm) compared to the distance between the adjacent points (~3e4 mm). 27 Equilibrium and dynamic moduli were determined from the biomechanical testing results.
After the MRI measurements, the samples were fixed in formalin and then decalcified in EDTA to soften them for histological sectioning. Sections were cut to include two biomechanical testing locations with matching distances from the edges of the block to ensure reliable location matching for microscopical imaging. Sections were stained using safranin-O for digital densitometry (DD) analysis to reveal proteoglycan content. 28 The staining was done in one process, using the same batch of dye for all the slices to ensure consistency of the stain concentration. Unstained sections were digested with hyaluronidase to remove proteoglycans to prepare them for polarized light microscopy (PLM), which was used to reveal collagen fiber angles and anisotropy. 29 DD was performed using a light microscope equipped with a monochromatic light source and a CCD-camera. The settings of the microscope and the camera were the following: objective 1Â, zoom 2Â, binning 1 and calibration with neutral density filters of 0e3.0 optical density. The optical density of each sample was calculated based on the calibration set. To calculate the average depthdependent profile for each measurement location, the subchondral bone was manually segmented from the images, followed by automatic extraction of profiles perpendicular to the cartilage surface and interpolation of profiles to 100 points. The analysis was performed with a custom Matlab algorithm (Matlab 2016b, Mathworks, Inc.).
PLM was performed using an Abrio PLM imaging system (CRi Inc., Woburn, MA, USA) which was mounted on a light microscope (Nikon Diaphot TMD, Nikon Inc., Shinagawa, Tokyo, Japan). The pixel size for both PLM and DD images was 3.5 mm Â 3.5 mm. The fiber angle anisotropy was calculated from fiber angle images using the following equation: where εðrÞ is the pixel-wise local entropy of a 5-by-5 pixel region in the collagen fiber orientation image. This method was used since it provides a good approximation of the true collagen fiber anisotropy 29 , which could not be calculated directly due to lack of access to the raw data acquired with the Abrio system. The fiber angle and anisotropy profiles were gathered from microscope images using a similar procedure as in optical density measurements.

MRI
MRI was performed in a 9.4 T vertical bore small animal scanner using a 19-mm-diameter quadrature RF volume transceiver (Rapid Biomedical, Rimpar, Germany). The samples were immersed in 1HMRI-signal-free perfluoropolyether oil (Galden HS 240, Solvay Solexis, Brussels, Belgium) inside a thin latex holder. During MRI, the cartilage surface was oriented approximately perpendicular to the main magnetic field of the scanner with extreme care, as orientation-anisotropy has been reported for both T2* and QSM of articular cartilage. 15,16,18 Imaging was conducted using a multi-echo gradient echo sequence with six echoes. The first echo time of the sequence was 2.0 ms and echo spacing was 3.05 ms. Isotropic voxel size of 100 Â 100 Â 100 mm 3 was used and the matrix size of the image was 200 Â 256 Â 200 voxels. The longest dimension, which was also the readout direction, was in line with the main magnetic field of the scanner.
After collecting the raw MRI data, the T2* relaxation time as well as QS-maps were reconstructed. In QSM post-processing, the "complex fitting"emethod was used to combine multi-echo data 30 , followed by Laplacian unwrapping to resolve residual wraps in echo combination data. 31 The unwrapped field map was then masked using cartilage segmentation. The susceptibility maps were calculated from masked unwrapped phase data using a total field inversion (TFI) method to limit the need of region of interest (ROI) erosion in the QSM dipole inversion. 32 The erosion of the ROI in conventional QSM methods such as projection onto dipole fields (PDF) þ morphology enabled dipole inversion (MEDI) 33 is depicted in supplemental material ( Supplementary Fig. S1). The susceptibility values were referenced by setting the mean over an individual sample to zero susceptibility. T2* relaxation time maps were calculated using 2-parameter linearized fitting of the data. In the T2* fitting procedure, voxels with T2* higher than 150 ms were considered erroneous and their T2* value was set to 150 ms. Three- dimensional surface maps were gathered from MRI parameters to visualize differences between the experimental and control samples.

Data analysis
Depth-wise profiles from the 3-D T2*-and QS-maps were obtained using cylindrical 3-D ROIs of 1-mm diameter and carefully matched with the biomechanical testing points based on mCT measurements and photographs of the samples. 25 Prior to the analysis, these profiles were interpolated into 100 depth-wise points. The profiles were then used to predict reference variables (i.e., equilibrium and dynamic moduli, proteoglycan content, collagen fiber angle and collagen anisotropy) using PLSR. Prior to the PLSR-analysis, depth profiles were normalized using standard normal variate. To avoid overfitting of the models, PLSR-analysis was conducted using 10-fold cross-validation. PLSR analysis was performed using a custom-made algorithm in MATLAB. Predictions from PLSR-analysis were compared with reference parameter values using Spearman's rank correlation analysis. In addition to PLSR-analysis, direct correlations were calculated between means of full and superficial (25% depth) MRI profiles and reference parameter values. For QSM, the range of the susceptibility values within the profile was used instead of the mean due to lack of absolute QSM referencing in the imaging.
Statistical significance of differences in tissue properties between experimental and control groups was tested with the ManneWhitney U test in SPSS (Version 25, SPSS Inc., IBM Company, Armonk, NY, USA). Statistical significance of Spearman correlations was tested using exact permutation distributions in MATLAB. In both tests P < 0.05 was considered as the limit for statistical significance. Non-parametric tests were utilized in this study due to the non-normality of the reference data (ShapiroeWilk test P < 0.0001).
Surface map analysis was conducted through the following steps. First, a tetrahedral grid was created in the cartilage mask (every 20th point of the mask was used as nodes). Surface triangles were then sought and their surface normals were calculated. After this, the 3-D QSM and T2* maps were sampled from within the volume along the surface normals. From these sampled profiles, the range for QSM and mean for T2* was calculated. Finally, the QSM range and T2* mean values were displayed on the surface triangles of the mesh to generate the parametric surface maps.

Data availability
All raw data and documentation, as well as key analysis codes used in this study, are available for download at Zenodo (https:// doi.org/10.5281/zenodo.2558172).

Biomechanical testing and histological analysis
Differences between experimental and control specimens were noted with all reference methods (Figs. 2e4). In biomechanical testing, both equilibrium and dynamic moduli were decreased in the experimental specimens compared to moduli of control samples at the same location (P < 0.05, for all distal and most of the central locations) and the differences became clearer when moving towards more distal aspect of the specimens (referring to the position in the joint) (Fig. 2). It is also worth noting that when moving towards the lesion sites (i.e., towards the medial aspect of the femoral groove), also the moduli of the control specimens decreased (Fig. 2) (see Fig. 1 for terminology). Further results of the biomechanical testing have been reported in a previous study. 25 Overall, lower proteoglycan contents were observed in experimental specimens compared to controls. The difference between control and experimental groups was largest at the regions nearby the lesions and decreased further away from the lesions for proximal and distal locations (Figs. 3 and 4). Conversely, for central locations differences in PG content became clearer while moving away from the lesions. The collagen fiber angles in the radial zone of the experimental group were lower than in the control group and the difference between groups was highest at the locations nearby the lesions (Figs. 3 and 4), except for the proximal regions, for which no significant differences in fiber angle were observed. There were also differences between the collagen fiber anisotropy of the experimental and control groups, but contrary to the optical density and fiber angle, the differences became clearer when moving away from the lesion site and notable changes were only observed on the distal side of the specimens (Figs. 3 and 4).

MRI
QSM exhibited differences between the control and experimental groups especially nearby the lesion site in the central locations and further away from the lesion site at proximal and distal locations (Figs. 5e7). Nearby the lesion at central locations, QSM profiles had a smaller range in the experimental samples when compared to the controls (Fig. 5). When moving further away from the lesion site at distal and proximal locations, the range of the profiles increased in experimental samples when compared to the control samples (Fig. 5). Overall, QSM seemed to have more spatial variation in the experimental samples (Fig. 7).
The T2* profiles had qualitatively slightly lower maximum values nearby the lesion and higher values further away from the lesion site (Fig. 5). However, the mean value of T2* was increased nearby the lesion sites (i.e., towards the abaxial aspect of the femoral groove) in both experimental and control groups but slightly more so in the experimental specimens (Fig. 7). Cartilage thickness seemed to have no appreciable effect on the T2* or QSM values (Fig. 7, Supplementary Figs. S2e4).

Correlation and PLSR analysis
The absolute value of Spearman's rank correlation between MRI and reference parameters was between r ¼ 0.03 and r ¼ 0.46, with the highest correlations observed between bulk T2* and collagen anisotropy (r ¼ À0.46) and between superficial QSM and bulk fiber angle (r ¼ 0.40) (Table I). In general, the T2* bulk value had a stronger correlation with the reference parameters, whereas for QSM the strongest correlations were observed between the superficial QSM and the reference parameters (Table I). Most of the correlations in the direct correlation analysis were statistically significant, and only correlations near r ¼ 0.00 were nonsignificant.
PLSR-modelling results showed that the estimation of the reference parameters from MRI data was more successful using T2* or a combination of T2* and QSM (apart from the superficial collagen fiber angle and anisotropy, the correlations were between r ¼ 0.49 and r ¼ 0.68), and less successful using QSM results only (correlations between r ¼ 0.41 and r ¼ 0.55) ( Table I, Supplementary Fig. S5). For combined T2* and QSM, mean PLM anisotropy showed the highest (r ¼ 0.68) and mean superficial PLM anisotropy showed the lowest (r ¼ 0.18) correlation between predicted and measured values (Table I). All correlations from PLSRanalysis were statistically significant.

Discussion and conclusions
In this study, we examined the potential of 3-D quantitative MRI to detect post-traumatic OA (PTOA) changes in cartilage properties. More specifically, QSM and T2* relaxation time mapping were investigated and correlated with biomechanical parameters, proteoglycan content and collagen network properties. The analysis revealed weak-to-moderate direct correlations   between MRI and reference parameters. PLSR-modelling of depth-wise MRI profiles increased the correlations between MRI and reference parameters, demonstrating that PLSR outperforms univariate analysis in quantitative MRI evaluation of articular cartilage.
Especially biomechanical testing, but also quantitative microscopic analyses (DD and PLM) revealed differences between the experimental and control samples (Figs. 2e4). Interestingly, large intrasample differences in both DD and PLM were observed within both experimental and control groups (Figs. 2e4). This could indicate that the cartilage nearby the trochlear ridges of femur is different from the cartilage in more central locations of the knee joint (Figs. 2e4). Particularly interesting were the PLM results, showing that the cartilage seemed to lose its usual trilaminar appearance nearby the lesion sites in both experimental and control groups and instead had fibers directed parallel to the cartilage surface throughout the cartilage thickness (Figs. 3 and 4). Marked structural changes in the properties of the collagen fiber network have been demonstrated earlier at the edges of canine humeral cartilage. 34 T2* relaxation time mapping was done using linear fitting through magnitude images of each individual echo of multi-echo gradient echo sequences. This procedure is well-established in cartilage MRI. T2*-mapping revealed the usual trilaminar appearance in both experimental and control groups further away from the trochlear ridge (Fig. 6). However, at locations nearby the ridge or lesion site, this trilaminar appearance was less visible, especially in the experimental group (Fig. 6). This supports the idea that cartilage near the lesion area exhibits degenerative changes. It is worth noting that for T2* results, the same main relaxation contributions should apply as for the more commonly used T2 relaxation time, because cartilage does not contain high susceptibility differences or other field homogeneity issues that would differentiate T2* from T2.
Both T2* and QSM displayed visual differences between the experimental and control samples (Figs. 5e7, Supplementary  Figs. S1e3). The surface visualization revealed more distinct changes in QS-than in T2*-maps. However, these changes were difficult to quantify, which may be reflected in the poorer PLSRmodelling results. One notable feature in QSM was the flattening of QSM depth profiles at central measurement points in the experimental group when compared to the control samples. This finding is in line with a recent study, statingtural that depth-wise standard deviation of QSM is decreased in OA. 19 The direct correlations between T2* and biomechanical and structural properties of cartilage were similar to those that have previously been reported between T2 and cartilage properties. 35,36 Furthermore, the PLSR-modelling yielded higher correlations between MRI and other cartilage properties than direct univariate analysis. Thus, in the future, modelling approaches such as artificial neural networks and PLSR would be preferred in the analysis of quantitative MRI of cartilage. While the PLSR-modelling could have certainly benefitted from an even larger data set, it outperformed direct correlation (r bulk ¼ 0e0.46) already in this limited data set and yielded high correlations (r PLS ¼ 0.41e0.68) between the predicted and measured reference parameters in analyses which took full-depth profiles into account.
In this study, QSM-post processing was based on the "complex fitting", Laplacian unwrapping and Total Field Inversion steps. 30e32 Several studies have been dedicated to finding the most appropriate processing pathways and methods, especially for brain imaging, which has been the main target for QSM. 9,37 For cartilage, QSM has not been fully established for imaging small ex vivo cartilage specimens, a concern that was also raised in a previous study. 16 The fundamental difficulty in QSM imaging of small ex vivo cartilage samples is that regardless of the methodology, the accuracy of QSM is sub-optimal nearby the boundaries of the ROI. 38,39 Due to the thinness of cartilage, the entire ROI comprising cartilage is always near to a tissue boundary, making QSM difficult. Here we utilized TFI because it doesn't necessitate erosion of the alreadythin cartilage mask. Supplementary material ( Supplementary  Fig. S1) provides a brief comparison of different methods and the implications with respect to small ROI. Encouragingly however, differences between the experimental and control groups were  Fig. 1 for sampling locations). Lines at the top of the image indicate regions, where the difference between control and experimental groups is statistically significant (P < 0.05, ManneWhitney U-test). detected in this study even though the processing is not fully optimized for such thin tissue structures of interest. In vivo imaging should at least partially overcome this problem, as there will be lots of other nearby signal sources in addition to cartilage, improving the susceptibility processing. Another fundamental difficulty (or perhaps a property) with the current QSM processing methods is that the susceptibility values are relative to the imaged target, i.e., if susceptibility is increased by 0.05 ppm in the whole target, this change may not be detectable. Thus referencing (e.g., adding a capsule of oil near the target tissue) would be beneficial in the QSM studies of cartilage or of any other tissue. 9,37 The main limitation of this study is the lack of referencing in QSM reconstructions, which was partially dealt with by using the range instead of the mean in the direct analysis. QSM could also have benefitted from imaging in signal-producing media such as saline instead of signal-free perfluoropolyether. On the other hand, this would have led to lower resolution due to potential wrapping in the phase encoding direction. Furthermore, the measurement of the T2* relaxation time may have been suboptimal since the last echo time (17.25 ms) was considerably shorter than the longest measured T2* relaxation times (~70 ms). However, this is considered to have only a small effect on the analysis of the data, since two-parameter linear fitting of T2* relaxation time is a relatively robust method. The robustness of the fit was further enhanced by limiting the maximum allowed value of T2* to 150 ms, which is a sufficient limit for the T2* of articular cartilage. The sequence parameters were optimized for QSM measurement, relieving the strain on the gradient system by reducing the number of echoes; the longest echo times were limited, thus briefly compromising the accuracy of the T2* relaxation time mapping. For PLSR-analysis, the dataset could have been larger, since intrasample differences were higher than expected (especially for PLM), which may affect the PLSR results due to the relatively high number of potential outliers in the data. While the number of animals as such may be considered low for the PLSR-analysis, the limitation was mitigated by performing multiple measurements at different locations in each sample, leading to a sufficiently high number of data points for the analysis. Furthermore, other qMRI parameters could have been included in the study, especially T1rho, which has been suggested as a good parameter in the evaluation of cartilage. 3,40,41 However, due to logistic reasons, MR-imaging time for the specimens was restricted and measuring additional parameters was not possible. In PLM, the lack of true anisotropy might have caused biased results when comparing MRI to anisotropy, since the anisotropy was calculated from the fiber angle images using image entropy.
This study showed that evaluation of the degradative state of articular cartilage using qMRI is improved when using PLSR-  modelling instead of direct correlations. Based on the present findings, PTOA processes change QSM-contrast in cartilage; while the exact reason remains unknown, the collagen fiber network structure appears to change as well, a finding that has been previously linked to changes in QSM contrast. 15,16 Moreover, T2* relaxation time was shown to be a more useful parameter than QSM in the evaluation of cartilage degeneration. Combining the results from both modalities resulted in slightly higher correlations between the measured and PLSR-predicted reference parameters, suggesting that QSM adds to qMRI evaluation of cartilage.
Tiitu: Histological sample preparation, manuscript editing. Mancini: animal experiments, design of the study, manuscript editing.
Visser: animal experiments, design of the study, manuscript editing.
Brommer: animal experiments, design of the study, manuscript editing.
van Weeren: animal experiments, design of the study, manuscript editing.
Malda: animal experiments, design of the study, manuscript editing.
T€ oyr€ as: design of the study, manuscript editing. Nissi: design of the study, MRI measurements, data-analysis, manuscript drafting and editing.
All authors: revision for intellectual content, final approval.

Conflict of interest
Authors declare no conflicts of interest.

Funding sources
Academy of Finland, Finnish Cultural Foundation, Northern Savo Regional Fund of Finnish Cultural Foundation, Dutch Arthritis Foundation.

Role of the funding sources
Funding sources had no influence on study design or interpretation.