2.1.

Sample Preparation

The samples of this study were originally collected and prepared for two other studies.^16,17 Intact bovine patellae ($n=14$) were obtained from a local abattoir (Atria Oyj, Kuopio, Finland). Ethical approval was not required, as the animals were not slaughtered for the purposes of this study. The bovine knees were stored at room temperature after slaughtering. The knee joints were delivered within a few hours of slaughter and were prepared immediately after delivery. Four osteochondral samples ($\text{diameter}=6\mathrm{mm}$) were prepared from lateral upper quadrants of each patella except for one, from which only two samples were prepared. Taken together, a total of 54 samples were used in this study. Half of the samples of the patella (altogether $n=27$) were treated with threose to increase the collagen cross-linking while the other half (altogether $n=27$) were used as control samples. Control samples were incubated for 7 days at 37°C in cell culture medium (DMEM low glucose $1\mathrm{g}/\mathrm{l}$, Lonza Cologne AG, Belgium) with $100\text{unit}/\mathrm{ml}$ penicillin, $100\mu \mathrm{g}/\mathrm{ml}$ streptomycin, 10 mM HEPES buffer solution (HEPES, EuroClone S.p.A., Italy), 1 mM L-glutamine and 10 mM vitamin C. Cross-link formation was induced by adding 100 mM of threose (Sigma Aldrich Co., St. Louis, Missouri) into the incubation media of the threose group. Finally, the samples were split in half, and one half was fixed in 10% formalin for 48 h, decalcified with ethylene diamine tetraacetic acid, dehydrated and immersed in liquid paraffin. Subsequently, the paraffin was hardened by cooling. The other half was used for biochemical reference analyses.

2.2.

Infrared Microspectroscopic Measurements

For FTIR microspectroscopy, $5\text{-}\mu \mathrm{m}$-thick sections were cut with a microtome. For dewaxing, the sections were immersed in series of xylene baths ($3\times 5\min$) and subsequently washed with ethanol. The sections were placed onto 2-mm-thick ZnSe windows and were let dry in room air. The sections were measured using the Perkin Elmer Spotlight 300 FTIR imaging system (Perkin Elmer, Shelton, Colorado) with spectral resolution of $4{\mathrm{cm}}^{-1}$, pixel size of $25\mu \mathrm{m}$, and using four scans per pixel. A rectangular region of interest (ROI) extending from cartilage surface to cartilage–bone junction was imaged from each section (Fig. 1). The ROI width was set to $400\mu \mathrm{m}$, whereas the ROI height, which varied depending on the cartilage thickness, was on average $(\pm \mathrm{SD})1693\pm 484\mu \mathrm{m}$. On average, 1088 spectra were collected per sample. A dry air purge (Parker Balston, Haverhill, Massachusetts) was used to minimize the variation in the measurement conditions. The spectra of each section were first averaged to obtain one mean spectrum. Thereafter, resonant Mie scattering correction (RMieSC) algorithm^18,19 was used to remove scattering-related effects from the spectra. After RMieSC, the spectra were vector normalized. Second derivative spectra were calculated using Savitzky–Golay algorithm with 11 smoothing points. In this paper, we use the following definitions for spectral regions: the amide I (1590 to $1720{\mathrm{cm}}^{-1}$), the amide II (1500 to $1590{\mathrm{cm}}^{-1}$), the mixed region (1200 to $1500{\mathrm{cm}}^{-1}$), and the carbohydrate region (800 to $1200{\mathrm{cm}}^{-1}$).²⁰ All data analysis was done in MATLAB (R2015a, MathWorks, Inc., Natick, Massachusetts).

Fig. 1

A rectangular ROI extending from cartilage surface to cartilage–bone junction was imaged from each section. The ROI width was set to $400\mu \mathrm{m}$, whereas the ROI height, which varied according to the cartilage thickness, was on average $(\pm \mathrm{SD})1693\pm 484\mu \mathrm{m}$. Therefore, the average size of the ROI was $16\times 68\text{pixels}$, which results in $\sim 1088\text{spectra}/\text{section}$. These spectra were averaged to obtain an average spectrum for the section.

2.3.

Univariate Analysis

To evaluate the overall changes in the spectra due to the increase in cross-linking, the difference spectrum (mean spectrum of threose-treated samples minus mean spectrum of control samples) was calculated. In addition, Pearson’s correlation coefficients between each variable (i.e., absorbance at each wavenumber) in the FTIR spectra of the samples and the concentrations of HP, LP, and Pent were calculated to find the relationship between the variables of the FTIR spectrum of AC and the cross-link concentrations.

2.4.

Partial Least Squares Regression

