Micro-CT based finite element models of cancellous bone predict accurately displacement once the boundary condition is well replicated_ A validation study

Non-destructive 3D micro-computed tomography (microCT) based finite element (microFE) models are used to estimate bone mechanical properties at tissue level. However, their validation remains challenging. Recent improvements in the quantification of displacements in bone tissue biopsies subjected to staged compression, using refined Digital Volume Correlation (DVC) techniques, now provide a full field displacement information accurate enough to be used for microFE validation. In this study, three specimens (two humans and one bovine) were tested with two different experimental set-ups, and the resulting data processed with the same DVC algorithm. The resulting displacement vector field was compared to that predicted by microFE models solved with three different boundary conditions (BC): nominal force resultant, nominal displacement resultant, distributed displacement. The first two conditions were obtained directly from the measurements provided by the experimental jigs, whereas in the third case the displacement field measured by the DVC in the top and bottom layer of the specimen was applied. Results show excellent relationship between the numerical predictions (x) and the experiments (y) when using BC derived from the DVC measurements (UX: y=1.07x−0.002, RMSE: 0.001 mm; UY: y=1.03x−0.001, RMSE: 0.001 mm; UZ: y=x+0.0002, RMSE: 0.001 mm for bovine specimen), whereas only poor correlation was found using BCs according to experiment set-ups. In conclusion, microFE models were found to predict accurately the vectorial displacement field using interpolated displacement boundary condition from DVC measurement.


Introduction
Bone tissue is a complex hierarchical material (Cowin, 2001). In order to address clinical and preclinical problems, it is important to study bone at the spatial scale that allows the most appropriate characterization of its mechanical behaviour. At the tissue scale, the interaction between bone mechanical stimuli and the biological function driven by the cell activity becomes more evident (Viceconti, 2012). The microCT based finite element (microFE) method has become a popular tool for non-destructive structural analysis of cancellous bone tissue (Hollister et al., 1994;van Rietbergen et al., 1995;Verhulp et al., 2008). The method involves the direct conversion of the 3D voxels of micro-computed tomography (microCT) images of the bone tissue (Feldkamp et al., 1989) into equally shaped and sized hexahedral elements. As microCT imaging has the ability to accurately resolve bone morphology in great detail (Bouxsein et al., 2010), specimenspecific microFE models that represent the structure of the specimen can be generated (Ulrich et al., 1998).
Every modelling method requires, before it can be considered reliable, a complete verification, validation, and uncertainty quantification assessment (Anderson et al., 2007). A systematic verification analysis of microFE models of bone tissue was recently published (Chen et al., 2014). However, for validation, the number of published reports is limited. While the predicted apparent properties (e.g. stiffness, strength) of each specimen can be compared with accurate experimental measurements (Christen et al., 2013;Pistoia et al., 2002;Wolfram et al., 2010;Yeni and Fyhrie, 2001), the validation of such models for local predictions is not trivial.
One possible approach is to use a full-field method, such as Digital Volume Correlation (DVC) techniques, to extract displacement or strain fields from repeated microCT scans performed during stepwise-compression experiments. In every validation study the boundary conditions (BCs) imposed in the model should be the same as in the experiments. However, even if we can measure accurately the resultant force applied, or the total displacement imposed during the stepwise compression test, the aspect ratio of the specimens typically used in these tests might be too small to assume valid Saint-Venant‫׳‬s Principle. If this is the case, it is not enough to reproduce in the model the loading resultant, but we need to consider also how such forces are locally distributed. To the authors' knowledge there are two studies in the literature that used DVC measurements to validate microFE models displacement predictions on cancellous bone specimens, and their results are somehow inconclusive. Recently, Zhu et al. (2015) compared the predictions of microCT based tetrahedral homogeneous models to DVC measurements for bovine bone interdigitated with acrylic cement and a cellular foam samples. Only qualitative comparison between models and DVC displacement only along the loading direction was reported. Zauel et al. (2006) was the first to use a DVC approach based on the one reported in Bay et al. (1999) to quantitatively validate a linear elastic microFE model of cancellous bone. They found very good correlation in displacement measured along the major loading direction (R 2 from 0.91 to 0.97, slopes between 0.93 and 0.98), but only poor correlation for transverse displacements (R 2 from 0.29 to 0.60, slopes between 0.33 and 0.88). This result is surprising as the precision error of their DVC method is isotropic (Dall‫׳‬Ara et al., 2014;Zauel et al., 2006) and the predictions of the microFE models should not be affected by the loading direction unless strong local anisotropy needs to be included in the models. The findings reported by Zauel et al. (2006) suggest that homogeneous isotropic microFE models is not reliable in predicting transverse displacement. Therefore, we need to further explore the ability of predicting local displacement from microFE models, widely used in the research community to estimate bone properties at the tissue levels.
The aim of this study was to determine the accuracy of microFE models, generated from microCT data, in predicting the displacements of cancellous bone specimens subjected to compression, when sufficiently accurate full-field displacement measurements are used for the validation, and appropriate boundary conditions are simulated.

