Spatial mapping of proteoglycan content in articular cartilage using near-infrared (NIR) spectroscopy

: Diagnosis of articular cartilage pathology in the early disease stages using current clinical diagnostic imaging modalities is challenging, particularly because there is often no visible change in the tissue surface and matrix content, such as proteoglycans (PG). In this study, we propose the use of near infrared (NIR) spectroscopy to spatially map PG content in articular cartilage. The relationship between NIR spectra and reference data (PG content) obtained from histology of normal and artificially induced PG-depleted cartilage samples was investigated using principal component (PC) and partial least squares (PLS) regression analyses. Significant correlation was obtained between both data (R 2 = 91.40%, p<0.0001). The resulting correlation was used to predict PG content from spectra acquired from whole joint sample, this was then employed to spatially map this component of cartilage across the intact sample. We conclude that NIR spectroscopy is a feasible tool for evaluating cartilage contents and mapping their distribution across mammalian joint.


Introduction
Articular cartilage is a connective tissue covering the ends of diarthrodial joints with specialized biomechanical functions of load-bearing and stress transmission to the underlying bone. The process of articular cartilage degeneration in the development of osteoarthritis (OA) is generally characterized by superficial collagen disruption and proteoglycan (PG) loss [1][2][3]. Several regenerative and repair techniques [4] have been proposed for treating this condition; however, the choice of treatment and effectiveness of surgical outcome depends strongly upon the availability of accurate sub-surface matrix information that transcends visual appraisal. This level of information is necessary to determine the extent of tissue degeneration from the lesion sites in order to optimize the amount of tissue removed during repair surgery. This is further compounded by the subjective nature of cartilage evaluation during surgery as evident in the poor interobserver reliability of articular cartilage grading during arthroscopic surgery [5,6].
To address these issues, previous research has focused on imaging cartilage to determine the content and distribution of the major matrix components using established methods such as quantitative magnetic resonance imaging (MRI) [7][8][9] and mid-infrared (mid-IR) spectroscopy [10,11]. MRI is non-rapid and relatively expensive; while mid-IR spectroscopy has limited depth of penetration (in the microns), making it impractical for deep tissue probing. Furthermore, these methods cannot provide the level of accuracy and reliability possible with histological evaluation of cartilage biopsies, a method also limited due to its destructive nature. Therefore, a new approach is needed to objectively assess and visualize articular cartilage properties non-destructively during surgery. Optical imaging approaches are emerging as promising non-destructive, real-time and high resolution modalities for cartilage evaluation, and prominent among these techniques are near infrared spectroscopy (NIR) and optical coherence tomography (OCT). This study focuses on the potential of NIR spectroscopy to provide quantitative information on articular cartilage structure, with regards to proteoglycan content, and how this could facilitate spatial imaging of the tissue during arthroscopy.
Near infrared absorption characteristics of biological tissues arise predominantly from CH, NH, OH and SH bonds [12], which reflect both micro-and macroscopic properties of the tissue. Therefore, the spectral response of articular cartilage incorporates latent information on structural and functional characteristics of its matrix. Furthermore, NIR is non-destructive and rapid, penetrates deep into biological tissues [13] (about 8.5 mm into neonatal head), beyond the thickness of articular cartilage [14,15], and needs no contrast agents or sample preparation. The capacity of this method to characterize cartilage integrity and properties has been demonstrated [14-22].
The primary purpose of this study is to demonstrate the potential of NIR spectroscopy to facilitate fast quantitative evaluation and spatial mapping of articular cartilage PG content. Hence, we hypothesized that there is a relationship between the near infrared (NIR) absorption spectra of articular cartilage and its PG content, and this correlation can be employed to quantitatively image the tissue. To test this hypothesis, we investigated and validated the relationship between articular cartilage NIR absorption spectrum and its PG content, the resulting relationship was then adopted to quantitatively map the PG content in whole patella sample.

