In vitro method for 3D morphometry of human articular cartilage chondrons based on micro-computed tomography

y Research Unit of Medical Imaging, Physics and Technology, Faculty of Medicine, University of Oulu, Oulu, Finland z Infotech Oulu, University of Oulu, Finland x Department of Applied Physics, University of Eastern Finland, Kuopio, Finland k Medical Research Center, University of Oulu, Oulu, Finland ¶ Institute of Biomedical Engineering, Ecole Polytechnique de Montreal, P.O. Box 6079, Station Centre-Ville, Montreal, Quebec H3C 3A7, Canada # Biomomentum Inc., 970 Michelin St., Suite 200, Laval, Quebec H7L 5C1, Canada yy Department of Diagnostic Radiology, Oulu University Hospital, Oulu, Finland zz Department of Laboratory Medicine and Pathobiology, University of Toronto, Toronto, Ontario, Canada xx Mount Sinai Hospital, Toronto, Ontario, Canada kk Groupe de Recherche en Sciences et Technologies Biom edicales, Polytechnique Montreal, P.O. Box 6079, Station Centre-Ville, Montreal, Quebec H3C 3A7, Canada ¶¶ Department of Physics, University of Helsinki, Helsinki, Finland ## Department of Neuroscience and Biomedical Engineering, Aalto University, Espoo, Finland


Introduction
With osteoarthritis (OA), chondrocytes are known to become hypertrophic and to form clusters, indicating that they may play a critical role in the increased metabolism observed in this condition 1,2 . The mechanisms that initiate these changes, however, are not fully understood. One explanation is that the chondrocytes produce more extracellular matrix (ECM) macromolecules trying to compensate the increased cartilage degeneration, which then leads to increased cellular clustering and hypertrophy 1 . Previous studies have also shown that not only the anabolic factors, but also the catabolic factors that contribute to ECM degradation are expressed in OA chondrocytes 3,4 . At the level of the chondron e the chondrocyte and its complex pericellular microenvironment 5e7 e previous studies have suggested that the cleavage of fibrillary collagens initiates chondron enlargement, which continues due to pericellular matrix deposition 8,9 . Both the chondrocyte and its pericellular matrix (PCM) thus appear to be modulated by OA 8,10 . While combining chondrocyte/chondron morphology with pathological pathways may yield new cues for understanding OA and cartilage metabolism, a lack of high-throughput tools currently exist for 3D chondron morphometry.
Some of the currently available methods for the quantification of 3D chondrocyte/chondron morphology include serial sectioning 11,12 , confocal 13 and multiphoton microscopy 14 , and synchrotron-based micro-computed tomography (mCT) 15 . Serial sectioning is a method in which multiple consecutive sections are cut from one sample and imaged; a 3D image is then reconstructed from the collected images 16 . The main problem with this method e besides the destruction of the sample e is that the cutting could damage the articular cartilage (AC) structures, including those of chondrons. While optical methods such as confocal and multiphoton microscopy 17,18 , can provide 3D images without destroying the specimen, their drawbacks are the limited light penetration in opaque tissues and the spherical aberration that results from refractive index mismatch 19 . It has been shown that synchrotronbased mCT can provide 3D information from chondrocytes 15 , although limited access to synchrotron facilities restricts the use of this approach. Other modalities, such as focused ion beam scanning electron microscopy and transmission electron microscopy, are also available for tissue 3D imaging within limited sample sizes 15 . Different mCT techniques are capable of providing volumetric information from the bone microarchitecture and cartilage structure of osteochondral samples 20,21 . Conventional desktop mCT imaging of AC requires contrast agents to provide distinguishable contrast between the structures 22e25 . However, current contrast enhanced mCT (CEmCT) methods are based on the contrast agent distribution in the tissue, which is affected by the electrostatic repulsion or attraction between the tissue constituents and the contrast agent molecule. Thus, the CEmCT can only provide indirect information about the composition of the AC, which may limit its use as a quantitative tool for morphometric chondrocyte analysis.
When imaging dried samples, the contrast arises from the tissue itself, not from external contrast agent distribution. And because native chondrocytes contain mostly water, which is removed during the drying process, the contrast between the chondrons and the ECM is enhanced. Hexamethyldisilazane (HMDS)-based air-drying was first proposed in 1983 as an alternative method for criticalpoint drying 26 . The surface detail in insect tissues and cells has been shown to be well maintained in dried HMDS-treated specimens 26e29 . Because HMDS has lower surface tension than water and is able to cross-link proteins, the sample will likely not fracture and collapse during the drying process 26,29 . HMDS drying has never been used in mCT cartilage imaging, however.
In this study, we present a novel HMDS-based sample processing protocol to enable in vitro chondron imaging from human AC samples using desktop mCT. The objectives were to develop a new sample processing protocol for mCT imaging and a semi-automatic algorithm to select and segment chondrons in 3D. The developed methodology was applied to human osteochondral samples in order to compare the chondron morphology within intact and osteoarthritic AC at different tissue depths.

