Hole in One: an element reduction approach to modeling bone porosity in finite element analysis

Introduction

Finite element analysis (FEA) is a continuum mechanics-based technique originally conceived and used in the engineering design process to predict the behavior (i.e., response) of structures to prescribed loading conditions. This technique uses discretized representations of real-world structures, thereby enabling the design of these systems to be optimized mathematically with minimum physical prototyping and testing (Dumont, Grosse & Slater, 2009; Zienkiewicz & Taylor, 2000). With advances in computer software packages that allow a seamless connection of FEA to CAD and image data based modeling, the simulation method has also been applied to functional morphological research in organismal biology, including extinct organisms (as reviewed in Ross, 2005; Rayfield, 2007; Tseng & Wang, 2010; Bright, 2014). FEA of feeding mechanics of living and extinct vertebrates have been used in comparative functional morphology for more than a decade (Rayfield, 2005; Alexander, 2006; Barrett & Rayfield, 2006; McHenry et al., 2006; Thomassen et al., 2007), and the method also has been applied in studies in other organismal systems such as insect flight and mechanoreception (Combes & Daniel, 2003; Dechant et al., 2006; Wootton, 2003), and plant biomechanics (Fourcaud & Lac, 2003; Niklas, 1999).

For the last decade or so, the boundaries of FEA have been pushed towards more accurate modeling of bone structures to better understand skeletal form and function (Rayfield, 2007; Bourke et al., 2008; Wroe et al., 2008; Strait et al., 2010). Still, porous structures like trabecular bone and other complex biological geometries remain problematic in FE modeling given their internal complexity, and the conversion from 2D to 3D of intricate structures that frequently generate errors in elemental overlaps and highly skewed elemental shapes in small anatomical regions. Based on our experience working with bone meshes, biological structures with a high amount of trabecular bone or porous components have higher chances of meshing errors in the FE solid meshing process (but see Fagan et al., 2007 for an alternative, albeit computationally more intensive, and resolution limiting, voxel-based modeling approach). When modeling this type of porous structure, it is common to avoid the complexity of creating a detailed trabecular network by modeling entire models as homogeneous cortical bone and ignoring trabecular geometry, and/or changing the material properties in different element groups within a model to represent cortical versus trabecular bones (Strait et al., 2005; Strait et al., 2009; Wroe, 2008; Attard et al., 2011; Chamoli & Wroe, 2011). This general simplification approach is used in most comparative studies using FEA that incorporate trabecular morphology, even though it has been demonstrated that trabecular structures can play a very important role in the performance of a mesh when using FEA (Parr et al., 2013).

Our objective in this study is to test an alternative, mechanical approach to trabecular bone modeling as a viable solution in addition to mathematical approaches (i.e., changing the material properties of solid models). Potential solutions to accommodate trabecular morphology in finite element modeling that can bypass time-consuming and scan resolution-dependent micro-modeling of trabecular structures are desired. We aim to test the hypothesis that percentage porosity adjustments in solid finite element meshes will generate simulation results comparable or closer to those using actual trabecular morphology, compared to solid models using only modified material property parameter values to simulate trabecular bone behavior.

Materials and Methods

We used five species samples to test a finite element reduction approach to trabecular bone modeling relative to actual trabecular structural models. Each species-specific test sample is represented by three types of experimental cylindrical models: one control cylinder (“CC”); one physically modified cylinder (“PC”); and one material-modified cylinder (“MC”). Definitions of each cylinder model are given below.

Script analyses: does random element reduction deliver consistent results?

Prior to comparing PC models to the CC group or MC group, we tested an additional set of 5 models to ascertain the internal consistency of the script (whether random element deletion delivers consistent results).

We applied the same script, set at 16.5% volume deletion (we chose 16.5% deletion as a mid-range value through our tested range), to five otherwise identical solid cylinder models. The remaining parameter values, such as material properties (Young’s Modulus: 20 GPa and Poisson’s Ratio: 0.3), the amount of force applied (1,000 N), nodes retrained (four nodes, at the end of a cross-section, at the bottom of the cylinder), and the area of application all remained identical (see Model Simulation Parameters). All the points sampled were identical through all of the five cylinders (Fig. 1).

If large differences in magnitude of the stress values are present in script-generated models across different replicates, the script would not represent a true randomized approach to element reduction. If the effects of the script are random, the variability in the results for all 5 additional models should be within comparable ranges of variation. Some variability is expected because the script is based on a random pattern. As a consequence, some arbitrary associations that affect stress values may occur. Overall, our assumption is that replication of porosity in trabecular structures by random reduction of solid element would result in replication of overall trabecular mechanical behavior.

