Imaging of Osteoarthritic Human Articular Cartilage using Fourier Transform Infrared Microspectroscopy Combined with Multivariate and Univariate Analysis

The changes in chemical composition of human articular cartilage (AC) caused by osteoarthritis (OA) were investigated using Fourier transform infrared microspectroscopy (FTIR-MS). We demonstrate the sensitivity of FTIR-MS for monitoring compositional changes that occur with OA progression. Twenty-eight AC samples from tibial plateaus were imaged with FTIR-MS. Hyperspectral images of all samples were combined for K-means clustering. Partial least squares regression (PLSR) analysis was used to compare the spectra with the OARSI grade (histopathological grading of OA). Furthermore, the amide I and the carbohydrate regions were used to estimate collagen and proteoglycan contents, respectively. Spectral peak at 1338 cm−1 was used to estimate the integrity of the collagen network. The layered structure of AC was revealed using the carbohydrate region for clustering. Statistically significant correlation was observed between the OARSI grade and the collagen integrity in the superficial (r = −0.55) and the deep (r = −0.41) zones. Furthermore, PLSR models predicted the OARSI grade from the superficial (r = 0.94) and the deep (r = 0.77) regions of the AC with high accuracy. Obtained results suggest that quantitative and qualitative changes occur in the AC composition during OA progression, and these can be monitored by the use of FTIR-MS.


Results
Cluster analysis. After careful preliminary analyses of the large data matrix (including all the samples), five clusters were found to be optimal for revealing the layered structure of AC. Different spectral regions were evaluated by clustering to expose the layered structure of AC tissue and to explore the changes in the resulting cluster images during OA progression. The layered structure of AC in each sample was observed in the cluster images using the carbohydrate region for clustering (Fig. 1A). The lowest intensity was observed in the average raw spectra of the green cluster in the carbohydrate region, which indicates the lowest relative PG content (Fig. 1B). Major spectral differences between the clusters were observed in the spectral region from 1058 cm −1 to 1080 cm −1 in the second derivative spectra when the carbohydrate region was used for cluster analysis (Fig. 1C). In Fig. 2, Safranin O-stained sections, cluster images that were obtained using the carbohydrate region and PLM images of three representative samples with different OA grades are presented (Fig. 2). The corresponding images for all 28 samples can be found in the Supplementary Material. Similar features can be observed in the images. For example, the superficial, the middle and the deep layers of the AC are visible in PLM and cluster images. Furthermore, the green cluster in the cluster images complies with the Safranin-O image. On the surface of the Safranin-O images, the amount of staining is low. Correspondingly, the green cluster in the cluster images is located in the superficial layer of the AC (Fig. 2). A visual correlation between the OARSI grade and the cluster image obtained from the carbohydrate region was not observed. Similarly, when using the extended carbohydrate region, the combination of the amide I and amide II regions, or the combination of the amide I, amide II and carbohydrate regions, there was no clear correspondence between the cluster images and the OARSI grade.
Partial least squares regression. PLSR predicted the OARSI grade more accurately from the second derivative spectra (Fig. 3A,B) (r Surface = 0.94, p < 0.0001, r Deep = 0.77, p < 0.001) than from the raw spectra (Fig. 3C,D) (r Surface = 0.78, p < 0.0001, r Deep = 0.66, p < 0.001) in both the surface and the deep layers of AC (Fig. 3). The percent errors for the second derivative spectrum models were also lower (E surface = 11.51%, E deep = 17.66%) than for the raw spectrum models (E surface = 18.30%, E deep = 21.25%). Overall, the OARSI grade was more accurately predicted from the surface layer than from the deep layer of AC. The number of latent variables (LVs) varied from two to six between the models. With the spectra from the surface, six and five LVs were used in the models built from the second derivative and the raw spectra, respectively. Correspondingly, five and two LVs were used in the models in the deep layer. The CARS algorithm selected eight optimal variables (wavenumbers) for the PLSR models in both the surface and the deep layer when the raw spectra were used. When the second derivative spectra were used, 21 and 20 variables were selected to the models for the superficial and the deep layer, respectively. The variables are marked with black dots in the spectra of corresponding model (Fig. 3A-D).
Univariate analysis. Statistically significant correlations were not observed between the OARSI grade and the amide I region, the carbohydrate region, the carbohydrate region normalized with the amide I region or the second derivative peak 1064 cm −1 (PG) in either of the analyzed zones (surface and deep). However, statistically significant negative correlations were observed between the OARSI grade and the collagen integrity parameter in both the surface and the deep zones (surface: r = − 0.55, p < 0.003, deep: r = -0.41, p = 0.03) (Fig. 4A,B). A statistically significant negative correlation (r = −0.43, p = 0.02) was also found between the OARSI grade and the second derivative peak at 1202 cm −1 (collagen) in the surface layer of AC (Fig. 4C). Statistically significant differences were found between the OA groups (group 1 = early OA, group 2 = intermediate OA and group 3 = severe OA) in all univariate parameters (amide I region, carbohydrate region, collagen integrity and the second derivative peak 1202 cm −1 ) (Fig. 5). In the amide I and carbohydrate regions, statistically significant differences were observed between intermediate OA and severe OA in the surface and middle layers (Fig. 5A,B). Collagen integrity was significantly different between early OA and severe OA in the surface layer and between early OA and intermediate OA in the deep layer (Fig. 5C). The peak at 1202 cm −1 exhibited significant differences between intermediate OA and severe OA in the surface layer and between early OA and severe OA in the middle layer (Fig. 5D).