Materials and methods
Two independent testing procedures were used in order to assess the sensitivity of the validation approach for two different experimental protocols and input images and to extend the validity of the results. In both cases similar workflows ( Fig. 1) were used.

Specimen preparation and scanning
All procedures on human tissue were performed with approval from the Research Ethics Committee for use of discarded bone material (LREC 2002/1/22).
Two cylindrical specimens (Specimen1: height equal to 13.2 mm, diameter equal to 10.6 mm; Specimen2: height equal to 11.5 mm, diameter equal to 10.6 mm) were extracted from the central part of two human femoral heads from patients who underwent total hip replacement using diamond-tipped cores (Starlite Indistries, Rosemount PA, USA). Specimen 1 was from an osteoarthritic male aged 68 and Specimen 2 from a 94 years old male. The ends of the core samples were cut parallel using a Buehler Isomet low speed saw (Buehler, Illinois, USA).
The third cylindrical specimen (Specimen3: height equal to 11.88 mm, diameter equal to 7.89 mm) was drilled (diamond core drill with nominal internal diameter equal to 8 mm mounted on a pillar drilling machine, GDM50B, Sealey, UK) from a bone slice cut (0.2 mm diamond bandsaw mounted on a 300 CP, Exakt Gmbh, Germany) from a bovine femoral greater trochanter (female, 18 months old). All operations were performed under constant water irrigation in order to reduce potential damage to be bone specimen. Animal tissue was extracted from a bovine femur, collected from an animal culled for alimentary purposes.
The specimens were then microCT scanned (Skyscan 1172; Specimen1 and Specimen2: voxel size 17.22 µm, 54 kV, 185 µA, 0.5 mm aluminium filter, exposure time 885 ms, no averaging; Specimen3: voxel size 9.92 µm, 59 kV, 169 µA, 1 mm aluminium filter, exposure time 1180 ms, averaged by two frames). Each image was cropped in order to include only the bone specimens and datasets were subsampled by a factor of two (ImageJ, V1.50a), resulting in a new voxel size equal to 34.44 µm and 19.84 µm for human and bovine specimens, respectively. For the Specimen1 and Specimen2 top and bottom slices with partial bone and air were removed, while for Specimen3 slices in the embedding material were removed. Bone volume fraction (BV/TV), trabecular thickness (Tb.Th), trabecular spacing (Tb.Sp), degree of anisotropy (DA) and angle between the main trabecular direction and the loading axis (α.Z) were computed with the ImageJ plugin BoneJ (Doube et al., 2010). Information for all specimens is summarized in Table 1.

