Recent studies demonstrated an association between atypical femoral fracture (AFF) and long-term bisphosphonate (BP) use for osteoporosis treatment. Due to BP treatment, bone undergoes alterations including increased microcrack density and reduced tissue compositional heterogeneity. However, the effect of these changes on the fracture response of bone is not well understood. As a result, the goal of the current study is to evaluate the individual and combined effects of microcracks and tissue compositional heterogeneity on fracture resistance of cortical bone using finite element modeling (FEM) of compact tension (CT) specimen tests with varying microcrack density, location, and clustering, and material heterogeneity in three different bone samples. The simulation results showed that an increase in microcrack density improved the fracture resistance irrespective of the local material property heterogeneity and microcrack distribution. A reduction in material property heterogeneity adversely affected the fracture resistance in models both with and without microcracks. When the combined changes in microcrack density and tissue material property heterogeneity representing BP treatment were evaluated, the models corresponding to BP-treated bone demonstrated reduced fracture resistance. The simulation results also showed that although microcrack location and clustering, and microstructure significantly influenced fracture resistance, the trends observed on the effect of microcrack density and tissue material property heterogeneity did not change. In summary, these results provide new information on the interaction of microcracks, tissue material property heterogeneity, and fracture resistance and may improve the understanding of the influence of mechanical changes due to prolonged BP use on the fracture behavior of cortical bone.

Introduction

Atypical femoral fracture (AFF) is a rare catastrophic fracture of femur which results in high rate of morbidity and mortality. Several recent studies demonstrated a relationship between AFF and long-term bisphosphonate (BP) use for osteoporosis treatment [16]. This possible relationship indicates that there are tissue level changes that take place as a result of BP treatment that may lead to adverse mechanical outcomes. However, the changes in the mechanical response of bone due to BP treatment are not well understood.

The recent studies have shown that the long-term use of BP treatment may result in a number of mechanical alterations in the bone tissue including a reduction in compositional heterogeneity and an increase in microcrack density. BP treatment suppresses bone turnover and reduces the removal of microdamage which consequently increases the microcrack density [7]. The studies on human and animal subjects showed that prolonged BP treatment increases microcrack density in bone [8,9]. The increase in microcrack density is shown to depend on the duration and dosage of BP treatment where more microcracks were observed with longer and higher dosage of BP treatment [10,11]. A recent case study conducted on a breast cancer patient who experienced an AFF after 9 years of BP treatment reported very high levels of microcrack density compared to the microcrack levels reported in other BP treatment studies [12]. On the other hand, some studies suggest that the microcrack density does not depend on BP treatment and is related to the patient's daily activities, i.e., mechanical loading of the bone [13,14]. In addition, the existing literature shows that BP treatment does not change the length of microcracks [10,11,15,16].

The studies conducted on bone biopsies from patients under BP treatment and patients with fracture showed that the tissue compositional heterogeneity was reduced compared to healthy bone [1724]. In addition, reduced heterogeneity in elastic modulus [25], a decreasing trend in energy to fracture [26], and reduced fracture toughness [27] were reported in experiments performed on long-term BP-treated bone. A deterioration in microindentation properties was also observed in AFF patients [28]. However, some studies found no relation between BP use and tissue compositional heterogeneity [2932]. These results indicate that there is no consensus on the effects of BP use on tissue compositional heterogeneity.

In the literature, microcrack formation and growth has been investigated using finite element modeling (FEM) both in trabecular and cortical bone. The existing models simulated microcrack formation and propagation mainly in small sections of cortical bone under tension or compression and represented the cortical bone as defect free at the beginning of the simulations (i.e., no preexisting microcracks in the models prior to the simulations) [3339]. As a result, these studies only focused on the influence of damage and crack formation that were the outcome of the simulations. Other studies developed a multiscale approach in which microcrack accumulation in trabecular bone was simulated and this information was transferred to macroscale femur simulations to assess fatigue damage [40,41]. In addition, at the macroscale, femur fracture was modeled through damage accumulation in cortical and trabecular bone [42]. Furthermore, some studies investigated the interaction of a single preexisting microcrack with osteons [43,44]. However, these studies did not incorporate the influence of preexisting microcrack distribution on fracture resistance. In addition, the influence of tissue heterogeneity on mechanical properties of bone was investigated by only a few studies. Two of these studies focused on the assessment of elastic material property distribution on energy dissipation [45,46]. Several studies also accounted for different material properties of osteons, interstitial bone, and cement lines [34,3739]. Our previous study evaluated the influence of heterogeneity of spatial distribution of material properties on fracture resistance of human cortical bone [47]. None of these studies implemented both preexisting microcrack distribution and compositional tissue heterogeneity and evaluated their combined effects on fracture resistance of cortical bone.

In summary, the previous experimental studies demonstrate the possible link between BP treatment and increased preexisting microcrack density and reduced tissue compositional heterogeneity. However, the effect of these factors on fracture resistance of bone has not been clearly identified either experimentally or computationally. The goal of the current study is to evaluate the interaction of preexisting microcracks and tissue compositional heterogeneity and quantify their influence on crack growth in cortical bone using computational modeling. The FEM approach developed in this study makes it possible to identify the individual and combined effects of preexisting microcracks and tissue compositional heterogeneity. The results provide a better understanding of the influence of BP treatment on the mechanical response of bone and may provide additional insight into the possible association of long-term BP treatment and AFF.

Materials and Methods

The effects of preexisting microcracks and tissue composition heterogeneity on crack propagation behavior in human cortical bone were investigated by modeling crack growth in compact tension (CT) specimen tests. The models represented varying (a) preexisting microcrack density, (b) preexisting microcrack location and clustering, and (c) material heterogeneity on three bone samples with different microstructures.

Finite Element Models

Model Generation.

In this study, the finite element models were generated based on the transverse microscopy images of three male cortical bones from the mid-diaphysis of tibia denoted as Models A, B, and C, respectively (Figs. 1(a)1(c)). The properties of each bone microstructure are reported in Table 1.

The 3D cortical bone models were generated by using abaqus finite element software (version 6.13, Simulia, Providence, RI) as outlined in detail previously [37,47]. In this approach, two-dimensional sketches of the images (including cement lines, osteons, Haversian canals, and interstitial bone) were extruded by 1 mm in the third dimension to obtain 3D models which have an approximate volume of 1 mm3 (Figs. 1(d)1(f)). After generation of the detailed microstructure region, the finite element model was assembled by placing the microstructure into the CT specimen with a 1 mm3 cutout section (Fig. 1(g)) such that osteons were oriented perpendicular to the crack growth direction. An initial crack with a depth of 0.1 mm was introduced to the midplane of the microstructure in order to reduce computational time without altering the outcome of the simulations. Mode I incremental displacement boundary condition was applied to the CT specimen at the loading pin holes (Fig. 1(g)). CT specimen was also fixed in all directions at the midplane at the far end from the loading point (Fig. 1(g)). All the models were evaluated at the same loading level which corresponded to the step at which the damage reached the detailed microstructure section boundary in the model with the highest crack growth amount. The microstructure and the CT specimen were meshed with tetrahedral elements except for the cement lines which were represented with zero thickness triangular interface elements. The total number of elements ranged between 1.6 and 1.9 million. The simulations were run using parallel computing systems utilizing 20 cores with a GPU accelerator or 48 cores without a GPU accelerator.

