Differences in tibial subchondral bone structure evaluated using plain radiographs between knees with and without cartilage damage or bone marrow lesions - the Oulu Knee Osteoarthritis study

Objectives To investigate whether subchondral bone structure from plain radiographs is different between subjects with and without articular cartilage damage or bone marrow lesions (BMLs). Methods Radiography-based bone structure was assessed from 80 subjects with different stages of knee osteoarthritis using entropy of Laplacian-based image (ELap) and local binary patterns (ELBP), homogeneity index of local angles (HIAngles,mean), and horizontal (FDHor) and vertical fractal dimensions (FDVer). Medial tibial articular cartilage damage and BMLs were scored using the magnetic resonance imaging osteoarthritis knee score. Level of statistical significance was set to p < 0.05. Results Subjects with medial tibial cartilage damage had significantly higher FDVer and ELBP as well as lower ELap and HIAngles,mean in the medial tibial subchondral bone region than subjects without damage. FDHor, FDVer, and ELBP were significantly higher, whereas ELap and HIAngles,mean were lower in the medial trabecular bone region. Subjects with medial tibial BMLs had significantly higher FDVer and ELBP as well as lower ELap and HIAngles,mean in medial tibial subchondral bone. FDHor, FDVer, and ELBP were higher, whereas ELap and HIAngles,mean were lower in medial trabecular bone. Conclusions Our results support the use of bone structural analysis from radiographs when examining subjects with osteoarthritis or at risk of having it. Key points • Knee osteoarthritis causes changes in articular cartilage and subchondral bone • Magnetic resonance imaging is a comprehensive imaging modality for knee osteoarthritis • Radiography-based bone structure analysis can provide additional information of osteoarthritic subjects


Introduction
Bony changes, including osteophytes or subchondral cyst formation, are clearly seen on plain radiographs and are providing useful morphologic information in diseases affecting bone density and structure, such as osteoarthritis (OA) or osteoporosis. Although, the plain two-dimensional (2-D) radiograph is a projection (summation) through the actual threedimensional (3-D) structure, bone density and bone structure as depicted by plain radiographs is significantly related with the actual 3-D structure of bone [1][2][3][4][5].
Diagnosis of OA is based on a subject's history and symptoms, physical findings, and characteristic changes on plain radiographs. Typically, the severity of OA is evaluated from radiographs using the Kellgren-Lawrence (KL) grading scale, which is based on the visual evaluation of joint space narrowing, subchondral bone sclerosis, presence of osteophytes, and deformation of bone ends [6]. As ordinal grading using the KL scale gives only a summary score of overall disease severity with varying intra-and inter-rater reliability [7][8][9], development of quantitative and userindependent image analysis algorithms that exploit additional radiographic information is important to potentially enhance the clinical value of plain radiographs in OA diagnostics.
Joint space width is the most common parameter measured quantitatively from plain knee radiographs [10,11]. Other parameters related to subchondral bone structure have also potential to be used as an additional measure in OA diagnostics and characterization with potential relevance for prediction of disease progression. Fractal analysis is the most popular method to assess bone structure from radiographs in OA research and a method called fractal signature analysis (FSA) has been shown to predict disease progression [12,13]. Furthermore, it has been reported that bone structure assessed from plain radiographs using Laplacian-based method, local binary pattern (LBP)-based methods, and FSA is significantly related with the 3-D microstructure of bone [5]. Recently, subchondral and trabecular bone structures evaluated using LBP-based and Laplacian-based methods have shown to differ between subjects with different KL grades [14]. In that study, the KL grading and structure analysis of bone was made for the same images making the measurements dependent on each other to some extent, since features evaluated in the KL grading, for instance, bone sclerosis, affect the structural parameters as well. In order to study further the potential relevance of the radiography-based bone structural analysis methods, these should be compared with independent reference methods.
Magnetic resonance imaging (MRI) is considered the most comprehensive imaging modality for assessment of knee OA in a research context [15]. Semi-quantitative scoring systems that evaluate features related with or altered in the knee OA process have been developed and used for the assessment of structural deterioration of tissues within the knee joint [16]. Among many different features, MRI enables direct evaluation of cartilage damage and subchondral bone marrow lesions (BMLs) that are known to be related with OA incidence and progression [17][18][19][20][21]. However, the differences in subchondral bone structure from radiographs among subjects with and without morphological changes of articular cartilage or BMLs has not been thoroughly clarified yet. Given the fact that the local biomechanics are altered as a result of structural damage, radiographic bone structural changes are expected with prevalent cartilage damage and BMLs or both. Therefore, the aim of our study was to investigate whether subchondral bone structure from plain radiographs is different between subjects with and without articular cartilage damage and between subjects with and without BMLs focusing on the medial tibia.

