Detection of Early-Stage Degeneration in Human Articular Cartilage by Multiparametric MR Imaging Mapping of Tissue Functionality

To assess human articular cartilage tissue functionality by serial multiparametric quantitative MRI (qMRI) mapping as a function of histological degeneration. Forty-nine cartilage samples obtained during total knee replacement surgeries were placed in a standardized artificial knee joint within an MRI-compatible compressive loading device and imaged in situ and at three loading positions, i.e. unloaded, at 2.5 mm displacement (20% body weight [BW]) and at 5 mm displacement (110% BW). Using a clinical 3.0 T MRI system (Achieva, Philips), serial T1, T1ρ, T2 and T2* maps were generated for each sample and loading position. Histology (Mankin scoring) and biomechanics (Young’s modulus) served as references. Samples were dichotomized as intact (int, n = 27) or early degenerative (deg, n = 22) based on histology and analyzed using repeated-measures ANOVA and unpaired Student’s t-tests after log-transformation. For T1ρ, T2 and T2*, significant loading-induced differences were found in deg (in contrast to int) samples, while for T1 significant decreases in all zones were observed, irrespective of degeneration. In conclusion, cartilage functionality may be visualized using serial qMRI parameter mapping and the response-to-loading patterns are associated with histological degeneration. Hence, loading-induced changes in qMRI parameter maps provide promising surrogate parameters of tissue functionality and status in health and disease.


Results
All 49 samples underwent complete MR imaging and biomechanical as well as histological assessment. Figure 1 gives a graphical illustration of the different Mankin sum scores (MSS). After sample dichotomization, 27 samples constituted the intact (int) group and 22 samples the degenerative (deg) group. Consequently, significant group-wise differences were found in MSS (int vs. deg: 2.4 ± 1.1 vs. 6.0 ± 0.9; p < 0.001) and Young's Modulus (YM) (0.6 ± 0.3 vs. 0.4 ± 0.3 [MPa]; p = 0.003) ( Table 1). Additional details of the histological, biomechanical and unloaded qMRI parameter values are given in Table 1. In the unloaded configuration, significant differences between int and deg samples were only found for T1 and T2, where int samples displayed significantly lower values than deg samples, however, only in the deep tissue zones (T1 [dp]: p = 0.029; T2 [dp]: p = 0.008).
Distinct qualitative and quantitative changes were observed with loading: First, mean sample height was significantly reduced (Δ 2.5 : −6.1 ± 22.0%; Δ 5.0 : −15.6 ± 25.8%; p < 0.001), while mean sample width remained about constant (Δ 2.5 : 1.9 ± 7.8%; Δ 5.0 : 3.1 ± 11.1%; p = 0.127). Correspondingly, the number of detected pixels significantly decreased in all zones ( 5.0 ], p < 0.001). Second, distinct loading-induced changes were found in absolute qMRI parameter values (Table 2 and Supplementary Table S1): For T1, consistent significant decreases were found in all tissue zones, irrespective of degeneration. For T1ρ, significant increases were observed in deg samples' deep zones only. This increase in deg samples rendered T1ρ changes over all samples significant, too. In int samples, however, T1ρ changes were inconsistent. For T2, significant decreases were found in the superficial zone and corresponding increases in the deep zone. Although these changes were significant in all regions-of-interest (ROIs) of int samples and only in the deep zone of deg samples, they were more pronounced in deg samples. For T2*, changes were equally more pronounced and significant in deg samples only.
Best sensitivities to differentiate int from deg samples were found for T1ρ and T2 (unloaded) and for T2* Δ 5.0 (loaded). Highest combined sensitivities were found for T2 -T2*Δ 5.0 and T1ρ -T2*Δ 5.0 , while highest combined specificities were found for T2 -T2Δ 2.5 and T1 -T2Δ 2.5 (Table 3). Group-wise comparisons of relative changes in qMRI parameters did not reveal significant differences (Supplementary Table S2). Similarly, no significant correlations were found between relative changes in qMRI parameters and YM or MSS (Supplementary Table S3).
Morphologically, sample width and signal intensity in PDW images remained grossly unaltered, while sample height consistently decreased (Fig. 2a-c). Serial qMRI parameter maps were reflective of the quantitative changes outlined above: While in int samples, relatively homogeneous parameter distributions at δ 0 were found that remained largely unaltered at δ 2.5 and δ 5.0 (Fig. 3), deg samples displayed inconsistent loading-induced changes: In some samples, diffuse signal alterations became more widespread to eventually involve the majority of the sample cross-sectional area (Fig. 4), while in other samples focal signal alterations were less distinct to nearly disappear altogether ( Supplementary Fig. S1).

