Detecting human articular cartilage degeneration in its early stage with polarization-sensitive optical coherence tomography.

Detecting articular cartilage (AC) degeneration in its early stage plays a critical role in the diagnosis and treatment of osteoarthritis (OA). Polarization-sensitive optical coherence tomography (PS-OCT) is sensitive to the alteration and disruption of collagen organization that happens during OA progression. This study proposes an effective OA evaluating method based on PS-OCT imaging. A slope-based analysis is applied on the phase retardation images to segment articular cartilage into three zones along the depth direction. The boundaries and birefringence coefficients (BRCs) of each zone are quantified. Two parameters, namely phase homogeneity index (PHI) and zonal distinguishability (Dz), are further developed to quantify the fluctuation within each zone and the zone-to-zone variation of the tissue birefringence properties. The PS-OCT based evaluating method then combines PHI and Dz to provide a G PS score for the severity of OA. The proposed method is applied to human hip joint samples and the results are compared with the grading by histology images. The G PS score shows very strong statistical significance in differentiating different stages of OA. Compared to using the BRC of each zone or a single BRC for the entire depth, the G PS score shows great improvement in differentiating early-stage OA. The proposed method is shown to have great potential to be developed as a clinical tool for detecting OA.


Introduction
Osteoarthritis (OA) is the most common degenerative joint disease [1]. Symptoms of OA include debilitating pain, tenderness, and stiffness of joint caused by gradual break down of articular cartilage. Those symptoms usually occur when the disease is already in its late-stage so that the patients have no other options than pain management or surgical intervention (e.g. total joint replacement) [2][3][4]. If OA can be identified in early stage, therapeutic strategies such as lifestyle modification and pharmaceutical drugs could be applied in time, which can delay or even prevent the irreversible damage [3]. However, currently there is no effective tool to diagnose or predict early-stage OA. Traditional X-ray detects cartilage loss by measuring the narrowing of joint space [1,2], which lacks the sensitivity to detect localized cartilage damage and thus is not capable of detecting early-stage OA. At the beginning of cartilage degeneration, changes only happen at the micro-scale, such as alteration in collagen fibril network and depletion of proteoglycan content [5]. Those changes contribute accumulatively to the gradual destruction of cartilage extracellular matrix (ECM) [6]. Such subtle alternations in ECM structure are not detectable by imaging the entire joint, such as with X-Ray. Therefore, a high-resolution imaging technique is needed for detecting early-stage OA.
Optical coherence tomography (OCT) is a non-invasive, high-resolution imaging technique based on the backscattering of light [7]. OCT can characterize the morphology and layered structures of tissue. Using OCT imaging, Nebelung et al. (2014) developed an optical irregularity index (OII) for assessing the tissue surface roughness, and an optical homogeneity index (OHI) for the tissue inhomogeneity [8]. They found higher surface irregularity and lower tissue homogeneity in moderate OA than in healthy articular cartilage but failed in detecting the surface or sub-surface disruption in early-stage OA. In a later study,  (from the same group) developed several parameters by characterizing the features such as ordinates parameters and amplitude distribution of the surface profile, which showed strong and highly significant correlation with histological OA grades [9]. Their work further confirmed the close association between the loss of surface integrity and the OA grades. However, it failed to differentiate between samples with close differences of damage, e.g. OA grade 1 vs. grade 2, or grade 2 vs. grade 3. Those studies were merely relying on the back-scattering intensity signal, which had limitations in assessing collagen organizations. Thus, the detection of early-stage OA was still a challenge. While OCT imaging can detect the morphological changes in cartilage degeneration, it lacks the sensitivity for subtle alterations in collagen organization due to ECM disruption.
Polarization sensitive OCT (PS-OCT) is a functional extension of OCT based on the polarization properties of light. PS-OCT can quantify tissue birefringence which is sensitive to the orientation and organization of collagen fibers [10]. Reduced birefringence was found to be associated with articular cartilage degeneration [11][12][13]. In a series of studies Xie et al. ( , 2008 showed that birefringence was not only associated with the healthy states [14], but also sensitive to the light illumination direction [15] and tissue topographical location [16]. Contemporarily, Shyu et al. (2008) quantified a birefringence coefficient (BRC), which was defined as the slope of the phase retardation curve in PS-OCT image [17]. They found that different OA states showed different BRC distributions. However, it was still difficult to differentiate OA severities. In a later study, Brill et al. (2016) also quantified BRC and other features from the banding patterns of PS-OCT image [18], such as the width and height of the bands and the total number of bands. They compared those features in 13 samples from 4 patients at different OA stages. However, they found the difference to be non-significant, which they claimed to be due to limited samples and insufficient number of diseased cases. ECM disruption has been observed in early cartilage degeneration under mechanical compression [19]. Recently Goodwin et al. (2018), applied BRC analysis on articular cartilage with and without mechanical compression on 26 bovine samples classified into three grades (10 healthy, 8 mild, and 8 moderate degeneration cases) [20]. Their results showed that the differentiation between healthy and early degeneration was statistically significant only when the imaging was conducted in a compressed state. However, applying mechanical indention is not practical in clinical diagnosis. In all the above-mentioned studies, the variation of collagen organization in the different structural zones was not considered.
Articular cartilage is composed of four zones, which are superficial, transitional, deep, and calcified zones, respectively. The collagen organization and orientation vary from zone to zone. For example, collagen fibers are aligned parallel to the tissue surface in the superficial zone but perpendicularly in the deep zone. Tissue birefringence is highly dependent on the fiber orientation. In our previous publication, we showed that PS-OCT can detect the change in birefringence at different depths in healthy swine samples [21]. We also developed a slope-based analysis to segment articular cartilage into three regions within the penetration depth of the PS-OCT, based on the slopes at different depths of the phase retardation curve.
In this study, we apply the slope-based analysis on human OA samples to study the zone disruption in the progression of OA. While the parameters of OII and OHI was defined based on the OCT intensity image, we extend similar analysis to the phase retardation image in PS-OCT. With the slope-based analysis, we segment the articular cartilage into three structural regions that have different BRCs. As OA progresses, the structural zones become less distinguishable. Thus, a segmentation irregularity index (SII) is defined to quantify zone destruction, and a phase homogeneity index (PHI) is defined to quantify localized extracellular matrix disruption. A PS-OCT based grade (G PS ) is developed based on the phase retardation image and the slope-based analysis. It is applied to ex-vivo study on human OA samples. We classify the OA severities by use of histological observation and then use the histological results to validate G PS . The slope-based analysis offers strong power to differentiate early-stage OA from healthy tissue. The PS-OCT based evaluating method is shown to have great potential in clinical diagnosis and related pathological investigations.

Sample preparation
Six patients (four males, two females; mean age 73 years old [range from 64 to 86 years old]) who underwent total hip joint replacement were consented for this study. Three patients had severe primary OA diagnosed clinically and radiographically, and the other three patients suffered non-OA diseases such as intracapsular femoral neck fracture. All procedures were performed under approval of the ethics board of the University of British Columbia.
Cartilage-bone samples were harvested from five sampling locations, the superior (S), posterior (P), anterior (A), superior anterior (SA), and superior posterior medial (SPM). The locations were aligned as much as possible to be the central part of each divided condyles to guarantee anatomic consistency. Cylindrical sample blocks of 18 mm diameter were extracted with a drill. The tissues were kept frozen at − 20°C until preparation for imaging. Prior to imaging, the samples were defrosted and allowed to equilibrate in a saline bath for up to 2 hours. For PS-OCT imaging, the laser beam shines perpendicularly to the cartilage surface. A total of 30 sample blocks were extracted and imaged by PS-OCT. Afterwards, the cartilage blocks were processed for histology imaging. A classification of OA severities based on histology images is carried out and used as the gold standard.

Histological classification of OA severities
Following PS-OCT imaging, the specimens were processed for histological observation and imaging. The specimens were formalin-fixed for 24 hours and demineralize in 10% formic acid for 4 days. The specimens were then dehydrated in ascending concentrations of ethanol, cleared in xylene, and paraffin-embedded. After deparaffinization, sections of 6-microns thickness were stained with fast green and safranin-O [22]. The tissue slides were observed under a light microscope (Carl Zeiss, Germany) and photographed by using a cooled CCD camera (QImaging, Canada). Images were processed by Zen software (Carl Zeiss, Germany).
Analysis of the samples was performed following the Osteoarthritis Research Society International (OARSI) osteoarthritis cartilage histopathology assessment system. In this system, a "grade" between 0-6 is assigned to evaluate the index of severity or biologic progression of OA in cartilage along the depth direction, and a "stage" between 0-4 is assigned to describe the horizontal extent of the disease. OARSI score is obtained by the product of "grade" and "stage" [23]. Semi-points are used to allow the score coverage over the scoring spectrum from 0 to 24. OA classification was conducted by three independent observers, according to the OARSI scoring system [23]. One of the observers is an experienced researcher in histopathology of bone and connective tissues. The other two observers have experience in medical and biological engineering for soft-tissue related applications. The three observers scored all the images for three repeated times, where each time the order of the images was randomized. The Intra-Class-Coefficient (ICC) values were over 0.75 for both intra-and inter-validations, which showed high reproducibility [24]. Furthermore, to ensure the reliability of the severity judgment, only the samples with all the scores within a 2-SD interval (twice the standard deviation (SD) to and from the average score) were used.
The OARSI based scoring results are summarized in Table 1. Four groups were classified based on their score distribution: Score 0−4.99 was defined as healthy (group S_1), 5−9.99 as mild OA (group S_2), 10−17.99 as moderate OA (group S_3), and ≥18 as late-stage (group S_4). The grade, stage, and score (grade × stage) were averaged over all the samples in each group. In this study, the late-stage group was not considered due to the limited cases and those samples were also excluded from the statistic study. After those filtering processes, the selected number of samples in this study became 23, including 7 in S_1 group, 8 in S_2 group, and 8 in S_3 group. Noticeably, the samples labeled as mild OA (group S_2) was defined as early-stage in this study. Among the 23 samples, there are 6 in location "S", 4 in "P", 5 in "A", 3 in "SA", and 5 in "SPM". There is no obvious relationship between the OA severity and location distribution.

PS-OCT image analysis
All the samples were imaged by a custom-built Jones-matrix-based PS-OCT system [21]. The sensitivity of the PS-OCT system is 92 dB. The axial and lateral resolutions are ∼8.1 µm and ∼19.2 µm, respectively. Using the PS-OCT system, co-registered OCT intensity images and phase retardation images were acquired from 3-D volumes. Effective signal-to-noise ratio (ESNR) was used to select the imaging depth for reliable phase retardation calculation with a threshold of 10 dB [21].

Slope-based analysis of phase retardation
Tissue birefringence can be quantified by BRC, which is defined as the linear fitting slope of the accumulative phase retardation A-line. A larger slope indicates higher birefringence. Previously, a single slope has been used to fit the entire depth range of the phase retardation A-line [20]. However, as the collagen organization in articular cartilage varies in the different zones along depth, a single slope is not sufficient to represent the depth dependent characteristics [17]. In the progression of OA, the disruption of ECM progresses from the surface to the deeper zones. Thus, characterizing and differentiating tissue birefringence in different structural zones will be important.
In our previous work, a slope-based analysis has been developed to segment the accumulative phase retardation data volume into three regions, corresponding to the superficial, transitional, and deep zones [21]. The procedure is summarized briefly here. First, an A-line is averaged with its surrounding A-lines over a window size of ∼100 µm ×100 µm, and then smoothed by a running median filter. Second, a surface detection algorithm and a thresholding are performed to determine the tissue surface and the maximum reliable imaging depth, respectively. Third, a three-segment linear fitting based on least squares regression is applied to segment the A-line into three regions with different slopes. Finally, the above steps are applied on each A-line in a 3D volume. From the slope-based analysis, the boundaries and the respective slopes of the three segmented regions can be obtained. The segmentation boundaries are written as b 0 (x, y), , which represent the tissue surface, the boundary between the first and second zones, and that between the second and third zones, respectively. The slopes of the three segments are sl 1 (x, y), sl 2 (x, y), and sl 3 (x, y), which represent the BRCs of the first, second, and third zones, respectively. By using the slope-based analysis, the depth dependent collagen organization in articular cartilage can be characterized. The slope-based analysis has been demonstrated with good performance in finding the structural zones in healthy articular cartilage [21], where the three segments represent superficial, transitional, and deep zones, respectively. When applying the analysis to OA samples [25], the three-segmented regions may not represent the actual structural zones, because the regular collagen ECM structure has been disrupted. More discussions will be provided in Section 3.1.

Phase accumulative changes (PAC) and phase homogeneity index (PHI)
In our study, we found that using the boundaries and slopes of the three segmentation regions alone is still not sufficient to classify the different stages of OA. In the progression of OA, the disruption of ECM increases the heterogeneity of collagen organization compared to normal articular cartilage. Higher heterogeneity in collagen organization will result in higher fluctuations in the phase changes along the A-line. Therefore, two new parameters, phase accumulative change (PAC) and phase homogeneity index (PHI), are defined to quantify such heterogeneity properties.
PAC quantifies the accumulative phase change in phase retardation image in PS-OCT. It is similar to OHI that quantifies the intensity change in the OCT intensity image [8,18]. By calculating the number of edges along the depth profile, OHI quantifies the inhomogeneity in the intensity of the scattering signal. Higher OHI value indicates higher fluctuation along the tissue depth. Similarly, PAC quantifies the phase change along the phase retardation A-line. Phase fluctuation peaks are the peaks which have relatively high phase change, defined as higher than 20% of the total phase accumulation within a segmentation zone. A standard edge detection function in MATLAB is applied on each phase retardation A-line to detect those phase fluctuation peaks. The PAC value for each zone is the summation of the phase fluctuation peaks. The analysis is applied on each A-line in a 3D volume, and thus PAC values for each A-line in the three zones can be obtained as PAC 1 (x, y), PAC 2 (x, y), and PAC 3 (x, y), respectively. In this calculation, the A-lines are not averaged or smoothed.
PHI is a weighed sum of the PACs from the three zones. Since higher phase fluctuation is usually associated with higher averaged phase accumulation, PHI is normalized to the averaged phase accumulation. The PHI value for each A-line can be obtained by: The averaged phase accumulation is calculated as y) is the zonal BRC, and Th i (x, y) is the zone thickness obtained from the boundaries b i (x, y).
Weighting factors w 1 , w 2 , and w 3 are assigned for the first, second, and third segmentation zone, respectively. The weighting factor is proportional to the averaged ESNR such that w i = Average(E s (x, y, z))/40, where E s (x, y, z) is the 3D ESNR matrix and the average is performed over all the voxels within a segmentation zone. The weighting factor w i has important meaning in penalizing the regions with low ESNR. The scaling by 40 in the denominator is because a ESNR>40 dB is considered to provide highly reliable phase calculation. In general, the PHI value is low in healthy articular cartilage and it increases as the severity of OA progresses. Figure 1 shows representative phase-retardation A-lines from a healthy (group S_1) and a moderate OA (group S_3) sample. Figures 1(a) and 1(b) are the original phase retardation A-lines before averaging. Figures 1(c) and 1(d) are the corresponding averaged phase retardation A-lines, where it is shown that the overall birefringence is reduced in the OA sample compared to the heathy sample. In Fig. 1(a), phase wrapping is observed when the phase accumulation is near π. The segmentation boundaries obtained from the slope-based analysis are shown as the vertical lines. The averaged phase accumulation (averaged over the 3D volume) drops from 4.477 rad to 1.320 rad. Other than the reduction of birefringence, the increase of phase fluctuation can also be observed in the un-averaged A-line [Figs. 1(a) and 1(b)]. Only those peaks above the 20% threshold are detected and used for PAC and PHI calculation. Notice that there are no fluctuation peaks detected in the second region (R-2) of the healthy sample. In the OA sample, the number of peaks and the absolute phase change of the fluctuation peaks are larger, which results in an increase of the PAC values in the three zones. Thus, a large difference in the PHI value is achieved, where the PHI for the healthy and OA samples are 1.189 and 16.628, respectively.