PLS regression, a method based on latent variables, was used to predict cross-link concentrations from FTIR spectra. Spectral region of 800 to $1800{\mathrm{cm}}^{-1}$ was used for PLS regression models. Root-mean-square error of cross validation (RMSECV) was used as the term to be minimized in PLS regression. Leave-one-out cross-validation was used when deciding the number of PLS components for the full spectrum PLS regression model. All models were validated using leave-one-out cross-validation.

2.5.

Backward Iterative Partial Least Squares Regression

Backward iterative partial least squares (biPLS) is a variable selection algorithm that evaluates the importance of different spectral regions in PLS regression problems.²¹ The spectra were divided into 25 equal-sized spectral windows ($20{\mathrm{cm}}^{-1}$ window size). The effect of each spectral window to the model was tested by building the model without the window. The window whose removal improved the model the most was removed in each cycle. The biPLS procedure was continued until there was only one window remaining. Finally, the combination of spectral windows that produced the lowest RMSECV was selected as the best biPLS model. The whole procedure was conducted for PLS components from 1 to 10. The optimal number of components for biPLS model was chosen based on the minimum of RMSECV.

2.6.

Backward Iterative Partial Least Squares Discriminant Analysis

PLS is most often used for regression analysis. However, PLS can also be used for classification tasks. For this purpose, a dummy variable describing the categories is used as the response variable. biPLS-DA algorithm was applied to second derivative spectra to classify the samples into control and threose-treated groups. The parameters of biPLS-DA were set to be the same as for biPLS. A binary dummy variable with a value of zero and one for control and threose-treated groups, respectively, was used as the response variable. The samples with predicted class values of under 0.5 were assigned to group 0 (control group) and with values of over 0.5 to group 1 (threose-treated group). The optimal number of PLS components was found experimentally. The model was validated using leave-one-out cross-validation.

2.7.

Genetic Algorithm for Wavenumber Selection

Genetic algorithms (GAs) are optimization methods that mimic the process of natural selection. GA was used in two ways: directly to FTIR spectra and in combination with biPLS (biPLS-GA). In biPLS-GA, biPLS was first used to remove 18 out of 25 spectral windows to reduce the number of variables before running GA. GA was then applied to the variables that were left after biPLS.

In GA, a chromosome is a binary vector with a size that corresponds to the number of variables in the spectrum. Variables of the chromosome, genes, indicate whether (1) or not (0) the corresponding variable in the spectrum is selected to the model. The GA is run for a predefined number of generations. In the beginning of GA, a population of randomly generated chromosomes is created. A PLS regression model is built using each of the chromosomes of the population, and RMSECV is used to evaluate the performance of the variables selected by each chromosome. The best chromosome is included in the next population, whereas the other chromosomes of the population are obtained by recombining the chromosomes of the earlier population by using cross-over of chromosomes and mutation of genes. Final generation contains the best chromosome of a single GA run. GA is run multiple times, and the importance of individual variables is evaluated based on their selection frequency in the best chromosome. In this study, GA was run for 50 times. Finally, the variables are added to the final model according to their selection frequency. The variable combination that results in the smallest RMSECV is chosen as the final model. The maximum number of PLS components was limited to 10. Other parameters of GA are summarized in Table 1. A more complete description of the algorithms can be found in earlier studies.^22,23

Table 1

Parameters used in the GA for selecting spectral variables.

2.8.

Competitive Adaptive Reweighted Sampling Partial Least Squares Regression

Competitive adaptive reweighted sampling (CARS)²⁴ is a variable selection that uses the absolute values of regression coefficients as estimates for the importance of each variable. An exponentially decreasing function is used to define how many variables are retained in each of the sampling runs. The parameters of the exponentially decreasing function are selected so that all variables are retained in the first sampling run, whereas only two variables are left in the final sampling run. In every iteration, the least important variables are removed (variables with the smallest absolute regression coefficient values). A new calibration model is then built using the remaining variables, and the regression coefficients of the new model are used as new indicators of their importance. Removal of variables is continued for a predefined number of sampling runs. In this study, the number of sampling runs was set to 100. Finally, the best variable combination is searched and used for the final model, which is validated by leave-one-out cross-validation. A more detailed description of the algorithm can be found in an earlier study.²⁴ The RMSECV was first monitored as a function of the number of PLS components in CARS-PLS. However, a clear minimum was not found, as the RMSECV steadily decreased as the number of PLS components was increased. Therefore, the maximum number of PLS components was limited to 10. CARS-PLS analyses were conducted in MATLAB using libPLS toolbox.²⁵

2.9.

Biochemical Analysis