Material Properties.

In this study, isotropic elastic properties were used for the CT specimen and the detailed microstructure region (Table 2) [48]. Poisson's ratio of 0.3 was assigned to all models. The crack behavior in the detailed microstructure region and the CT specimen was modeled by cohesive extended finite element method (XFEM) with the exception of cement lines which were represented by cohesive interface elements [47]. The details of the computational modeling technique are provided in the following section (Sec. 2.2). The microstructure of the CT specimen was not modeled in detail; therefore, the homogenized material properties were assigned to it for all simulations (Table 2) [4851]. The detailed microstructure region was assigned both homogeneous and heterogeneous material properties. The homogeneous material properties were represented by mean values in Table 2 based on the experimental data reported in the literature [48,50,52,53]. The models with compositional heterogeneity in the detailed microstructure region were assigned ten different material property values in the cement lines, osteons, and the interstitial bone. This number was chosen based on our previous study [47] and corresponds to the highest fracture resistance among other heterogeneity levels investigated. The heterogeneous models were created with ten different material property sets composed of different elastic moduli, ultimate strengths, and fracture toughness values based on the mean values and standard deviations reported in Table 2 [48,50,52,53]. Cement lines, osteons, and the interstitial bone each have their own sets of materials. The heterogeneous material properties were created and assigned randomly to the elements in the detailed microstructure region. All models were assigned the same material property sets with different distributions. The random material property generation and assignment were done by a matlab (MathWorks, Natick, MA) script.

Microcracks.

The number and the length of the preexisting microcracks were determined based on the previously reported values in the literature for healthy bone [13,54]. The distribution of the preexisting microcracks in the models is shown in Fig. 2. A microcrack density of 5 microcracks/mm2 was selected for all the models representing healthy bone based on the curve fit for age-related variation in microcrack density provided in the literature [54]. This value corresponds to 81-year-old bone which is the highest age group among the models in the current study. The number of preexisting microcracks was reported to increase between 1.3- and 2-fold in BP-treated bone compared to healthy bone [16,55]. Based on this information, the number of preexisting microcracks in the healthy bone models was doubled to 10 microcracks/mm2 for models representing long-term BP treatment. In order to create preexisting microcracks in the models, random coordinates were selected from the midplane region of the interstitial bone within the detailed microstructure region. The regions around the preexisting microcracks were meshed finer compared to the rest of the model to capture the initial area of influence of the preexisting microcrack accurately. The studies in the literature show that the BP treatment does not alter the preexisting microcrack length [10,15,16], thus the length of preexisting microcracks was selected randomly based on the reported range of preexisting microcrack lengths for healthy bone in the literature [10,13,56,57]. The preexisting microcrack lengths in this study ranged from 35 μm to 52 μm.

Finite Element Simulations Run.

A total of 32 simulations were performed in this study. The simulations investigated the influence of variation in preexisting microcrack density, location, and clustering as well as material heterogeneity.

The preexisting microcrack density was evaluated by simulating two preexisting microcrack density levels for all models which were 5 microcracks/mm2 representing healthy bone and 10 microcracks/mm2 representing BP-treated bone. For comparison purposes, simulations with no preexisting microcracks were also performed for all models. The effect of preexisting microcrack location was assessed by generating two preexisting microcrack distributions for each model composed of different and identical size and distribution of preexisting microcracks in each model (Fig. 2). The influence of compositional heterogeneity was evaluated by generating models both with homogeneous and heterogeneous material properties as outlined in Sec. 2.1.2.

To identify the simulations, labels “A,” “B,” “C” for model name, “5,” “10” for preexisting microcrack density, “HM,” “HT” for material composition, and “I,” “D” for preexisting microcrack distribution were used. For example, A5HMD represents model A with homogeneous material properties and has five preexisting microcracks with different preexisting microcrack sizes and distributions compared to other models. The full naming scheme for simulations is shown in Table 3. Additionally, two simulations for model A with 10 clustered preexisting microcracks were performed in order to understand how clustered preexisting microcracks influence the crack propagation. The size of the preexisting microcracks in the clustered models was the same as the identical preexisting microcrack models but the preexisting microcracks clustered in a certain region of the interstitial bone instead of being randomly distributed.

Cohesive Finite Element Modeling.

The cohesive finite element modeling was implemented to simulate crack propagation in the models. This modeling technique captures the nonlinear behavior of the process zone around a crack. The incorporation of the cohesive modeling approach in the models was done by using the XFEM [58] in the osteons and the interstitial bone and cohesive interface elements [59] at the cement lines.

In the current study, a previously developed cohesive traction-crack opening displacement relationship for bone was implemented (Fig. 3) [47,60]. This relationship is expressed as a function of the ultimate strength (σc), the critical energy release rate (Gc), and the critical crack opening displacement (δu). A quadratic equation was used to define the damage initiation (Tnnc)2 + (Tssc)2 + (Tttc)2 = 1, where Tn, Ts, and Tt are the current stress values in normal and shear directions and σnc, σsc, and σtc are critical normal and shear strengths, respectively [61]. The damage accumulation starts once the quadratic equation equals to one. Damage evolution follows the traction-crack opening displacement curve and the crack follows a mixed-mode fracture criterion (Gn/Gnc) + (Gs/Gsc) + (Gt/Gtc) = 1, where Gn, Gs, and Gt are the current energy release rates in normal and shear directions and Gnc, Gsc, and Gtc are critical energy release rates in normal and shear, respectively [61]. Full fracture occurs once the fracture criterion is satisfied. For damage initiation and evaluation, the values for shear in both directions were assumed to be the same (Gsc = Gtc, σsc = σtc). The undamaged response of the interface elements is characterized by an initial slope that is a penalty stiffness in the numerical formulation [61], whereas the initial slope is not required for the XFEM formulation. The details of the cohesive modeling technique for bone were outlined in Refs. [37], [47], and [60].

The cohesive fracture behavior of the cement lines was modeled as zero thickness interface elements. The rest of the model was formulated by XFEM which is defined by local enrichment functions to represent discontinuities created within the elements by the formation of cracks [58]. In this method, the crack does not need to align with the element boundaries and is independent of the mesh.

Analysis of Results.