Discussion
The most important finding of this study is that loading-induced changes in some serial qMRI parameter maps are related to histological degeneration and improve their diagnostic accuracy, which provides a solid scientific framework to assess cartilage functionality in future applications.
As demonstrated previously 22 , the MRI-compatible whole-knee joint loading device delivers standardized and reproducible compressive loading of chondral samples. Distinct patterns of intra-tissue changes were detectable in morphology, sample height and pixel numbers, indicating effective sample pressurization. Loading-induced changes in qMRI parameters are reflective of sufficient sample pressurization. Consistent decreases in T1 were found in all samples, irrespective of degeneration, while significant increases in T1ρ and T2* were observed in deg samples only. For T2, changes were related to the tissue zone with decreases (superficial zone) and increases (deep zone) observed alike.
Loading-induced decreases in T1 have been reported before 13,19 , even though loading protocols differed. Xia et al. observed consistent reductions in T1 and a clear relation to the applied strain intensity 13 , which is reflected by our data as greater changes in T1 were found at δ 5.0 than at δ 2.5 . Although reductions in T1 were observed throughout the entire tissue depth, superficial-zone changes were larger than deep-zone changes, which is due to marked structural and compositional differences in cartilage 23 : The superficial zone is softest because of its limited fixed charge density and water content 24 . As T1 is widely considered a marker of tissue hydration 25 , considerable water redistribution within and -most likely-out of the tissue happens with loading. Saarakkala et al. found that the water content hardly changes in early degeneration 9 ; hence, similar loading-induced changes in T1, irrespective of degeneration, become plausible in light of the dominating effect of water on T1 characteristics 26 . Correspondingly, the diagnostic profile of T1 in differentiating int from deg cartilage is relatively weak as is demonstrated by only moderate sensitivity.
Disparate depth-and degeneration-dependent loading patterns were observed for T1ρ. Significant increases were found in the deep zone of deg samples only, although changes were similar, yet non-significant, in int samples. Opposite, yet non-significant changes were found in the superficial zone with T1ρ decreases in int samples and increases in deg samples. In absolute terms, changes in T1ρ were more pronounced in deg than int samples. These findings are in line with recent in-vivo data, as Souza et al. reported that changes in T1ρ are considerably larger in OA patients than controls 16 . However, they also observed significant decreases in the superficial and increases in the deep zone, irrespective of degeneration. Most likely, this discrepancy is secondary to doubtful patient allocation procedures based on radiographic evaluation, which is coarse and questionable when applied as www.nature.com/scientificreports www.nature.com/scientificreports/ reference measure 2 . Additionally, T1ρ changes of the medial compartment are not linearly corelated with overall OA severity 27 , thereby further compromising this reference measure.
Nonetheless, T1ρ is a promising indicator of biologically meaningful intra-tissue adaptive processes 12,15,16,18 that may be visualized using serial T1ρ mapping, even though the exact determinant of T1ρ remains to be defined 2,4,28 .
Physiologically, these adaptive intra-tissue processes involve water redistribution and reductions in tissue thickness, thereby increasing the relative proteoglycan concentration. Additionally, the ECM is condensed, deformed and altered in its orientation so that, in conclusion, the complex interplay of the solid and fluid cartilage phases determines the tissue's loading response in healthy and diseased cartilage 18 . The extent of these adaptive intra-tissue changes, however, is related to degeneration with larger changes noted in more degenerative tissue 9,10 .
Similar observations as for T1ρ were made for T2 with significant decreases in the superficial and increases in the deep zone. As T2 is an indicator of water content, collagen composition and collagen anisotropy 29 , these changes are also reflective of intra-tissue adaptations outlined above. The diagnostic profile of T2 is equally promising; yet, the discriminatory power of T2 and T1ρ in static and functional contexts is still discussed 12,17,30 and remains to be defined in future in-vivo studies.
For T2*, significant loading-induced increases were only observed in deg samples. The exact structural and compositional correlate of T2* is still debated 31,32 ; hence, the target structure and/or mechanism thus measured in its functional contribution still needs to be identified. Even though reports on whether biologically meaningful changes are detectable by T2* mapping are conflicting 7,31-33 , our study suggests that T2* can be successfully applied to assess degeneration-dependent adaptive intra-tissue changes in response to loading, in particular at  Segmentation included the entire cartilage sample (ECS) as well as superficial (sf) and deep (dp) zones. Upon log-transformation, repeated measures ANOVA was used to detect differences between δ0 (unloaded), δ2.5 (2.5 mm displacement) and δ5.0 (5.0 mm displacement). Data are mean ± standard deviation [ms] and p-value. Significant differences are displayed in bold-type followed by consecutive numbers indicating post-test details (see Supplementary Table S1). Abbreviations as in Table 1. −16.7; 20.7 n/a n/a n/a n/a n/a n/a n/a n/a <−16.7; >20.7 n/a n/a n/a n/a n/a n/a n/a n/a www.nature.com/scientificreports www.nature.com/scientificreports/ large strains, which is in line with earlier findings 34 . However, multi-gradient echo sequences (used for T2* mapping) are more prone to susceptibility artefacts secondary to local magnetic field inhomogeneities [31][32][33] , specifically when interfaces (as in cartilage) are present. Also, T2* measurements are affected by refocusing pulses, gradient spoilers and echo spaces 32 , which limits inter-study comparability and questions the reliable inter-patient and inter-study quantification of T2* in cartilage, regardless of loading. Against this background, intra-patient referencing to the unloaded configuration as in our study seems to provide a scientifically sound and clinically feasible approach to eliminate inter-patient and inter-study concerns pertaining to variability in imaging protocols, unloading times, analysis routines, scanner and coil configurations, and methods of segmentation and registration 8 .
Nonetheless, this study suggests that cartilage functionality assessment ought to be multiparametric and serial T1ρ, T2 and T2* mapping techniques seem to be most promising to differentiate the tissue's status in health and disease, not least due to their distinctly different sensitivity profiles 26 . If structural integrity of cartilage needs to be confirmed based on qMRI parameter values and their response-to-loading patterns, sensitivity needs to be highest and stand-alone T1ρ, T2 and the loading-induced changes in T2* should be assessed. Correspondingly, if degeneration needs to be diagnosed with a high degree of confidence, stand-alone T2 and the loading-induced changes in T2 should be studied.
Surprisingly, no significant correlations were found between relative changes in qMRI parameters and biomechanical or histological properties. Although significant differences in YM were found as a function of degeneration, the biomechanical properties are largely determined by ECM integrity rather than composition 35 . Therefore, compositional changes (as primarily assessed by MRI techniques) are not necessarily reflective of structural tissue properties.
Considerable standard deviations were observed throughout this study, which diminished statistical power. Biologically, the substantial inter-individual variability observed at δ 0 is reflected by the samples' variable loading responses. For the future in-vivo translation the predictive ability and diagnostic accuracy of these metrics have to be thoroughly addressed, especially in view of additional complexities such as joint positioning and compression www.nature.com/scientificreports www.nature.com/scientificreports/ status. Nonetheless, tissue degeneration is one contributory factor to cartilage functionality (among others) that needs to be considered alongside other patient-and joint-level factors such as age, sex, constitution and sports activities.
This study has limitations that involve technical, engineering and biological aspects. Considerable stress relaxation (of up to 48%) 22 was observed during loading, which is not surprising, given the displacement-controlled loading setup. However, as the same order of sequence acquisitions (i.e. T2*-T2-T1ρ-T1) was maintained throughout the study, samples have experienced different loading conditions during T2* than T1 measurements. For improved standardisation, pressure-controlled loading configurations should be implemented. However, in this regard, further research is necessary to determine how exactly loading-induced changes in cartilage (as  For all qMRI parameters, loading-induced increases were observed. Irregular signal hyperintensities (at δ 0 ) became more diffuse and widespread to eventually involve the majority of the sample's cross-sectional area (in particular T1ρ and T2). Histologically, superficial clefts (indicated by single arrow in e 1 ) and pannus formation (indicated by double arrows in e 1 ) were found alongside diffuse hypercellularity (only visible at higher magnification [not shown]) and severe discoloration on proteoglycan staining. MSS 8. Image details as in Fig. 3. www.nature.com/scientificreports www.nature.com/scientificreports/ assessed by qMRI parameters) are related to stress in comparison to strain. Also, this setup's realization of truly physiological exam conditions was limited as meniscus and ligaments were not included, which increases cartilage stress levels upon loading 36 . Additionally, any other physiological motion beyond uni-axial compression along the mechanical leg axis was precluded (i.e. the 'screw-home' mechanism at full extension). Loading therefore equalled a well-conformed, yet confined, compression test. Moreover, wirosil ® (i.e. the artificial cartilage material) is slightly less stiff than human cartilage 22 and future studies need to determine whether the observed response-to-loading patterns are similar to physiological cartilage-cartilage bearings. Yet, recent in-vivo data indicated that cartilage experiences a combination of compression and shearing when bearing weight 37 . As, physiologically, cartilage is relatively compliant, yet incompressible, and constrained by the underlying subchondral bone, compressive loading induces lateral expansion and secondary shearing, which is more relevant to the loading response than compression itself 37 .
Another important aspect is sample standardization. Chondral samples were cut to 3 mm thickness, which led to variability in the deeper tissue zones to be included per sample as cartilage thickness is variable. This is functionally relevant as deep tissue zones are essential for load-bearing due to their role in fluid convection and pressurization 38 . Hence, human cadaver or in-patient studies (with subsequent tissue harvest, e.g. response-to-loading assessment by MRI prior to total joint replacement) need to address these aspects to pave the technique's clinical translation. In this context, standardization of loading of a given joint compartment or region and exact matching to reference measures are perspective challenges to tackle. Our sample source (i.e. joint replacement material) is problematic as even grossly intact samples exhibited signs of early degeneration such as matrix discoloration or hypercellularity, thereby rendering the dichotomization of int vs. deg samples somewhat arbitrary in view of the continuous degenerative changes. Additional tissue sources (i.e. organ donor networks or amputations) may help overcome this issue.
In conclusion, distinct patterns of qMRI parameter changes were found in this in-situ study and in response to loading. Changes in T1ρ, T2 and T2* were different in deg as compared to int cartilage, while changes in T1 were consistent in all samples irrespective of degeneration. Even though these metrics' diagnostic accuracy and predictive ability needs to be better defined in entire joint configurations, the non-invasive assessment of cartilage functionality based on serial qMRI mapping provides an exciting framework to further stratify cartilage degeneration beyond mere static analysis.

Methods
Industry support. This study was supported by Philips Healthcare (Hamburg, Germany) by providing the T1ρ sequence. The authors had and have full control over the data and information submitted for publication. study design. This study was designed as a prospective comparative ex-vivo imaging study of cartilage samples that were obtained from total knee replacements at our institution between 10/2015 and 10/2016. Local Institutional Review Board approval of all experimental protocols and the use of the cartilage-bone material (Ethical Committee, RWTH Aachen University, Germany, AZ-EK157/13) was obtained beforehand as well as individual informed patient consent. The methods were carried out in accordance with the relevant guidelines and regulations. This study's datasets are available on sensible request by contacting the corresponding author.
Compressive loading device. The MRI-compatible whole-joint compressive loading device was described and validated in an earlier study 22 . Briefly, a right formalin-fixed human knee had been scanned by computed tomography in its native configuration and digitally processed to create standardized femoral and tibial bone models. The bone models had been covered by cartilage-mimicking polyvinyl-siloxane (Wirosil) as artificial femoral and tibial cartilage layers in their native configuration. A standardized defect at the central medial femoral condyle (8-mm diameter, 3-mm depth) had been created at the site of the initial contact into which the native www.nature.com/scientificreports www.nature.com/scientificreports/ chondral samples of corresponding dimensions were placed. Care was taken to precisely align the samples to avoid any step-offs (Fig. 5a). The standardized knee joint was surrounded by an artificial joint capsule filled with PBS buffer (Gibco-BRL) (Fig. 5b). Uni-axial compressive loading was performed by displacing the mobile tibia versus the immobile femur to 2.5 mm (δ 2.5 ) and 5.0 mm (δ 5.0 ) displacement (as measured from the initial contact point). Displacement-controlled loading resulted in mean forces on the entire joint as determined by use of a digital hydraulic force gauge (#HKMD29D, Induk, Wuppertal, Germany) of 141 ± 8 N (δ 2.5 ) and 906 ± 38 N (δ 5.0 ), corresponding to 20% and 110% of standard body weight. Mean local pressures on the chondral samples as determined by use of digital electronic pressure-sensitive sensors (K-Scan 4000, 10.000 psi, Tekscan, Boston, MA, USA) were 0.7 ± 0.1 MPa (δ 2.5 ) and 1.1 ± 0.1 MPa (δ 5.0 ). Of note, loading-induced strains within the chondral samples had not been determined. Compressive loading induced significant decreases in sample height (δ 0 : 2.86 ± 0.25 mm; δ 2.5 : 2.56 ± 0.25 mm; δ 5.0 : 2.02 ± 0.16 mm; p < 0.001 [repeated-measures ANOVA]), while sample circularity remained unchanged. Similarly, T2 signal intensity was decreased as a sign of sufficient tissue pressurization during loading. The interested reader may be referred to 22 for additional information on the reference and validity measurements. MRI compatibility and the absence of significant B 0 -inhomogeneity was confirmed by B 0 -mapping.
Cartilage sample preparation. Cartilage-bone material was harvested from 49 patients undergoing total knee arthroplasty at our institution (26 men, 23 women; 23 right and 26 left knees; mean age, 67.7 years [range, 46-93 years]) as before 19 . Primary OA of at least one joint compartment as determined radiographically (i.e. Kellgren-Lawrence grades ≥2 39 ) was defined as the inclusion criterion, while all forms of secondary OA or other bone and joint disorders (e.g. avascular necrosis or rheumatoid arthritis) and a history of previous trauma or surgery to the knee were defined as exclusion criteria. After excision, the cartilage-bone material was collected in sterile Dulbecco's modified Eagle's medium (DMEM) containing 100 U/ml penicillin, 100 µg/ml gentamycin and 1.25 U/ml amphotericin-B (Gibco-BRL, Gaithersburg, US) and prepared according to standard: First, the medial femoral condyle was identified and only cartilage from that joint region was included for reasons of topoanatomic consistency. An 8-mm diameter skin biopsy punch and #13-scalpel (both from PFM-Medical, Cologne, Germany) were used to create cylindrical chondral samples by removing the subchondral bone. A dedicated metallic cutting device (i.e. metal block with circular moulds of 8-mm diameter and 3-mm depth) was used to cut chondral samples to standard thickness of 3 mm, which was confirmed using a standard digital micrometre (Mitutoyo-293-521, Tokyo, Japan). Second, samples were evaluated macroscopically by the first author (SN) according to the International Cartilage Repair Society (ICRS) classification 40 . Due to this study's focus on differentiating intact from degenerative cartilage, only ICRS grades 0 (normal), 1 (superficial lesions) and 2 (fraying and lesions extending < 50% depth) were included, while more degenerative cartilage was discarded. Third, chondral samples were placed into the standard defect created within the artificial cartilage layer of the medial femoral condyle of the compressive loading device.
Based on earlier comparable studies 16,19 , sample size was estimated using a dedicated online tool (http://www. statstodo.com [Sample Size for Differences in Measurements Between Unpaired Groups Tables]). Minimum sample size of 34 was determined with the following parameters: statistical power 0.9, type-I-error 0.05, and effect size [mean of paired difference (to be detected)/expected standard deviation of paired difference] 0.8. To avoid sample pooling, one sample from each patient was harvested, i.e. a total of 49 samples.

MRI measurements.
After positioning the sample-loaded device centrally within a clinical 3.0 T MRI scanner (Achieva, Philips, Best, Netherlands), imaging was performed using standard general-purpose coils (Sense-Flex M, Philips) attached to the device's upper and lower surfaces (Fig. 5c). The coils' position was centered around the medial femoral condyle for maximized signal-to-noise ratio. MRI measurements were performed at three displacement positions as described above: a) unloaded (δ 0 ), b) at δ 2.5 , and c) at δ 5.0 . Measurements at δ 0 confirmed proper sample positioning and were used to determine reference qMRI parameter values. At each displacement position, stable sample position and proper displacement to δ 2.5 and δ 5.0 were confirmed using proton density-weighted (PDW) sequences acquired in the coronal, sagittal and axial plane (Fig. 2a-c). Sagittal and axial views were used to guide sections along the mid-coronal plane, thus creating a centrally bisecting section through the sample. Sequentially, T2*, T2, T1ρ and T1 sequences were obtained in this order; the sequence details are given in Table 4. Once the desired displacement position was set, an equilibration period of 5 min was observed prior to scanning. Using the inbuilt digital caliper tool of the picture archiving and communication system (PACS, Philips), the chondral sample's height and width were determined at the sample center on mid-coronal PDW images. Measurements were undertaken at room temperature, which was monitored during the measurements (19.3 ± 0.4 °C).