Zonal distinguishability
As OA progresses, the ECM of articular cartilage is disrupted gradually, during which both the zonal division and the BRCs are affected. For healthy articular cartilage, the zone division shows clear and regular boundaries [21], and the corresponding BRCs usually have high distinguishability.
However, in the diseased cases, the zone division will be disrupted, and the BRCs of the different segmented regions could be hard to distinguish [25]. Those features can also be observed in Figs. 1(c) and 1(d). As OA progresses, the zone thickness may first increase due to swelling, then decrease due to ECM erosion, and finally the zones may disappear due to loss of articular cartilage layers. Similarly, the BRCs may decrease due to ECM disruption, but may also increase in some specific conditions [14]. Thus, the disease severity is not directly associated with the absolute values of the boundaries and slopes, but more related to the "clearness" of zonal division, or the zonal distinguishability.
To quantify the zonal distinguishability, two parameters are introduced based on the results from the slope-based analysis. The segmentation irregularity index (SII) is based on the boundaries b i (x, y), and the slope dis-similarity (SDS) is based on the slopes sl i (x, y). SII quantifies the fluctuation of the segmented boundaries of phase retardation image in PS-OCT. It is similar as OII that measures the irregularity of the boundary profile of intensity image in OCT [8,18]. The irregularity of the boundary profile is calculated as the difference between the original boundary and its smoothed regression. The irregularity is calculated for the boundaries b 0 (x, y), b 1 (x, y), and b 2 (x, y), respectively. Finally, SII is a weighted sum of the irregularity of the three boundaries: Here function g OI ( ) represents the OII analysis [8] adapted for the PS-OCT image. Since each boundary b i (x, y) is associated with two segmentation zones, their weighting factors w i and w i+1 are both included in Eq. (2). For the surface boundary b 0 (x, y), the weighting factor w 0 is associated with the air layer above the tissue. It can be assumed that w 0 =1 because the stochastic phase in air is stable and reliable after averaging. Here α is a normalization constant such that i=0,1,2 αw i w i+1 = 1. Another parameter, SDS, is developed to represent how clear the zone division is at the boundary. Two adjacent zones are considered to be mixed together when the BRCs from them are similar to each other. A two-sample Z test is used to estimate the global difference between the BRCs from two adjacent zones. Positive Z value means a decrease of the BRCs of the two adjacent zones along depth, and vice versa. Higher value of |Z| occurs if the BRCs of the two zones are very different and distinguishable. The SDS is calculated as: where Z 1 is obtained between sl 1 (x, y) and sl 2 (x, y), and Z 2 is obtained between sl 2 (x, y) and sl 3 (x, y).
Since Z 1 and Z 2 are associated with all the three segmentation zones, the weighting factors w 1 , w 2 , w 3 are compounded as w z = (w 2 w 3 )/(w 1 w 2 ) = w 3 /w 1 . Here (1-C X ) is a normalization factor based on the localized dis-similarity. The localized dis-similarity is calculated as C X = C 1 + w z * C 2 , where C 1 and C 2 are the zero-normalized cross-correlation (ZNCC) [26] between sl 1 (x, y) and sl 2 (x, y), and between sl 2 (x, y) and sl 3 (x, y), respectively. In this study, we use C X to characterize the sub-surface homogeneity from the zonal BRCs. Higher C X indicates the situation that zonal BRCs are more homogeneously distributed over the region of interest (ROI). Whether the normalization by (1-C X ) is necessary or not depends on the sign of Z 1 Z 2 . When Z 1 Z 2 >0, it means the BRCs over the three zones change monotonically. This is usually the case for healthy articular cartilage, where the BRCs change regularly and consistently. In this case, C X is close to 1 and normalization to (1-C X ) could introduce additional noise, and thus it is eliminated; When Z 1 Z 2 <0, it means the BRCs first increase and then decrease, or vice versa. This is usually the case in OA, where the BRCs change more irregularly. In this case, normalization by (1-C X ) is necessary. The zonal distinguishability D z is defined as D z = SDS/SII. Higher slope dis-similarity will contribute to higher zonal distinguishability, while irregular segmentation boundary (high SII) is closely associated with mixed zones (low D z ). For the examples shown in Fig. 1, D z value of 9.290 and 0.448 are obtained for the healthy and OA cases, respectively, which shows a large difference.