Study protocol
This study includes two objectives: i. investigating the relationship between articular cartilage NIR absorption spectra and its PG content, and ii. employing the resulting correlation to spatially map this matrix component on whole patella sample. The methodological approach involved in undertaking these objectives is summarized in the flow diagram of Fig. 1.

Sample preparation
Bovine patellae samples (N = 6, 5 with visually normal cartilage surface and 1 with a visible defect), harvested from prime oxen within 24 hrs of slaughter, were used in this study. The samples were wrapped in 0.15 M saline soaked towels and stored at −20°C until required for testing. Prior to testing, the visually normal patellae were thawed in 0.15 M saline at room temperature for 6 hrs, then osteochondral specimens (n = 20, l x b x h = 15 x 15 x 10 mm) were extracted from the samples. The remaining patella sample (with defect) was set aside for rapid imaging of PG distribution. All tests were conducted with the specimens fully hydrated in 0.15 M saline. NIR spectral data were acquired from the normal samples before enzymatic degradation to remove PGs.

Proteoglycan depletion protocol
Prior to PG depletion, thin lateral surface-to-bone osteochondral slices (for histological evaluation) were extracted from the normal samples. The samples were then incubated in trypsin (T4667, Sigma Aldrich, Sydney, Australia) of 1mg/ml concentration in 0.15 M phosphate buffered saline solution at 37°C to remove PGs. The samples were briefly removed from the solution at intervals of 1hr and NIR spectral data were acquired from them, followed by further extraction of thin osteochondral slices for histological evaluation. This protocol was performed for a period of 4hrs and 7µm osteochondral sections were prepared from the slices for histological evaluation of PG content.

Near Infrared spectroscopy -instrumentation and data acquisition
Near infrared (NIR) spectroscopy was performed in the full NIR range 12500 -4000 cm −1 wavenumber (800-2500 nm wavelength) using a Bruker MPA FT (Fourier Transform) NIR spectrometer (Bruker Optics, Germany) fitted with a diffuse reflectance fibre optic probe with 5 mm diameter window. Instrument trigger and data acquisition was controlled via OPUS 6.5 software (Bruker Optics, Germany). Precautions from preliminary experiments and previous studies [14] were observed to ensure accuracy and repeatability of acquired data. Prior to data acquisition, a reference spectrum was taken from a reflectance standard (99% Spectralon ® , Labsphere Inc., North Sutton, USA) embedded in the spectrometer, for converting raw data into absorbance spectra. Spectral data were obtained over the full range at 16 cm -1 resolution, with the output of each scan consisting of 64 co-added spectra. A single scan was taken from the center of each specimen. It was also ensured that further analyses were conducted on tissue extracted from the same region as that exposed to NIR.
To obtain spectral data for mapping the distribution of PG content across whole joint sample, the excluded patella sample was placed on a 14 x 14 grid (grid size = 5 x 5 mm). Each point on the patella corresponding to a grid location was probed. The acquired spectral data were reserved for analyses to determine the PG content at each point using the optimal correlation model.

Histological staining protocol and proteoglycan grading
The sections were fixed on microscopic slides and stained with safranin-O using a previously published staining regime [23]. Optical absorbance provided visualization and qualitative indicator of the PG content and distribution from surface to bone. This method was based on the stoichiometric binding of safranin-O to chondroitin 6-sulphate and keratan sulphate that has previously been validated by biochemical comparison in solution [24] and in cartilage tissue sections [25,26]. Linearity of this relationship between dye and PG was maintained through alcohol fixation to transform safranin-O to orthochromatic form [24,26]. This was followed by absorbance profiling under monochromatic light source using a Nikon Labo-Phot light microscope to obtain images of each sections before and after staining.
Image analysis to generate the absorbance-depth profile, i.e. the intensity and distribution of PG in the tissue, was performed on the images using custom-written macro  Table 1. Since samples exposed to trypsin would experience different amounts of PG loss depending on their initial PG concentration, the samples were grouped according to the amount of PG lost (or left) instead of the time of exposure. Grade 0 represents normal and fully stained micrographs, while grade 4 represents severe PG loss with minimal stain.

