Mechanics of osteoporotic trabecular bone

– Osteoporosis is a disease characterized by a loss of bone density and an altered bone architecture. These modiﬁcations lead to an increased risk factor for bone fracture, particularly of the femoral neck. This disease can be explained by a disorder in the bone remodeling process which is triggered by the apparition of micro-cracks within the bone. According to Frost’s theory [1], these micro-cracks appear for a speciﬁc local strain threshold. Thus, the knowledge of the microarchitecture and quality of trabecular bone is essential to determine this local strain threshold. This paper studied the mechanical trabecular bone behavior of 43 patients diagnosed as osteoporotic whose femoral heads were replaced by hip prosthesis. From each patient, a cylinder-shaped of trabecular bone samples was cored. Each sample was scanned by X-ray micro-tomography before a compression test in order to reconstruct a reliable Finite-Element (FE) model of the bone architecture in Abaqus. The force-displacement curves were recorded for all the samples and calibrated by the experimental responses. The force-displacement numerical curves were adjusted to the experimental ones, by modifying the tissue microscopic mechanical behavior. This process leads to the determination of the local strain threshold responsible for triggering the bone remodeling process.


Introduction
Osteoporosis is a disease of bone metabolism which decreases the bone density and alters the architecture of the trabecular bone. With an estimation of 6.3 million femoral neck fractures in 2050 according to the WHO, osteoporosis is a major health issue. The determination of the mechanisms involved in osteoporosis is an essential point to better predict the fracture risk. In this context, a biomechanical approach is strongly relevant to determine the link between the mechanical properties and the microarchitecture of the trabecular bone.
Knowledge of the bone remodeling process is of importance to gain understanding about the origin of osteoporosis. This process is equilibrium between destruction of the "old bone" by osteoclasts and formation of a "new one" by osteoblasts. One of the theories is that this mechanism is triggered by micro-cracks appearing in bone. Frost et al. [1] have given a range of the local strain threshold which varies, in the physiological range, from 400 to 3000 μ def (i.e. 10 −6 def). Patients submitted to 400 μ def or less (i.e. under the minimum in vivo strain) are considered "under-stimulated" and patients submitted to over 3000 μ def (i.e. over the maximum in vivo strain) are considered "over-stimulated".
The determination of these thresholds requires a thorough knowledge of the local mechanical behavior of trabeculae. This determination was the subject of many studies in the literature. For this purpose, several experimental techniques are commonly used to determine the mechanical properties of both cortical [2] and trabecular bone such as: (i) micro/nanoindentation [3][4][5][6][7][8] (ii) three-point bending tests [9,10] (iii) single trabeculae tension [11,12] (iv) ultrasound measurements [13][14][15]. Many papers made a review on the different values of mechanical properties [16]. The microscopic Young's modulus (i.e. tissular) for trabecular bone given in the literature ranges from 8 to 20 GPa. The measured value depends on the region where measurement is performed, on the wet or dry sample state [17] or on the experimental technique used for the measurement. By considering only the values obtained for human patient trabecular bone in the top region of the femur, the microscopic Young's modulus value ranges from 11.4 to 18.14 GPa [18]. This value can even reach 21.8 GPa in a more recent study [19].
However, these experimental methods are more suitable to measure the properties of cortical bone where the Haversian micro-architecture is revealed. The numerical techniques give another way to measure the local mechanical parameters of the trabecular bone while having information on its architecture (bone density, anisotropy. . . ). Several studies [20][21][22][23][24][25][26] simulate the mechanical trabecular bone behavior by taking into account the geometry. In all these cases, the treatment of images [27,28] provided by micro-computed tomography (μ CT scan) is one of the key point to ensure a reliable reconstruction of the trabecular bone sample.
The Finite Element method, used as an identification tool, leads to several drawbacks: (i) underestimation of the microscopic Young's Modulus E [26,27] (ii) high values dispersion and CPU resources [21,29] (iii) lack of representativeness of the micro-architecture [21,[30][31][32][33] (iv) use of nanoindentation measurement to build the finite element model [34]. To remove these drawbacks, a numerical method is proposed in the present study which reproduces the architecture of a fresh human bone sample while using reasonable CPU resources. In the present approach, the local mechanical behavior of the trabecular bone was determined by a manual identification procedure using an experimental compression test as a reference. The reliability of this method was analysed by comparison of the microscopic Young's modulus to literature values. This approach made it possible to numerically assess the local strain threshold involved in Frost's theory.
The micro compression tests performed on human trabecular bones, the microarchitecture analysis and the identification procedure are first described in Section 2. The determination method of the strain threshold used in Frost's theory is also presented in this section. Section 3 details the results obtained at macroscopic and microscopic scale. All mechanical parameters are compared to bone porosity and patient's age. The strain threshold responsible for triggering bone remodeling is given. In Section 4, the reliability of the numerical method and the relevance of the local strain threshold involved in Frost's theory are discussed. Conclusions are provided in Section 5.