PS-OCT based OA grade
Compared to healthy articular cartilage, the heterogeneity of collagen organization in OA cases typically results in higher PHI value, which means higher phase irregularity. The D z value is usually lower, which means that the zone-distinguishability is lower, and the zones are mixed together due to ECM disruption. Therefore, a PS-OCT based OA grade G PS is defined by combining PHI and D z as: Due to the different magnitudes of PHI and 1/D z , a scaling factor F s is applied to balance those two terms. In this study, F s = 2 is selected by trial and error. Here, PHI (x, y) is a 2D matrix in the ROI used for evaluation and the operation · indicates averaging over the corresponding ROI. Meanwhile, SII and SDS are single values.

Statistical analysis
Statistical analysis was conducted through MATLAB. One-way ANOVA was used to identify if any relationships existed between the three groups (S_1, S_2, S_3) followed by Tukey's post-hoc test using the multi-compare function to quantify the significance of the relationship. P-values ≤ 0.05 were considered statistically significant. "Strong" and "very strong" significance were defined by p-values ≤ 0.01, and ≤ 0.001, respectively.

Phase retardation slope-based analysis
PS-OCT imaging and slope-based analysis are carried out on the human cartilage specimens. Figure 2 shows representative results from a healthy (S_1 group, left column), mild OA (S_2 group, middle column), and moderate OA (S_3 group, right column) samples, respectively. Figures 2(a)-2(c) show the photograph of the three specimens, respectively. Only the moderate OA sample appears a distinct difference, while the healthy and mild OA samples look similar. The corresponding OCT intensity images are shown in Figs. 2(d)-2(f). From the intensity images, it is observed that the tissue surface changes from smooth to uneven, and to swelling, as the disease progresses. It is noticed that the mild OA sample show minor changes near the surface, e.g. some superficial abrasion and slightly higher irregularity, which may imply early disruptions to the ECM structure. The disruptions of collagen organization can be revealed by PS-OCT phase retardation images, which are shown in Figs. 2(g)-2(i), respectively. The color mapping corresponds to phase change over 0-π. The band-like patterns indicate how the accumulated phase retardation changes, which is related to the level of tissue birefringence. In Fig. 2(g), the sample shows relatively uniform structure, even surface, and strong tissue birefringence, which represent features of healthy cartilage. In Figs. 2(h) and 2(i), the samples show uneven surface and reduced birefringence, which indicate OA cases. A representative A-line from each sample is shown in Figs. 2(j)-2(l), respectively. From the selected phase retardation A-lines, one can further confirm that the healthy cartilage (S_1 group) exhibits strong birefringence, while the diseased cartilage samples (S_2 and S_3 groups) show reduced birefringence. The features of the phase retardation A-line are sensitive to the collagen organizations. In healthy articular cartilage, the collagen fibers are aligned parallel to the tissue surface in the superficial zone, bent and aligned obliquely in the transitional zone, and aligned perpendicularly to the tissue surface in the deep zone. The strong birefringence indicates that the collagen fibers in the healthy sample are highly organized and aligned. However, the accumulation of phase retardation becomes weak with very flat slopes in Figs. 2(k) and 2(l), which indicate low tissue birefringence. This is likely because the OA progression has caused disruption of the collagen fiber organization which resulted in the reduction of birefringence. In our study, the OA specimens usually show reduced tissue birefringence than healthy ones. In a few cases, regions of high birefringence have also been observed in OA samples and those samples show large fluctuations of tissue birefringence at different local regions. Contradictory findings were also reported in literature by statistical studies with a large number of samples [14,18]. Many factors can affect the measurement of tissue birefringence, including topographical locations and illumination angle of light, etc. Therefore, further quantification of the distinguishability of the zones by our PS-OCT based OA grading is carried out. One important observation is that the phase retardation images show increased fluctuations along both axial and lateral directions as the disease progresses, which implies an increased level of PHI and SII. Moreover, the zones become mixed together, indicating a reduction of SDS. By combining those features, the PS-OCT based grade G PS is expected to increase as the disease progresses.