Study subjects
This cross-sectional (level 3) study is part of the Oulu Knee OA (OKOA) study [22] and includes 80 subjects (49 women, 31 men) with different stages of symptomatic knee OA (KL grades from 0 to 4; Table 1). Written informed consent was obtained from each participant. The study was carried out in accordance with the Declaration of Helsinki and approved by the Ethical Committee of Northern Ostrobothnia Hospital District, Oulu University Hospital (number 108/2010).

Acquisition of radiographs
Weight-bearing posterior-anterior fixed-flexion radiographs of both knees with the X-ray beam angle at 10 degrees were obtained from all subjects. Since the original pixel size varied between images (mean: 138 μm, standard deviation: 20 μm, range: 100-168 μm), all images were resampled to the pixel size of 143 × 143 μm 2 using bicubic interpolation to ensure comparability of the structural parameters, without producing as much artifacts as bilinear or nearest neighbor interpolation algorithms.

Selection of regions of interests
Two rectangle-shaped regions of interests (ROIs) were defined (Fig. 1). The locations of the ROIs were based on previous literature [5,14,23,24]. One ROI (size: 14 mm × 6 mm) was placed into the subchondral trabecular bone in the middle of the medial tibial plateau immediately below the cartilagebone interface. This ROI is referred to as the subchondral bone ROI in this study, although different bone types are mixed in the ROI [25,26]. Furthermore, another ROI (14 mm × 14 mm), referred to as trabecular bone ROI, was placed immediately below the dense subchondral trabecular bone. Trabecular bone ROIs were aligned horizontally with the subchondral bone ROIs. Anatomical landmarks for the ROIs were tibial spine, subchondral bone plate, and outer borders of the proximal tibia. A custom-made MATLAB software (version R2014b, The MathWorks, Inc., Natick, MA, USA) was used for the manual placement (JH) of the ROIs. Reproducibility of the texture parameters from the aforementioned locations has been shown to be high [14].

Structural analysis of bone
Prior to the bone structural analysis, radiographs were medianfiltered (3 × 3 pixels) to remove high-frequency noise and grayscale values of the image were expanded to full dynamic range. Bone structure was evaluated from the radiographs using Laplacian-based methods [5,14,27], LBP-based methods [5,28] and FSA [5,23,29].

Laplacian-based analysis
The Laplacian-based method [5,27] enhances the appearance of bone trabeculae and quantifies the variation in the grayscale values of the Laplacian-based image. Laplacians were calculated in the horizontal and vertical directions and summed into the one matrix. Subsequently, the unprocessed ROI was multiplied by the square root of the Laplacian matrix to enhance the bone, and grayscale values were expanded to the full dynamic range to obtain the final Laplacian-based image. To measure the randomness of the grayscale values in the Laplacian-based image, entropy of the image (E Lap ) was calculated using the following equation:

