Participants
Ten participants with RA in physician-classified remission were recruited using the Rheum4U platform (19) from the Division of Rheumatology at the University of Calgary. In order to be included in the study, participants had to be over the age of 18 and meet the American College of Rheumatology/European League against Rheumatism 2010 Classification criteria for rheumatoid arthritis diagnosis (20). Participants with a history of MCP injury, replacement, or whose conventional radiographs showed no joint space in the MCPs were excluded.
HR-pQCT Image Acquisition and Analysis
HR-pQCT scans of the 2nd and 3rd MCP joints of the participants’ hand were acquired (XtremeCTII, Scanco Medical, Brüttisellen, Switzerland) at baseline and at a 6-month follow-up. Quality control scans were performed daily using a phantom with rods of known density. A same day repeat scan after repositioning was obtained at the follow-up visit. The hand was secured in a custom positioning device and used an adaptation of an established protocol (21) to reduce the impact of motion on image processing. The scan was acquired with a nominal isotropic resolution of 60.7 mm using the manufacturer standard settings (68 kVp, 1470 mA, 43 ms integration time). We obtained a reference x-ray in the coronal plane, and a reference line was placed at the distal cortical surface of the 2nd metacarpal head. Three 10.2 mm sections or “stacks” were acquired separately with a 25% (2.55 mm) overlap of the previous stack resulting in a total scan length of 25.5 mm. A 25% overlap was selected after testing because it allowed for consistent registration results, while maintaining an adequate scan distance to include both the proximal phalange and distal metacarpal bones across both the 2nd and 3rd MCPs. Each stack was evaluated for motion using the manufacturer’s standard scoring system from 1-5 (22). Data that included a motion score of 4 or 5 were excluded from the longitudinal analysis and scan-rescan analysis. Erosions were identified by a trained reader. Image processing was completed using Image Processing Language (IPL v5.42, Scanco Medical). The analysis scripts generated and used during this study will be made available through the Manske Lab GitHub repository.
Multi-Stack Registration
The proximal, middle, and distal stacks of the 2nd and 3rd MCP joints were registered in order to achieve optimal alignment. Each stack was registered using the 25% overlap of the previous stack. The registration was initialized using the mass centers and moment of inertia to give an estimate of the new orientation. Once the registration was initialized, the proximal stack was registered to the middle stack using a simplex optimization to determine rigid transformation parameters by maximizing the correlation coefficient. The proximal stack was then transformed to the middle stack image space using a cubic interpolator. The same procedure was followed to register the distal stack to the middle stack. Masks of the overlapped regions were created and used to remove the common regions with the middle stack from the distal and proximal stacks. Finally, the transformed proximal and distal stacks (excluding the overlapped region) were concatenated to the middle stack (including the overlapped region) creating an image of the three stacks optimally aligned. The final image of the joint includes the proximal and distal stacks registered, transformed, and aligned with the middle stack to give the appearance of contiguous bone, with the stack shift artifact eliminated (Figure 1).
Longitudinal Image Registration
In order to investigate bone damage and healing phenomena in RA, longitudinal 3D rigid body image registration was implemented to align and superimpose the follow-up image with the baseline image. The first step in the image registration workflow was to determine a transformation matrix to rotate and translate the follow-up image to the baseline image space. For the purpose of this registration, the periosteal mask of either the metacarpal or phalange bone was dilated and then used to crop the desired bone, with a 5-voxel border of background. From there, the center of mass and moment of inertia was used to initialize the registration. This rigid registration used a simplex optimizer approach and a correlation object function to iteratively find the optimal fit between the follow-up and the baseline image. The transformation matrix was then applied to the follow-up image to create a translated and rotated image that was resampled using cubic interpolation. Finally, the images were masked within the largest common volume (LCV) to exclude any voxels outside the common region.
Bone Remodeling Analysis
Once the baseline and follow-up images were aligned in the same image space, they were analyzed to find regions of bone formation and bone loss. This was accomplished by subtracting the grey-scale density images from one another voxel-by-voxel within the LCV based on a method reported previously (18). The resulting image therefore represents the differences in densities between the two images, which can represent areas of local bone remodeling. The follow-up image was subtracted from the baseline image and then segmented with a threshold value of 125 mgHA/cm3 so that any voxels above this threshold would represent bone resorption. The baseline image was subtracted from the follow-up image and then segmented so that any voxels greater than the threshold value of 125 mgHA/cm3 would represent bone formation. To eliminate noise from the difference image, a connected components filter was then applied so that clusters of voxels connected by fewer than 5 voxels were removed from the image. These segmented regions of bone loss and bone gain were then overlaid with the segmented baseline scan to create a comprehensive image depicting bone loss and bone gain (Figure 2). The threshold and connected components cluster size was determined based on previous work (17,18), and adapted for the second generation HR-pQCT scanner and MCP joints to provide a plausible depiction of bone formation and bone resorption based on visual inspection of the grey-scale images at baseline and follow-up as well as the final segmented image. Bone formation and resorption fractions were then quantified as a percentage of the total voxels labelled as formation or resorption from the respective difference image. These formation and resorption fractions were also calculated using the same bone remodeling script for the same day scan-rescan after repositioning. This was used as a negative control to ensure that any changes observed over the 6-months were due to bone change rather than repositioning in the scanner or image artifact.
Effects of Image Interpolation on Bone Remodeling Assessment
After registering the follow-up (moving) image to the baseline (fixed) image, the moving image was transformed to the fixed image space. Following this transformation, cubic interpolation was performed to realign the moving image to the fixed image’s coordinate system. Image interpolation can be seen as a linear convolution with a low-pass filter, and therefore, may have an effect on bone morphometric analysis (23,24). To analyze the effects of interpolation on the bone remodeling analysis, the bone remodeling values were compared for registration of the 6-month follow-up scans to their respective same-day rescans where both scans underwent a cubic interpolation and where only the moving image underwent interpolation.
Statistical Analysis
Descriptive statistics were used to observe if there were any improvements in the number of erosions that could be analyzed with multi-stack registration, compared to scans without multi-stack registration. Kruskal-Wallis tests were used as a non-parametric method to compare the differences in bone formation and resorption fraction across different motion grades against an alpha value of 0.05. A Wilcoxon signed rank test was used to compare the difference in bone formation and resorption fraction between registrations with only the moving image interpolated compared to registrations where both scans undergo interpolation against an alpha value of 0.05. Results are reported as median (interquartile range, IQR) unless otherwise stated. Analysis was completed using R (v3.4.3) in RStudio (v1.1.423).