Finite element modeling of multiple density materials of bone specimens for biomechanical behavior evaluation

The finite elements method allied with the computerized axial tomography (CT) is a mathematical modeling technique that allows constructing computational models for bone specimens from CT data. The objective of this work was to compare the experimental biomechanical behavior by three-point bending tests of porcine femur specimens with different types of computational models generated through the finite elements’ method and a multiple density materials assignation scheme. Using five femur specimens, 25 scenarios were created with differing quantities of materials. This latter was applied to computational models and in bone specimens subjected to failure. Among the three main highlights found, first, the results evidenced high precision in predicting experimental reaction force versus displacement in the models with larger number of assigned materials, with maximal results being an R2 of 0.99 and a minimum root-mean-square error of 3.29%. Secondly, measured and computed elastic stiffness values follow same trend with regard to specimen mass, and the latter underestimates stiffness values a 6% in average. Third and final highlight, this model can precisely and non-invasively assess bone tissue mechanical resistance based on subject-specific CT data, particularly if specimen deformation values at fracture are considered as part of the assessment procedure.


Introduction
Bone tissue is formed by tissue components with varying mechanical properties. This results in complex bone geometry as cortical and trabecular tissues have different mechanical traits and are distributed heterogeneously. As such, the creation of computational models that facilitate accurate assessments of osseous behavior first depends on developing a non-invasive technique able to analyze the geometry and the spatial distribution of tissues, as well as respective mechanical properties.
Currently available non-invasive clinical methods for estimating osseous resistance in a specific bone include bone densitometry (i.e., dual-energy X-ray absorptiometry) and peripheral quantitative computed tomography. However, these techniques are limited to information on regional and geometric density and do not consider factors such as three-dimensional structure, materials distribution within a structure, mechanical properties, or load condition [1].
The finite elements method (FEM) based on computerized axial tomography (CT) involves applying mathematical techniques to create a bone tissue model that incorporates information related to three-dimensional architecture and the distribution of osseous density, as provided via a CT scan of the bone under study [2,3]. Greater precision predicting structural behavior occurs as this mesh becomes more refined, as well as by increasing the quantity of materials in the model.
Several bone models using the FEM already exist [4][5][6]. These models could potentially serve to non-invasively evaluate a bone through the application of distinct mechanical force tests. Furthermore, these models can reveal compression, tension, fatigue, and fracture zones, as well as the respective force limits needed. Importantly, these models can be subject specific; in other words, they can be characterized for each specimen under investigation [7,8]. Complementary to what was previously published, predictive models of the highest precision are required, which would be capable of non-invasively and individually assess mechanical properties of bone in patients with pathological bone structure.
The use of a multiple materials FEM scheme [6,11] renders an enhancement of the FEM technique as it can realistically include the distribution of the mechanical properties of the subject-specific bone specimen. When loads and boundary conditions are well determined spatially and do not change significantly over time, deformations are bounded and materials properties are homogeneous; then, the FEM scheme in general will render results, from which output practical decisions could be made. In recent years, there has been a common tendency to include non-linear behavior of materials as part of the FEM scheme. This in order to achieve a higher fidelity of the expected mechanical behavior of the material modelled system. On one hand, non-linear FEM schemes are computationally costly [9,10] in terms of implementation and of the CPU time required to solve the simulation, as commonly an iteration procedure is needed to achieve convergence. While on the other hands, it means to include additional unknown parameters (depending upon the constitutive model used) or an adjustment procedure to represent the non-linear behavior of the involved materials. The latter is otherwise difficult to find or fit even from phenomenological data [10]. In despite of this fact, whenever a linear model is used, overestimation of the onset of damage stress is a known resulting feature of this simplification. However, the clinical praxis may overcome this overestimation by considering the experience of the specialist, who then can from the linear model prediction, correct the diagnosis using, for example, the data associated with the deformation values occurring at the onset of fracture, plus of course his or her expert knowledge. The tradeoff between this a posteriori correction and the additional costs of a more sophisticated non-linear model implementation is yet a matter of assessment and discussion. However, linear model implementation would be more feasible to transfer initially into the clinical realm in the early stages of a computational-assisted diagnosis initiative.
Regarding clinical application, the subject-specific prediction of mechanical bone properties is a highly valuable tool as this technique allows to precisely replicate morphology and structure, as shown in recent studies (Ramezanzadehkoldeh and Skallerud [12], Bahia [13] and Väänänen [14]). Particularly, if one considers this FEM analysis to enter the clinical daily praxis, where more concern is centered toward estimating critical loads at which fracture could occur rather than precisely knowing the load level at which transition to damage occurs. For ex vivo mechanical force tests in small animals, the femur is commonly used because of its size, length-to-width ratio, and consistent cross-sectional shape along the entire length [15]. Other studies have shown that the three-point bending test is a relevant biomechanical test that can provide information about the structural and material properties of the bone [15][16][17].
The objective of this study was to carry out subjectspecific simulations for the effects of three-point bend tests in porcine femur models. The models were developed using the finite element method and individual CT scan data. Distinct quantities of materials in the models were also analyzed. In parallel, the biomechanical properties of osseous tissue were characterized for the assessed animal model, and a comparison was made between computersimulated fractures and real-life ex vivo biomechanical testing.

