The development of a multi-axial failure criterion for trabecular skull bone has many clinical and biological implications. This failure criterion would allow for modeling of bone under daily loading scenarios that typically are multi-axial in nature. Some yield criteria have been developed to evaluate the failure of trabecular bone, but there is a little consensus among them. To help gain deeper understanding of multi-axial failure response of trabecular skull bone, we developed 30 microstructural finite element models of porous porcine skull bone and subjected them to multi-axial displacement loading simulations that spanned three-dimensional (3D) stress and strain space. High-resolution microcomputed tomography (microCT) scans of porcine trabecular bone were obtained and used to develop the meshes used for finite element simulations. In total, 376 unique multi-axial loading cases were simulated for each of the 30 microstructure models. Then, results from the total of 11,280 simulations (approximately 135,360 central processing unit-hours) were used to develop a mathematical expression, which describes the average three-dimensional yield surface in strain space. Our results indicate that the yield strain of porcine trabecular bone under multi-axial loading is nearly isotropic and despite a spread of yielding points between the 30 different microstructures, no significant relationship between the yield strain and bone volume fraction is observed. The proposed yield equation has simple format and it can be implemented into a macroscopic model for the prediction of failure of whole bones.

## Introduction

An understanding of the failure response in trabecular skull bone due to multi-axial loads is important for developing an in-depth characterization of its failure behavior. Bones in the body, including the skull, experience multidirectional loads on a daily basis [1,2]. In certain scenarios, extreme loads can lead to undesirable bone damage and fracture—resulting in injury and complications. Studying the multi-axial failure response of microscopic samples of bone when it is subjected to complex loads can give information necessary for the development of a failure criterion. This criterion could be incorporated into constitutive models at the macroscopic level that would enable simulations to predict the mechanical response and resulting fracture of the bone under extreme loads, such as ballistic impacts.

There has been extensive research on trabecular bone in order to understand its failure behavior. These studies analyzed bones such as the spinal vertebrae [3,4], radius [5], tibia [6–8], and femur [9–11] to get a general understanding of the fracture mechanics using both experimental and computational methods [12–14]. These studies typically analyzed the response under uniaxial tension or compression loading [7,9,15] and sometimes biaxial loading [6]. Although many in-depth studies have been conducted on uniaxial and biaxial analyses, triaxial/multi-axial studies have been limited [16,17]. In terms of failure, it has been found that, like most materials, trabecular bone is stronger in compression than in tension, which a von Mises criterion (von Mises stress or equivalent stress is a scalar stress value that predicts yielding of materials under multi-axial loading conditions using results from simple uniaxial tensile tests, so the yield envelop is symmetric in compression and tension) does not account for Refs. [6], [7], [9], and [18]. Mechanical tests show that for trabecular bone, on average, the ratio between tensile to compressive yield strains is about 0.6 and this strength asymmetry must be included for accurate modeling [13,15]. In addition, there is a lack of available data for tissue-level failure properties for trabecular bone [15]. Due to technical challenges of conducting experimental mechanical tests, substantial discrepancies remain in the published data that makes it hard to make trusted conclusions [7,13,15,19,20]. For example, many studies present conflicting results on the relationship between trabecular bone porosity and yielding properties. Some studies report that yield strains are independent of bone porosity [21–26], while others claim that the bone is stronger when it is more dense [27–31]. Developing a multi-axial yield criterion, excluding variability, requires repeated failure testing on a single sample, which is impossible to achieve with mechanical tests. As a result, experimental multi-axial yield data based on repeated failure testing on a single sample do not exist. Thus, finite element simulations must be used as “virtual” experiments, in place of physical experiments, to obtain the thousands of failure points necessary to create the continuous statistically based yield surface. These simulations may provide a means to understand and connect mechanisms of failure, such as trabecular ligament fracture, with the resulting failure envelope. Insight into these physical mechanisms will elucidate the mechanical response of trabecular bone.