The simulation results were assessed to compare the fracture behavior of cortical bone under varying preexisting microcrack density, location, and material heterogeneity. A previously developed Python script [47] was used for extracting the simulation results. The fracture resistance was assessed by quantifying the total volume of the elements which were fully fractured as well as the damage energy density that represents the total dissipated energy per unit volume associated with damage and crack formation. An element starts accumulating damage when the damage initiation criterion outlined in Sec. 2.2 is satisfied. An element is considered fully cracked when the fracture criterion outlined in Sec. 2.2 is attained. The crack volume was reported instead of using the macrocrack length due to the out-of-plane crack deflection. The crack volume information captures the amount of crack formation. A higher crack volume indicates that it is easier to form cracks in a model and corresponds to a lower fracture resistance compared to a model with less crack volume. The damage energy density is calculated by dividing the total energy dissipated due to damage and crack formation by the total volume of elements that are damaged or fractured. This parameter reflects the intensity of damage in the tissue. In addition, the crack formation pattern in the models was determined by identifying the elements which have undergone full fracture and have accumulated damage but have not formed a complete crack.

Results

The simulation results showed that an increase in preexisting microcrack density led to improved fracture resistance. The crack volume (Fig. 4) and the damage energy density (Fig. 5) were the highest when there were no preexisting microcracks in all three models for both homogeneous and heterogeneous material properties and for both preexisting microcrack distributions. Increasing the preexisting microcrack density from 0 to 5 and then to 10 microcracks/mm2 resulted in lower crack volume and damage energy density in all models irrespective of the material heterogeneity and preexisting microcrack distribution (Figs. 4 and 5). Among the models with different preexisting microcracks (5HMD, 10HMD, 5HTD, and 10HTD), model A had the highest crack volume reduction rate. The highest crack volume reduction rate was seen in model C in models with identical preexisting microcracks (5HMI, 10HMI, 5HTI, and 10HTI). The lowest crack volume reduction rates were observed in model B with the exception of 5HTD and 10HTD models that had the lowest crack volume reduction in model C. The reduction in crack volume with increasing number of preexisting microcracks was higher in HM models for model A and model C, and in HT models for model B (Fig. 4). Introducing preexisting microcracks to both homogenous and heterogeneous models reduced the crack volume mainly in interstitial bone and to a limited degree in osteons (Fig. 6).

The heterogeneous material distribution resulted in lower crack volume compared to the homogeneous material properties in all models with and without preexisting microcracks and for both preexisting microcrack distributions (Fig. 4). The damage energy density was also lower in models with heterogeneous material properties than the homogeneous ones with the exception of B5HTD and B10HTD (Fig. 5). This difference in behavior was caused by a large amount of crack volume with low amount of damage that reduced the per volume damage energy amount in these two models.

The effect of BP treatment was evaluated by comparing the models with homogeneous material properties and 10 microcracks/mm2 to models with heterogeneous material properties and 5 microcracks/mm2 representing BP-treated and healthy bone, respectively. The results showed that the crack volume was higher in models representing BP-treated bone than models representing bone with no BP-treatment with the exception of Model A with different preexisting microcrack distributions (Fig. 7(a)). Similar trends were observed in damage energy density with the exception of model B with different preexisting microcrack distributions (Fig. 7(b)).

The distribution of preexisting microcracks had a significant influence on crack propagation. Two different preexisting microcrack distributions simulated for each model resulted in different levels of crack volume and damage energy density (Figs. 4 and 5). The relative values of crack volume and damage energy density obtained in the models varied based on the microstructure.

Clustered preexisting microcrack models (AHMC1 and AHMC2) resulted in lower crack volume compared to their zero crack HM model counterpart (AHM) but the reduction was much lower compared to the distributed preexisting microcrack models resulting in higher crack volume than A10HMI and A10HMD (Fig. 4). Even though the crack volumes of AHMC1 and AHMC2 were close, the crack propagation patterns were completely different from each other (Fig. 6(d)). The crack tended to propagate away from the clustered regions and demonstrated a larger amount of crack growth away from these regions. Damage energy density was higher in AHMC2 but lower in AHMC1 compared to the distributed preexisting microcrack models (A10HMI and A10HMD). The difference in damage energy density behavior was caused by larger low damage regions leading to lower damage energy density in AHMC1 compared to AHMC2.

Cortical bone microstructure had a significant effect on the crack propagation behavior. When preexisting microcracks with identical size and location were incorporated in all three models, the extent of crack propagation attained in each model varied (Fig. 4). However, the same trends were observed in all microstructures in terms of the influence of preexisting microcrack density, location, and material heterogeneity. The highest crack volume was observed in model B followed by models A and C (Fig. 4). Damage energy density was also predominantly the highest in model B with the exception of B5HMD (Fig. 5).

Discussion

This study presented a new approach to assess the combined influence of preexisting microcracks and tissue composition heterogeneity on fracture resistance of cortical bone. The computational studies in the literature did not simultaneously incorporate the influence of preexisting microcrack distribution and spatial variation of material properties on fracture resistance. The results presented here address the knowledge gap on the interaction of preexisting microcracks and tissue material property heterogeneity in determining the fracture resistance in cortical bone.

The simulation results showed that an increase in preexisting microcrack density within the preexisting microcrack density range up to 10 microcrack/mm2 improved the fracture resistance evidenced by the reduction in crack volume and damage energy density irrespective of the local material distribution, preexisting microcrack location, and microstructure. This finding is in line with the previous experimental studies that demonstrated the contribution of preexisting microcracks to fracture toughening in cortical bone [62]. Limited microcracking provides a protective mechanism against crack growth by promoting energy dissipation and activating other toughening mechanisms such as crack deflection [63]. There is conflicting information in the literature about the influence of BP on preexisting microcrack formation. Some studies reported an increase in preexisting microcrack density with BP treatment in animals and humans [8,9], whereas other studies reported no significant difference in preexisting microcrack density [13,14]. In this study, a twofold increase in preexisting microcrack density was incorporated in the models representing the BP-treated bone to account for a possible increase in preexisting microcrack density with treatment. Without any other changes in the tissue, such as alterations in material property heterogeneity, our results indicate that only a possible twofold increase in preexisting microcrack density compared to healthy bone does not lead to an impairment in the fracture response of cortical bone. However, the simulations demonstrated that the heterogeneous material property distribution had a pronounced influence on the crack resistance of cortical bone by reducing the crack volume and damage energy density in models both with and without preexisting microcracks. This indicates that although limited number of additional preexisting microcracks that may form as a result of BP treatment does not individually impair fracture resistance, further changes in the tissue due to treatment such as reduction in material heterogeneity may lead to an overall decrease in fracture resistance.