Selection of bone specimens
Five femur specimens from 3-month-old Landrace pigs were used for analyses. Macroscopic examination and X-rays were utilized to rule out any bone abnormalities in the specimens. This study was approved by the Animal Welfare Ethics Committee, School of Medicine, Pontificia Universidad Católica de Chile, and it complies with the requirements regarding the allowed number of tested specimens. Research Article the femurs were supported by resin molds with the exact shape of the bone to restrict any movement. For each specimen, vertical force was applied to the midpoint of the diaphysis at 0.5 mm/min, with registries every second until bone failure (Fig. 1). The exact location of bone failure and the failure pattern were recorded for later comparison with the finite elements model. These coincide with forces suffered in trauma-induced fractures.

Boundary conditions
Molds from the epiphyses of each femur (distal and proximal sides) were resin casted and used to rigidly hold the femur horizontally in place inside the testing machine. This prevented sliding of the bone epiphyses during the test as the middle point of the diaphysis was pressed down by the punch rod. Therefore, boundary conditions implemented in the model consisted in fixing these femur lateral regions in two nodes at each side. These nodes were arbitrarily selected within the neighborhood of the distal and proximal epiphyses, restricting their motion along the x-, y-and z-directions, nonetheless allowing for numerical convergence of the numerical solution to result. The midpoint of the diaphysis was calculated as the semi-length between each pair of nodes.

Numerical modeling
Numerical models were created through the following steps sequence (Fig. 2a), which ranged from radiological assessments to model creation through the FEM: 1. Acquiring CT Information: A helical CT scan was conducted using a GE LightSpeed Ultra CT Scanner (GE Medical Systems, Chicago, IL, the USA) with a 0.625mm slice separation. The information was digitized   [6], resulting in the assignment of specific mechanical properties to the models. 5. Gray Calibration Values: This process was used to establish density values from CT grayscale volume values, as already reported by several other authors. [2,5,20]. A linear regression curve was then constructed that correlated each scalar value with a corresponding tissue density. This process was confirmed by calculating the total bone mass and comparing it against the mass of the FEM mesh. The equation relating density and CT numbers is linear and can be written according to Taddei et al. [6], as follows: where n is the uniform density assigned to the n element of the mesh, CT n is the uniform CT number, and α and β are the calibration coefficients provided by the user. An iterative process was used to calculate the minimum and maximum bone densities from the CT scan of each specimen (see Fig. 2b). 6. Relations of Density and Elasticity: The mechanical properties of each defined material were estimated based on density values in the form of power laws, as proposed in other studies [2,3,21]. For this, techniques proposed in the literature were applied [22][23][24].
A distinction was made between the properties of the trabecular and cortical bones to create a robust model. The correlation between apparent density (grams/ cm 3 ) to cortical and trabecular tissue elastic modulus (MPa) is shown in Eqs. 2 and 3, respectively. The latter study did not provide a correlation linking Poisson ratio with apparent density; however, an average value for Poisson ratio equal to 0.3 was then adopted. Quantity of Materials in the Generated Models: Five models with increasing quantities of materials were created for each specimen. The reference model (r) was created with two materials, while the following models (a thru d) were created with 20, 40, 75, and 140 material types, respectively, where each material was defined as having a specific tissue density. The material distribution is illustrated in Figure 4. The difference in the number of materials to be incorporated into each model is carried out by defining a given spacing between consecutive elastic modulus values. Results from the reference model were compared with regard to results obtained from the models that incorporate variable material properties. Table 1 illustrates the number of materials of each model assigned at Specimen 5, both at the cortical (c) and trabecular regions (t), and the value range of the mass density and elastic modulus. 7. Finite Elements Analysis: This analysis was performed using the ANSYS v11.0 software (ANSYS Inc.), using tetrahedrons with ten nodes for each element [2,20,21,25]. Approximately 20,000 elements considering approximately 35 different material groups were chosen for the calibration procedure. Three to four iterations were conducted to reach a total mass value closed enough to the experimentally measured value for each bone. In all cases, values with less than 1% of difference were achieved. Once the iterative process was finished, the density value that made it possible to calculate a total bone mass closer to the obtained in the experimental results was chosen. Regarding the mesh sensitivity analysis, significant differences across the parameters under study (longitudinal component of elements stress and longitudinal component of nodal stress) were not observed. The measured differences in no case exceed 6%. The mesh was chosen for sensitivity analysis then consisted of 34,000 nodes and 26,000 elements. These amounts varied depending on the geometric characteristics of each bone´s anatomical features.

