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.
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 [1–6]. 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 . 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 . 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 [17–24]. In addition, reduced heterogeneity in elastic modulus , a decreasing trend in energy to fracture , and reduced fracture toughness  were reported in experiments performed on long-term BP-treated bone. A deterioration in microindentation properties was also observed in AFF patients . However, some studies found no relation between BP use and tissue compositional heterogeneity [29–32]. 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) [33–39]. 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 . 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,37–39]. Our previous study evaluated the influence of heterogeneity of spatial distribution of material properties on fracture resistance of human cortical bone . 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
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.
In this study, isotropic elastic properties were used for the CT specimen and the detailed microstructure region (Table 2) . 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 . 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) [48–51]. 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  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.
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 . 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  in the osteons and the interstitial bone and cohesive interface elements  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 (Tn/σnc)2 + (Ts/σsc)2 + (Tt/σtc)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 . 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 . 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 , whereas the initial slope is not required for the XFEM formulation. The details of the cohesive modeling technique for bone were outlined in Refs. , , and .
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 . 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  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.
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).
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 . Limited microcracking provides a protective mechanism against crack growth by promoting energy dissipation and activating other toughening mechanisms such as crack deflection . 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 . 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  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 . The same computational study also demonstrated a reduction in fracture strength with increasing porosity . 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  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 . 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 . 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.
National Science Foundation (Award No. CMMI-1434412).
The computational resources for this study were supported by XSEDE (Award No. MSS-160001).