Experimental protocol
Trabecular bone samples were retrieved from the femoral head of 43 patients aged 80 ± 11.6. All patients had a femoral neck fracture with a hip joint replaced by a prosthetic implant. Cylinder-shaped samples were removed from the femoral head (see STEP 1 Fig. 1) thanks to a trephine in a water bath. The two extreme surfaces were polished (grit #1200) to obtain a fresh standard sample with a mean radius of 3.5 mm and a mean length of 8.7 mm.
The bone samples were placed in Eppendorf tubes, filled with PBS Solutions to prevent desiccation and analysed with a Skyscan 1172 micro-CT (Skyscan, Aartselaar, Belgium) operating in the cone-beam method (see STEP 2 Fig. 1). The tubes were fixed on brass tubes with plasticine, analysed at a magnification of ×35 (pixel size = 8 microns in all three spatial directions) and the images were stored in the .tif format. Projection images on a CCD camera were obtained at 59 kV and 100 μA with a 0.5 mm aluminum filter. A rotation of 0.25 • was used between each image acquisition, providing a series of 720 projection images. From this series of projection images, a stack of 2D sections was reconstructed for each bone sample and stored in the .bmp format, with indexed grey levels ranging from 0 (black) to 255 (white). This imaging step was required to build a FE mesh of the actual architecture of each sample. For each sample, a compression test (see STEP 3 Fig. 1) was performed along the axis of the cylinder followed by an unloading step. The mechanical tests were performed on a specific micro device [35] at a constant displacement velocity of 1 μm.s −1 until reaching the maximum force of the first plateau. The force-displacement curves were recorded and used as references for the numerical identification process.
Finally a last micro-CT acquisition of the sample was performed after the compression test with the same acquisition parameters than in STEP 2 (see STEP 4 Fig. 1). The real residual plastic strain after the compression test is determined using the length of the sample measured on the reconstructed CT scans image with a pixel precision. This residual plastic strain allowed a correction of the experimental force-displacement curve by a coefficient equal to the real residual plastic strain divided by the residual plastic strain measured on the curve. This procedure allowed to access the real displacement of the sample in absence of extensometer.

Numerical method
To identify the mechanical behavior of each bone sample, a special treatment of CT images was required to build an accurate FE mesh. The choice of a tetrahedral mesh was adopted [36]. In a first step, the microCT images were converted to an exploitable form (8 bit images in axial slices i.e. 256 grey-level images) and were smoothed using a 3D median filter to reduce noise coming from the experimental acquisition (see Fig. 2a). 256 grey-level images were then binarised using a threshold calculated from the minimum of the %pixel/grey-level histogram of all slices. All details of this procedure can be found in [37] (see Fig. 2b).
All voxels with grey level higher than the threshold were considered as bone. Bone density (BV/TV: Bone Volume/Tissue Volume) was determined using CTan software from Skyscan. This bone density corresponds to the ratio of white voxels (i.e. 3D pixels representing the bone volume) on the total number of voxels (black and white) included in the cylinder-shaped ROI (Region Of Interest).
A 3D surface contour was built thanks to the software Avizo fire r based on the shape of the white voxels of the binarised images (see Fig. 2c). A simplified version of the surface-based contour (see Fig. 2d) enabling to create a tetrahedral FE mesh was then built (see Fig. 2e). The final FE mesh faithfully reproduced the original micro architecture of the considered sample. For each patient, the FE model of the compression test was reliably designed to the experimental one.
FE mesh obtained through this image treatment was composed by tetrahedral 4-noded elements with approximately 250 000 nodes and 800 000 elements for each patient. All the islands (i.e. unconnected regions of the mesh) were removed prior to the mesh generation to be considered as a single part. The FE model was then designed using the software Abaqus 6.9-2. In a first approach, the bone mechanical behavior of the samples was considered linear isotropic elastic. The macroscopic strain was applied through a prescribed displacement at the nodes located on the upper surface of the sample while the displacements at the nodes on the bottom surface were fixed in the compressive loading direction. Once the FE model was implemented, the manual identification process can be achieved.

Determination of microscopic mechanical parameters
In a first approach, the identification process only involved the elastic part of the compression. As there was only one parameter to identify (i.e. the Young modulus as the Poisson's ratio was assumed to 0.3) a simple identification procedure was used. The linear elastic part of the numerical force-displacement curve was fitted to the corrected experimental one to identify the local Young's modulus of each patient. The difference between the two curves was minimized using the least squares approach. Due to the weak number of parameters, i.e. only the Young's modulus, this step was done manually.
In a second step, the local strain threshold that initiates the apparition of micro-cracks triggering the remodeling process was determined. The micro-cracks are assumed to appear at the bone yield strength which corresponds to the first deviation of the linear part of the force/displacement curve. A simulation of the compression test, using the numerical model and the local Young's modulus identified, was performed for each patient. The local strain state of each element of the mesh was computed through the equivalent deformation ε eq at the yield strength step.
The local strain threshold should correspond to the highest value of ε eq in the mesh. However, due to numerical boundary conditions, several local strain values were overestimated and had no physical meaning. To remove these overestimated values, the different values of ε eq obtained in all elements of the mesh were plotted for all patients in decreasing order. As an example, the first hundred maximum values of ε eq were plotted for one patient in Figure 3. This typical curve is composed of two parts: a linear part for weak values, followed by a strong increase for the highest values. These highest values came from the numerical boundary conditions and were not considered.

Macroscopic results
The macroscopic Young's modulus E app was measured from the experimental data (slope of the linear part of the stress-strain curve) for the 43 patients. This macroscopic Young's modulus ranged from 427 MPa to 3621 MPa with an average value of 1538 MPa (i.e. 1538 ± 674 MPa). It is worth to note that this parameter plays a significant role in the identification method that requires a good degree of confident on its value. These results should be compared to 2730±1060 MPa measured byÖhman et al. [38], 3230 ± 936 MPa measured by Morgan et al. [39] and the results of Chevalier et al. [19] that ranged from 632 MPa to 2483 MPa. All these results were obtained in the case of human femoral neck in compression [38,39] and in the case of a numerical study using results from experimental tests [19].
In fact, as all biomechanics experiments, the results are strongly dependent of the patient and of the anatomic area tested. However, the E app obtained for the 43 patients of the present study are of the same order of magnitude as the results found in literature. The difference could be assigned to the direction of loading in regard of the microstructure organization or to the bone density of the samples. The BV/TV mean value was 21.8% for the 43 patients analysed in the present study as against 26.5% forÖhman et al. Figure 4 shows the relationship between the macroscopic yield stress parameter σ and the patient's age or the bone porosity. The first graph between the age of the patient and the macroscopic yield stress parameter (see Fig. 4a) showed no correlation that which was expected in the literature.
By contrast, a correlation is obtained between the macroscopic yield stress parameter and the BV/TV (bone volume/tissue volume) with a square correlation coefficient (R 2 ) close to 50%. As expected, this result was close to those obtained in [38] that compare the ultimate stress σ max to the BV/TV values and obtained a R 2 = 68%.
As the measurements of the macroscopic mechanical parameters were in an acceptable range of value compared to values found in literature, the reference used in the identification process is reliable.

Microscopic results from the identification process
Results obtained for the microscopic Young modulus (E loc ) of the 43 patients ranged from 5.25 GPa to 29 GPa with a mean value of 18.6 GPa (i.e. 18.8 ± 5.2 GPa) which is closer to the higher value of the values reported in the introduction section and to the results obtained using nanoindentation technique [3,19].
Based on the microscopic Young's modulus E loc the strain threshold ε loc was determined for each patient. For this purpose, the finite element simulation was conducted to reach the macroscopic yield strength. The results obtained ranged from 266 to 5400 μ def with a mean value of 1821 μ def (i.e. 1821 ± 976 μ def). These values are in the same range than the values expressed in the Frost's theory. It is worth to note that 2 patients are below the minimum in vivo strain and 2 others are above the maximum in vivo strain (see Fig. 5).
All these results were obtained by assuming that the microscopic Young's modulus is constant all over the microarchitecture. To investigate the effect of the local tissue heterogeneity on the value of ε loc , a random dispersion of the microscopic Young's modulus within the microarchitecture was tested. The microscopic Young's modulus E loc for the homogeneous model was split into five classes of values: E loc − 4 GPa, E loc − 2 GPa, E loc , E loc + 2 GPa, E loc + 4 GPa. This heterogeneous simulation was performed on three extreme samples in term of bone density. The results showed that the value of ε loc is weakly modified (less than 2% error between the two heterogeneous and the homogeneous models) by introducing this dispersion.

Discussion
In the present study, we assumed that bone tissue was homogeneous and isotropic all over the trabecular bone microarchitecture. The mechanical behavior of the bone tissue is modeled by only one parameter, the microscopic Young's modulus. However, introduction of heterogeneous local bone behavior didn't change strongly the results. Both experimental and numerical procedures were repeated with the same protocol for all the samples which ensured that all the results obtained afterwards can be compared to each other.
The mean values of the macroscopic Young's modulus E app and the maximum stress σ max measured experimentally for the 43 patients are in good agreement with literature data (cf. Fig. 5). The good range of mechanical parameters value obtained at macroscopic scale shows that the studied patients are representative. In the same way, the microscopic Young moduli obtained with the present identification procedure are close to literature values. It is worth to note that these microscopic Young's moduli from the literature were obtained by various experimental technics.
Local strain thresholds obtained from the identified microscopic numerical model are in the range of values (even considering four cases which are out of the physiological range of Frost et al.) described by Frost et al. [1]. As the local deformation thresholds ε loc were determined at the yield strength of the trabecular sample they correspond to the maximum of the local strain threshold. In fact, in the case of a bone physiological loading the yield strength which corresponds to the beginning of the macroscopic damage is never reached. Bone remodeling and the microcrack process associated involve fatigue loading. However, it is now well known that microarchitecture of bone is adapted or optimize to resist in in vivo loading. Then the maximum value ε loc obtained in the present model was in the same order of magnitude than the real strain threshold.

Conclusion
In the present study, three objectives were achieved: (i) developing a numerical FE model based on the microarchitecture of the trabecular bone to reproduce the experimental compression test performed on the real trabecular bone sample (ii) using an identification procedure to determine a reliable and personalized microscopic Young's modulus (iii) finding a local deformation threshold ε loc responsible for triggering the bone remodeling. All the mechanical parameters identified by this process were compared to the bone porosity and the age of patient. This comparison allowed on one hand to pinpoint which parameters play an important role in the osteoporosis and on the other hand to find relations between microscopic and macroscopic behavior.
The relationship analysis between the microarchitecture parameters (SMI, Tb.th. . . ), the patient features (Height, weight, disease. . . ) and the mechanical parameters is in progress. The aim of this future work will be to determine the main parameters governing the fracture risk prediction. The four patients whose local strain thresholds are not in the physiological range determined by Frost et al. should give new insight in regard to their biomedical past.