Sample preparation
Tibial plateaus of two cadaveric human donors asymptomatic of cartilage-related diseases (ages 26 and 49, body mass indices (BMIs) 18.4 and 30.6, one male and one female; RTI Surgical tissue bank, FL, USA) and four patients who had undergone total knee replacement surgery (ages 51e67, mean BMI 29.8, two males and two females; Maisonneuve-Rosemont Hospital, Montreal, Canada) were used in this study. The study was conducted under institutionally approved ethic committee certificates (C ER-13/14-30 for Polytechnique Montreal and C ER 14060 for Maisonneuve-Rosemont Hospital). The tibiae were preserved at À80 C, thawed once, transported at À20 to 0 C, and preserved again at À80 C before thawing for sample preparation. Osteochondral cores (diameter 4 mm) were prepared from both medial and lateral sides of the tibial plateaus and cut in half. One half was subjected to histological analysis and the other to HDMS-based mCT imaging. A pipeline figure visualizing the workflow of this study is shown in Fig. 1.

Histological analyses
Osteoarthritis Research Society International (OARSI) histological grading 30 was used to evaluate the OA progression of the samples (first half of the osteochondral core). The grading was performed for three consecutive 3-mm-thick Safranin Oestained histological sections. Sections were imaged with a virtual light microscope (Aperio AT2, Leica Biosystems, Wetzlar, Germany) using 40Â magnification and 0.25 mm pixel size. The OARSI grading was first performed independently by three graders [SSK, LR, IK; inter-observer reliability: intraclass correlation coefficient (ICC) ¼ 0.93, 95% confidence interval (CI) ¼ (0.84; 0.98)], and their consensus grade was given as a final grade for each sample. Finally, twelve samples were selected into two groups: intact cartilage (n ¼ 6, OARSI grades 0 and 1) and degenerated cartilage (n ¼ 6, OARSI grades 3 and 3.5).

HMDS-based mCT imaging
After preparation, the mCT samples (second half of the osteochondral block) were fixed in 4% saline-buffered formaldehyde for at least 5 days. Fixed samples were then dehydrated in ascending ethanol concentrations (30%e50%e70%e80%e90%e96%e100%) for a minimum of 3 h in each step, treated with HMDS for 3 h, and finally air-dried in a fume hood at room temperature overnight