Multivariate statistical approach
To investigate the influence PG content on articular cartilage spectrum, principal component analysis (PCA) [29] was employed. PCA reduces data dimensionality to a set of linearly uncorrelated composite variables (principal components-PC) without losing too much information, and the first few PCs explain most of the variation in the original data set. PC scores, which are linear combinations of the original variables in the principal component space, were obtained and used to characterize the samples relative to their PG content. The spectral data (predictor variables) for all samples were correlated with their corresponding PG content data (response variables) using partial least squares (PLS) regression [30]. To optimize multivariate correlation between predictor and response variables, preprocessing algorithms including straight line subtraction, constant offset elimination, and derivative pretreatment, to correct spectral non-linearities and baseline offsets resulting from light scattering variations in reflectance spectroscopy [31], were adopted. Analyses were performed with and without pre-processing for comparison.
In addition to pre-processing, correlation analyses was performed on sub-regions of the spectrum (with and without pre-processing (Table 2)) using the single y-variable PLS (PLS1) algorithm. Leave-one-out (LOO) cross-validation method was used to estimate prediction error and coefficient of determination for model selection, and a maximum of 10 PLS components were tested. The left out sample is predicted with the model and the procedure repeated with each sample being left out of the calibration set. The correlation coefficients between the predicted and measured values are then calculated during calibration and validation. To avoid under-or over-fitting, while accounting for as much variation in the predictor variables, optimal model selection was based on a balance between the highest coefficient of determination (R 2 ), lowest root mean square error of cross-validation (RMSECV), and minimum PLS factor (to avoid over-fitting) ( Table 2). PCA and PLS spectral analyses were performed using OPUS Quant2 software (Bruker Optics, Germany).

Histological grading of proteoglycan content via depth loss of SO staining
Using the grading description in Table 1, the samples were categorized into groups based on the absorbance-depth profile of their safranin-O stained sections obtained from image analysis. Typical profiles for the different groups are presented in Fig. 2(a). The number of samples in the different PG content group is presented in Table 1.

Spectral analyses: classification (PCA) and prediction (PLS)
The effect of PG content on the NIR spectra of the cartilage samples ( Fig. 2(b)) can be observed from the PCA scores plot of the first and second principal components, i.e., PC1 vs. PC2 (Fig. 3). The plot of these two components provided the best distinction between group memberships. It revealed sample clustering into two distinct classes along the 1st PC axis: "class 1", on the right hand side of the plot, consisting of samples with zero to mild PG loss (groups 0-1); and "class 2", on the left hand side of the plot, consisting of samples with moderate to severe PG loss (groups 2-4). A sub-group within "class 2" consisting mostly of groups 2 and 3 PG content samples (demarcated by the ellipse) can also be observed in the scores plot. The samples within the same class are distributed along the 2nd PC. The results of PLS regression analysis on sub-regions of the NIR spectral data, with and without pre-processing, are presented in Table 2. Only the best two pre-processing methods are presented in the table. Interestingly, optimal correlation was obtained without pre-processing of the spectral data in the combined sub-regions of 5,440-6,540 and 9,959-12,489 cm −1 .

Table 2. Assessment statistics of articular cartilage NIR to PG content grade correlation based on specific regions of the spectrum with and without preprocessing. RMSECV = root mean square error of cross-validation [ † indicates the region comprising the water peak of articular cartilage]
Regions Range (cm -1 )

Straight Line Subtraction No preprocessing
Calibration The significantly high correlation between the predictor (spectral data) and response variables (PG content data) can be observed in the calibration and validation plots based on the optimal model ( Fig. 4(a) and 4(b), respectively). A comparison of the PG content grading data from assessor 1 (AS1) and assessor 2 (AS2) shows good agreement (Fig.  4(c)). In addition, high correlation was obtained between the grading results of AS2 and NIR predicted PG content data modeled on the results of AS1 (Fig. 4(d)). Fig. 4. Partial least squares calibration (a) and validation (b) plots for correlation between near infrared (NIR) spectral data and PG content for normal and artificially degraded cartilage samples. Agreement between AS1 and AS2 grading results (c), and validation of NIR-predicted PG content data against results of AS2 (d).