Discussion
In the present study, human AC samples with various histological grades of OA were studied using FTIR-MS. The aims of the study were to detect quantitative and qualitative osteoarthritic changes in the FTIR spectra of AC. We The second derivative spectra of each cluster. Major differences between the second derivative spectra can be seen around the peak at 1060 cm −1 . The colors of the spectra correspond to the cluster colors.
hypothesized that the layered structure of AC could be observed from the cluster analysis images and that the particular cluster will concentrate on a certain OARSI grade. The results obtained indicate that the structural and compositional changes caused by OA can be observed with both the uni-and multivariate analyses of spectral data. The present study is the first to conduct cluster analysis on FTIR images of multiple human AC samples simultaneously, i.e., by combining hyperspectral data from several individual samples into one large data matrix. This is advantageous because it allows direct comparison of the clusters between the samples. Unfortunately, compared to our original hypothesis, cluster analysis did not reveal any OA-related clusters. However, the PLSR models demonstrated that OA-related changes occur in both the surface and the deep layers of AC.
The layered structure of AC in cluster images is similar to the structure observed in PLM orientation images. Based on the visual evaluation of cluster images from all the spectral regions, the best correspondence with PLM images was observed when the carbohydrate region was used. It can be assumed that the relative amounts of collagen and PGs are the main differences between the clusters in the carbohydrate region. For example, the green cluster is found in the superficial layer of AC in the cluster images (Fig. 1A). The AC surface contains low amounts of PGs. Furthermore, the mean spectra of clusters (Fig. 1B) demonstrate that the overall absorbance in the carbohydrate region is at its lowest in the green cluster. Therefore, it can be assumed that the green cluster indicated the  regions with the lowest relative amounts of PGs in the sample sections. In contrast, by examining the overall absorbance of the mean spectra in the carbohydrate region (Fig. 1B), the relative PG content is at its highest in the grey cluster that is found in the deep zone. The other three clusters lie between the grey and green clusters in terms of their relative PG content (Fig. 1B). In the second derivative spectra, the main differences are observed in the spectral region around the peak at 1060 cm −1 (Fig. 1C). The vibrations from this region have been shown to arise from the C-O stretching vibrations of carbohydrate residues in collagens and in PGs and from − SO 3 stretching vibration in sulfated glycosaminoglycans [46][47][48] . Because Safranin-O is a cationic dye, it binds to the negatively charged PGs, and therefore, it can be used for observing the distribution of PGs 49 . When the cluster image is compared with the Safranin-O stained sections, the green cluster can be observed to concentrate in the regions where the amount of the Safranin-O stain is low. This supports the idea that the clustering result follows the relative amount of the PGs.
Because the layered structure of AC is traditionally defined by the organization of the collagen network, the amide region was presumed to reveal the layered structure more clearly than the other spectral regions that were investigated. The characteristic layered structure of healthy rabbit and bovine AC could be revealed with clustering algorithms using the amide I, amide II or carbohydrate region 40 . In earlier studies, the clustering results followed the concentration of AC constituents and the organization of the collagen network. In our study with human AC, the layered structure was not evident in all samples when using the amide I or amide II regions alone. One factor that affects the amide I and amide II regions is the traces of water vapor in the measurement chamber. Water vapor displays several sharp bands that overlap the amide I and amide II regions 51 . A second factor that may explain the difference is the different species investigated (rabbit vs. human). Finally, in the study by Kobrina et al. 40 , the clustering was performed on one sample at a time, whereas in this study clustering was performed simultaneously on all samples.
Compared to our original hypothesis, cluster analysis did not detect a cluster that is directly related to osteoarthritic changes in AC. K-means clustering is an unsupervised algorithm, i.e., it was not specifically taught to find osteoarthritic changes, but it groups the spectra according to their similarities based on Euclidean distance. The clusters that are detected seem to be related to variations in the relative PG and collagen contents in the AC. It is possible that clusters associated with OA would be found if the number of clusters would be significantly increased. However, the interpretation of all clusters would become more difficult as their number increases; therefore, based on our preliminary analyses it was not seen feasible in this study.
Although our clustering results did not correlate strongly with the OARSI grade, the PLSR analysis demonstrated that the OARSI grade can be predicted from the FTIR spectra. It is often suggested that the first signs of OA in AC are observed in the superficial layer 5,6 . A competing view suggests that OA could be initiated in the subchondral bone, as during OA changes also occur in the subchondral bone 52 . In our study, the OARSI grade was predicted more accurately from the superficial layer than the deep layer (as indicated by both the correlation and the percentage error analyses). This supports the idea that the most significant biochemical changes occurred in the surface of the AC in this relatively heterogeneous sample set. The OARSI grading system does not directly take into account the biochemical changes that occur in AC. However, as the PLSR satisfactorily predicted the OARSI grade from the deep AC layer, biochemical changes may also occur in the deep layers of AC during OA progression. Furthermore, as expected, the results of the second derivative spectra, which enhance the separation of the overlapping peaks 53 , resulted in a more accurate prediction of the OARSI grade.
With regard to univariate analysis of the FTIR spectra, no significant correlations were observed between the OARSI grade and the integrated absorbances of the amide I and carbohydrate regions, even though an earlier study suggested that quantitative correlations may exist 54 . However, a significant correlation was observed in the superficial layer between the OARSI grade and the second derivative peak at 1202 cm −1 , which can be used to monitor the collagen content 53 . The difference in the results between the amide I and the second derivative peak might be due to the better separation of overlapping peaks obtained using second derivative spectroscopy 53 . However, the lack of a correlation between the collagen content and the OARSI grade would not be that surprising, as some earlier studies did not detect a difference in the collagen content between healthy and OA human AC 55 . This suggests that the biochemical changes in AC are not linearly related to the OARSI grade. Group comparisons revealed that changes are observed primarily in the first 40% of the thickness from the cartilage surface: the collagen content, indicated by both the amide I region and the second derivative peak at 1202 cm −1 , and the PG content are significantly lower in severe OA compared to intermediate OA. Surprisingly, the collagen and PG contents appear lower in early OA compared to intermediate OA, although the differences were not statistically significant. However, if this finding were true, it may be explained by the increased levels of collagen and PG expression in chondrocytes after initial cartilage damage 56,57 .
Our results clearly indicated that the collagen integrity parameter decreases while the OARSI grade increases in both the superficial and deep zones of AC. This finding is in agreement with an earlier study, which demonstrated a statistically significant difference in the collagen integrity parameter between the Collins grade 1 (nearly normal) and grade 3 (degraded) human AC samples 29 . The collagen integrity parameter is believed to be sensitive to the triple helical structure of type II collagen. Therefore, it can be used to assess the degradation of collagen 29,58 . The enzyme collagenase degrades the collagen network during the development of OA, especially near the focal lesions 56,57 . The collagen integrity parameter has been validated earlier specifically with collagenase treated samples 59 . Therefore, the collagen integrity parameter, derived from the FTIR spectral data, may provide a measure of the denaturation level of the collagen in AC. In earlier studies, only group comparisons between the collagen integrity and OA stage have been reported 4,29,59 . In this study, a correlation between the collagen integrity and the OARSI grade was observed, which supports the sensitivity of this univariate parameter for OA changes. The collagen integrity parameter is at a similar level in early and intermediate OA in the superficial layer. However, the collagen integrity in the intermediate OA group decreases quickly to the same level with the severe OA group in the deeper layers of AC. The differences in the collagen integrity between the groups were mostly statistically insignificant probably due to the low number of samples that were evaluated. In the OARSI system, the condition of the AC surface mainly determines the OA grade of the sample. The surface layer of AC is already completely degraded in the OARSI grade of 4.0 10 . This may partially explain the lack of correlation between the OARSI grade and the amide I or carbohydrate regions in the superficial layer. Furthermore, the OARSI system does not directly take into account the changes in the collagen or PG content in AC during OA progression. Therefore, significant correlations may not be expected especially in a small sample set. In addition, our sample set consisted of samples with OARSI grades higher than 1.0, which is the threshold for OA in the OARSI grading system 10 . Therefore, we did not have completely healthy specimens in our sample set, and we could not study the very early OA changes in the AC composition. It is possible that most systematic changes in the biochemical content occur during the early stages of OA. Another limitation of the study is the relatively low number of samples (N = 28) that were evaluated. Further studies with a larger sample size are required, as multivariate methods benefit greatly from a larger sample population.