Among the current models, homogeneous material properties with 10 microcracks/mm2 and heterogeneous material properties with 5 microcracks/mm2 can be assumed to represent BP-treated and healthy bone, respectively, based on experimental observations summarized earlier. Comparison of the results of these models showed that the crack volume and the damage energy density were predominantly lower in models representing no BP treatment compared to models representing BP-treated bone for both preexisting microcrack distributions. These results further demonstrate that the adverse effects of reduced tissue compositional heterogeneity on fracture resistance may overshadow the improvement in fracture resistance due to preexisting microcracks. In addition, it was observed that the beneficial influence of the preexisting microcracks may be more pronounced in the absence of local tissue material property heterogeneity as demonstrated by the results obtained for models A and C. This result may indicate that the loss of protective mechanism that is provided by tissue heterogeneity may be mitigated to a degree by a possible increase in preexisting microcracks.

As previously postulated, the occurrence of AFF may be due to the transformation of one of the microcracks into a macrocrack due to the lack of material heterogeneity rather than due to the existence of many preexisting microcracks [63]. This is supported by our results where an increase in the number of preexisting microcracks did not have a detrimental influence on the fracture response. In our study, individual preexisting microcracks did not attain a state that will transform them into macrocracks, therefore, this aspect of AFF was not evaluated in the current study.

The simulation results showed that the fracture resistance of bone was not only influenced by the preexisting microcrack density but also by the location of preexisting microcracks. Two different preexisting microcrack distributions incorporated in the same models led to different crack propagation amounts and patterns highlighting the effect of preexisting microcrack location on fracture resistance. The results demonstrated that the crack growth was impeded in the vicinity of the preexisting microcracks. For example, comparison of A10HMD and A10HMI (Fig. 6(a)) showed that the three preexisting microcracks located on upper right region of A10HMD provided a higher resistance to crack growth compared to A10HMI which had only one preexisting microcrack in the same region. Furthermore, the clustered preexisting microcrack models showed that the crack volume increased significantly compared to distributed preexisting microcrack models with similar preexisting microcrack size and density. There was also a significant difference in the crack propagation patterns between the two clustered models due to the location of the clustered preexisting microcracks. These results show that a region with a high density of preexisting microcracks may not provide the same benefit as distributed preexisting microcracks as it will influence only a restricted area instead of the entire bone area.

The microstructure also played a significant role in the crack growth patterns. This was highlighted by different amounts of crack growth and damage energy density in each model when preexisting microcracks with identical location and size were incorporated in the models. The important contribution of microstructure to crack growth was established in our previous study [47] as well as in other computational [38,39] and experimental studies in the literature [49,64]. Although, the number of samples investigated in the current study was not sufficient to establish a correlation between the microstructural parameters and the simulation results, a qualitative evaluation indicated that models A and C had similar osteonal area and density had similar amount of crack propagation. The higher amount of crack propagation observed in model B is most likely due to the extensive interstitial bone region in the initial crack growth region where there are no osteons to impede crack growth (Fig. 2(b)). This observation is supported by the previous computational studies that showed reduced fracture strength with reduced osteon density [38]. The same computational study also demonstrated a reduction in fracture strength with increasing porosity [38]. In the bone specimens that were investigated in the current study, the porosity amount was similar between the models, therefore, the variation in crack growth behavior was not related to porosity. Although the level of crack propagation was different in each microstructure, the trends in the effect of preexisting microcrack density, material heterogeneity, and preexisting microcrack location were the same in all microstructures. Therefore, although percent variations may vary with microstructure, the main conclusion of the study is not dependent on microstructure.

In this study, we simulated a fracture toughness experiment at the bulk tissue level via a CT specimen to measure the fracture resistance of bone at the millimeter-scale. The purpose of these experiments is to determine the local fracture response of bone which contributes to the fracture behavior at the whole bone scale. The typical outcome of these tests is the stress intensity factor which combines the influence of load and crack size. In the current study, due to the complexity of crack propagation perpendicular to the osteons, including crack deflection at the cement lines, the equations developed for calculating stress intensity factors for the CT specimen [65] are no longer valid. As a result, the amount of crack formation (represented by crack volume) in each model at the same load level and the damage energy density were reported. Comparing the crack volume at the same load level provides information on the fracture resistance where a high crack volume indicates that crack can easily propagate. If the amount of crack formation is less in a model, this indicates an increased resistance to fracture in the model (Fig. 6). In addition, damage energy density was used as a surrogate measure of fracture resistance similar to the previous studies in composite materials [66]. A high damage energy density indicates large amount of energy dissipation to form a small crack representing increased fracture resistance. On the other hand, a low damage energy density corresponds to low energy dissipation to generate a large crack representing reduced fracture resistance.

The strength of the current study is its capability to evaluate the independent and combined effects of tissue material property heterogeneity and preexisting microcrack on fracture resistance of cortical bone. This approach isolated the influence of each factor as well as assessed the interaction among these factors in determining the fracture behavior of bone which cannot be attained in experiments due to other confounding factors that may be present in specimens.

One of the limitations of the study is that the age-related changes in the material properties and the preexisting microcrack densities were not incorporated in the models. All models were assigned the same properties and preexisting microcrack densities due to the lack of complete set of data presented in Table 2 for all age groups. Furthermore, modifications in the material properties due to BP treatment besides the spatial variation were not included in the models. In the current study, the material properties utilized in the models (Table 2) represent material properties of healthy bone. Due to BP treatment, the degree of mineralization is expected to increase which may modify these material properties. However, a complete set of experimental data for BP-treated bone that reflects these possible modifications is not currently available in the literature. As a result, it was not possible to implement these changes in the current study. The measurements performed on canine bone have shown that there is no significant difference between elastic modulus and ultimate strength but there is a reduction in bone toughness in BP-treated bone compared to healthy bone at the bulk tissue scale [10,67]. Incorporation of a reduced toughness in models representing BP-treated bone will further decrease the fracture resistance observed in our simulations and lead to more extensive crack propagation. Therefore, if these material changes are incorporated in the models, the trends reported in this study in Fig. 7 between models representing BP treatment and no BP treatment will become even more pronounced.

In the current study, the material heterogeneity was represented by experimentally measured average material properties and their standard deviations instead of locally measured bone material properties specifically for each bone sample. This approach was chosen due to the lack of experimental measurements of local specimen specific elastic and fracture material properties of bone at the length scale of material property assignment in the current models. However, the computational modeling approach utilized in this study can readily be combined with direct experimental local measurements of tissue material properties on bone biopsies when such measurements become available.

Another limitation is that isotropic elastic properties were assigned to the bone instead of anisotropic elastic properties. A full set of anisotropic elastic engineering constants has not been determined in the literature at the scale of material assignment in the current models. Therefore, elastic anisotropy was not implemented in the current study. On the other hand, anisotropy in the fracture toughness of bone was captured in this study by the detailed cortical bone microstructure in the models. It was shown experimentally that the fracture toughness is approximately the same in transverse and longitudinal crack growth in cortical bone at very small crack lengths [53]. It starts to diverge as the crack grows due to the influence of microstructure which leads to the anisotropy in bone fracture toughness. As a result, elastic anisotropy may affect the stiffness of the model but is not expected to influence the crack growth behavior significantly which is governed by the crack growth direction with respect to the microstructure and the material distribution both of which were already included in our models.