Chondron selection and segmentation algorithm
Volumes of interest (VOIs) with sizes of 300 Â 300 Â Z (480 Â 480 Â Z h mm 3 , Z h being the height of the sample) were chosen for analysis, as middle of the image stacks as possible, but avoiding the imaging artifacts. A custom-made algorithm developed in Matlab (v.8.5, Mathworks, Natick, MA, USA) was applied to automatically select and segment the chondrons (see Supplementary material). Briefly, chondron selection was performed by assessing the amount of the volumetrically connected (within the 3D vicinity) voxels as possible chondron. A sub-VOI was then generated for each potential chondron to be volumetrically segmented. For each sub-VOI, 3D histogram equalization was applied to enhance the contrast between the chondron and the ECM. Then, a multiscale 3D local binary patternsebased method, adapted from 31,32 , was used to segment the chondron itself. This method assesses the volumetric continuity of the voxels that are considered part of the identified chondron by fitting multiple spheres with different radii and evaluating the gray-level distribution calculated at their surface. As a preliminary criterion, a minimum threshold of segmented volume (400 mm 3 ) was used to remove potential segmented artifacts.

Manual verification of the segmentation
A second custom-made Matlab algorithm was used to manually verify the automatic segmentations. It contains a graphical user interface that visualizes the superposition of the original image and the segmented mask in a 3D orthogonal view from the center of the segmented volume (see Fig. S4 in the Supplementary material) allowing the user to determine the accuracy of the segmentations. The algorithm also enables the manual annotation and separation between segmented chondrons containing either a single cell or multiple cells (cluster).

Algorithm validation
The intra-repeatability of the manual verification of the chondron selection was evaluated from 1000 segmentations from one sample (at two time points). To validate the performance of the automatic 3D segmentation script, 20 chondrons containing a single cell and 20 chondrons containing a cluster were randomly selected from different AC depths (range: 2e93% from the AC surface) from the manually verified automatic segmentations (22 chondrons from the healthy group, and 18 chondrons from the OA group). They were then segmented manually with MIMICS software (v.17.0.0.435, Materialise NV, Leuven, Belgium) by two independent users; the similarity of the manual and automatic segmentations was then evaluated by calculating a Dice similarity coefficient (DSC) for each chondron: where A is the sum of the common segmented voxels between the two compared segmentations, B is the sum of the voxels from the first segmentation, and C the sum of the voxels from the second segmentation. Finally, the average DSC [standard deviation (SD)] was calculated between the compared segmentations. Furthermore, the linear relationships between the automatic and manual segmentations were calculated (Fig. S5). The repeatability of the developed methodology and the manual segmentations are reported in the Supplementary material.

Volumetric parameters
Chondron density was calculated from the number of chondrons selected by the script. For each sample, script-identified chondrons The pipeline of the methodology. First, the osteochondral cores were prepared from the tibiae and cut in half. The first half was subjected to histology (OARSI grading), and the second half to mCT. Based on the OARSI grades, 12 samples were selected for this study: six samples with OARSI grades 0e1.0, and six samples with OARSI grades 3e3.5. After the HMDS-based sample processing protocol, the samples were imaged with the mCT device, followed by image reconstructions and VOI selection. The chondron selection and segmentation algorithm was then applied to the VOIs, and finally, chondron analysis were performed.
were checked to discover any false detection or eventual duplicates. Of all the segmentations, 500 chondrons/sample were randomly picked throughout the AC thickness for further investigation. After manual approval of the segmentation accuracy, the volume and sphericity 33 of the accurately segmented chondrons were calculated as well as their depth (%) from the AC surface. These parameters were calculated using CTAn software's (v.1.15.4.0, Bruker microCT) individual object analysis. The sphericity, first introduced by Wadell 33 , was calculated as follows: where V is the object volume and S its surface area. A sphericity value of 1 refers to a sphere, and values smaller than 1 refer to nonspherical and complex 3D objects. The depth from the AC surface was calculated as the distance between the AC surface and the centroid z c coordinate of each chondron obtained with CTAn. The parameters were then divided into three AC depth zones (zone 1: 0e10%; zone 2: 10e40%; zone 3: 40e100% from the AC surface, similarly to previous literature 34,35 ) and grouped by the OARSI histological grades (OARSI grades 0e1.0, n ¼ 6; OARSI grades 3.0e3.5, n ¼ 6).