Materials and Methods
Samples. The ethical Committee of the Northern Ostrobothnia Hospital District, Oulu, Finland (191/2000 and 14-291) approved all experimental procedures, and the study was performed in accordance with the approved guidelines. Human AC samples (n = 28) were collected from patients (N = 26) undergoing total knee arthroplasty at Oulu University Hospital. Informed consent was obtained from all patients. Cylindrical osteochondral blocks were drilled from human tibial plateaus with different visual grades of OA. Subsequently, the blocks were fixed in formalin and decalcified in ethylenediaminetetraacetic acid (EDTA). After decalcification, the blocks were dehydrated and embedded in paraffin. For the FTIR-MS measurements, 5 μ m thick sections were cut using a microtome (Thermo Scientific, Microm HM 355S, Waltham, MA, USA) and, after removal of paraffin with xylene, placed on 0.5 mm thick Zinc-Selenide (ZnSe) windows (Crystran Ltd., Poole, United Kingdom). In addition, three adjacent sections were stained with Safranin-O and evaluated according to the OARSI histopathology grading system 10 . Briefly, the OARSI system divides the OA into seven main grades (0-6), where healthy cartilage is graded as 0, and the cartilage is denuded at grade 5. The highest grade 6 includes bone remodeling. In In the plots, black stars indicate the difference between the groups 2 and 3, red stars indicate the difference between the groups 1 and 3, and green stars indicate the difference between groups 1 and 2. * p < 0.05, * * p < 0.01. addition, subgrades of 0.5 were also used in the grading. Samples with an OARSI grade greater than 4 were not used in this study because the AC is already completely worn off at grade 5.