Mechanical study of computational models
Force values applied to real specimens were used as boundary conditions in the generated computational models. The latter was run under a load equivalent to 1500 N in the femoral diaphysis (midpoint) as schematically indicated in Fig. 5. Moreover, Table 2 indicates the exact node location and value of the applied vertical deformation (i.e., boundary condition) for each of the five femur models, as well as the corresponding value of the application point.
The similarities between reaction and tension forces were compared through linear regression of the force values measured from each of the five specimen models as presented in the following graph (Fig. 6). From this correlation, one can obtain the coefficients of linear regression, the intercept, and curves by using the method proposed in Helgason et al. [2] and prior studies [18,25]. These results show how the linear model can be used to predict the onset of damage in the bone. Fracture locations in the models were inferred through the von Mises yield criterion according to Keyak et al. [26].  Specimen fracture was evaluated considering the surface area containing elements whose stress values exceeded the strength of the bone tissue. Thus, it was necessary to choose a criterion to assess the computed stresses. According to the literature, there are several fracture criteria that can be considered for bone material: the von Mises stress, maximum principal stress and maximum principal strain. However, none of these criteria has been used to evaluate juvenile porcine femur fracture. The approach adopted in this study was the von Mises stress criteria, which has proven to be useful for studying other long bones (e.g., tibiae) that are anatomically similar to the femur (Keyak et al. 1994). On the other hands, this criterion yields values which fall in between other relationships such as the one of Schileo et al. 2008 (maximum principal strain) and Keyak et al. 1994 (maximum principal stress), so its use would avoid over-or underestimation. In summary, to evaluate bone fracture at local levels, stresses were calculated according to the von Mises criterion and their values compared to the estimate given by the relationship proposed by Keyak et al. (1994): where sigma_lim is the maximum stress (MPa) value the bone can withstand before fracture and corresponds to the mineral density tissue value (g/cm3). Mineral density is calculated from the apparent density according to the ratio as suggested by Schileo et al. (2007). Table 3 shows maximum stress values calculated from the average density of each specimen model.

Statistical analysis
We performed statistical analysis using Microsoft Excel (Microsoft 2017. Microsoft Excel version 15.39) and Stata 13 (StataCorp. 2013. Stata Statistical Software: Release 13). Data distribution assessed using the Shapiro-Wilk test. Descriptive statistics were used to report specimen data and results. Correlations between the experimental model and the finite elements model were determined using Wilcoxon signedrank test for paired data and Bland-Altman analysis. Linear regression analysis of the experimental model was also conducted, the results of which were compared with the FEM. Statistical significance was set at p value < 0.05.