Statistical analyses
Linear mixed models (SPSS, v.22.0, IBM SPSS, Armonk, NY, USA) were used to compare chondron density, volume, and sphericity between the different OARSI grade groups, separately in the three different zones. Chondron density, volume, and sphericity were treated separately as dependent variables, patient number was set as a subject, and OARSI grade group as a fixed variable. Furthermore, in volume and sphericity analyses, different models were used for all well-segmented chondrons, chondrons containing single cells, and chondrons containing clusters. A P-value smaller than 0.05 was considered statistically significant. The data in the figures are presented as means, with 95% CIs. The raw data {means (SD) and medians [interquartile range (IQR)]}, together with sample-wise visualizations of the morphological results, are presented in the Supplementary material.

HMDS-based mCT imaging
Volumetric visualizations of HMDS-dried osteochondral samples from OARSI grades 0e4 imaged with desktop mCT are shown in Fig. 2. When compared to conventional Safranin Oestained histological sections (Fig. 3), the features of OA in different OARSI grades are similarly visualized in both the 2D mCT slices and the histological sections. The cartilage thickness, however, appears slightly smaller in most of the mCT images than in the histological sections, and some inconsistencies in image quality occurred, primarily in the AC surface and at the cartilage-calcified cartilage interface.

Automatic selection and segmentation performance
After the manual validation of the selected objects, a total of 12,804 chondrons were automatically segmented from all 12 samples. Out of 1000 segmented objects from one sample, only five selections were excluded during the first manual verification round, whereas during the second round, those same five were excluded along with three other selections, thus suggesting that the algorithm was able to accurately identify chondrons. Furthermore, 25% of the 500 chondrons/sample e the segmentation accuracy of which was determined using the manual verification algorithm e were approved for further analyses. The user input time for the verification of the 6000 chondrons (500 chondrons/sample) when using this new algorithm was roughly 48 h (on average 30 s per inspection), thus significantly decreasing the estimated time of 600 hours that a fully manual segmentation would require. The average DSC (SD) between the automatic and the two different manual segmentations were 0.85 (0.07) and 0.80 (0.07), respectively, and between the two manual segmentations, 0.78 (0.11).

Chondron morphology vs OARSI grades
The differences in chondron density between the two OARSI grade groups (Table I, Fig. 4) were not statistically significant in any depth zone. In contrast, the chondron volumes were significantly larger (zone 1: P < 0.05; zone 2: P < 0.001; zone 3: P < 0.001) in the OARSI grade 3.0e3.5 group compared to the OARSI grade 0e1.0 group [ Fig. 5(A)]. Similar results were observed when considering only the chondrons containing clusters, but the differences were significant only in zones 2 and 3 (P < 0.001) for chondrons containing only single cells [ Table II, Fig. 5(B) and (C)]. The chondron sphericity in zones 2 and 3 was significantly smaller (P < 0.001) in the OARSI grade 3.0e3.5 group compared to the OARSI grade 0e1.0