The overarching aim of this study is to expand on previous work conducted by others concerning the uniaxial and biaxial response of trabecular bone [6,7,9,12], to include a multi-axial characterization of the failure response. Our objective is to develop a statistically based multi-axial yield surface that characterizes the mechanical response of the bone under various loading cases. To achieve this, we used finite element analysis (FEA) to examine hundreds of multi-axial loading scenarios each applied to 30 different microcomputed tomography (microCT)-derived trabecular bone samples. The final result from the simulations is a three-dimensional (3D) yield surface plotted in principal strain space. We used an approach similar to one that has been published previously by other researchers [32], but instead, considered the porcine skull bone rather than the human femur, analyzed additional microstructure models (30 rather than 3), and completed a statistical analysis of yield points in relation to the porosity of the bone models. The resulting compilation of yield points from thousands of simulations gives an overall characterization of the yielding behavior for trabecular bone and describes this behavior for uniaxial, biaxial, and triaxial strain states. In addition, by compiling results from 30 different microstructure samples, the surface captures differences in response as a result of varying bone architecture in the models.

## Methods

### Image Analysis and Model Development.

For this project, the skull of a 6 month old Gottingen mini pig was obtained and scanned at a spatial resolution of 50 *μ*m using a microCT scanner. This resolution falls within the range of resolutions used in other studies for similar research to minimize the error and capture important characteristics of the bone [15,33,34]. Since the main focus for this study is on the trabecular portion of the bone, microstructure samples were extracted from the full-skull dataset by means of cropping it to smaller cubic subsections. Cropped sections were roughly 3.15 × 3.15 × 3.15 mm^{3} in volume and taken from various locations throughout the front portion (anterior) of the skull as shown in Fig. 1 on the left. These cropped sections were analyzed in BoneJ [35] to determine the degree of anisotropy (DA). The images were then rotated so that their faces were aligned with the principal trabecular orientation. The *X*, *Y,* and *Z* axes are aligned to the fabric tensor of the sample's anisotropic structure. BoneJ automatically aligns the image to the fabric tensor. The largest eigenvalue is associated with the principal trabecular orientation becomes our *Y-*axis in the simulations. This was done consistently for all microstructures. Most models had a DA = 0.6. A total of 30 of these sections were cropped from various locations throughout the skull volume. Each cropped section thus had a different level of porosity and bone architecture. The range of bone volume to total volume, namely bone volume fraction (BV/TV), for each of the cropped sections was 0.392–0.481.

Based on these microCT scans, the three-dimensional surface geometry of the trabecular bone was generated computationally using Avizo (Visualization Sciences Group) by segmenting the bone to create a binary dataset of voxels with values of one for the labeled bone, and zero for anything other than bone. Figure 1(b) shows a sample of some of the cropped binary slices from the microstructure datasets. A three-dimensional triangulated surface of each of the microstructures was then generated from the segmentation using a modified form of the marching cubes algorithm available in Avizo. Smoothing of the surface model was performed in Avizo to eliminate surface staircasing (or castellation) using the smoothing algorithm of Taubin [36], which smooths while nearly preserving the original internal volume of the original surface. Figure 2 shows the 3D rendered surfaces of the 30 microstructures.

Volume meshes were created for each of the 30 microstructure samples in icemcfd (Ansys, Inc. 2600 ANSYS Drive, Canonsburg, PA) using four-noded tetrahedral elements. A mesh sensitivity and convergence study was completed to determine the balance between mesh density and reasonable compute time for the simulations and to determine the mesh density at which results started to converge. Meshes of three different densities were analyzed. The coarse mesh had around 50,000 elements, the medium had around 500,000, and the fine had around 1 million. The front surfaces of each mesh are shown in Fig. 3. We used two metrics for determining the mesh density, the first was that the stress response seemed to converge and the second was that the failure surfaces had similar or convergent topology. After running simulations on each mesh, the medium mesh, with an average edge length of approximately 0.05 mm, was chosen as it had a fine enough resolution to resolve the stress response and the failure surface in the bone without sacrificing quality of the bone's representation, but still had a reasonable run time.

### Boundary Condition Setup.