Local binary patterns (LBP)-based methods
LBP-based methods were used to measure the randomness of local patterns and the variations in the orientation of adjacent local patterns. First, the image was divided into bone and non-bone regions by determining a local threshold for every image pixel using the Otsu method [30] with a 9 × 9-pixel window size. Next, the LBP operator (eight-neighborhood on a circle with a radius of 1) was applied in the bone regions and the bone edges (i.e., nonbone regions next to the bone). The pixel was considered to be an edge pixel if at least one of the eight neighbors of the center pixel was a bone pixel. To reduce the number of irrelevant patterns, grouping of patterns was carried out by determining the main orientation and the number of valid neighbors (i.e., markers) for each pattern. The main orientation angle was calculated using principal component analysis. The angles (0°, 45°, 90°, and 135°) were calculated only for the patterns consisting of two to five consecutive markers; otherwise, the patterns were assigned as non-uniform.
To measure the randomness of the patterns occurring in the image, entropy of the grouped patterns (E LBP ) was determined using the Eq. 1. If E LBP = 0, there is only a single pattern occurring in the image.
The homogeneity index for the orientation of the valid patterns was derived from the co-occurrence matrix of the angles. The co-occurrence matrices were calculated in 0°, 45°, 90°, and 135°directions with one pixel distance. The non-uniform and non-bone areas were excluded from the cooccurrence matrices. To take into account the orientation of bone trabeculae in the analysis, co-occurrence matrices of 0°a nd 135°directions were combined together as well as 45°a nd 90°directions to calculate the homogeneity index perpendicularly to the bone trabeculae (HI Angles,Perp ) and along the trabeculae (HI Angles,Paral ), respectively. Furthermore, the mean homogeneity index (HI Angles,mean ) was calculated from the cooccurrence matrix as the sum of the four possible directions. The interpretation of the HI Angles parameters used in this study is the following: If all adjacent patterns have similar orientation, HI Angles is equal to one, while a large variation in the orientation of local patterns results in a low HI Angles value.

Fractal signature analysis (FSA)
To estimate fractal dimension, related to roughness and complexity of an image, the FSA method was used [23,29]. In brief, the image was dilated and eroded in horizontal and vertical directions with a rod-shaped one-pixel-wide structuring element and the volume, V, between dilated and eroded images was calculated. Calculations were repeated by varying the element length r from two to four pixels. The surface area, A(r), was calculated using the Eq. 2: Subsequently, a log-log plot was constructed by plotting log of A(r) against log of r. Finally, the fractal dimension was estimated using a regression line to points between 2 and 4. We restricted our analyses only to the small-scale fractal dimensions in order to keep the results compact. When the structuring element is pointing in the horizontal direction, fractal dimension of vertical structures (FD Ver ) is produced and vice versa. High fractal dimension values are associated with high complexity of the image, whereas low complexity results in low fractal dimension values. acquisition matrix: 320 × 320, slice thickness: 3 mm). A musculoskeletal radiologist (AG, 15 years of experience in semiquantitative MRI analysis of knee OA) scored the medial articular cartilage damage and BMLs using the MRI OA knee score (MOAKS) [16]. A subject was included in the medial tibial cartilage damage group if he/she had any cartilage loss [defined as MOAKS grades 1-3 for the size of any cartilage loss (including partial and full-thickness loss)] in the medial anterior, central, or posterior tibia ( Fig. 2A). Similarly, a subject was included in the medial tibial BML group if he/she had any BML (MOAKS grades 1-3 for the size of BML by volume; including ill-defined lesions, bone marrow edema and subchondral cysts) in the medial anterior, central, or posterior part of tibia (Fig. 2B).

Statistical analyses
Analysis of covariance was used to compare the bone structure between groups with and without cartilage damage or BMLs. Bone structural parameters were adjusted with gender, age and body mass index. Unadjusted mean values (standard deviations) are shown in the result tables. Mean KL grade and its standard deviation are reported in Tables 2 and 3 for clarity although the KL grade is an ordinal scale variable. Difference in the KL grade between groups was tested using the Mann-Whitney U test. Statistical analyses were conducted using IBM SPSS Statistics for Windows (Version 22.0, IBM Corp., Armonk, NY, USA).

