3D morphometric analysis of calcified cartilage properties using micro-computed tomography

Powered by TCPDF (www.tcpdf.org) This material is protected by copyright and other intellectual property rights, and duplication or sale of all or part of any of the repository collections is not permitted, except that material may be duplicated by you for your research use or educational purposes in electronic or print form. You must obtain permission for any other use. Electronic or print copies may not be offered, whether for sale or otherwise to anyone who is not an authorised user. Kauppinen, S.; Karhula, S. S.; Thevenot, J.; Ylitalo, T.; Rieppo, L.; Kestilä, I.; Haapea, M.; Hadjab, I.; Finnilä, M. A.; Quenneville, E.; Garon, M.; Gahunia, H. K.; Pritzker, K. P.H.; Buschmann, M. D.; Nieminen, H. J. 3D morphometric analysis of calcified cartilage properties using micro-computed tomography


Introduction
Osteoarthritis (OA) is typically considered as a progressive joint disease that produces several pathological changes, such as cartilage matrix erosion, tidemark duplication, bone sclerosis and osteophyte formation 1,2 . However, the pathogenesis of these features is still poorly understood. Many hypotheses have been introduced about the pathogenesis of OA, and some studies focus more on the changes in cartilage matrix and chondrocytes as the source of disease progression 3e6 , while some focus on the subchondral bone (SCB) as the main site of OA initiation 7e9 . However, there is evidence that in OA, the morphological changes in periarticular bone predate the changes in articular cartilage 10,11 . This promotes the importance of studying mineralized regions such as calcified cartilage (CC) and SCB. Furthermore, current understanding of OA has also led to the phenotype theory, where initiation site, pathogenesis and underlying other diseases differ, and, therefore, lead to different manifestations of OA. 12 Noncalcified cartilage (NCC) and SCB are commonly separated by a layer of CC. This layer is important as morphological changes in it can affect the crosstalk 13 , biomechanical force distribution 14 and solute transportation 15,16 between the NCC and SCB. The transition boundary from NCC to CC is called tidemark 17 , where collagen fibrils of NCC are anchored to the CC 18 . The tidemark is metabolically active structure, representing the calcification front where the NCC becomes the CC possibly associated with matrix vesicles released from the surfaces of deep chondrocytes 19,20 . This process is inhibited by proteoglycans 21 and matrix Gla protein (K-vitamin dependent protein that binds free calcium ions) 22 . Therefore, there is at least some form of metabolical homeostasis in normal tidemark, which is disrupted in OA.
During the progression of OA, the tidemark advances into the NCC, slowly shifting towards superficial cartilage, but typically not exceeding the deep zone cartilage 23 . This is known as tidemark duplication, i.e., more than one tidemark may be present. Simultaneously, the CC is replaced by bone via chondroclastic resorption 24 leading to endochondral ossification 25 . Furthermore, hypocellularity and the presence of hypertrophic chondrocytes are linked to OA progression indicating a diminished number of chondrocytes in OA cartilage 26 . It has also been hypothesized that OA hypertrophic chondrocytes do not inhibit calcification, but act in the same manner as hypertrophic chondrocytes in endochondral ossification leading to degradation of the surrounding extracellular matrix 3 . In previous literature, chondrocytes located at the tidemark can be seen as dimples in scanning electron microscopy (SEM) images 27 . As the chondrocyte size and number affect CC topography and metabolic homeostasis in the tidemark is disrupted in OA, the CC surface topography would be expected to be changed with disease progression.
The CC contains vascular and avascular channels 28 . The subchondral plate, consisting of CC and SCB 29 , has three types of channels: (1) open channels encapsulated by CC, (2) closed channels encapsulated by bone and (3) open channels, which reach and penetrate the tidemark 28 . The latter ones are of special interest since they are open channels for nutrition and fluid exchange. The avascular channels can contain soft tissues e.g., adipocytes and bone marrow 28 . Thus, distinguishing avascular from vascular channels with micro-computed tomography (mCT) would require a specific contrast agent to differentiate between vascular tissue from other soft tissues. In this study, we use mCT without contrast agent and therefore we refer to these channels as open subchondral channels (OSCC).
To date, most studies investigating the tidemark or channels through the calcified layer are done from 5 mm thick histological sections 13,30,31 . For instance, studies have been conducted to assess how the tidemark roughness (TMR) is altered as a function of OA 30,31 . Roughness of a 2D surface can be easily calculated as a ratio between the true length of the surface and the projected length of that surface, where the minimum value of one defines a perfectly smooth surface 32 . Vessels in CC are also often studied from histological sections 33e35 . Using histological sections for quantifying CC morphometrics is prone to poor or irregular staining and cutting artifacts, and consequently often requires manual segmentation. A complex 3D structure can be derived from stacked sections of a volume, but the workload of serial sectioning is high and the coregistration of the consecutive histological slices is limited by the sectioning artifacts.
The mCT provides volumetric imaging of tissue blocks, thus overcoming the limitations of conventional histology. High resolution mCT imaging and user-independent segmentation of the tidemark could potentially capture the smallest variations in the tidemark and also the calcified content accurately; importantly, this could be achieved in 3D. Therefore, a mCT approach with software algorithms has the potential to quantify the true TMR with minimal operator input and bias. However, 3D methods to quantify the OA related changes in CC topography and vessel perforations are limited. Therefore, this study aims to establish a mCT based method to quantify the CC topography and OSCC perforations through CC in 3D. Furthermore, we demonstrate that the methods can be used to investigate the relations between OA, tidemark surface roughness and OSCCs.