To develop a complete multi-axial yield criterion, loads that span all areas of strain space (uniaxial, biaxial, triaxial, and tension and compression) must be applied. To do this, we used a systematic approach to determine the required displacement boundary conditions.

Displacement loads were applied uniaxially, biaxially, or triaxially to the model in combinations of tension and compression. These displacements were applied in directions normal to the face of the microstructure in various ratios in the *x-*, *y-,* and *z*-directions. In order to determine these ratios, the equation of a sphere was spanned to obtain different points along the sphere's surface, where the *x*, *y*, and *z* coordinates of each point correspond to the applied displacements in the respective *x-*, *y-,* or *z*-direction. The sphere equation was spanned at increments of *θ* and Φ to give 376 unique combinations of *x*, *y*, and *z* coordinates. For each given value of Φ, a circle with a fixed *z*-coordinate and radius equal to sin(Φ) is created (shown by each individual circle in Fig. 4(a)). To determine the *x* and *y* coordinates, the circle with a fixed *z*-coordinate is spanned between 0 deg and 360 deg to get different *x* and *y* coordinates along the circle's perimeter. A graphical explanation of this method is outlined in Fig. 4 depicting the various circles with fixed *z*-coordinates and fixed radii along with the *x* and *y* coordinates spanned throughout each individual circle.

### Simulation Setup and Result Analysis Algorithm.

In order to accelerate and standardize this process, an algorithm was written in C to automatically setup, submit, and analyze all simulations for each model. The source code is available at the website.^{2} The algorithm works as follows: first, a file hierarchy is created with individual directories for each unique simulation. The code then writes ls-dyna (7374 Las Positas Road, Livermore, CA) input files with unique nodes, elements, and boundary conditions specific to the microstructure model. Next, each simulation is submitted and run on Penn State's high-performance computing cluster [34] using the portable batch system queuing system. Simulations were run in parallel on 16 processors and each took between 22 and 45 min to run depending on availability of computing resources during run time. To complete simulations on all 30 microstructures, a total of 11,280 simulations were executed (for about 135,360 central processing unit-hours). After the simulations were complete, another algorithm was developed for postprocessing the results. The stress and strain results were analyzed and averaged for each time-step in the specified node set of the microstructure model. To determine the yield point, a slope analysis was used to obtain the 0.2% offset method. First, the Young's modulus was determined by calculating the slope of the initial linear portion of the stress versus strain relationship, then, the yield point was determined using the 0.2% offset method. A combination of the stress–strain response in the *X*, *Y,* and *Z* directions were used to determine the location of the yield point in stress or strain space. The 0.2% offset method is used in all three directions, and the yield point in each direction is calculated independently of yielding in the other two directions. Then each of those individual stress–strain values in the *X*, *Y,* and *Z* directions are combined together to create a single coordinate that can be plotted in either 3D stress or strain space.

### Finite Element Simulations.

Nonlinear finite element analysis in ls-dyna was used for all simulations, and dynamic equilibrium was enforced by means of the weak form of the equation of motion, using explicit dynamics [37]. Displacement boundary conditions were used for all tension and compression simulations. An example of one of the applied boundary conditions can be seen in Fig. 4(c). Simulations were run at a nominal strain rate of 60 s^{−1} that has been confirmed to be quasi-static. For example, below a strain rate of 100 s^{−1}, rate effects are negligible [38].

A rate-independent, strength-asymmetric material model (*MAT_124 in ls-dyna) was used for all simulations, as symmetric yield properties for the tissue fail to capture the asymmetric behavior at the apparent level [7,15]. This allowed for an initial linear elastic (*E* = 13.8 GPa) relation between stress and strain up until the yield stress, followed by a nonlinear plastic portion where elements yielded plastically up until a certain plastic strain threshold. Elements reaching this threshold value of plastic strain were then removed to keep the simulation stable. Material properties were taken from the literature [15] as a best approximation for our given sample. The tissue modulus used was 13.8 GPa based on a sample of papers that offered a range of values. For example, 14 GPa reported by Ashman and Rho [39], 7 GPa by Zysset et al. [13], and 13.5 by Rho et al. [40]. The nonlinear plastic portion was perfectly plastic such that there was no hardening in the response. The nonlinear plastic portion was perfectly plastic such that there was no hardening in the response. An outline of additional material parameters can be seen in Table 1.