Results
Basic specimen measurements are summarized in Table 4. The median length of specimens was 130 mm (128-135 mm). The median diaphysis diameter was 19 mm    (18-20 mm), and the median total mass was 94.78 g (91.9-106.8 g). The median reaction force at fracture, measured experimentally, that occurred in the specimens was 1005 N (968-1857 N). Maximum displacement ranged between 2.62 and 4 mm, with the median being 3 mm (Fig. 7). The sigmoidal shape of the five curves is due to several factors, namely an initial sliding and accommodation of the bone epiphysis within the molded fixture at the initiation of the bending test. This is followed by a decreasing increment behavior of the ultimate force at the end of the test due to the non-linear elastic response of the bone material. The middle regions show clearly a linear behavior as expected from elasticity theory at low elongations, with quasi-parallel slopes. Elastic stiffness obtained as the slope in the linear zone account to m 1 = 714 N/mm (Femur 5), m 2 = 600 N/mm (femurs 1 and 3), m 3 = 576 N/mm (femurs 2 and 4). Femur 5 has the largest mass (106.8 g) thus highest stiffness, while Femur 2 has the lowest mass (91.9 g) and lowest stiffness. It can be also observed that fractures occur just prior to peak forces in all specimens, suggesting that the material has not yielded before break and that it has rather followed a non-linear elastic behavior, particularly at the last stage of deformation.
On the other hands, from the computational models, the median maximum reaction force prior to fracturing was 1177 N (904-2087 N). When comparing maximum force, models a (20 materials) and b (40 materials) significantly differed from the experimental values obtained (p = 0.043 and p = 0.044, respectively); while models c (75 materials) and d (140 materials) were not significantly different than the specimen values (p = 0.078 and p = 0.225, respectively). Although reaction force versus displacement showed a tendency of linear behavior (Fig. 8), the median displacement was also 3 mm (2.60-4 mm). It can be noticed the similarity of the simulated reaction force magnitude to that from the experimental model; however, a stronger linear tendency is observable. The linearity of the latter five curves is expected as the implemented FEM model is based on a linear elastic model. In this case elastic stiffness obtained as the slope in the linear zone accounts to m 1 = 667 N/mm (Femur 5), m 2 = 633 N/mm (Femur 3), m 3 = 546 N/mm (femurs 1, 2 and 4). Here again, femur 5 has the highest stiffness value, while femur 2 has the and lowest value. However, it is worth noticing that predicted stiffness values slightly fall within a range of 5-7% lower values than real stiffness values.
Moreover, from Figs. 7 and 8, the % error associated with the underestimation of the specimen deformation at fracture is less than the % error associated with the overestimation of the force at fracture. In the case of Femur 5, the % error for the maximum force accounts to (2076-1848 N)/1848 N = 12.4% versus the % error for the deformation at fracture which comes to Finally, linear regression of the mechanical behavior of real specimens versus the computational models revealed that the model with the greatest quantity of materials (model d) provided more precise results. Specifically, an R 2 and R2adj ≥ 0.97 and a root-mean-square error range of 3.29-5.29% were observed ( Fig. 9 and Table 5). However, use of less materials also allows acceptable simulation precision of the experimental model.
Furthermore, visual inspection of the locations where damage accumulation has occurred in the five specimens (i.e., signaling fracture onset) can be obtained from Fig. 10. Comparison between the experimental specimens and the ones obtained from the computational models coincides significantly and closely resembles fractures in the bone specimens.
Finally, Bland-Altman analysis comparing real specimen deformation to model deformation showed low bias across all five specimens (Fig. 11). Bias did not change significantly secondary to increase in deformation averages, with a median bias of 0.67% (0.01 mm), with a range from 0.35-0.89%. This analysis illustrates a high degree of correlation between the experimental measurements and the finite elements model (d) regarding specimen deformation values.