Materials and methods
Sample preparation and mCT imaging Samples were harvested from both knees of six patients (age 49e67) (Maisonneuve-Rosemont Hospital, Montreal, Canada) undergoing total knee replacement (TKR) surgery and two asymptomatic cadavers (age 26 and 49) donated by RTI Surgical tissue bank (FL, USA), resulting in 16 tibial plateaus. Osteochondral cores (n ¼ 15, Ø ¼ 4 mm) were drilled from the weight bearing area of lateral tibial plateaus. One sample was excluded as calcified tissue was missing from the lateral weight bearing area. Samples were fixed in buffered 4% formaldehyde and imaged in the fixation media after minimum of 3 days of fixation. Imaging was conducted with the desktop mCT (Skyscan 1272, Brüker microCT, Kontich, Belgium). The following settings were used for projection image acquisition: X-ray tube voltage ¼ 50 kV, X-ray tube current ¼ 200 mA, exposure time ¼ 3200 ms, filter ¼ 0.5 mm Al, frame averaging ¼ 3, amount of projection images ¼ 1200, isotropic voxel size ¼ 2.8 mm. Projections were reconstructed with Nrecon software (v1.6.9.8) with beam hardening and ring artifact corrections applied. Sample preparation procedure is shown in Fig. 1 (Fig. 2). Histological slices were imaged with an automatic slide scanner (Aperio AT2, Leica Biosystems, Wetzlar, Germany) using 40Â  magnification and a 0.25 mm pixel size. The OARSI grading was performed by three individual graders, (SK, SSK, LR) who were blinded to sample quality information during grading. In case of disagreement, a unanimous consensus grade, reached after discussion, was used as the final OARSI grade. Quantitative histological image analyses (Fig. 2) were performed on three slices per sample and average of the three measurements was used as the final parameter. All quantitative histological analyses were done in Fiji 36 using a custom-made macro. First, CC and SCB plate were segmented manually, and masks of the segmentations were used for calculating the cross-sectional areas and feret diameters (longest distance parallel to the joint surface, red line in Fig. 2) as described previously 30 . Local thicknesses of CC and SCB plate were calculated within the macro using BoneJ plug-in 37 as described previously 38 . The calculated parameters were: number of tidemarks (TM.number), CC cross-sectional area per feret diameter (CC.T mean ), SCB plate cross-sectional area per feret diameter (SCB.T mean ), local thickness of CC (CC.T local ), and local thickness of SCB plate (SCB.T local ) (see Table I for list of parameters).