In the current study, all preexisting microcracks were introduced in the interstitial bone and did not incorporate preexisting microcracks that have penetrated into osteons. The incorporation of all preexisting microcracks in the interstitial bone was implemented based on the previous studies which reported that the majority of the preexisting microcracks are located in the interstitial bone [54,57]. Therefore, the models are expected to closely reflect in vivo conditions. Another limitation is that the models were assumed to be uniform in the third dimension which is not the case in an actual bone sample. However, this assumption is not expected to significantly affect the results due to the very small size of the detailed microstructure region.

In conclusion, this study provided a new computational modeling approach to evaluate the individual and combined effects of preexisting microcracks and tissue compositional heterogeneity on fracture resistance of cortical bone. These factors reflect the tissue level changes that may occur as a result of long-term BP treatment. The new approach outlined here provides a tool that can be used to evaluate the mechanical response of BP-treated bone and may also bring additional insight into the possible association of long-term BP treatment with AFF.

Funding Data

  • National Science Foundation (Award No. CMMI-1434412).

  • The computational resources for this study were supported by XSEDE (Award No. MSS-160001).

References

References
1.
El Rachkidi
,
R.
,
Sari-Leret
,
M. L.
, and
Wolff
,
S.
,
2011
, “
Atypical Bilateral Pedicle Fracture in Long-Term Bisphosphonate Therapy
,”
Spine (Phila Pa 1976)
,
36
(
26
), pp.
E1769
–E
1773
.
2.
Dell
,
R. M.
,
Adams
,
A. L.
,
Greene
,
D. F.
,
Funahashi
,
T. T.
,
Silverman
,
S. L.
,
Eisemon
,
E. O.
,
Zhou
,
H.
,
Burchette
,
R. J.
, and
Ott
,
S. M.
,
2012
, “
Incidence of Atypical Nontraumatic Diaphyseal Fractures of the Femur
,”
J. Bone Miner. Res.
,
27
(
12
), pp.
2544
2550
.
3.
Gedmintas
,
L.
,
Solomon
,
D. H.
, and
Kim
,
S. C.
,
2013
, “
Bisphosphonates and Risk of Subtrochanteric, Femoral Shaft, and Atypical Femur Fracture: A Systematic Review and Meta-Analysis
,”
J. Bone Miner. Res.
,
28
(
8
), pp.
1729
1737
.
4.
Meier
,
R. H.
,
Perneger
,
T. V.
,
Stern
,
R.
,
Rizzoli
,
R.
, and
Peter
,
R. E.
,
2012
, “
Increasing Occurrence of Atypical Femoral Fractures Associated With Bisphosphonate Use
,”
Arch. Intern. Med.
,
172
(
12
), pp.
930
936
.
5.
Schilcher
,
J.
,
Koeppen
,
V.
,
Aspenberg
,
P.
, and
Michaëlsson
,
K.
,
2015
, “
Risk of Atypical Femoral Fracture During and After Bisphosphonate Use: Full Report of a Nationwide Study
,”
Acta Orthop.
,
86
(
1
), pp.
100
107
.
6.
Shane
,
E.
,
Burr
,
D.
,
Abrahamsen
,
B.
,
Adler
,
R. A.
,
Brown
,
T. D.
,
Cheung
,
A. M.
,
Cosman
,
F.
,
Curtis
,
J. R.
,
Dell
,
R.
,
Dempster
,
D. W.
,
Ebeling
,
P. R.
,
Einhorn
,
T. A.
,
Genant
,
H. K.
,
Geusens
,
P.
,
Klaushofer
,
K.
,
Lane
,
J. M.
,
McKiernan
,
F.
,
McKinney
,
R.
,
Ng
,
A.
,
Nieves
,
J.
,
O'Keefe
,
R.
,
Papapoulos
,
S.
,
Howe
,
T. S.
,
van der Meulen
,
M. C. H.
,
Weinstein
,
R. S.
, and
Whyte
,
M. P.
,
2014
, “
Atypical Subtrochanteric and Diaphyseal Femoral Fractures: Second Report of a Task Force of the American Society for Bone and Mineral Research
,”
J. Bone Miner. Res.
,
29
(
1
), pp.
1
24
.
7.
Hirano
,
T.
,
Turner
,
C. H.
,
Forwood
,
M. R.
,
Johnston
,
C. C.
, and
Burr
,
D. B.
,
2000
, “
Does Suppression of Bone Turnover Impair Mechanical Properties by Allowing Microdamage Accumulation?
,”
Bone
,
27
(
1
), pp.
13
20
.
8.
Mashiba
,
T.
,
Hirano
,
T.
,
Turner
,
C. H.
,
Forwood
,
M. R.
,
Johnston
,
C. C.
, and
Burr
,
D. B.
,
2000
, “
Suppressed Bone Turnover by Bisphosphonates Increases Microdamage Accumulation and Reduces Some Biomechanical Properties in Dog Rib
,”
J. Bone Miner. Res.
,
15
(
4
), pp.
613
620
.
9.
Komatsubara
,
S.
,
Mori
,
S.
,
Mashiba
,
T.
,
Li
,
J.
,
Nonaka
,
K.
,
Kaji
,
Y.
,
Akiyama
,
T.
,
Miyamoto
,
K.
,
Cao
,
Y.
,
Kawanishi
,
J.
, and
Norimatsu
,
H.
,
2004
, “
Suppressed Bone Turnover by Long-Term Bisphosphonate Treatment Accumulates Microdamage but Maintains Intrinsic Material Properties in Cortical Bone of Dog Rib
,”
J. Bone Miner Res.
,
19
(
6
), pp.
999
1005
.
10.
Mashiba
,
T.
,
Turner
,
C. H.
,
Hirano
,
T.
,
Forwood
,
M. R.
,
Johnston
,
C. C.
, and
Burr
,
D. B.
,
2001
, “
Effects of Suppressed Bone Turnover by Bisphosphonates on Microdamage Accumulation and Biomechanical Properties in Clinically Relevant Skeletal Sites in Beagles
,”
Bone
,
28
(
5
), pp.
524
531
.
11.
Allen
,
M. R.
,
Iwata
,
K.
,
Phipps
,
R.
, and
Burr
,
D. B.
,
2006
, “
Alterations in Canine Vertebral Bone Turnover, Microdamage Accumulation, and Biomechanical Properties Following 1-Year Treatment With Clinical Treatment Doses of Risedronate or Alendronate
,”
Bone
,
39
(
4
), pp.
872
879
.
12.
Iwata
,
K.
,
Mashiba
,
T.
,
Hitora
,
T.
,
Yamagami
,
Y.
, and
Yamamoto
,
T.
,
2014
, “
A Large Amount of Microdamages in the Cortical Bone around Fracture Site in a Patient of Atypical Femoral Fracture After Long-Term Bisphosphonate Therapy
,”
Bone
,
64
, pp.
183
186
.
13.
Chapurlat
,
R. D.
,
Arlot
,
M.
,
Burt-Pichat
,
B.
,
Chavassieux
,
P.
,
Roux
,
J. P.
,
Portero-Muzy
,
N.
, and
Delmas
,
P. D.
,
2007
, “
Microcrack Frequency and Bone Remodeling in Postmenopausal Osteoporotic Women on Long-Term Bisphosphonates: A Bone Biopsy Study
,”
J. Bone Miner Res.
,
22
(
10
), pp.
1502
1509
.
14.
Forwood
,
M. R.
,
Burr
,
D. B.
,
Takano
,
Y.
,
Eastman
,
D. F.
,
Smith
,
P. N.
, and
Schwardt
,
J. D.
,
1995
, “
Risedronate Treatment Does Not Increase Microdamage in the Canine Femoral Neck
,”
Bone
,
16
(
6
), pp.
643
650
.
15.
Chapurlat
,
R. D.
, and
Delmas
,
P. D.
,
2009
, “
Bone Microdamage: A Clinical Perspective
,”
Osteoporosis Int.
,
20
(
8
), pp.
1299
1308
.
16.
Stepan
,
J. J.
,
Burr
,
D. B.
,
Pavo
,
I.
,
Sipos
,
A.
,
Michalska
,
D.
,
Li
,
J.
,
Fahrleitner-Pammer
,
A.
,
Petto
,
H.
,
Westmore
,
M.
,
Michalsky
,
D.
,
Sato
,
M.
, and
Dobnig
,
H.
,
2007
, “
Low Bone Mineral Density is Associated With Bone Microdamage Accumulation in Postmenopausal Women With Osteoporosis
,”
Bone
,
41
(
3
), pp.
378
385
.
17.
Roschger
,
P.
,
Rinnerthaler
,
S.
,
Yates
,
J.
,
Rodan
,
G. A.
,
Fratzl
,
P.
, and
Klaushofer
,
K.
,
2001
, “
Alendronate Increases Degree and Uniformity of Mineralization in Cancellous Bone and Decreases the Porosity in Cortical Bone of Osteoporotic Women
,”
Bone
,
29
(
2
), pp.
185
191
.
18.
Boskey
,
A. L.
,
2013
, “
Bone Composition: Relationship to Bone Fragility and Antiosteoporotic Drug Effects
,”
BoneKEy Rep.
,
2
, pp. 1–11.
19.
Bala
,
Y.
,
Depalle
,
B.
,
Farlay
,
D.
,
Douillard
,
T.
,
Meille
,
S.
,
Follet
,
H.
,
Chapurlat
,
R.
,
Chevalier
,
J.
, and
Boivin
,
G.
,
2012
, “
Bone Micromechanical Properties are Compromised During Long-Term Alendronate Therapy Independently of Mineralization
,”
J. Bone Miner. Res.
,
27
(
4
), pp.
825
834
.
20.
Monier-Faugere
,
M.-C.
,
Geng
,
Z.
,
Paschalis
,
E. P.
,
Qi
,
Q.
,
Arnala
,
I.
,
Bauss
,
F.
,
Boskey
,
A. L.
, and
Malluche
,
H. H.
,
1999
, “
Intermittent and Continuous Administration of the Bisphosphonate Ibandronate in Ovariohysterectomized Beagle Dogs: Effects on Bone Morphometry and Mineral Properties
,”
J. Bone Miner. Res.
,
14
(
10
), pp.
1768
1778
.
21.
Donnelly
,
E.
,
Meredith
,
D. S.
,
Nguyen
,
J. T.
,
Gladnick
,
B. P.
,
Rebolledo
,
B. J.
,
Shaffer
,
A. D.
,
Lorich
,
D. G.
,
Lane
,
J. M.
, and
Boskey
,
A. L.
,
2012
, “
Reduced Cortical Bone Compositional Heterogeneity With Bisphosphonate Treatment in Postmenopausal Women With Intertrochanteric and Subtrochanteric Fractures
,”
J. Bone Miner. Res.
,
27
(
3
), pp.
672
678
.
22.
Boskey
,
A. L.
,
Spevak
,
L.
, and
Weinstein
,
R. S.
,
2009
, “
Spectroscopic Markers of Bone Quality in Alendronate-Treated Postmenopausal Women
,”
Osteoporosis Int.
,
20
(
5
), pp.
793
800
.
23.
Gourion‐Arsiquaud
,
S.
,
Lukashova
,
L.
,
Power
,
J.
,
Loveridge
,
N.
,
Reeve
,
J.
, and
Boskey
,
A. L.
,
2013
, “
Fourier Transform Infrared Imaging of Femoral Neck Bone: Reduced Heterogeneity of Mineral‐to‐Matrix and Carbonate‐to‐Phosphate and More Variable Crystallinity in Treatment‐Naive Fracture Cases Compared With Fracture‐Free Controls
,”
J. Bone Miner. Res.
,
28
(
1
), pp.
150
161
.
24.
Ciarelli
,
T. E.
,
Tjhia
,
C.
,
Rao
,
D. S.
,
Qiu
,
S.
,
Parfitt
,
A. M.
, and
Fyhrie
,
D. P.
,
2009
, “
Trabecular Packet-Level Lamellar Density Patterns Differ by Fracture Status and Bone Formation Rate in White Females
,”
Bone
,
45
(
5
), pp.
903
908
.
25.
Tjhia
,
C. K.
,
Odvina
,
C. V.
,
Rao
,
D. S.
,
Stover
,
S. M.
,
Wang
,
X.
, and
Fyhrie
,
D. P.
,
2011
, “
Mechanical Property and Tissue Mineral Density Differences Among Severely Suppressed Bone Turnover (SSBT) Patients, Osteoporotic Patients, and Normal Subjects
,”
Bone
,
49
(
6
), pp.
1279
1289
.
26.
Acevedo
,
C.
,
Bale
,
H.
,
Gludovatz
,
B.
,
Wat
,
A.
,
Tang
,
S. Y.
,
Wang
,
M.
,
Busse
,
B.
,
Zimmermann
,
E. A.
,
Schaible
,
E.
, and
Allen
,
M. R.
,
2015
, “
Alendronate Treatment Alters Bone Tissues at Multiple Structural Levels in Healthy Canine Cortical Bone
,”
Bone
,
81
, pp.
352
363
.
27.
Allen
,
M. R.
,
Reinwald
,
S.
, and
Burr
,
D. B.
,
2008
, “
Alendronate Reduces Bone Toughness of Ribs Without Significantly Increasing Microdamage Accumulation in Dogs Following 3 Years of Daily Treatment
,”
Calcified Tissue Int.
,
82
(
5
), pp.
354
360
.
28.
Güerri-Fernández
,
R. C.
,
Nogués
,
X.
,
Quesada Gómez
,
J. M.
,
Torres del Pliego
,
E.
,
Puig
,
L.
,
García-Giralt
,
N.
,
Yoskovitz
,
G.
,
Mellibovsky
,
L.
,
Hansma
,
P. K.
, and
Díez-Pérez
,
A.
,
2013
, “
Microindentation for In Vivo Measurement of Bone Tissue Material Properties in Atypical Femoral Fracture Patients and Controls
,”
J. Bone Miner. Res.
,
28
(
1
), pp.
162
168
.
29.
Roschger
,
P.
,
Misof
,
B.
,
Paschalis
,
E.
,
Fratzl
,
P.
, and
Klaushofer
,
K.
,
2014
, “
Changes in the Degree of Mineralization With Osteoporosis and Its Treatment
,”
Curr. Osteoporosis Rep.
,
12
(
3
), pp.
338
350
.
30.
Tamminen
,
I. S.
,
Misof
,
B. M.
,
Roschger
,
P.
,
Mäyränpää
,
M. K.
,
Turunen
,
M. J.
,
Isaksson
,
H.
,
Kröger
,
H.
,
Mäkitie
,
O.
, and
Klaushofer
,
K.
,
2014
, “
Increased Heterogeneity of Bone Matrix Mineralization in Pediatric Patients Prone to Fractures: A Biopsy Study
,”
J. Bone Miner. Res.
,
29
(
5
), pp.
1110
1117
.
31.
Bousson
,
V.
,
Bergot
,
C.
,
Wu
,
Y.
,
Jolivet
,
E.
,
Zhou
,
L. Q.
, and
Laredo
,
J.-D.
,
2011
, “
Greater Tissue Mineralization Heterogeneity in Femoral Neck Cortex From Hip-Fractured Females Than Controls. A Microradiographic Study
,”
Bone
,
48
(
6
), pp.
1252
1259
.
32.
Boskey
,
A. L.
,
Donnelly
,
E.
,
Boskey
,
E.
,
Spevak
,
L.
,
Ma
,
Y.
,
Zhang
,
W.
,
Lappe
,
J.
, and
Recker
,
R. R.
,
2016
, “
Examining the Relationships Between Bone Tissue Composition, Compositional Heterogeneity, and Fragility Fracture: A Matched Case-Controlled FTIRI Study
,”
J. Bone Miner. Res.
,
31
(
5
), pp.
1070
1081
.
33.
Gao
,
X.
,
Li
,
S.
,
Adel-Wahab
,
A.
, and
Silberschmidt
,
V.
,
2013
, “
Effect of Random Microstructure on Crack Propagation in Cortical Bone Tissue Under Dynamic Loading
,”
J. Phys.: Conf. Ser.
,
451
(
1
), p.
012033
.
34.
Budyn
,
E.
, and
Hoc
,
T.
,
2010
, “
Analysis of Micro Fracture in Human Haversian Cortical Bone Under Transverse Tension Using Extended Physical Imaging
,”
Int. J. Numer. Methods Eng.
,
82
(
8
), pp.
940
965
.
35.
Li
,
S.
,
Abdel-Wahab
,
A.
,
Demirci
,
E.
, and
Silberschmidt
,
V.
,
2013
, “
Fracture Process in Cortical Bone: X-FEM Analysis of Microstructured Models
,”
Int. J. Fract.
,
184
(1–2), pp. 43–55.
36.
Mischinski
,
S.
, and
Ural
,
A.
,
2011
, “
Finite Element Modeling of Microcrack Growth in Cortical Bone
,”
ASME J. Appl. Mech.
,
78
(
4
), p.
041016
.
37.
Mischinski
,
S.
, and
Ural
,
A.
,
2013
, “
Interaction of Microstructure and Microcrack Growth in Cortical Bone: A Finite Element Study
,”
Comput. Methods Biomech. Biomed. Eng.
,
16
(
1
), pp.
81
94
.
38.
Budyn
,
E.
,
Hoc
,
T.
, and
Jonvaux
,
J.
,
2008
, “
Fracture Strength Assessment and Aging Signs Detection in Human Cortical Bone Using an X-FEM Multiple Scale Approach
,”
Comput. Mech.
,
42
(
4
), pp.
579
591
.
39.
Jonvaux
,
J.
,
Hoc
,
T.
, and
Budyn
,
É.
,
2012
, “
Analysis of Micro Fracture in Human Haversian Cortical Bone Under Compression
,”
Int. J. Numer. Methods Biomed. Eng.
,
28
(
9
), pp.
974
998
.
40.
Hambli
,
R.
,
2011
, “
Apparent Damage Accumulation in Cancellous Bone Using Neural Networks
,”
J. Mech. Behav. Biomed. Mater.
,
4
(
6
), pp.
868
878
.
41.
Hambli
,
R.
,
2011
, “
Multiscale Prediction of Crack Density and Crack Length Accumulation in Trabecular Bone Based on Neural Networks and Finite Element Simulation
,”
Int. J. Numer. Methods Biomed. Eng.
,
27
(
4
), pp.
461
475
.
42.
Hambli
,
R.
,
Lespessailles
,
E.
, and
Benhamou
,
C. L.
,
2013
, “
Integrated Remodeling-to-Fracture Finite Element Model of Human Proximal Femur Behavior
,”
J. Mech. Behav. Biomed. Mater.
,
17
, pp.
89
106
.
43.
Najafi
,
A. R.
,
Arshi
,
A. R.
,
Eslami
,
M. R.
,
Fariborz
,
S.
, and
Moeinzadeh
,
M. H.
,
2007
, “
Micromechanics Fracture in Osteonal Cortical Bone: A Study of the Interactions Between Microcrack Propagation, Microstructure and the Material Properties
,”
J. Biomech.
,
40
(
12
), pp.
2788
2795
.
44.
Raeisi Najafi
,
A.
,
Arshi
,
A. R.
,
Saffar
,
K. P.
,
Eslami
,
M. R.
,
Fariborz
,
S.
, and
Moeinzadeh
,
M. H.
,
2009
, “
A Fiber-Ceramic Matrix Composite Material Model for Osteonal Cortical Bone Fracture Micromechanics: Solution of Arbitrary Microcracks Interaction
,”
J. Mech. Behav. Biomed. Mater.
,
2
(
3
), pp.
217
223
.
45.
Tai
,
K.
,
Dao
,
M.
,
Suresh
,
S.
,
Palazoglu
,
A.
, and
Ortiz
,
C.
,
2007
, “
Nanoscale Heterogeneity Promotes Energy Dissipation in Bone
,”
Nat. Mater.
,
6
(
6
), pp.
454
462
.
46.
Yao
,
H. M.
,
Dao
,
M.
,
Carnelli
,
D.
,
Tai
,
K. S.
, and
Ortiz
,
C.
,
2011
, “
Size-Dependent Heterogeneity Benefits the Mechanical Performance of Bone
,”
J. Mech. Phys. Solids
,
59
(
1
), pp.
64
74
.
47.
Demirtas
,
A.
,
Curran
,
E.
, and
Ural
,
A.
,
2016
, “
Assessment of the Effect of Reduced Compositional Heterogeneity on Fracture Resistance of Human Cortical Bone Using Finite Element Modeling
,”
Bone
,
91
, pp.
92
101
.
48.
Rho
,
J. Y.
,
Tsui
,
T. Y.
, and
Pharr
,
G. M.
,
1997
, “
Elastic Properties of Human Cortical and Trabecular Lamellar Bone Measured by Nanoindentation
,”
Biomaterials
,
18
(
20
), pp.
1325
1330
.
49.
McCalden
,
R. W.
,
McGeough
,
J. A.
,
Barker
,
M. B.
, and
Court-Brown
,
C. M.
,
1993
, “
Age-Related Changes in the Tensile Properties of Cortical Bone. The Relative Importance of Changes in Porosity, Mineralization, and Microstructure
,”
J. Bone Jt. Surg. Am.
,
75
(
8
), pp.
1193
1205
.
50.
Reilly
,
D. T.
, and
Burstein
,
A. H.
,
1975
, “
The Elastic and Ultimate Properties of Compact Bone Tissue
,”
J. Biomech.
,
8
(
6
), pp.
393
405
.
51.
Zioupos
,
P.
, and
Currey
,
J. D.
,
1998
, “
Changes in the Stiffness, Strength, and Toughness of Human Cortical Bone With Age
,”
Bone
,
22
(
1
), pp.
57
66
.
52.
Dong
,
X. N.
,
Zhang
,
X.
, and
Guo
,
X. E.
,
2005
, “
Interfacial Strength of Cement Lines in Human Cortical Bone
,”
Mech. Chem. Biosyst.
,
2
(
2
), pp.
63
68
.
53.
Koester
,
K. J.
,
Ager
,
J. W.
, III
, and
Ritchie
,
R. O.
,
2008
, “
The True Toughness of Human Cortical Bone Measured With Realistically Short Cracks
,”
Nat. Mater.
,
7
(
8
), pp.
672
677
.
54.
Schaffler
,
M. B.
,
Choi
,
K.
, and
Milgrom
,
C.
,
1995
, “
Aging and Matrix Microdamage Accumulation in Human Compact Bone
,”
Bone
,
17
(
6
), pp.
521
525
.
55.
Allen
,
M. R.
, and
Burr
,
D. B.
,
2007
, “
Mineralization, Microdamage, and Matrix: How Bisphosphonates Influence Material Properties of Bone
,”
BoneKEy-Osteovision
,
4
(
2
), pp.
49
60
.
56.
Akkus
,
O.
, and
Rimnac
,
C. M.
,
2001
, “
Cortical Bone Tissue Resists Fatigue Fracture by Deceleration and Arrest of Microcrack Growth
,”
J. Biomech.
,
34
(
6
), pp.
757
764
.
57.
O'Brien
,
F. J.
,
Taylor
,
D.
, and
Lee
,
T. C.
,
2003
, “
Microcrack Accumulation at Different Intervals During Fatigue Testing of Compact Bone
,”
J. Biomech.
,
36
(
7
), pp.
973
980
.
58.
Moës
,
N.
, and
Belytschko
,
T.
,
2002
, “
Extended Finite Element Method for Cohesive Crack Growth
,”
Eng. Fract. Mech.
,
69
(
7
), pp.
813
833
.
59.
Ortiz
,
M.
, and
Pandolfi
,
A.
,
1999
, “
Finite-Deformation Irreversible Cohesive Elements for Three-Dimensional Crack-Propagation Analysis
,”
Int. J. Numer. Methods Eng.
,
44
(
9
), pp.
1267
1282
.
60.
Ural
,
A.
, and
Vashishth
,
D.
,
2006
, “
Cohesive Finite Element Modeling of Age-Related Toughness Loss in Human Cortical Bone
,”
J. Biomech.
,
39
(
16
), pp.
2974
2982
.
61.
Camanho
,
P. P.
,
Davila
,
C. G.
, and
de Moura
,
M. F.
,
2003
, “
Numerical Simulation of Mixed-Mode Progressive Delamination in Composite Materials
,”
J. Compos. Mater.
,
37
(
16
), pp.
1415
1438
.
62.
Vashishth
,
D.
,
Behiri
,
J. C.
, and
Bonfield
,
W.
,
1997
, “
Crack Growth Resistance in Cortical Bone: Concept of Microcrack Toughening
,”
J. Biomech.
,
30
(
8
), pp.
763
769
.
63.
Ettinger
,
B.
,
Burr
,
D. B.
, and
Ritchie
,
R. O.
,
2013
, “
Proposed Pathogenesis for Atypical Femoral Fractures: Lessons From Materials Research
,”
Bone
,
55
(
2
), pp.
495
500
.
64.
Yeni
,
Y. N.
,
Brown
,
C. U.
,
Wang
,
Z.
, and
Norman
,
T. L.
,
1997
, “
The Influence of Bone Morphology on Fracture Toughness of the Human Femur and Tibia
,”
Bone
,
21
(
5
), pp.
453
459
.
65.
ASTM
,
2012
, “
ASTM Standard E399 Standard Test Method for Linear-Elastic Plane-Strain Fracture Toughness KIc of Metallic Materials
,” ASTM International, West Conshohocken, PA, Standard No. ASTM E399.
66.
Mast
,
P.
,
Nash
,
G.
,
Michopoulos
,
J.
,
Thomas
,
R.
,
Badaliance
,
R.
, and
Wolock
,
I.
,
1995
, “
Characterization of Strain-Induced Damage in Composites Based on the Dissipated Energy Density—Part I: Basic Scheme and Formulation
,”
Theor. Appl. Fract. Mech.
,
22
(
2
), pp.
71
96
.
67.
Allen
,
M. R.
, and
Burr
,
D. B.
,
2011
, “
Bisphosphonate Effects on Bone Turnover, Microdamage, and Mechanical Properties: What We Think We Know and What We Know That We Don't Know
,”
Bone
,
49
(
1
), pp.
56
65
.