Discussion
This study has presented a novel HMDS-based sample processing protocol for osteochondral samples, which permits highresolution imaging of chondrons at full AC thickness using conventional desktop mCT. We also developed a semi-automatic 3D selection and segmentation algorithm and applied this protocol to quantify the morphology of intact and osteoarthritic human AC chondrons in 3D. The proposed method is time-efficient and provides a high-throughput approach to investigate changes that occur in AC chondrons in OA with minimum user input.
Unlike previous methods for 3D chondrocyte/chondron imaging, such as confocal microscopy 1,13 and synchrotron-based mCT 15 , our protocol is relatively fast and easy to apply, and it does not require any extra equipment or specific stains. Our protocol also allows for imaging considerably larger volumes than e.g., with confocal microscopy. HMDS was selected as a drying solution because it allows the complete drying of the cartilage, thus enabling contrast in the mCT images to arise from the natural X-ray attenuation of the tissue itself, rather than from external contrast agent distribution.
The second objective of this study was to develop and validate an automated 3D selection and segmentation algorithm for chondrons imaged with mCT. The main advantage of this algorithm is that it allows rapid 3D analyses of multiple chondrons and, unlike threshold-based algorithms, it takes into account the volumetric information of the chondron. This approach is also easy to implement in other laboratories and for other segmentation purposes. When compared to manual segmentation, this new algorithm is roughly 12 times faster, the only user task being to check the performance of the segmentation provided by the script. From the 500 chondrons/sample selected for manual inspection, 25% were approved and further analyzed. This seemingly low percentage was primarily the result of the cartilage surface area, where the segmentation was fairly difficult e even manually e due to inconsistencies in image contrast and the chaotic-like cartilage surface morphology (especially in the case of OA samples). The area near the interface between the non-calcified and calcified cartilage is also challenging to segment because of beam-hardening which results in streaking artifacts and incorrect attenuation values near the interface. In this protocol, these artifacts occurred because of large differences in the X-ray attenuation between the noncalcified and calcified tissue, and could be prevented by decalcification of the sample. While the percentage of approved chondrons easily could have been increased by adding conditions related to locations in the algorithm, we decided not to adjust the script to avoid missing any potential chondrons in these challenging areas. The second reason for some of the inaccurate segmentations was due to brighter areas (the remnants of chondrocytes) inside or at the border of chondrons, which affected the automatic process. The accuracy of the automatic segmentations was visually evaluated from three orthogonal 2D views rather than from the full 3D volume. While the latter method would be more accurate, it would also be more time-consuming because of higher computational costs related to the generation and rotation of 3D models.
The average DSC calculations were used to evaluate the performance of the automatic 3D segmentation algorithm, as in a previous publication 36 . The average DSC (SD) between the automatic and manual segmentations for the first and second segmenters were 0.85 (0.07) and 0.80 (0.07), respectively, which suggests that both segmenters agreed well with the automatic segmentation; we can therefore conclude that the script for automatic segmentation was accurate. The average DSC between the two segmenters, however, was slightly lower [0.78 (0.11)], which may have indicated that the manual segmentation was not completely inter-reproducible.
The developed segmentation script was able to identify a total of 12,804 chondrons from all 12 samples. Even though the segmentation of chondron borders occasionally failed, the script was able to automatically identify a large number of chondrons. The chondron density in zones 1 and 2 was observed to be slightly lower in the OARSI grade 3e3.5 group compared to the OARSI grade 0e1 group, but not at a statistically significant level. In zone 3, chondron density was similar between the two OARSI grade groups. The absolute chondron density results in the middle and deep AC were in the same range as previously reported 11,37,38 , although our density values in zone 1 were considerably lower than what previous studies have found. As mentioned above, zone 1 was problematic for the segmentation script, which could explain the high variation and seemingly low density values in that zone. Previous literature shows that in healthy AC, chondrocyte density seems to decrease at greater tissue depths 11,37,38 , which we also observed in our results when not considering zone 1.
From the morphological analysis, we found that the chondrons were significantly larger in the OARSI grade 3e3.5 group. These results concur with previous studies that have reported increasing chondron volumes with OA 8,9 . The observed hypertrophic morphological changes in chondrons could be explained by their upregulated metabolism induced by OA activity 7,39 . Horikawa et al.
(2004) suggested another explanation for chondron hypertrophy; the authors showed that the volume ratio of the pericellular microenvironment to chondrocyte increases with OA 40 . We observed in our study that the chondrons were less spherical in the OARSI grade 3e3.5 group in zones 2 and 3, which suggests that the increase in chondron size did not occur evenly in all directions, but the hypertrophic changes made the chondrons more cylindrical. Alexopoulos et al. (2005) have previously shown that the Young's modulus of PCM is significantly lower than that of the surrounding ECM 41 . This finding indicates that the PCM shape might be modulated by the properties of the surrounding tissue and could partly explain the morphological changes in the chondrons observed in our study.
In zone 1, statistically significant differences in chondron volume and sphericity were not systematically observed, which indicates chondron hypertrophy to be less common in the superficial AC layer, although this observation could also be explained by either a deficient segmentation of smaller chondrons or as a sample processing artifact in the superficial AC layer. The superficial layer remains a challenge for chondron segmentation, mainly due to the shrinkage of the cartilage surface during the HMDS drying protocol. Poole et al. (1987) have shown that surface chondrons do not have the pericellular capsule that surrounds the pericellular matrix of the chondrons in the middle and deep AC 6 . This could explain why surface chondrons are more prone to shrinkage during the drying. Another previous study has also shown that in OA, type VI collagen expression and distribution is increased in the lower-middle and upperdeep zones but not in the upper zones 42 . This finding suggests that the protective barrier of the chondrons in the middle and deep AC is stronger than in surface chondrons, thus making the surface chondrons more vulnerable to changes caused by the sample processing. Previously reported chondron/chondrocyte/PCM volumes have shown relatively large variation 1,11,13,15,37,38,40,43 , most likely due to the different species, sample processing, and imaging modalities used. In a comparison of our volumetric results of chondrons containing single cells to those of previous human studies 1,11,37,38,40 , our results were more or less in the same range, although our observed values were slightly larger. The lack of studies that have quantitatively investigated the morphology of clusters prevents a thorough comparison of our observations with the previous studies. Finally, it is important to note that chondrons are dynamic 6 (i.e., they respond to changes in their environment), and therefore their morphological properties may differ depending on the OA stage; a comparison of the results obtained from different modalities thus may not be directly relevant.
One major limitation of this study was the unfortunate freeze/ thaw cycles before the sample processing. This could not be prevented because of the distant origin of the samples (Montreal, Canada) and because they had been used in other studies as well. All the samples did go through the same processes, however, and therefore the results obtained in this study should be comparable to each other. Although HMDS has relatively low surface tension and the ability to cross-link proteins 26 , some shrinking was observed in the cartilage because of the drying, with zone 1 being affected the most. The sample processing protocol could thus be developed further to minimize (and potentially even prevent) the cartilage shrinking. With the segmentation verification system used in this study, the differentiation between the chondrons containing a single cell and those with a cluster was not completely accurate, as some clusters may not have been visible in the three orthogonal 2D views, and thus they could have been labeled as chondrons containing a single cell. To avoid any misinterpretation in the verification process arising from the current system, verification of the automatic segmentations could be conducted from the full 3D volume instead of from the orthogonal views. Another limitation was that our sample size was relatively small, which was partly due to the novel method we were testing. Thus, a larger number of samples (and more thorough set of OARSI grades) should be included in the future studies. Because the selection of cores was based on their histological grades, some variation in core locations occurred. Separate datasets were not used in the validation of the methodology and the chondron morphometry assessment, which could also be considered a minor weakness. Finally, the zonedivision approach we used is accurate only for samples with remaining cartilage surface (OARSI grades 0e3.5). Other zonedivision methods (e.g., those based on polarized light microscopy) would be necessary if investigating more advanced OA grades.
To conclude, we have developed a novel sample processing protocol for imaging AC chondrons in 3D with a conventional desktop mCT, and a validated semi-automatic algorithm for  chondron selection and segmentation for morphological 3D chondron analysis. The protocol and algorithm were applied to quantify the morphology of human AC chondrons in 3D and to assess the differences between intact and OA samples. We did not find a statistical difference in chondron density between the intact and OA samples but chondrons did become larger and more elongated in OA, which indicates that the catabolic events in the surrounding ECM had an effect on chondron morphology. The HMDSbased sample processing protocol also has the potential for more than chondron imaging: it has been suggested that this method could be useful for 3D quantification of the collagen distribution in AC as well 44 . This new mCT protocol provides a new approach to spatially evaluate not only chondron/chondrocyte properties but also ECM macromolecule distribution in 3D during different phases of OA.