## Results

### Analysis of Yield Points.

From all of the simulations, yield points were combined together for analysis. Each unique loading condition resulted in a unique set of yield points. Not only are the yield points different for different loading conditions, but for the same loading condition, results also vary between the microstructures. Due to the differences in applied load, and the bone architecture of the microstructures, the resulting failure modes, fracture surfaces, and yield points were different. Yield strains from uniaxial and biaxial simulations are plotted in Fig. 5. Open and filled circles represent points from each of the 30 microstructure models, while the filled-in blue triangles are an average between all microstructures. Figure 5(a) shows yield points from simulations with loading in the *X-* and *Y*-directions. Figure 5(b) shows yield points from simulations with loading in the *X-* and *Z*-directions. Figure 5(c) shows yield points from simulations with loading in the *Y-* and *Z*-directions.

Bone architecture plays a large role in the failure of the bone and is responsible for the noticeable spread in yield points between the different microstructure models. This spread is attributed not only to a range of BV/TV in the microstructures but also to different bone architecture. Even for models with the same, or very similar, BV/TV, the results were quite different. Pore size and orientation and bone strut thickness and shape all play a role in the failure modes of the bone. Although an in-depth analysis of these specific characteristics beyond BV/TV was not conducted, it is evident that their differences are a contributing factor to yielding properties [41], as samples with similar BV/TV did not behave identically.

### Band of Failure Area.

To combine results together to incorporate yielding for all of the microstructure models in a more condensed fashion, the yield points were averaged between microstructures for each loading condition to get an averaged overall response. However, in order to fully characterize the response, the spread in the data must also be described. To do this, two additional sets of points were generated with a distance of one standard deviation (SD) from the mean (i.e., mean + standard deviation, and mean—standard deviation) shown in Fig. 6. Figure 6 shows the yield strains for biaxial simulations. Blue points represent average of points between all microstructure models under the same loading conditions. Yellow points are −1 SD from the mean, while orange points are +1 SD from the mean. Standard deviations were taken, exclusively between yield points from a certain loading scenario, and as such, each set three points (a yellow, blue and orange), characterizes the response for all 30 microstructures under one of the unique loading conditions. Including the standard deviation points in addition to the average gives a better representation of the overall response among all microstructures. Looking at the standard deviation for each point, the SD values are quite variable between all of the points and they tend to be larger in the compressive quadrant than that in the tension quadrant. For example, in the *XY* biaxial simulations, the average SD in the compressive quadrant is 0.29% while the average SD in the tension quadrant is 0.22%. From the spread of results shown in Fig. 6, the dependence on bone architecture is clear for the same loading condition, each microstructure with completely different bone architecture produces a different yield point.

### Yield Surface Optimization.

*x*,

*y,*and

*z*directions; $c1$, $c2$, and $c3$ are the shift in the center coordinates with respect to the origin; $b1$, $b2$, and $b3$ are the offsets of the surfaces from the center; and $k$ is slope of the surfaces. In the case when $c1=c2=c3$ = $c$ and $b1=b2=b3=b$, the yield surface is identical in three biaxial principal strain planes. After this simplification, an isotropic version of the yield surface with three parameters can be determined as follows:

*X*,

*Y,*and

*Z*directions when yielding happens. The magnitude of the three components (square root of the sum of the square of the three components) is the $||\epsilon FEA||$.

The optimized values can be found in Table 2. The negative center coordinates should be noted, as a strength asymmetric material model was used for the bone making it stronger in compression than in tension, thus off centering the surface in the purely negative quadrant. Although the yield surface is not identical in the three principal directions based on the seven-parameter model, the difference is very small. This indicates that the yield behavior under multi-axial loading is nearly isotropic.