Results
Subjects with medial tibial cartilage damage had significantly different bone structure in medial tibial subchondral and trabecular bone ROIs (Table 2). In the medial subchondral bone region, FD Ver (Fig. 3A) and E LBP were higher, whereas E Lap , HI Angles,mean (Fig. 3B), HI Angles,Perp , and HI Angles,Paral were lower among subjects with medial tibial cartilage damage ( Table 2). In the medial trabecular bone region, FD Hor , FD Ver (Fig. 3A), and E LBP were higher, whereas E Lap , HI Angles,mean (Fig. 3B), HI Angles,Perp , and HI Angles,Paral were lower among subjects with medial tibial cartilage damage ( Table 2).
Subjects with medial tibial BMLs had also significantly different bone structure in medial tibial subchondral and trabecular bone ROIs (Table 3). In the medial subchondral bone region, FD Ver (Fig. 3C) and E LBP were higher, whereas E Lap , HI Angles,mean (Fig. 3D), HI Angles,Perp , and HI Angles,Paral were lower among subjects with medial tibial BMLs (Table 3). In the medial trabecular bone region, FD Hor , FD Ver (Fig. 3C), and E LBP were higher, whereas E Lap , HI Angles,mean (Fig. 3D), HI Angles,Perp , and HI Angles,Paral were lower among subjects with medial tibial BMLs (Table 3).
When tibial cartilage damage analyses were adjusted for medial tibial BMLs along with gender, age, and body mass index, FD Ver was higher (p = 0.003) and E Lap was lower (p = 0.023) among subjects with medial cartilage damage in the medial subchondral bone region. In the medial trabecular bone region, FD Hor (p = 0.016) and FD Ver (p = 0.042) were higher, whereas E Lap was lower (p = 0.023). When BML analyses were adjusted for medial cartilage damage along with gender, age, and body mass index, significantly higher HI Angles,mean values (p = 0.036) between the control and medial tibial BML groups was detected.
When tibial cartilage damage analyses were adjusted for KL grade, gender, age, and body mass index, FD Ver was higher (p = 0.001) and E Lap was lower (p = 0.008) among subjects with medial cartilage damage in the medial subchondral bone region. In the medial trabecular bone region, FD Hor (p = 0.020), FD Ver (p = 0.005) and E LBP (p = 0.030) were higher, whereas E Lap was lower (p = 0.009). When BML analyses were adjusted for the KL grade gender, age, and body mass index, only HI Angles,Perp was significantly lower (p = 0.041) in medial subchondral bone ROIs, and FD Ver was significantly higher (p = 0.006) in medial trabecular bone ROIs among subjects with medial tibial BMLs.