Results

Our results show that physically modified cylinder replicates, assigned the same specific settings, have largely uniform outputs (Fig. 2, Tables S2–S3). There was only a small problematic region, located at the bottom (points 8 to 10) of one of the transects in cylinder IV (Fig. 2A, Table S2). Because there are no differences between the cylinders beside the random arrangements that the script may have produced, the higher stress values on the nodes correspond to a more localized deletion at the sampled area. The higher deletion around that area would affect how the applied force is transmitted and distributed in that location, and thereby extend influence to contiguous areas (as subsequent points show higher stress values). This inconsistency effect is diluted by average stress values across analogous locations of the four sampling transects on each cylinder (Fig. 2B, Table S3). Therefore, all subsequent simulated stress values are reported as mean values.

There is a better overall performance of the physically modified cylinders in comparison with materially modified cylinders when referring to the control cylinders. In the experimental groups for 26.1% and 46.6% porosity (Figs. 3C and 3E, Tables S6–S8), we see a consistent performance of the PC. We can see a slightly more accurate overall trend in physically modified cylinders (it underestimates in certain regions and it is not able to imitate the peaks of CC, but replicates the general trend). The particular, the bottom section (sampled nodes 8 to 10) of the PC cylinders has a more accurate performance (relative to the control) than the materially modified cylinders. MCs in both 26.1% and 46.6% models exhibit a linear transect trend with relatively low stress changes.

On the other hand, in the experimental group of 16.5% porosity (Fig. 3B, Table S5), PC seemed to be unable to correctly replicate both trend and stress values of the control group. For the experimental group of 7.8% porosity (Fig. 3A, Table S4), PC and MC seem to perform equally well in most of the sampling points (same stress values or off by less than 10 MPa). Except at the beginning and the end (where higher variability may be present, close to the area of force application and nodal restraints), there are minimal differences in stress values across all modeling approaches tested at the lowest porosity level.

In the experimental group of 35.8% porosity (Fig. 3D, Table S7), the differences in stress values seem to be consistent with what we observe in groups with 26.1% and 46.6% porosity (Figs. 3C and 3E, Tables S6–S8). PC replicates the overall CC trend but its values are offset by 60 to 80 MPa, especially at the center core region. MC shows a less accurate trend, with a more linear pattern, and even less resemblance to the CC trend. As seen in all experimental groups (Figs. 3A–3E, Tables S4–S8) the combined PC+MC approach presents the same stress values as the PC group results. The differences are indistinguishable between PC and PC+MC results. The overall stress trends relative to porosity changed are plotted separately for the control versus element-reduction models in Fig. 4.

Strain energy values, used here as an overall measure of the stiffness or degree of deformation experienced by the different models, show that when porosity is low, modifying material properties of solid models provides the closest approximation to trabecular geometry models (Fig. 5; Table 1). However, as porosity increases to 16.5% and above, the element reduction models provide the closest approximation to trabecular geometry model strain energy values. Furthermore, a combined modeling approach of element reduction plus material modification generates close approximations of strain energy values at lower porosities, but is less accurate than either material or element reduction approaches alone, in models with higher porosities (Fig. 5).

Discussion

Element reduction is potentially a more accurate approach for modeling trabecular stress and strain energy than modification of regional material properties. We tested the hypothesis that, even if they are not 100% replicates of trabecular bone models, porous FE models can at least behave in a comparable way, and provide a closer approximation of mechanical behavior than only modifying overall material property parameters of solid models. Our results indicate that an element reduction approach to modeling bone porosity produced stress magnitudes that are generally closer to values generated from models containing actual trabecular bone geometry, compared to only modifying material properties to simulate bone porosity (Figs. 3 and 4). Furthermore, element reduction produces models with strain energy values that are comparable to those estimated by trabecular bone and material property modified models at low porosity values, and values that best approximate trabecular bone model outcomes at higher porosity values out of all the modeling approaches tested (Fig. 5).

Overall, the peaks and valleys of stress in the trabecular geometry models are not well-replicated by alternative modeling approaches. Stress peaks in the transect plot for the control group might be explained by how close the sampled node was to a physical hole or opening on the model surface (in other words, adjacent to an internal porous network) (Fig. 4A). The nodal values may be influenced by elevated stress values associated with such porosity. Thus, creating a cover layer of plate elements, then sampling from that surface, could be one modeling solution to account for the source of that possible noise. This could be considered in further studies, but our goal for this first study was to compare relative performances between the mechanical approach and the mathematical approach (PC vs MC); rather than specifically creating a protocol to mimic actual bone. As apparent from the results, the overall stress magnitude changes, but not the precise stress peaks, are replicated using the element-reduction approach proposed in this study (Fig. 4B).