The proposed parallelepiped envelope matches the FEA simulation results well as shown in Figs. 7 and 8. The yield points from FEA simulations are distributed closely to the parallelepiped surfaces. When using the three-parameter model, the *XY*, *XZ,* and *YZ* plane shapes are the same and Fig. 9 compares this shape to all the biaxial yield strain points from FEA simulations. Yield strains appear to be similar in the three principal directions.

If we plot the error distributions for the seven-parameter model and the three-parameter model (Fig. 10), we can see similar distribution of error and the accuracy of the seven-parameter model is just slightly higher. The mean $\xb1$ SD error norm is 3.77 $\xb1$ 3.89% for the seven-parameter model and 5.30 $\xb1$ 3.93% for the three-parameter model.

In addition to yield strains, we also investigated the yield stress of those microstructures. Figure 11 plots all the biaxial yield stress points from FEA simulations (points are averaged across microstructures).

## Discussion

The primary objective of this research was to develop a three-dimensional yield surface in strain space of trabecular bone. The equation for the three-dimensional yield surface derived from the resulting yield points can be used in new material models to govern the failure behavior of trabecular bone under a certain strain state, and this model must account for bone architecture. As shown by the spread of yield points in Fig. 6, the dependence on bone architecture is clear. In light of this, bone architecture must be included in models in order to obtain an accurate response. The bone architectural properties, such as porosity, how the trabeculae are oriented, their size and thickness, and their shape, may all have significant effect on the yield strain. As for porosity, some studies report that yield strains are independent of bone porosity [21–26], while others claim that the bone is stronger when it is more dense [27–29]. Our research result matches the former as we did not see a relationship between the yield strain and the BV/TV in our simulations (a plot of yield strain versus BV/TV is provided in Appendix B).

The yield strain of porcine trabecular bone under multi-axial loading is nearly isotropic while the yield stress shows more anisotropy, primarily in the *XY* plane. There is no significant difference between the biaxial yield strain in *XY*, *XZ,* and *YZ* planes as shown in Fig. 9. However, more variability has been seen in the biaxial yield stress among the principal planes (Fig. 11). Previous research on human trabecular bone [29], cancellous bone [25], and bovine trabecular bone [39] has also shown that the yield behavior is isotropic in strain space and this is consistent with our research. The anisotropy in stress as observed in Fig. 11 indicates that the trabecular bone has different strength in different directions. This is due to the underlying orientation of trabeculae with respect to the loading direction. Note that microstructures were aligned so the principal direction of trabeculae was aligned in the *Z*-axis.

In a previous work, a “modified super-ellipsoid” mathematical model was used to fit to the yield surface data; however, we found that this description did not match well with the current results. This could be explained by the different architectures of human trabecular bone and porcine trabecular bone. First, the human trabecular bone is more porous. The BV/TV for human trabecular bone used in the modified super-ellipsoid yield model [29] is 0.28–0.38 while the BV/TV for porcine trabecular bone used in our simulation is 0.392–0.481. Second, from the scan images, we can see human trabecular bone has sharp-edge pores while porcine trabecular bone has round-edge pores. The different pore shapes will affect strain distribution inside the bone and thus result in different overall yield strains. Thirdly, even within human trabecular bone, experimental data have shown that the modified super-ellipsoid yield model cannot fit well due to the broad range of architectural properties of trabecular bone [13].

There are a number of reasons why the presented computational micromechanical simulations are advantageous. First, 376 unique multi-axial loading cases were simulated on each of the 30 microstructure models and this is impossible to achieve with mechanical tests due to the destructive nature of mechanical test. In addition, mechanical tests are costly in both time and money. Second, the proposed yield function (1) for porcine trabecular bone has very simple format and it can be implemented into a macroscopic model for the prediction of failure of whole bones. The simplified function (2) only has three coefficients and this further facilitates the implementation of this model. Third, 11,280 simulations provides us with useful statistical data such as the SD of the yield strains among 30 samples. The spread of data between each of the microstructures is important to note as it appears to be characteristic of the microstructural response. A yield surface that captures this response could incorporate this spread into the model. We calculated the average yield strains for the 30 samples (i.e., the yield envelop in Fig. 7) and we also calculated the SD of the yield strains (SD = 0.25%). When we use the envelop for prediction, we know that it is quite possible that the actual yield strain will fall between the predicted value ±0.25%. Or we could run the macroscopic model with a variation of the critical yield stress by the SD. The standard deviation of points in our analysis also aids with this development. The inclusion of the standard deviation gives much more conclusive information about the envelop that points generally fall within, than just the average alone. Finally, the yield surface is formulated in strain space and the dependency on the apparent density of each sample could be eliminated.