Validation with histology images
To validate the PS-OCT results, histology images are taken from the same samples and used for comparison. The locations of PS-OCT and histology images are matched by marking the imaging location with notches on the specimen. The PS-OCT and histology images are compared in Fig. 3 for one specimen from S_1 group (1st column), one from S_2 group (2 nd column), and two from S_3 group (3 rd and 4 th columns). Their averaged OARSI scores are 4.38 ± 1.17, 5.85 ± 1.09, 15.94 ± 1.23, and 17.31 ± 1.31, respectively. The higher OARSI score indicates a more severe case (e.g. 17.31 indicates very close to late-stage). The OCT intensity images are shown in Figs. 3(a)-3(d).  For the healthy sample, the intensity image [ Fig. 3(a)] shows a smooth surface and homogeneous signal below the surface. Relatively high birefringence is observed from the phase retardation [ Fig. 3(e)], indicated by the strong rainbow-like patterns. The segmentation boundaries are regular and distinct. It is commonly observed from healthy and relatively healthy samples that the regions of the three structural zones can be differentiated by the slope-based analysis as R-1 for superficial zone, R-2 for transitional zone, and R-3 for deep zone, respectively. Features indicating healthy articular cartilage can be noticed in the histology image [ Fig. 3(i)], where intact surface, normal architecture, and appropriate chondrocyte orientation can be observed. In the mild OA sample, some minor superficial abrasions can be observed from the intensity image [ Fig. 3(b)]. In its phase retardation image [ Fig. 3(f)], a slight reduction in birefringence accumulation is indicated by the missing of the second blue band; some localized bumps appear at the boundary between R-1 and R-2 (the second dashed segmenting curve). In its histology image [ Fig. 3(j)], cellular proliferation and chondrocytes clustering can be revealed, which is a sign of early-stage OA. From the sample with moderate OA (3 rd column), uneven surface and inhomogeneous sub-surface signal can be observed from the intensity image [ Fig. 3(c)], with reduction in birefringence being observed from the phase retardation image [ Fig. 3(g)].
The segmented regions become irregular and, to some extent, disrupted. The corresponding histology image shows that the surface becomes disrupted and uneven; however, the cartilage maintains its affinity for Safranin-O stain (red color) which is indicative of the persistence of proteoglycans in the tissue [ Fig. 3(k)]. In the more severe case of moderate OA (4 th column), the top layers of the articular cartilage are eroded and a thinner cartilage layer is observed in the intensity image, where cartilage-bone boundary is indicated by a yellow dashed line [ Fig. 3(d)]. Further reduction in tissue birefringence is observed in the phase retardation image [ Fig. 3(h)]. Correspondingly, the histological image reveals that the cartilage has relatively severe OA, where the cartilage layer becomes very thin and shows no affinity for Safranin-O, indicating low concentration of proteoglycans [ Fig. 3(l)]. Thus, the results from PS-OCT imaging match well with the observations from the histological imaging. PS-OCT imaging and analysis can be used for distinguishing OA severities.