A more detailed description of the biochemical analysis is found in previous studies.^16,17,26 Briefly, the amount of the collagen specific amino acid, hydroxyproline, was measured spectrophotometrically to estimate the total amount of collagen.²⁷ HPLC was used to separate hydroxylysyl pyridinoline (HP), lysyl pyridinoline (LP), and pentosidine (Pent).²⁸ Pure compounds of HP, LP, and Pent with known concentrations were used as standards and for converting the measurement results into concentration values. The results are expressed as mol of cross-link per mol of collagen. The HPLC system used consisted of a Quaternary Gradient Pump unit, PU-2089 Plus, Intelligent AutosamplerAS-2057 Plus, and Intelligent Fluorescence Detector, FP-2020 by Jasco (Jasco Scandinavia AB, Mölndal, Sweden). Data were analyzed using Jasco Chrompass software (Jasco, Sweden). The LiChroCART^® 125-4 column was from Merck Hitachi (Merck KGaA, Darmstadt, Germany)

2.10.

Statistical Analysis

The values predicted by the PLS models were compared to the reference information by using the Pearson’s correlation analysis. Furthermore, the Spearman’s correlation analysis was also conducted as the Anderson–Darling test revealed that LP and Pent values did not follow the normal distribution. The limit for statistically significant correlation was set to $p<0.05$. All statistical analyses were conducted using MATLAB (R2015a, MathWorks, Inc.).

3.1.

Difference Spectrum and Univariate Analysis

Mean absorbance spectra of control and threose-treated samples are shown in Fig. 2(a). The difference spectrum, i.e., mean spectrum of threose-treated samples ($n=27$) minus mean spectrum of control samples ($n=27$) shows a negative peak in the amide I (1700 to $1590{\mathrm{cm}}^{-1}$) and a positive peak in the carbohydrate region (1000 to $1100{\mathrm{cm}}^{-1}$). Pearson’s correlation coefficients between the variables in the FTIR spectra of the samples and the concentrations of HP, LP, and Pent were calculated to find the best variables for univariate analysis. For HP, the highest absolute correlation coefficient is found in the carbohydrate region ($\sim 1120{\mathrm{cm}}^{-1}$) ($r=-0.31$, $p<0.05$) [Fig. 2(b)]. Some samples were excluded from the LP ($n=14$) and Pent ($n=16$) correlation coefficient analysis as the biochemical analysis did not detect any LP or Pent cross-links in some of the control samples due to the limitations of the HPLC method. The correlation coefficients of LP [Fig. 2(c)] and Pent [Fig. 2(d)] follow the shape of the difference spectrum. The highest absolute correlation coefficient values are $r=-0.45$ ($p<0.001$) in the amide II region and $r=0.50$ ($p<0.001$) in the amide I and the carbohydrate regions for LP and Pent, respectively.

Fig. 2

(a) Mean absorbance spectra of control samples (solid black line) and threose-treated samples (dashed black line). There is an added offset between the mean spectra for better clarity. The difference spectrum (multiplied by a factor of five for better visualization) is shown as a solid gray line. Pearson’s correlation coefficients between the variables in the FTIR spectra of the samples and the cross-link concentrations are shown for (b) HP, (c) LP, and (d) Pent.

3.2.

Backward Iterative Partial Least Squares Discriminant Analysis

BiPLS-DA was applied to second derivative spectra for separating the threose-treated samples from intact ones. Two PLS components were used for classification. The sensitivity (true positive rate) and specificity (true negative rate) for threose-treated samples were 92.5% and 88.9%, respectively. The accuracy (percentage of correctly classified samples) was 90.7%. The class predictions are shown in Fig. 3.

Fig. 3

The predicted class memberships obtained from biPLS-DA. The control and the threose-treated samples are marked with circles and squares, respectively. Correctly and incorrectly classified samples are displayed with solid and open markers, respectively. Two PLS components were used in the final biPLS-DA model to separate the control and threose-treated groups from each other. The samples were classified with an accuracy of 90.7%.

3.3.

Partial Least Squares Regression

The biochemical analysis did not detect any LP or Pent cross-links in some of the control samples. These samples were excluded from the PLS regression models when predicting LP or Pent concentrations. The results of raw spectrum [Fig. 4(a)] models are summarized in Table 2 and Fig. 4. In general, equal or better results were obtained when using second derivative spectra [Fig. 5(a)] instead of raw spectra. Therefore, we focus on the results obtained using second derivative spectra.

Table 2

The results in predicting cross-link (HP, LP, and Pent) concentrations using PLS regression with different variable selection techniques from the raw spectrum. The number of PLS components (Comp), the number of used variables (Var), the Pearson’s correlation coefficient (r), and the Spearman’s correlation coefficient (rs) between the predicted values and the reference values, and RMSECV are shown for each PLS method.