MRI data extraction.
After importing the MR raw data into Matlab R2016a software (Matlab, Natick, MA, US), respective time constants for each pixel of the mid-coronal image were calculated using predefined mono-exponential fitting routines as before 7,19 . Spatially resolved quantitative T1, T1ρ, T2, and T2* maps were generated for each sample and displacement position. For fitting, all values were included for the T1 and T1ρ maps, while the first echo and echo times >60 ms were disregarded for T2 and T2* because of potential fitting inaccuracies 41 and too low signal-to-noise ratios. R² statistics adjusted to the degrees of freedom were used to check fit quality. Sample outlines were segmented manually based on the mid-coronal PDW image by choosing pixels that safely lay within the cartilage sample (Fig. 2d 1 ). Correspondingly, boundary pixels (at the surface or bottom) were excluded to prevent partial volume effects. Segmented outlines were validated against T1, T1ρ, T2 and T2* maps. For zonal analysis, sample outlines were partitioned into two equal zones, i.e. the superficial (sf) and deep zone (dp), which were defined as the sample halves at the cartilage surface or subchondral bone, respectively (Fig. 2d 2-3 ). Thus, the entire cartilage sample (ECS) with its cartilage zones (sf, dp) constituted the (2019) 9:5895 | https://doi.org/10.1038/s41598-019-42543-w www.nature.com/scientificreports www.nature.com/scientificreports/ ROIs for which mean qMRI parameter values were calculated. For each sample and displacement position, ROIs were individually defined.
Biomechanical analyses. After MRI measurements, the chondral sample was retrieved to undergo biomechanical testing as before 22 : unconfined compression, compression rate: 0.005%/s, maximum strain: 21%. A relatively low strain rate was chosen deliberately to assess the contribution of the ECM, which bears more than 80% of the applied load at a strain rate of 0.005%/s 42 . Load-displacement data were recorded and YM (defined as the stress-strain ratio) was determined by fitting a tangent to the strain range of 10-20% 43 . Throughout, samples were kept hydrated.
To assess the diagnostic performance of individual qMRI parameters, their relative changes and combinations in the differentiation of int and deg cartilage were classified as intact (i.e. belonging to the int group) or degenerative (i.e. belonging to the deg group) based on the respective qMRI parameter values. Threshold values for individual qMRI parameters were chosen based on the respective parameters' ranges determined for int samples (Table 1): lower limit, M int − SD int ; upper limit, M int + SD int . Scatter plots of data (i.e. qMRI parameter values vs. Mankin sum scores) were segregated into true positive, false positive, true negative, and false negative. The attribute true or false was determined by belonging to the int group as determined histologically, while the attribute positive or negative was determined by the respective qMRI parameter value. Once segregation was completed, the sensitivities and specificities of the respective qMRI parameters, their relative changes and several representative combinations was calculated. To maximize sensitivity, the "Believe the positive" rule was applied, i.e. a  Table 4. Acquisition Parameters of MR sequences. n/a -not applicable.