Quantification features
Using the slope-based analysis, three segmentation regions R-1, R-2, and R-3, can be obtained based on their different slopes. For each region, the BRC and the region thickness are quantified. Table 2 shows the averaged BRC and thickness over all the samples in the groups S_1, S_2, and S_3, respectively. For each sample, the BRC and thickness are analyzed over a ROI of 150 B-scans by 300 A-lines (∼3 mm by 3 mm). In S_1 group, the BRCs are relatively high and have a descending tendency from 9.190 ± 4.874, to 6.439 ± 1.839, and to 4.093 ± 4.069, in R-1, R-2, and R-3, respectively. In S_2 group, the BRCs are lower than that in S_1 group, and the difference among the different regions are not significant. Similar properties of the BRCs in S_2 group are also seen in S_3 group. It is observed that the ratio of standard deviation to the mean is higher in S_2 and S_3 groups than that in S_1 group, which implies the destruction of the ECM structure. In literature, a single slope (w-BRC) could be obtained by fitting the phase retardation A-line with one linear regression [20]. This single slope approach is also applied on our samples with the results shown in Table 2. As we can see, the single slope is not a distinct feature among the three groups. Therefore, a single slope cannot represent the complex situation in collagen degeneration or the zonal difference among the different stages of OA. Meanwhile, the thickness of each zone decreases from S_1, to S_2, and to S_3 as OA progresses, which may indicate the erosion of the ECM structure. To achieve better quantification, the PS-OCT based OA grade G PS defined in Section 2.3 is also applied to the samples, over the same ROIs (∼3 mm by 3 mm). The results are presented in the statistical analysis below.
The boxplots obtained from the statistical analysis conducted on the 23 samples (7 in S_1, 8 in S_2, and 8 in S_3) are shown in Fig. 4. The p-values between every two groups are calculated and marked by stars. The situation with p value < 0.05 is considered as statistical significance (*), p value < 0.01 as strong significance (**), and p value <0.001 as very strong significance (***).  Fig. 4(a)], there exists statistical significance between S_1 and S_2/S_3, but there is no significant difference between S_2 and S_3. From the BRC of zone R-2 [ Fig. 4(b)], there is a strong significance between S_1 and S_2/S_3, but S_2 and S_3 show no difference. For the BRC of zone R-3 [ Fig. 4(c)], there is a strong statistical significance between S_2 and S_3, and a statistical significance between S_1 and S_3, but there is no difference between S_1 and S_3. In Fig. 4(d), it shows that there exists a A single slope of the whole depth (w-BRC) is obtained by fitting the accumulative phase retardation A-line in the region between 15% -85% of the phase profile without unwrapping [20]. Due to the different implementations, the w-BRC is not comparable with the zonal BRCs.
no difference between the groups when using a single slope (w-BRC) for the full depth. This is consistent with the statistical studies reported in Ref. [20]. Figure 4(e) shows the statistical result for our PS-OCT based OA grade G PS value, which shows distinguishing power among the three groups. The difference between S_1 and S_2/S_3 is very strong (p<0.001), and there exists strong significance (p<0.01) between S_2 and S_3 as well. Compared with BRCs for each zone or for the full depth, the G PS value shows much stronger capability in distinguishing OA stages. From the boxplots in Fig. 4, G PS is the only parameter that shows consistent significance between the sample groups of healthy versus mild, and mild versus moderate. Thus, with the proposed G PS value, PS-OCT offers reliable capability in detecting early-stage OA. Green dots indicate the original data points. Statistical significance is indicated by three levels, * for significance (p value < 0.05), ** strong significance (p value < 0.01), and *** for very strong significance (p value < 0.001).