Discussion
As the number of materials increased, the spacing between the values of continuous elastic modules decreased, and clusters with lower density values were  created. Due to this, a progressive decrease was observed in the maximum reaction force value. Figure 12 shows the results of the maximum reaction force value calculated computationally for each specimen, using an increasing number of materials. The sharp decline recorded in the first points of the curve is associated with a notable difference between the reference model (homogeneous 2 materials model) and the simplest model that incorporates spatially variable material properties (model-a incorporating 20 materials). On the other hands, the computations involving models having spatially variable material properties show percentage variations near to 15% when comparing the simplest (model-a) with the most detailed model (model-d) that incorporates ca. 140 materials. The more precise models tested approximated a final reaction force very close to that found in the experimental models. In other words, by incorporating a high number of material types, certain elements of the finite mesh were given isotropic properties (i.e., independent of spatial direction). Values were assigned depending on bone density where the element was located. Therefore, an anisotropic problem, where properties depend on spatial orientation, would be approximated as the sum of distinct, discrete and spatially minute isotropic problems distributed over the space. The three-point bend test was chosen because of its simplicity and capacity to provide information about the structural and material properties of the bone [15][16][17][18]. Various other bone models using the FEM have been reported as shown in Table 6. Huiskes and Chao [27] provide a good review on the first research studies conducted in this area. Other references to initial models combined with the use of CT scans are given by Huang et al. [28] and Seitz and Rüegsegger [29], among others. Current models can predict forces of stress with high precision through this combined technique (R 2 > 0.9, root-mean-square error < 10%) [21,29]. Worth noting, the method presented in our research achieves the highest precision, and the applied mesh morphing has been used with a high degree of predictive precision in recent clinical applications [19,30,31]. Nonetheless, once a mathematical model is obtained, according to Choi, the finite element method can be implemented to solve the model; however, it is difficult for numerical simulations to consider every mechanical condition arising in real applications [32].
Despite the reduced quantity of specimens used in the present study, the fact that various scenarios were assessed for each specimen (25 in total) increased the amount of obtained information. Furthermore, the uniformity of results supports the validity of the employed method; especially considering the use of the von Mises criterion for determining fracture location. After the original article describing the technique [26], many articles Research Article compared the use of strain base methods (principal strain) and stress-based criteria (principal stress, von Mises) [33,34]. Although Schileo et al. [5] proposed the use of a strain criteria after the analysis of three cadaveric femurs, it continues to be a lack of consensus on the preferred criteria [35]. Furthermore, our current study supports the use of stress-based von Mises with a high correlation between the location of the specimen fracture and the predicted computational model. Finally, from the recent study of Fung et al. [36], Schileo et al. [37] and others [38,39], subject-specific prediction of mechanical bone properties is a useful tool, particularly in clinical applications, allowing morphology and structure replication. The potential applications of non-invasive prediction include predicting fracture risk in pathological bone, managing rehabilitation in complex skeletal reconstructions, and evaluating bone consolidation and the incorporation of bone grafts, among other applications. Furthermore, the herewith presented model evidences the possibility of obtaining high precision with a manageable (i.e., not so high) materials quantity, thus providing the advantage of decreasing data processing times without sacrificing model predictive precision. Particularly, when specimen deformation measurements at fracture are considered, as the % error is less compared to the force estimation % error. These latter could mean advancements in applicability in clinical contexts, where increasing age of the population is leading to increased patients in risk of pathological bone fractures that could benefit from a predictive model to guide prophylactic bone fixation surgery.

Conclusions
The present study compared the experimental biomechanical behavior of porcine femur specimens with distinct types of computational models generated through the finite elements' method based on computed axial tomography under a multiple density materials scheme.
One advantage of this method is the ability to use an increasing number of materials, (i.e., from 2 to ca. 140), thus, to explore the compromise between the overestimation of the ultimate stress from a single isotropic linear material model and the more accurate yet intricate and computationally costly non-linear elastic model.
Main highlight result obtained, in first place, shows a high precision prediction for measured reaction force versus displacement between the experimental three-point bending tests and simulations (i.e., R 2 of 0.99 and a minimum root-mean-square error of 3.29% in the model with greatest definition). Second highlight result indicates that elastic stiffness, both measured and computed, follows the trend, higher the mass-higher the stiffness; however, stiffness predictions are in average a 6% below real stiffness values. As third highlight, the model can precisely and non-invasively assess bone tissue resistance based on subject-specific CT data, particularly if specimen deformation values at fracture are taking into consideration.
The present study represents a contribution to the development of biomechanical bone behavior assessment based on computed axial tomography and the finite elements method.
Acknowledgements Work presented in this manuscript was supported by the Pontificia Universidad Católica Medicine School and   [42] 2013 Quantitative CT 0.81 Luisier et al. [43] 2014 Quantitative CT 0.84 Dall'Ara et al. [44] 2012 Quantitative CT 0.87 Koivumäki et al. [45] 2012 Standard CT 0.87 Anez-Bustillos et al. [46] 2013 Quantitative CT 0.89 Mirzaei et al. [47] 2015 Quantitative CT 0.89 Taddei et al. [25] 2006 Standard CT 0.91 Helgason et al. [2] 2008 Standard CT 0.92 Hambli and Allaoui [48] 2013 the Mechanical and Metallurgical Engineering Department from the same institution. The authors would also like to acknowledge FON-DEF Project D06I1026 for its collaboration in the accomplishment of this research work. We would also like to acknowledge University of Notre Dame, Aerospace and Mechanical Engineering Department, for granting a Melchor Visiting Scholarship to the corresponding author, time during which part of this manuscript was prepared.
Funding None.

Declarations
Conflict of interest Each author certifies that he or she has no commercial associations (e.g., consultancies, stock ownership, equity interest, patent/licensing arrangements, etc.) that might pose a conflict of interest in connection with the submitted article.
Ethical approval This study was approved by the Animal Welfare Ethics Committee, School of Medicine, Pontificia Universidad Católica de Chile.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.