Spatial mapping of proteoglycan content in whole patella
The model that optimizes the relationship between NIR absorption spectra and PG content data (based on region F in Table 1) was employed to predict PG data from the NIR spectra acquired from multiple locations on the whole patella sample. Predicted PG data for all points on the sample were converted to percentage PG loss and formatted as a 14 x 14 matrix with each data on the matrix corresponding to grid points on the sample. Using standard imaging function in MATLAB® (Mathworks, Natick, MA), a 14 x 14 pixel resolution map of the PG content distribution across the patella sample was constructed from the formatted data ( Fig. 5(b)). The low spatial resolution of the resulting image is evident due to the low number of pixels because spatial resolution depends on the number of pixels utilized in the construction of a digital image; higher spatial resolution would require greater number of pixels. Bilinear interpolation algorithm with 8 iterations was used to increase the number of pixels by a factor of 238 (yielding a 3329 x 3329 pixel resolution image), thus improving the spatial resolution of PG content map across the whole patella sample (Fig. 5(c)).

Discussion
In this study, we first investigate and validate the relationship between the PG content of articular cartilage and its NIR absorption spectrum. We then demonstrate, for the first time, the potential of NIR spectroscopy to spatially map (image) cartilage component in whole joints. The relationship between the spectral data and PG content is consistent with our previous results [17] where we demonstrated the capacity of NIR spectroscopy to accurately predict changes in the functional response of articular cartilage with gradual (but controlled) PG removal. Furthermore, in comparison to earlier studies in the literature where NIR spectroscopy was applied for evaluating articular cartilage integrity in both experimental [14-20] and clinical studies [21,22], the outcomes of this study constitute a significant improvement in the adaptation of this technique for cartilage assessment.
The clustering of samples into two classes observed in the PCA scores plot (Fig. 2), one class consisting of samples with grades 0-1 PG content, and the other class consisting of samples with grades 2-4 PG content, suggests that the 1st PC (PC1) is strongly associated with group membership, and hence PG content of the tissue. This suggests that the NIR spectral data is effective in distinguishing cartilage samples with mild PG loss from those with moderate to severe loss. At closer inspection of the scores plot, the subgroup within "class 2" (encompassed within the ellipse and mostly consisting of samples with grades 2-3 PG content) suggests that this optical technique may be able to distinguish between moderate (grades 2-3 PG content) and severe (grade 4 PG content) PG loss, although more data would be required to clearly achieve this.
The models developed from data from in regions A and B (Table 1) present significantly high correlations with moderate errors, with and without pre-processing. Absorption in the wavelength range of region A is associated with third overtone CH n and ROH vibrations which characterizes cartilage proteoglycans, and NIR light penetrate deeper into articular cartilage matrix and into the underlying bone [15], allowing fulldepth tissue probing. However, since absorption and scatter is generally lower in these regions, caution must be exercised when using data from this wavelength ranges for prediction analyses.
On the contrary, models based on spectral data from regions C and D offer poor correlation and high validation errors, particularly with pre-processing. The weak correlation observed in the wavelength range of region C is arguably due to saturation of the NIR light by OH bonds in water which is known to exhibit strong NIR absorptions in this range [32]. Furthermore, the high error values of models based on this region is likely to negatively influence the accuracy of any model developed using the whole spectrum for analysis. The weak correlation and high error observed in region D, which is characterized by first overtone CH n and SH bond vibrations arguably arising from cartilage PGs, is unclear. However, a lack of full depth light penetration into articular cartilage in this wavelength range, and significant absorption from the collagen fibres, are plausible reasons for this poor relationship.
The combination of spectral data in the wavelength range spanned by regions A and D yielded the optimal model with low validation errors, with and without pre-processing. This suggests that spectral data in the combined region accounts better for the effects of PG on cartilage optical response than data from the individual regions, as evident in the models developed based on the individual region. It is interesting to note that the optimal model was achieved without pre-processing (Table 2), and this is likely due to minor loss of information during pre-processing. In addition, predictor and response variable centring and standardization prior to PLS analysis available in some multivariate statistical packages, including the one used for spectral analyses in this study, seem to be sufficient pre-processing for optimizing. The significantly high correlations (Table 2; Fig.  4(a), 4(b)) demonstrates the capacity of this method to accurately estimate cartilage PG content, and the agreement between the predicted PG content and results from assessor 2 (AS2) (Fig. 4(d)) further validates this. Restricting the analyses to only the necessary wavelength ranges of the spectrum presents a more efficient option, both in terms of data acquisition and computation, to using the whole spectrum, since only sections of the full data is required for prediction. This translates to lower latency in converting data to useful information necessary for diagnosis and decision-making.
With the spatial mapping capability demonstrated in this study (Fig. 5(c)), NIR spectroscopy could be a potentially useful tool for in vivo mapping of the components and integrity of cartilage at desired portions of whole samples. In particular, this could be beneficial in demarcating areas of damaged tissue around cartilage lesions, thus providing information that can help optimize the amount of degenerated tissue removed during surgery. In addition, quantitative information on tissue surrounding the defect (Fig. 5(c)), which is not apparent from a normal arthroscopic view (Fig. 5(a)), can be obtained to enhance decision-making. More so, this could enable patient-specific treatment, indicating whether total-, partial-replacement, or regenerative medicine (such as localized cell-based treatment) is adopted, with consequence for cost and patients rehabilitation. In addition to PG content, parameters such as thickness [14], functional properties like reswelling [17], and even subchondral bone parameters like mineral density [15] can also be mapped using this protocol; however, the protocol would need to be optimized for real-time application, and this is easy to achieve. It is worth noting that unlike current non-invasive techniques such as X-ray and MRI, our research has so far indicated that direct access and placement of the diffuse reflectance NIR probe on the surface of cartilage is important for optimal spectral response and accurate tissue characterization. Nevertheless, due to the unique depth-penetrating capacity of NIR into soft biological tissues compared to other optical techniques, it can be inferred that the system can be used as a post-surgery tool and also during follow-up arthroscopy to evaluate the integrity/health of regenerating cartilage after repair surgery, using an approach similar to that adopted by Khalil et al. [33]. However, further studies to investigate the effect of other soft tissues, such as skin and adipose tissue, on the resulting spectral output would be required before this method can be reliably applied non-invasively. This would be necessary to isolate the spectral effect of overlying tissues from the optical response of cartilage.
Although the reference data used in developing the prediction models in this study were based on artificially degraded tissue samples, the performance of the optimal model in predicting the PG content for spatial mapping of its distribution on the whole patella sample was sufficiently high to demonstrate the concept. However, calibration model based on reference data from normal and naturally degenerated human samples, where changes in matrix components could be non-linear, would provide a first step towards optimizing and validating this technique for clinical application in humans. The limitation encountered in this study, i.e. probe vibration during data acquisition, is more pronounced when acquisition time is long and this is instrument-based. While this was controlled during experiment, hardware and software stabilization techniques are possible ways of addressing this issue. Furthermore, with recent advances in NIR spectroscopy developments, spectrometers equipped with passive linear variable filter (LVF) allow rapid data acquisition (in milli-seconds), faster than current spectrometer with moving parts, such as the Fourier Transform (FT) based device used in this study.

Conclusion
We conclude that NIR spectral data of articular cartilage correlate significantly with its proteoglycan content, and this relationship can be adapted for non-destructive and quantitative spatial imaging of the tissue properties. Thus, NIR spectroscopy is a practical tool for evaluating and mapping the properties, and hence integrity, of articular cartilage with the potential to augment traditional arthroscopy.