In situ mechanical testing
All three specimens were mechanically tested within the microCT system. Specimen1 and Specimen2 were fully hydrated, and tested with the compressive device provided by the manufacturer of the microCT with a 440 N load cell. The specimens were positioned in between two parallel loading plates, in the middle of the device. A first scan (undeformed) was performed with the specimens under a small preload of 7.22 N in order to avoid motion artefacts. Afterwards, a compressive step up to 1% apparent strain was applied without repositioning and the specimen was scanned in its deformed configuration ( Fig. 2, left).
Specimen3 was tested in a custom-made compressive device to be positioned within the same microCT model. The load was applied by a manual screw-ball joint mechanism and was measured with a 2 kN load cell (LPM530, Cooper Instruments & Systems, Warrenton, USA). The 1.5 mm external portions of the specimen were embedded in PMMA (Technovit 4071, Heraeus Kulzer Gmbh, Wehrheim, Germany) after proper alignment with the loading axis of the jig. A first scan (undeformed) was performed with the specimen under a preload of 2 N in order to avoid motion artefacts. Afterwards, a compressive step to 120 N was applied without repositioning and the specimen was scanned in its deformed configuration. A liquid-filled chamber was used in order to keep the specimen submerged in a 0.9% NaCl solution during the test (Fig. 2, right).

DVC measurement of displacement
The DVC method computes the field of displacements by registering elastically the couple of undeformed and deformed images for each specimen (34.44 µm and 19.84 µm). In the present study, we used an in-house elastic registration library (Sheffield Image Registration Toolkit, ShIRT) (Barber and Hose, 2005;Barber et al., 2007;Khodabakhshi et al., 2013). The registration equations are solved in the nodes of a grid superimposed to both images to be registered and with certain nodal spacing (NS), assuming a linear behaviour in displacement in between the nodes. In the current study, we used NS equal to 12 voxels (~413 µm) for human cancellous bone (Specimen1 and Specimen2) and NS 25 voxels (~496 µm) for bovine cancellous bone (Specimen3). With this NS the accuracy and precision in displacement is approximately 0.00016 ± 0.0034 µm (~400 µm) for Specimen1 and Specimen2 and 0.0000098 ± 0.00014 µm (~500 µm) for Specimen3 (Palanca et al., 2015). As for all other methods described in the literature, also the accuracy of our method was calculated by using virtually moved images; it should be noted that this test provides an estimate of the lower boundary of the error. Still, the method used here is significantly more accurate than any other described previously.

microFE models
To reduce the computational cost, the original image datasets were first subsampled by a factor of two as for the DVC measurements (34.44 µm for Specimen1 and Specimen2; 19.84 µm for Specimen3). All image datasets were filtered by a Gaussian filter with σ=1.2 and a support of 3 voxels to reduce the high frequency noise (Christen et al., 2014). Then the images were binarized using a single-level threshold by finding the mean value between two peaks (one representing the bone tissue, one representing the background) in the grayscale histograms. Voxels below the threshold value were deleted and for those above the threshold value, a connectivity filter (Matlab, R2014b, Mathworks, Inc.) was applied to remove the isolated voxels. In particular, only elements with surface (four nodes) connectivity were kept in the model. Finally, each remaining voxel in the image datasets was converted directly into equally sized 8-node hexahedral elements (Fig. 3). The material properties for all microFE models were assumed to be linear and isotropic, with a uniform Young‫׳‬s modulus of 17 GPa (Bayraktar et al., 2004;Morgan et al., 2003) and a Poisson‫׳‬s ratio of 0.3 (Pistoia et al., 2002). Following the procedure by Chen et al. (2014), a convergence study was performed on models with three BCs (mentioned in Section 2.5). All models (element size equal to 34.44 µm for Specimen1 and Specimen2 and 19.84 µm for Specimen3) reached a local convergence in displacement, with a difference of less than 1% compared to models of the most mesh refined microFEs with element size equal to the original microCT voxel size (17.22 µm for Specimen1 and Specimen2 and to 9.92 µm for Specimen3).