Discussion and conclusions
Previous studies based on PS-OCT have used a single slope (w-BRC) or related features to quantify cartilage degeneration, which did not consider the zone-based ECM distribution and disruption [14,17]. Therefore, those studies still found challenges in detecting OA severities. Both the ECM structure and the disease progression vary over articular cartilage depth. Therefore, by differentiating the structural zones in articular cartilage and utilizing the information based on the zonal segmentation, enhanced sensitivity has been achieved in distinguishing OA severities. With this method, a significant relationship between the PS-OCT based grade G PS and the histology scoring of the cartilage health states is achieved. The two important parameters involved in this study, PHI and D z , are developed from the morphological features (OHI and OII) of OCT intensity images. After a proper development to fit the context of phase retardation accumulation and careful modification with the prior knowledge of articular cartilage zonal architecture, PHI and D z perform very well in detecting the tissue inhomogeneity and the zonal distinguishability. This demonstrates the advantage of using PS-OCT with a proper quantification analysis in early-stage OA detection. In this study, accumulative phase retardation is used to quantify and segment the different structural zones in articular cartilage. By analyzing the heterogeneity and distinguishability of features in cartilage tissue, our method shows promising results in detecting early stage OA. Nevertheless, with accumulative phase retardation, the effect from local birefringence and optic axis cannot be separated. Recently, several groups have developed approaches to separate and quantify local birefringence and optic axis in tissues. Ravanfar et al. developed an optical polarization tractography based on PS-OCT, which revealed the local birefringence and fiber orientation, improving the visualization and characterization of the zonal structure in human cartilage [27]. Li et al. used PS-OCT to detect the local optical axis orientation in sheep cornea, which revealed cornea lamellar organization into thin sections [28]. Detecting the local birefringence and optic axis can provide more details about the collagen organization. However, in order to recover the local birefringence and optic axis orientation, there are also more challenges due to the effect of preceding tissue layers, transmission through fiber and system elements, and imperfect system alignment [29]. Such effects need to be compensated or calibrated through correction factors or special hardware. Nevertheless, it will be a promising future direction.
The light penetration depth is limited to <1 mm into the cartilage in this study, which is normally less than the typical articular cartilage thickness (1∼1.5 mm). For healthy OA samples, the PS-OCT imaging reported in this study typically reaches into the deep zone. With the progression of OA, the cartilage wears out gradually and becomes thinner, which makes it possible for the illumination to penetrate into the subchondral bone below the cartilage [ Fig. 3(c)]. All the parameters developed in this study have penalized the impact of the signals from the deeper regions by either assigning a weighting factor based on ESNR or avoiding the use of the third inner boundary. Since the disease usually starts from the tissue surface, the penetration depth will not have much impact on OA severity detection.
One limitation of this study is the relatively small sample size (30 samples from six patients). In this study, the OA severity is considered to be localized. Five sampling locations are extracted from each joint and they are independently graded by histology and PS-OCT based approaches. Different disease severities could be obtained from the different sampling locations from the same joint. The authors have validated that the degeneration levels could vary a lot topologically. This assumption is reasonable for the sake of developing a technical solution for OA detection and evaluating, which greatly releases the heavy burden for a large sample size. Nevertheless, such assumption should be avoided in future clinical studies. Pathologically, despite the existence of complex local lesions, OA is more believed to be a disease of whole joint, where simultaneous changes occur in cartilage, bone, bone marrow, synovium, menisci, ligaments and even neural tissues [30]. In this study, frozen articular cartilage samples are used. The freezing/thawing process has been reported to affect the viability and integrity of chondrocytes [31]. However, no obvious difference is observed in the PS-OCT images when comparing frozen versus fresh cartilage samples. This is likely because the PS-OCT contrast comes from the collagen structure which tends to maintain its features. In the future, fresh samples with a larger sample size should be imaged to further validate the approach.
OA treatment is as important as its diagnosis. For decades, cartilage regeneration has drawn attention as an important treatment form [1,[32][33][34]. Procedures such as microfracture, implantation with autologous chondrocyte and embryonic stem cells, scaffold supports, and biofilm culturing, would help recover the disease yet the mechanism and pathological utility remain elusive [35]. PS-OCT offers great potential in monitoring many of those processes with its powerful capability in assessing mesoscopic structure and its functional bio-sensitivities, especially focusing on the tissue surface and a few hundred micrometers below. For achieving in-vivo diagnosis and real-time treatment monitoring, an important future direction is to develop an arthroscopic setup for PS-OCT imaging. With this setup, the PS-OCT based OA evaluating method will be further validated and improved in more practical manners.
Apart from ECM disruption and remodeling, mineralization is another important process in both OA progression and cartilage regeneration [1,34]. Changes of subchondral bone (usually 1 mm below surface) with increased metabolism and sclerosis (mineralization levels), are also detectable alterations at the very beginning of the OA process [34]. In healthy articular cartilage, the thickness of the calcified zone is maintained by a balance between progression of the tidemark into the unmineralized cartilage and changing into bone by vascular invasion and bony remodeling. Such balance will be broken in OA where progressive mineralization occurs, which results in the thinning of cartilage and the thickening of the calcified area [36]. Rebuilding this balance is very critical in OA treatment. Such changes could possibly be assessed by the proposed PS-OCT imaging and the parameters (e.g. PHI) defined in this study, although the penetration may still be a limitation for in-vivo detection. PS-OCT offers additional bio-sensitivities than OCT intensity imaging. With the parameters developed in this study, PS-OCT would have more capability to be applied in many related areas. Thus, PS-OCT imaging, together with proper quantification strategies, has the potential to be applied to investigating the pathological development, regenerative recovering, and tissue-engineering in OA diagnosis and treatment.