mCT data analysis
From mCT datasets, volumes of interest (VOIs) were generated with Matlab (v. R2015a, Mathworks, Natick, MA, USA) by fitting a 600 Â 400 Â Z pixel VOI to the center of the tidemark surface. Z varied between 97 and 330 pixels (271.6e924 mm) depending on the subchondral plate thickness. The calcified content from the VOIs was binarized using fuzzy C-means clustering 39 in a separate custom-made Cþþ software. The raw data was divided into background and sample with a sample threshold probability of 0.64. Cmeans convergence was set to 0.0001 which means that iterations stop when the change in probability reaches this value. The continuity of the volume was checked and only the largest object was left as a result. All the subchondral channels were segmented with CTAn software (v1.14.4.1, Bruker microCT, Kontich Belgium) by using ROI-shrink wrap operator in 3D to fit a region of interest to the whole volume. The original volume was then removed from this volume, thus, leaving us with only the channels and some boundary artifacts. Most boundary artifacts were removed using common morphological operations such as closing and opening. Remaining artifacts were removed using despeckle techniques. Subsequently, the binary mask of all the channels within the calcified content was added to the binary mask of the calcified content resulting in a mask where no channels are present. The surface of this combined mask (defining the interface between NCC and CC) was segmented and only the channels in connection with this interface were considered as channels that perforate through bone and CC to the tidemark. As a result, only the channels connecting NCC to the epiphyseal compartment of the tibia were used in the analyses. TMR was calculated as an average roughness from 600 slices from the VOI by finding the shortest route through an automatically segmented surface using an A-star search algorithm (v1.0) 40 and fitting a simple linear regression to each route. The generated regression lines were used as reference lengths for each shortest route through the tidemark surface. If channels were found in any slice, a combined length of these channels was removed from the calculations in order to only analyze the CC surface. Visual explanation of TMR calculation is presented in Fig. 2(E). The equation used to calculate the TMR is: where n ¼ number of slices and i ¼ running slice number. Parameters for the channels were calculated from a projected channel mask and they were: channel area fraction (CAF), channel density (CD), and mean area of the channels (MC.area) [ Fig. 2

(F)].
A local binary patterns (LBP) -based 41 analysis was conducted for the CC surface to assess the distributions of the local surface orientations in 3D. The same VOIs were used for the LBP-based analysis as for TMR, and the channel mask was used to exclude the locations of the channels perforating to the surface. The LBPbased method used here was adapted from a previous study 42 . Briefly, a specific pattern is calculated for each voxel adjacent to the segmented surface, representing the local volumetric orientation of the surface. Each studied voxel is the center of a sphere with equidistant points fitted on it. A pattern is then calculated by assessing which points are within the surface. Only the patterns that have a clear orientation (uniform patterns) were selected for the analysis. Eventually, a pattern represents the local volumetric orientation of the surface within the neighborhood of the studied voxel. In our analyses, a radius of 2.8 mm was used for the neighborhood distance and 32 points were fitted on each sphere to define the amount of potential different patterns (here two 32 ). Amount of different patterns (ADP) gives count of all different local patterns in the sample surface. The higher the ADP, the more different local volumetric orientations are observed within the studied surface. The count was then normalized by the analysis area due to removal of channel locations. To quantify the consistence of each pattern along the sample, the entropy of the patterns distribution (EP) was calculated as follows: P i contains a count of a specific pattern i. EP describes the randomness of the local patterns, which represents how likely a specific pattern might occur within the sample 42 . The homogeneity index (HI) 42 of local angles was also calculated. This parameter describes the continuity of the angles calculated from the patterns by quantifying the overall local variations in orientations along the surface. An angles co-occurrence matrix ebased approach is used to assess how the orientations of the discrete areas of the surface differ from their local surroundings. The method is described in detail in Thevenot et al. (2017) 42 .