It is remarkable that even without a cover of cortical bone (or a thick layer that might homogenize the values at the nodal transect regions) the mechanical modeling approach still has a certain consistency (results are similar in all four experimental groups for PC+MC models). Based on our results, the ability of PC models to approximate stress values in the control group models is best in moderate density models. As shown in Fig. 3D, the overall curvature of the stress sampling transect in the control model are mimicked by PC, whereas materially modified cylinder trends show a low-sensitivity trajectory, indicating that the overall performance of materially modified cylinders is less accurate than observed for data in the PC group.

We note that the element reduction script generated models with holes in a random pattern, whereas the actual species trabecular geometries contain holes surrounding a network of bony struts. As a consequence, PC models are more homogeneous in how they distribute forces. In other words, when compared to the CC group, the PC models perform as a more rigid material. This is probably related to their lack of internal heterogeneity in arrangements or concentration of large pores/bony struts that may not be represented by the mechanical modeling approach. This is another key factor to consider in future research into improving accuracy of trabecular bone modeling in FE simulations.

It is also quite clear that material properties modified cylinders behave as an even more rigid material than the other two groups. The von Mises stress values, which reflect the likeliness of a certain structure to fail, are apparently lower in MC. This stiffness, or lack of it, may be related to the internal network influence on the overall performance (Parr et al., 2013). Our results point to a plausible explanation to why so many published FE analyses report high stiffness values compared to those known for actual bone: most modeling approaches use a “solid approximation” of trabecular bone, but differences in material distribution, which clearly influences stress estimates as per results in this study, are rarely taken into account.

Bone tissue can behave as a homogeneous material on a microscale (Müller, 2009) with both individual trabeculae and compact bone having similar material properties (Rho, Ashman & Turner, 1993). Therefore, changing material properties to differentiate compact versus trabecular bone may not adequately replicate bone behavior in FE simulations. The adjustment of bone porosity based on the internal density of corresponding trabecular geometry models did better replicating the stress and strain energy values of the control group than MC models (Figs. 3 and 5). Accordingly, overall, the change in the material properties is a less effective way to approximate model mechanical behavior than physically reducing the element density of solid mesh models via the randomization approach (Fig. 6).

In addition, models with both physically introduced porosity and material property changes combined behaved similarly to the models with only introduced porosity, suggesting the dominant role of element reduction in dictating surface mechanical behavior of the cylinder models in our study. This finding could be even more important when considering animals with different amounts of trabecular volume. When modeling low trabecular volumes, such as those in the bones of birds or pterosaurs, the considerations should differ from modeling animals with higher bone densities such as some mammals. Consequently, comparative analysis or paleontological reconstructions should consider the nature of the structure that is being analyzed. An adequate adjustment relative to cortical and trabecular bone ratios when modeling could produce more accurate models that reflect mechanical behavior of trabecular geometries.

All things considered together, we summarize an important take home message for modeling complex geometries of skeletal elements in general: without element reduction or specific modeling of trabecular morphology, comparative FEA studies that include taxa with a wide range of cortical-to-trabecular bone ratios (e.g., relatively low cortical ratio in vertebrates such as birds and extinct taxa such as pterosaurs) will exhibit larger deviations in both stress and strain energy values than those models that represent morphologies with high cortical ratios. This is an additional source of uncertainty in the modeling protocol that could strongly influence simulation outcomes and thereby, functional morphological inferences made from those simulation results. If close approximation of both stress and strain energy values in finite element simulations involving trabecular structures are of interest, findings from this study demonstrate that an element reduction approach would be preferred over existing material property modification protocols.

Conclusions

We demonstrated that an element reduction approach to modeling trabecular structure could more closely simulate behavior of trabecular geometry compared to changing material properties in solid models. Additionally, when element reduction and material property modification methods are combined, the effects of element reduction (i.e., generation of porosity) far outweighs the relative effects of material property modification. We suggest that, unless the complex geometry of trabecular bone is precisely accounted for during the model building process, researchers should first consider modeling the porosity of the material instead of changing material properties. This recommendation is supported by our findings that indicate physical internal porosity generation better approximates mechanical performance (stress and strain energy values) of trabecular structures compared to material property changes. Therefore, we recommend considering bone porosity in such a physical manner in biomechanical modeling of complex trabecular bone geometries in comparative functional morphological studies, as a fast and effective way to approximate trabecular geometry, and to alleviate potential biases in finite element modeling protocol towards taxa that exhibit high trabecular bone ratios in their morphology.

Supplemental Information