Fig. 4

(a) The mean raw spectrum of the dataset, and the variables selected from the spectra by the studied variable selection algorithms for the prediction of (b) HP, (c) LP, and (d) Pent concentrations, respectively.

Fig. 5

(a) The mean second derivative spectrum of the dataset, and the variables selected from the second derivative spectra by the studied variable selection algorithms for the prediction of (b) HP, (c) LP, and (d) Pent concentrations, respectively.

Four PLS components produced the best model ($r=0.41$, $p<0.01$) when the full spectrum (800 to $1800{\mathrm{cm}}^{-1}$, 500 variables) was used for prediction of the HP concentration, whereas the other models used 9 or 10 PLS components. Variable selection using biPLS improved the correlation coefficient to $r=0.70$ ($p<0.001$). The model utilized portions of the amide II region, the mixed region, and the carbohydrate region [Fig. 5(b)]. The correlation coefficient between the HP concentration and the concentration predicted by the PLS-GA regression model was $r=0.62$ ($p<0.001$). PLS-GA selected variables quite evenly throughout the spectral region of 800 to $1700{\mathrm{cm}}^{-1}$ [Fig. 5(b)]. BiPLS-GA utilized mainly the carbohydrate region and, to a smaller extent, the amide II and the mixed regions [Fig. 5(b)]. The correlation coefficient between the HP concentration and the concentration predicted by the biPLS-GA regression model was $r=0.83$ ($p<0.001$). CARS-PLS selected variables throughout the spectral region of 800 to $1600{\mathrm{cm}}^{-1}$ [Fig. 5(b)] and resulted in a correlation coefficient of $r=0.83$ ($p<0.001$) [Figs. 6(a) and 6(d)].

Fig. 6

Scatter plots between the biochemically determined concentrations and the concentrations predicted by CARS-PLS for (a) LP, (b) HP, and (c) Pent cross-links. (d–f) Corresponding Bland–Altman plots are shown below each correlation plot.

When predicting the LP concentration, two PLS components produced the best model ($r=0.20$, $p>0.05$) in the full spectrum (800 to $1800{\mathrm{cm}}^{-1}$, 500 variables) model. Also BiPLS, PLS-GA and biPLS-GA used only two PLS components. BiPLS improved the correlation coefficient to $r=0.47$ ($p<0.01$). The model relied on the carbohydrate and mixed regions [Fig. 5(c)]. The correlation coefficient between the LP concentration and the concentration predicted by the PLS-GA regression model was $r=0.52$ ($p<0.001$). PLS-GA selected variables from all regions [Fig. 5(c)]. In addition to the regions used in biPLS, BiPLS-GA used the additional spectral window of 800 to $840{\mathrm{cm}}^{-1}$ [Fig. 5(c)]. The correlation coefficient between the LP concentration and the concentration predicted by the biPLS-GA regression model was $r=0.51$ ($p<0.001$). CARS-PLS used nine PLS components for prediction of the LP concentration. CARS-PLS used variables mainly from the amide I, the amide II, and the carbohydrate regions [Fig. 5(c)]. The correlation coefficient between the LP concentration and the CARS-PLS model was $r=0.87$ ($p<0.001$) [Figs. 6(b) and 6(e)].

Nine PLS components produced the best model ($r=0.66$, $p<0.001$) in the full spectrum (800 to $1800{\mathrm{cm}}^{-1}$, 500 variables) model when predicting the Pent concentration. In cases of biPLS and biPLS-GA, the selected spectral regions were very similar [Fig. 5(d)]. PLS-GA selected a relatively high number of variables throughout the whole spectrum [Fig. 5(d)]. The correlation coefficients between the Pent concentration and biPLS, PLS-GA and biPLS-GA were $r=0.84$ ($p<0.001$), $r=0.90$ ($p<0.001$) and $r=0.93$ ($p<0.001$), respectively. CARS-PLS selected variables from all regions of the spectrum [Fig. 5(d)] and resulted in a correlation coefficient of $r=0.92$ ($p<0.001$) [Figs. 6(c) and 6(f)]. The results of all second derivative spectrum-based PLS regression models are summarized in Table 3.

Table 3

The results in predicting cross-link (HP, LP, and Pent) concentrations using PLS regression with different variable selection techniques from the second derivative spectrum. The number of PLS components (Comp), the number of used variables (Var), the Pearson’s correlation coefficient (r), and the Spearman’s correlation coefficient (rs) between the predicted values and the reference values, and RMSECV are shown for each PLS method.