Statistical analyses
Linear associations between the parameters were evaluated using Pearson's correlation coefficient, r p , for continuous variables and Spearman's rho, r s , if either of the tested variables were ordinal. Correlation coefficients of 0.40e0.59 were considered as moderate, 0.60e0.79 as strong, and above 0.80 as very strong 43 . Samples were also divided into three groups based on the OARSI grading (healthy (OARSI ¼ 0, n ¼ 5), early OA (OARSI ¼ 1e2, n ¼ 7) and advanced OA (OARSI >2, n ¼ 3)) in order to see if changes in CC can be seen in the earliest stages of OA, where there are minimal changes in the superficial layer of cartilage. Differences between the OARSI groups were tested using a linear mixed model which takes correlation between the tibias of the same subject into account. OARSI grade was set as a fixed and the subject as a random variable. Restricted maximum likelihood estimation was used in the model. A P-value smaller than 0.05 was considered statistically significant. Statistical analyses were performed using SPSS (v. 22.0, IBM SPSS, Armonk, NY, USA).
Compared to healthy samples, early OA had lower values of EP (P ¼ 0.005), ADP (P ¼ 0.012) and TMR (P ¼ 0.037), and advanced OA had lower values in EP (P ¼ 0.014). No statistically significant differences were found between healthy and advanced OA in ADP (P ¼ 0.064) and TMR (P ¼ 0.142), or between healthy and early OA, or advanced OA in HI (P ¼ 0.102 and P ¼ 0.055, respectively). Early OA and advanced OA did not differ in any of the test parameters. Box plots visualizing the group comparisons are shown in Fig. 3.
All the measured parameters and their standard deviations (SDs) in the three OARSI groups are presented in Table IV. Feret diameters that were used for normalizing the average thickness measurements, had a coefficient of variation (CV) of 8.14% for CC measurements and 8.29% for SCB measurements in the whole data set. Healthy, early and advanced groups of OARSI had a CV of 4.67%, 8.49% and 9.90% for CC measurements and 4.79%, 9.51% and 9.73% for SCB measurements, respectively.