Boundary conditions
Three different BCs for each microFE simulation were used. The first set of BCs (hereinafter referred to as "force BCs" was defined based on the experimentally measured force. The layer of nodes of the microFE belonging to the fixed surface of the bone were fixed in the loading direction and a vertical force was distributed equally on each node of the loaded surface so that their resultant was equal to the measured axial force (42N for Specimen1, 162N for Specimen2 and 120N for Specime3). The nodes of the both loaded and fixed surface were free to move in transverse direction for the Specimen1 and   Specimen2 (free boundary conditions) and fixed in the transverse directions for Specimen3 (simulation of embedding). The second set of BCs (later referred to as "displacement BCs") was defined based on the experimentally measured displacements. The layer of nodes of the microFE belonging to the fixed surface of the bone were fixed in only the loading direction and the experimentally measured vertical displacement was imposed to each node of the loaded surface (130 µm for Specimen1, 115 µm for Specimen2; for Specimen3 no displacement was measured and this type of BCs was not modelled), leaving free to move in the transverse directions.
The third set of BCs (later referred to as "interpolated BCs") was imposed from the DVC measurements. The top and bottom layers of nodes from the microFE were superimposed into the displacement field measured from the DVC approach. The displacements of the boundary nodes were assigned by linear interpolation using element shape functions.

Comparison between experimental and computational results
The comparison was limited to the middle 80% of the specimen to avoid boundary effect. The DVC procedure provides the displacements at the nodes of its grid (with NS 12 or 25 voxels according to the specimen) over the whole image. Due to way the DVC grid and of the microFE models are generated, each node of the DVC grid lays either in the centroid of a microFE element or in the marrow. The comparison of the results was performed only in the locations of the DVC nodes corresponding to microFE elements (i.e. in the bone). For those locations, the displacements measured with the DVC were compared to the interpolated value in the centroid of the element of the microFE. Following this procedure, we obtained the following number of comparison pairs: 4041 for Specimen1, 3671 for Specimen2 and 589 for Specimen3. In each specimen, we compared the displacements predicted by microFE under the different BCs to the DVC results.

Statistics
Any observation with Cook‫׳‬s distance (Fox and Long, 1990) larger than five times the mean Cook‫׳‬s distance was considered as outliers and removed from the analysis. This approach removed 1% to 4% points for each analysis. The comparison of displacement for microFE models and experiments was performed using linear regression, where the slope and intercept of the equation as well as coefficient of determination (R 2 ) were reported. For each comparison the Root mean square error (RMSE), the RMSE divided by the maximum experiment value (RMSE%), the largest difference between microFE prediction and DVC measurements (Max.error), the Max.error divided by the maximum experimental value (Max.error%), and the intra-class correlation (ICC) were computed. The number of comparison pairs remained in each analysis was reported in Table 2.

Results
All coefficients calculated from the correlations between predicted and measured displacements are reported in Table 2.
The displacement predicted by microFE models with three different BCs and DVC measurements were all significantly linearly correlated (p < 0.01). However, microFE models with "force BCs" and "displacement BCs" were far from the 1:1 relationship, underlined by the low ICCs (from 0.02 to 0.42). Conversely, microFE models with "interpolated BCs" lead to excellent correlations, with slope close to one (range: 0.98 to 1.07), intercept close to zero (range: −0.006 to 0.006 mm), high R 2 (range: 0.97 to 0.99) and high ICC (0.99). In that case similar results were found for the three specimens (Fig. 4), with RMSE% lower than 2.5% (with maximum equal to 2.4% for predictions of Uz for Specimen2 and of Uy for Specimen3) and Max.err% lower than 11% (with maximum equal to 10.7% for predictions of Ux for Specimen2). These models overall predicted better U Z (displacements along the major compression direction) with RMSE% from 1.1% to 2.4% and Max.err% from 3.5% to 5.6%) compared to the displacements along the transverse directions (RMSE% from 1.7% to 2.1% and Max.err% from 5.6% to 10.7% for U X and RMSE% from 1.7% to 2.4% and Max.err% from 5.2% to 7.1% for U Y ). The best correlation was found for predictions of U Z for Specimen3, with slope equal to 1 and intercept equal to 0.0002 mm. In Fig. 5 comparisons between the predicted and measured vertical displacement for the "interpolated BCs" are reported. From those graph it can be noted that all specimens have rotated to some extent during the experiments, especially Specimen1 and Specimen2.
The scatter plots for all three specimens using "force BCs" and "displacement BCs" can be found in the Supplementary materials (available through the FigShare on-line https://figshare.com/s/ 0fa91152b09dfe5c5142).

Discussion
The aim of the present study was to compare DVC measurements of displacement in cancellous bone samples with the predictions of microFE simulated with different BCs.
All microFE models using "interpolated BCs" accurately predicted the displacement field measured experimentally, as induced by loads small enough to assume linear elasticity, with similar accuracy in all spatial directions. It should be noted that similar results were found over experiments generated by two separate research groups, using different experimental protocols, and imaged at different resolutions, underlying the robustness of the method.
On the other hand, microFE models using "force BCs" and "displacement BCs" in general predicted displacement poorly, when compared to DVC measurements. These poor predictions could be due to the fact that the friction between loading plates and surface of the sample (which in some studies is artificially reduced, but that can never be nullified) was not modelled (for "force BCs"), or by the potential parallelism error between the two flat surfaces of the specimens. ("displacement BCs", for Specimen1 and Specimen2 which were not embedded) and/or by the potential micro-movements of the loading plate ("displacement BCs", all three specimens). This effect does not seem to be related to the different inclination between the main trabecular direction of each sample and the loading direction (from 9°to 60°), at least for relatively dense samples as those considered in this study (23% < BV/TV < 30%). The excellent correlations found for displacements along the  . However, in that case much worse predictions were found for the transverse directions (R 2 from 0.29 to 0.60; slope from 0.33 to 0.88; intercept from −954 to 40 µm) while in the present study they were excellent, even if with slightly larger scatter (R 2 from 0.97 to 0.99; slope from 1.00 to 1.07; intercept from −6 to 6 µm). The small differences between the predictions of the displacements along the different directions found in our study underline that the assumption of isotropic material property for every element in case of microFE models is well posed. The improvement in the predictive accuracy of the displacements along the transverse directions in this study compared to the previous report (Zauel et al., 2006) may be due to the improved accuracy of the DVC method used in this work (Dall‫׳‬Ara et al., 2014;Palanca et al., 2015). It remains to be investigated how sensitive the models are for strain predictions.
The main limitation of this study is that it validates microFE models only with respect to displacements, whereas in many studies these models are primarily used to predict strains. However, the accuracy and precision of the current experimental methods for strain measurement in each element (as small as 10-20 µm) when modelling whole bone biopsy (10-20 mm large) are too low (Grassi and Isaksson, 2015;Palanca et al., 2015). All available DVC methods can provide reasonably accurate quantifications of strains only over much coarser resolutions (i.e. 40-50 voxels); while in principle we could sub-sample the microFE model at the same lower resolution, this comparison would be meaningless, as we would need to assume bone as a continuum. Such validation has been already done at the organ level by using straingauges (Cristofolini and Viceconti, 2000;Cristofolini et al., 1996;Taddei et al., 2010;Viceconti et al., 2001) or digital image correlation (Grassi et al., 2016). Further research is required to produce a controlled experiment where the strain field can be measured with sufficient accuracy with a special resolution of 10-20 µm. Here we investigated the simplest type of microFE models, based on Cartesian homogenous elements. It might be interesting to explore if the predictive accuracy is further improved by adding local heterogeneity (Gross et al., 2012), or local smoothing (Muller and Ruegsegger, 1995) to overcome the inevitable partial volume effect induced by homogeneous models. Moreover, if larger compressive loads would be analysed, nonlinear models (Harrison et al., 2013) would become fundamental in order to simulate the local yielding of the trabecular bone structure. Finally, more specimens should be tested in order to increase the applicability of the model, possibly with increased aspect ratio in order to reduce the effect of the loading plates on the displacement field. All these improvements and further application to different bone microstructures are currently under development.
In conclusion, microFE models with "interpolated BCs" predict local displacements in cancellous bone samples with excellent accuracy in all spatial directions. Having provided consistent results across all three specimens tested by different experimental set-ups, we can conclude microFE models are robust once the BC assigned is well matched with experiments and, therefore, can be used as a modelling framework in the future to study bone mechanics at the tissue level, at least to predict displacement. It remains to be investigated whether this hold valid for the local strain and if improvement in the models (e.g. by adding local heterogeneity, smooth surfaces, nonlinearities) would lead to more accurate predictions.