Fourier Transform Infrared Microspectroscopy (FTIR-MS). FTIR-MS measurements were conducted
using an FTIR spectrometer (Tensor 27, Bruker Inc., Billerica, MA, USA) coupled with a Bruker Hyperion 3000 microscope (Bruker Inc., Billerica, MA, USA) equipped with a focal plane array (FPA) 64 × 64 detector. The system was controlled with the manufacturer's Opus 7.2 software. First, panoramic images were collected in visual imaging mode using the same microscope. From the visual images, 500 μ m wide areas from the cartilage surface to subchondral bone (500-2000 μ m) were selected for FTIR-MS. All measurements were conducted in transmission mode, where incident infrared light is beamed through the sample section, and an FTIR absorption spectrum is obtained for every pixel. The spectral region from 900 cm −1 to 3800 cm −1 (corresponding to wavelengths from 2.6 μ m to 11.1 μ m) was collected. The spectral resolution was set to 4 cm −1 , and each spectrum was averaged 16 times. Spatial resolution was set to 21.6 μ m using pixel binning. The sample chamber was purged with dry air to stabilize measurement conditions. Data preprocessing. All preprocessing and analysis of spectral data were performed with custom-made MATLAB (MathWorks Inc., MA, USA) scripts and the commercial MATLAB toolbox Cytospec 2.00.01 (Berlin, Germany). A Resonant Mie Scattering Correction (RMieSC) algorithm was used to remove the scattering effects from the spectra 60 . Predefined thresholds for peak intensity and signal-to-noise ratios (SNR) were used as quality criteria, and poor quality spectra were removed from the dataset using the quality test functions in Cytospec. Principal component analysis (PCA) was used to detect the cartilage-bone interface. PCA is a method that creates a new set of linearly uncorrelated variables named principal components (PCs) 61 . PCs are classified so that the first PC explains most of the variance in the data and the second PC the second most and so on. PCA was performed using the carbohydrate region of the data matrix that included all the samples. A gray scale image for discriminating the bone from the AC was created using score values of the first PC. The score image was visually inspected and the bone spectra were manually removed from the dataset. After the removal of bone spectra, the spectra were truncated to a fingerprint region (950-1800 cm −1 ) that contains detailed information about proteins and carbohydrates. When calculating the collagen integrity parameter, a linear two-point baseline correction was performed on the CH 2 side-chain and the amide II regions individually. Furthermore, for the cluster analysis, the spectra were vector normalized after RMieSC. Thus, the differences in the amounts of material, i.e., the total absorbance, are removed, and clustering is based only on qualitative properties, i.e., the shape of the sample's spectra. The second derivatives spectra for the univariate and PLSR analyses were calculated using the Savitzky-Golay algorithm with nine smoothing points 62 .
Cluster analysis. Cluster analysis methods are unsupervised, i.e., no a priori information about class memberships is given. The K-means clustering (KMC) algorithm was used in this study. The basic mathematical principle of the KMC algorithm has been presented by MacQueen 63 . The KMC algorithm is a "coarse" method for classifying the spectra because a single spectrum can be a member of only one class (cluster). The algorithm classifies the spectra so that the spectral differences within clusters are minimized and the spectral differences between clusters are maximized. The Euclidean distance is used to define the differences between the spectra 38 . First, the algorithm selects a predefined number (K) of random spectra from the data set as the centroid spectra. Then, the algorithm calculates the distance between every spectrum to each centroid spectrum and classifies the spectra in K classes according to the distances. New centroid spectra are created by calculating the average spectrum for each class. Again, the distances between each spectrum to each centroid spectrum are calculated and new K classes are created. In this study, the number of these iterations was set to 20, which was visually observed to be an adequate number (no notable differences were observed in the resulting cluster images when the number of iterations was increased). As a result, the algorithm generates an image in which one color represents the particular class (cluster) 38,39 . Black pixels in the cluster images are due to the removal of poor quality spectra.
The 28 individual hyperspectral data images that were measured were combined into one large data matrix for cluster analysis. Different spectral regions can be used in KMC classifications. In this study, the following spectral regions were evaluated for cluster analysis: the carbohydrate region (985-1140 cm −1 ), the extended carbohydrate region (950-1400 cm −1 ), the amide I region (1585-1720 cm −1 ), the amide II region (1485-1585 cm −1 ), and the combination of the carbohydrate, amide I and amide II regions. Spectral differences between the clusters were investigated by observing the average spectra of each cluster. In addition, second derivative spectra of the cluster averages were also calculated.
First, the zonal structure of AC in the clustering images that were obtained was visually compared with the collagen fibril orientation images obtained using polarized light microscopy (PLM). Second, the obtained cluster images were compared with the Safranin-O stained sections to obtain information about the PG contents of different clusters.
Partial least squares regression (PLSR). PLSR is a powerful multivariate analysis method, which is widely used in chemometric analysis. The spectral peaks from different molecular vibrations tend to overlap in FTIR spectra of biological tissues. PLSR is insensitive to such collinear variables and tolerates a large number of variables. PLSR projects the spectral data onto orthogonal latent variables (LVs) 30 . In the present study, the PLSR was used to predict the OARSI grade of each sample using the FTIR spectra as predictors. The models were created separately for the spectral data from the surface layer (0-15% of the thickness from the surface towards the tidemark) and the deep layer (65-100% of the thickness from the surface toward the tidemark). Using the second derivative and the raw spectra, four models were created. Additionally, the competitive adaptive reweighted sampling (CARS) algorithm 64 was used to select an optimal combination of variables (wavenumbers) for the PLSR Scientific RepoRts | 6:30008 | DOI: 10.1038/srep30008 model. The maximum number of LVs was limited to 20. Leave-one-out cross-validation was used for validation of the multivariate models. Spearman's correlation was used to measure the agreement between the reference values and the predicted OARSI grades. In addition, the mean percent error was calculated for each model. Univariate analysis. The integrated areas under the amide I and carbohydrate regions as well as the ratio of the integrated areas under the region around the 1338 cm −1 peak and the amide II region were calculated for each spectrum for univariate analysis. Additionally, to avoid sample thickness variations, the integrated area of the carbohydrate region was also normalized by the integrated area of the amide I region. Consequently, depth-wise average profiles were calculated for every parameter. The regions of interest for the superficial and the deep layers were selected similarly to those used in the PLSR analysis. From the second derivative spectra, the peaks at 1064 cm −1 and 1202 cm −1 were used to analyze PG and collagen contents, respectively 53 . The results were compared with the OARSI grade using correlation analysis.
Polarized light microscopy (PLM). PLM was used as a reference method when studying the layered structure of AC 65 . All 28 samples were imaged using an Abrio PLM system (CRi, Inc., Woburn, MA, USA) mounted on a conventional light microscope (Nikon Diaphot TMD, Nikon, Inc., Shinagawa, Tokyo, Japan). The Abrio system consists of a green bandpass filter, a circular polarizer, and a computer-controlled analyzer composed of two liquid crystal polarizers and a CCD camera. In a PLM orientation image, each pixel consists of information about the average orientation of the collagen fibrils. Statistical analysis. The Spearman's correlation test was used for all correlation analyses. For the group comparisons, the OARSI grades were pooled into three groups. Group 1 (early OA) included the grades from 1.0 to 2.5 (n = 10), group 2 (intermediate OA) included the grades from 3.0 to 3.5 (n = 9), and group 3 (severe OA) included the grades from 4.0 to 4.5 (n = 9). The Kruskall-Wallis test with pairwise comparison was used to analyze the statistical significance between the groups. The group comparisons were performed for the amide I region, carbohydrate region, collagen integrity and second derivative peak at 1202 cm −1 . The statistical analyses for the univariate parameters were performed using SPSS 22 software (IBM SPSS Statistics, SPSS Inc., Chicago, IL, USA).