In Fig. 11, we observed that the average biaxial yield stress in the *XZ* and *YZ* directions are similar in both the first and second principle directions. However, while the biaxial yield stress in the *XY* direction was also similar to the *XZ* and *YZ* directions in the second principle direction, we observed that it was slightly larger (about 36% difference) in the first principle direction. This suggests that the underlying microstructure of bone, including its directionality, may influence the macroscopic yielding response. Directionality of bone has been observed to play an important role in other bones, like the femur, where bone structure plays an important role in distributing load [42]. Our results show there is some slight anisotropic yield response due to the underlying geometric configuration of trabeculae, but the fundamental reasoning for this is not clear. Further research is needed to characterize the anisotropy of trabecular bone with respect to loading patterns the skull experiences, say from internal pressure or growth.

## Limitations

Despite many strengths of our approach, there are limitations. Material properties were taken from experiments run on human femoral trabecular [12] since there are very limited experiments on porcine trabecular skull bone. The biological difference could affect the simulated values. The spatial resolution is 50 *μ*m for the scan images and smoothing has been done to eliminate surface staircasing in Avizo. This smoothing nearly preserves the original internal volume but still has some effect on the sample geometry. All the 30 samples have the same dimensions of 3.15 × 3.15 × 3.15 mm^{3}. Previous research has shown that 5 mm length representative volume element is large enough to capture the features of trabecular bone [40–42]. Our simulation results with the sample size of 3.15 mm also captured many features such as isotropic yield behavior in strain space but it could be possible that some features are not fully represented by the relatively smaller sample size.

We conclude that yield strain is independent of BV/TV based on our simulations of the 30 samples, in agreement with other studies [27–29]. However, this conclusion is limited by the BV/TV span of our samples ranging from 0.392 to 0.481. The wide spread of yield strains of the 30 samples as shown in Fig. 5 indicates that the variability of the bone architecture affects the yield strain. The bone architecture is dependent on the sample location in the skull but we have not fully evaluated all locations of the skull. This research will be done when more samples are available to reach an additional conclusion. Additionally, the gender of the specimen was not provided, and failure results could differ between the two genders. Furthermore, although pore/bone strut size and shape most likely play a role in the yielding properties of the bone, this study only analyzed the effect of BV/TV and simply aimed to characterize the overall response mathematically. Further in-depth analysis would need to be conducted to systematical determine the effect on yielding properties. However, note that porosity is a variable that changes from location to location in the skull. Therefore, we believe our approach of considering a wide range of porosity and reporting the average response is appropriate.

## Conclusions

Developing a multi-axial failure criterion for porcine trabecular skull bone has many clinical and biological implications such as the use in head injury studies. Since there has been limited research conducted around multi-axial loading of trabecular bone, our aim was to use microstructure-level simulations of trabecular bone to develop a computationally derived yield surface. This was done by analyzing thousands of loading scenarios on 30 different microstructure models, and using the results from those to develop a yield surface plotted in three-dimensional strain space. A mathematical equation was optimized to fit to the resulting yield points and gave a good representation of the failure behavior in three-dimensional space.

## Acknowledgment

Thanks to ARL, HPCMP, and Penn State for using of their computing resources. Thanks to Stephen Alexander, Christopher Hoppel, and Michael Kleinberger from ARL for providing data, discussion, and support for the work.

## Funding Data

Army Research Laboratory through Massachusetts Institute of Technology, Institute for Soldier Nanotechnologies (W911NF-13-D-0001).

### Appendix A

Below are BV/TV for each microstructure.

### Appendix B

Below is a plot of yield strain versus BV/TV, where no relationship is observed.