Discussion
In this study, we established a new method to quantitatively evaluate the roughness, the local surface variations of the tidemark and perforations through CC in 3D. Furthermore, we demonstrated that the established parameters from the tidemark and CC could be used to quantitatively study relationships between pathological features relevant to OA by comparing the parameters with the histological OARSI grade (a measure for OA progression).
While we had samples from two asymptomatic cadavers, only the 26-year-olds samples were classified as healthy (OARSI grades: Manually segmented SCB plate cross-sectional area divided by its feret diameter 30 .
Average values were calculated from three slices spaced at a minimum of 30 mm.

CC.T local Calcified cartilage local thickness from histological sections [mm]
Local thickness of manually segmented calcified cartilage was calculated using Bonej plug-in 38 .
Average values were calculated from three slices spaced at a minimum of 30 mm.

SCB.T local SCB plate local thickness from histological sections [mm]
Local thickness of manually segmented SCB plate calculated using Bonej 38 .
Average values were calculated from three slices spaced at a minimum of 30 mm.

TM.number Number of tidemarks from histological sections
Maximum amount of separate tidemarks. Average values were calculated from three slices spaced at a minimum of 30 mm.  Our results demonstrate association between TMR and OA severity. An OA-related decrease in CC surface roughness was statistically significant in all the parameters describing CC topography. LBP-based parameters (ADP and EP) correlated stronger than TMR with OARSI. Conventional TMR analysis is done slice by slice in 2D. Therefore, complex structures in the surface might not be included in the 2D method. While the LBP-parameters are non-standard measurements, we decided to include them in this study to have a local 3D structural orientation assessment along the sample surface. The real 3D approach includes complex structures, which could explain the differences in the strength of the correlation. Furthermore, the direct relation between HI and OARSI grade was observed. An increase in HI indicates that the studied surface is in overall smoother with less local variations along the surface discrete curvatures. This suggests that, in the current sample set, the surface irregularities in the tidemark diminish with OA progression. Grouping of the samples into three groups of OARSI also revealed that changes of CC topography especially in early OA can be captured with this method, providing tools for understanding the pathophysiology of OA.
The behavior of roughness of the tidemark during OA could be dominated by small dimples of normal chondrocyte lacunae in healthy cartilage. Especially in healthy samples, chondrocyte lacunae can be seen from our 3D mCT reconstructions [ Fig. 4(A)], similarly as in previous SEM studies 27 . The loss of normal chondrocyte function could lead to increased calcification and entrapment of chondrocytes inside CC 44 , thus leading to smoothing of the tidemark. To simplify, there could be two kind of roughness components in the tidemark: "high frequency" roughness, which originates from the chondrocyte lacunae and calcification inhibitory components; and "low frequency" roughness where the more overall form of the tidemark starts to appear as a wavy/undulating shape [ Fig. 4(A)] which could lead to an increase in low resolution roughness measurements as can be seen in the work by Deng et al. (2016) 30 . Combined findings from these studies would lead to a decrease-increase-decrease response of the TMR in OA. First, the high frequency roughness disappears (entrapment of the chondrocytes), then the low frequency undulations start to appear (uneven increase in calcification) and in the end (after the cartilage is eroded), the tidemark is completely smoothed by physical force of the bones grinding together.   30 , as roughness was seen to increase with OA progression (e.g., OARSI grade). The same roughness parameter was used in these previous studies, but they did not present the resolution of the acquired images. As the resolution determines the minimum size of the measurable feature in the surface, the results might not be fully comparable.
We found no statistically significant association between the OARSI grade and OSCC area fraction, density or average size. However, there was a strong correlation between the OSCC parameters and the surface parameters (Table II). For the first time, these results demonstrate a relationship between the tidemark surface irregularities and the perforations through the subchondral plate. This indicates an increase in TMR e also in local variations e when more OSCC perforations are present. This warrants for further investigation of this matter using a larger sample size.
OARSI grade showed an indirect relationship with CC.T mean . This relation can be explained by the location of the samples as there is a thicker CC in load bearing areas and healthy samples have a thicker CC. However, this would mean that the SCB replaces the CC as OA progresses and a direct relation between OARSI and SCB.T mean would be expected but this relation was not seen in either thickness measurements of the SCB. We saw a positive correlation between CAF and CC.T mean . As we are using the samples from the load bearing area of the tibial plateau, where CC and SCB are believed to be thicker 46 , these results are in line with findings of Lane et al. (1977) 47 . They found that remodeling of the SCB and CC was increased at the load bearing areas of the human femur. Furthermore, they detected a higher number of vessels in the load bearing sites regardless of age of the patient.
Number of tidemarks showed a negative correlation with both thickness measurements from SCB and a positive correlation with CC.T local . These results indicate that with a thicker SCB, fewer tidemarks are present and a thicker CC has a greater number of tidemarks than thinner CC. This is in line with Doube et al. (2007) 24 , where they found that CC thickness in 18 month old equine distal metacarpal condyle is dominated by chondroclastic resorption (bone replacing CC) rather than advancing of the  tidemark, which leads to multiple tidemarks 24 . This could potentially mean that SCB plate has replaced enough CC to hide "older" tidemarks in the deep CC as seen in Fig. 4(F), where in advanced OA the bone has replaced the lowest tidemark in the right side of the image. The number of tidemarks was also seen to decrease when the amount of separate OSCCs increase, but no correlation was found when compared to CAF and MC.area. The effect that the OSCCs may play is most likely related to their function (avascular or vascular).
Studying CC surface with mCT is superior to histology as it negates the need of manual segmentation and histological staining, which both can affect the reliability of the results. The segmentation of calcified content from soft tissue from mCT images is relatively easy as the X-ray attenuation between these tissues is so distinctive. Therefore, automatic user independent thresholding, such as fuzzy C-means clustering, is a repeatable and reliable way to segment calcified content from soft tissue. Unfortunately, the cement line (CC-SCB-interface) is indistinguishable with basic desktop mCT. However, with different imaging modalities, such as synchrotron X-ray and phase contrast Xray imaging, this method could be used for quantifying the cement line topography also 48 . We have no means of separating different types of channels from each other in this study. Limitations in contrast of soft tissue and resolution also make this method incapable of answering the question whether or not blood vessels are present within the detected channels. To differentiate the vascular from avascular channels with mCT, further studies with vessel and soft tissue specific contrast agents would be required. While the developed methods reveal quantitative associations between the pathological parameters related to CC, the results may or may not be representative for large population. However, OARSI grading focuses mainly on changes in the NCC while our method focuses on changes in CC topography. Decrease in CC surface roughness co-occurs with increasing OARSI grade. This association highlights that changes in CC topography are relevant in the progression of OA. Most importantly, the results demonstrate the feasibility of the developed methods for investigating these associations. In conclusion, we have established a new 3D method for quantifying tidemark topography and perforations through CC. The established methods provide means to study the structure of CC quantitatively from mCT images.
Author contributions SK, HJN, and SS contributed to the conception and design of the study. SK, HJN, SSK, IK, IH, EQ, MB, MG, and SS. participated in acquisition of the data. SK, SSK, TY, JT, LR, MAF, MH, HJN, and SS participated in analysis of the data. All authors contributed to interpreting the data, drafting or revising the manuscript, and have approved the submitted version of the manuscript.
Competing interests HJN, TY, and SS are inventors in a patent application related to image analysis of AC. Other authors report no conflicts of interest.
Role of the funding sources Funding sources are not associated with the scientific contents of the study.