Discussion
Our results demonstrate that medial tibial subchondral and trabecular bone structures from 2-D plain radiographs were different between subjects with and without medial tibial articular cartilage damage and between subjects with and without medial tibial BMLs.
The finding that trabecular bone structure from radiographs is different between subjects with and without cartilage damage is in line with an earlier study [31]. In that study, the fractal dimension of horizontal structures in the lateral compartment and fractal dimension of vertical structures in both medial and lateral compartments were significantly higher among   subjects with cartilage defects (in medial, lateral, or both compartments) compared to subjects without cartilage defects [31]. As a complement to that study, the subchondral bone area was also analyzed in the current study. Furthermore, we used novel LBP-and Laplacian-based algorithms along with FSA to assess bone structure. As a novel finding, we were able to confirm that, in addition to trabecular bone structure, the structure of subchondral bone assessed from 2-D radiographs is also different between subjects with and without medial cartilage damage. Consequently, the assessment of the subchondral bone area immediately under the cartilage-bone interface should be considered in future studies using plain radiographs as that area has an important role in OA pathology [25,26]. Bone structural parameters in both subchondral and trabecular bone ROIs were different between subjects with and without BMLs in medial tibia. The higher fractal dimension and E LBP and lower HI Angles values suggest that the bone structure is more disoriented among subjects with medial tibial BMLs. In a previous MRI study, subjects with BMLs had higher apparent bone volume fraction (calculated from MRI data), higher trabecular thickness and number, and lower trabecular separation [32]. The current results are in line with that study since, for example, higher FD Ver and E LBP and lower HI Angles,mean values are known to be related with lower trabecular separation values [5]. It is known that BMLs are more common among subjects with advanced OA and this might explain partly the observed differences in the calculated bone structural parameters. It should also be noted that all subjects with BMLs had concurrent cartilage damage which further limits our BML-related analyses.
In this study, we showed that FSA, Laplacian-based and LBP-based methods are able to detect differences between subjects with and without cartilage damage or BMLs. One explanation for the differences in structural parameters of bone is that subchondral and trabecular bone in the proximal tibia is not as well-organized among subjects with cartilage damage or BMLs. This explanation is supported by the previous studies which have proposed that the trabecular bone is more disorganized and the orientation of trabecular bone is shifted, probably due to abnormal loading, in OA [33][34][35]. With the FSA method, fractal dimensions of vertical and horizontal structures are typically reported [12,29,[36][37][38]. However, bone trabeculae are not perfectly oriented vertically in the proximal tibia. To account for this, HI Angles was calculated approximately perpendicularly and parallel to bone trabeculae. It should be noted that due to differences in calculation of fractal dimension and HI Angles values, FD Ver and HI Angles,Perp (and FD Hor and HI Angles,Paral ) are related to each other, although these variables are not directly measuring the same phenomenon.
One limitation of our study is the difference in imaging conditions of radiographs. The conditions varied because the study subjects were imaged at different institutions on different imaging systems. Since texture analysis methods used in this study are affected by pixel size, the comparability of the calculated structural parameters between different images was ensured by resampling the pixel size of all images to 143 × 143 μm 2 . Furthermore, texture analysis methods are not very sensitive for the differences in imaging conditions (e.g., tube voltage and current), as the methods do not usually evaluate plain grayscale values directly. Another limitation is that many of the subjects had multiple concurrent tissue changes in MOAKS. However, due to limited sample size and restricted knowledge on the most important factors for the bone structure appearance as determined from radiographs, we decided not to adjust bone analyses with MOAKS features. Furthermore, in the OKOA data, radiographs were available only for symptomatic subjects and, thus, our Bcontrol^groups might have had some OA-related changes. On the other hand, it is likely that the differences would have been even higher if non-symptomatic true controls were used. Finally, since our study design was cross-sectional, we could not study the causality of the tissue changes; i.e., it is not possible to determine whether the cartilage damage or BML induces the alteration in the subchondral bone structure or vice versa, for example. The longitudinal relevance of the subchondral bone structure analysis from plain radiographs warrants further exploration.
Plain radiography is a cheap, fast, and widely available imaging method suitable for imaging of large subject cohorts. Analysis of subchondral bone structure can be performed relatively fast using the presented methods. The only manual user input needed is to place the ROIs in the correct predefined location, which takes about one minute depending on the experience of the user. The remaining analyses are performed by the software algorithm. Semi-automated and automated methods have also been developed for the ROI placement [39][40][41]. In this study, we confirmed that subchondral and trabecular bone structure from radiographs is different among subjects with cartilage damage or BMLs, which both are significantly related with OA incidence and progression [17][18][19][20][21]. Therefore, our results suggest that structural analysis of bone from radiographs can be used as a supplementary tool when evaluating subjects with OA or at risk of having it.

Compliance with ethical standards
Guarantor The scientific guarantor of this publication is Simo Saarakkala.

Conflict of interest
The authors of this manuscript declare relationships with the following companies: Ali Guermazi is the President of Boston Imaging Core Lab (BICL), LLC and a consultant to MerckSerono, Genzyme, OrthoTrophix and TissueGene. Frank W Roemer is the CMO and a shareholder of BICL. The other authors declare that they have no competing interests.
Funding This study has received funding by the Radiological Society of Finland, the Orion Corporation Research Foundation, Päivikki and Sakari Sohlberg Foundation, Academy of Finland (project 268378), University of Oulu (strategic funding), Sigrid Juselius Foundation, and European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement no. 336267. The funding sources had no role in the study design, data collection or analysis, interpretation of data, writing of the manuscript, or in the decision to submit the manuscript for publication.
Statistics and biometry Marianne Haapea, Ph.D. kindly provided statistical advice for this manuscript.
Ethical approval Institutional review board approval was obtained.
Informed consent Written informed consent was obtained from all subjects (patients) in this study.
Study subjects or cohorts overlap The study subjects have been previously reported in Scientific Reports by Podlipská et al. 2016 [22].

Methodology
• retrospective • cross sectional study • performed at one institution Open Access This article is distributed under the terms of the Creative Comm ons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.