Crack propagation in cortical bone is affected by the characteristics of the cement line: a parameter study using an XFEM interface damage model

Bulk properties of cortical bone have been well characterized experimentally, and potent toughening mechanisms, e.g., crack deflections, have been identified at the microscale. However, it is currently difficult to experimentally measure local damage properties and isolate their effect on the tissue fracture resistance. Instead, computer models can be used to analyze the impact of local characteristics and structures, but material parameters required in computer models are not well established. The aim of this study was therefore to identify the material parameters that are important for crack propagation in cortical bone and to elucidate what parameters need to be better defined experimentally. A comprehensive material parameter study was performed using an XFEM interface damage model in 2D to simulate crack propagation around an osteon at the microscale. The importance of 14 factors (material parameters) on four different outcome criteria (maximum force, fracture energy, crack length and crack trajectory) was evaluated using ANOVA for three different osteon orientations. The results identified factors related to the cement line to influence the crack propagation, where the interface strength was important for the ability to deflect cracks. Crack deflection was also favored by low interface stiffness. However, the cement line properties are not well determined experimentally and need to be better characterized. The matrix and osteon stiffness had no or low impact on the crack pattern. Furthermore, the results illustrated how reduced matrix toughness promoted crack penetration of the cement line. This effect is highly relevant for the understanding of the influence of aging on crack propagation and fracture resistance in cortical bone. Electronic supplementary material The online version of this article (10.1007/s10237-019-01142-4) contains supplementary material, which is available to authorized users.


Introduction
Bone tissue is a complex material that adapts during the course of life, and structural and material changes occur at all hierarchical length scales of bone Zimmermann et al. 2015). Old bone is characterized by reduced fracture toughness with compromised ability to slow down or deflect propagating cracks (Chan et al. 2009;Koester et al. 2011;Nalla et al. 2004Nalla et al. , 2006Zimmermann et al. 2011), increased porosity (Cooper et al. 2007;Mirzaali et al. 2016;Stein et al. 1999) and reduced plasticity at the nanoscale (Zimmermann et al. 2011). One detrimental effect of aging is the increased risk of fractures in patients with osteoporosis (Cummings and Melton 2002;Hernlund et al. 2013). Yet, it is very challenging to experimentally measure damage parameters locally (Kruzic et al. 2009) and to distinguish how the fracture resistance is affected by local material or structural alterations. Instead, existing methods measure average bulk behavior of cortical bone tissue (Koester et al. 2008(Koester et al. , 2011Nalla et al. 2005). Furthermore, crack deflections along cement lines are identified as important toughening mechanisms in cortical bone (Chan et al. 2009;Koester et al. 2008Koester et al. , 2011Nalla et al. 2006;Zimmermann et al. 2009Zimmermann et al. , 2010, but still there is no consensus regarding the mechanical properties of the interface (Burr et al. 1988;Milovanovic et al. 2018; Montalbano and Feng 2011;Skedros et al. 2005). At this point, computer models can improve the understanding of how microstructure contributes to tissue fracture toughness.
Traditionally, there have been two types of theoretical models for composite materials for analyzing crack propagation at an interface: strength-based models, as introduced by Cook and Gordon (1964), where the difference in strength between interface and substrate (i.e., osteon in the case of cortical bone) dictates whether the crack will be deflected by the interface, and energy-based models where the crack path is determined by the maximum energy release when comparing different possible crack paths (He and Hutchinson 1989). These ideas are combined in cohesive damage models that account for both strength and energy when modeling crack propagation. Parmigiani and Thouless (2006) used cohesive elements to show that both strength and energy can be important for the crack trajectory at an interface. Mischinski and Ural (2011) modeled crack propagation in cortical bone by outlining two possible crack paths with cohesive elements: one penetrating the osteon and the other deflecting along the cement line. They concluded that low cement line strength was the most critical factor for promoting crack deflection (Mischinski and Ural 2011).
Another option for modeling crack propagation is the extended finite element method (XFEM) which has the benefit, that is, it does not require a predefined crack path (Belytschko and Black 1999;Melenk and Babuska 1996). A handful of XFEM models have been used to model crack propagation in 2D at the microscale in cortical bone, and the maximum principal strain (MAXPE) criterion has been commonly used to model damage initiation with cohesive cracks (Abdel-Wahab et al. 2012;Budyn and Hoc 2007;Idkaidek and Jasiuk 2017;Li et al. 2013;Vergani et al. 2014). However, the MAXPE criterion cannot be used to model realistic crack trajectories around osteons, as it always predicts crack penetration of the cement line interface (Gustafsson et al. 2018a). To handle this problem, we recently proposed an XFEM interface damage model that is able to capture both crack deflections along cement line interfaces and crack propagation through osteons (Gustafsson et al. 2018a). Additionally, Marco et al. (2018a) proposed heterogeneous maximum tangential stress criteria to model crack propagation along cement lines. However, the lack of experimental data is reflected in the computational models, where material parameters governing crack propagation are loosely defined and based on assumptions. As an effect, one set of damage parameters (e.g., parameters used in Abdel-Wahab et al. 2012) has been used in almost all XFEM models simulating crack propagation in cortical bone. Furthermore, there is a large uncertainty regarding the cement line stiffness where values used in recent modeling studies range from more than 160 times lower (Marco et al. 2018a) to 20% higher (Gustafsson et al. 2018a) than the matrix stiffness. Hence, there is a need to thoroughly evaluate the parameters used to model crack propagation in cortical bone.
The aim of this study was therefore to identify the most important material parameters for crack propagation in cortical bone by evaluating their effect on the maximum force, fracture energy, crack length and crack trajectory. By identifying important material parameters numerically, this study will reveal what parameters need to be studied or measured with higher accuracy experimentally. As experimental studies show that the crack pattern at the microscale affects the fracture resistance of cortical bone tissue, e.g., (Chan et al. 2009;Koester et al. 2008), this study focused on the interactions between a propagating crack and an osteon. Crack propagation was simulated in 2D using simplified microstructural geometries depicting an osteon surrounded by a cement line interface in different orientations (Gustafsson et al. 2018a).

Methods
The study design involved a parameter study in three steps, where important parameters were identified in each step and selected for further analysis. An overview of the procedure is given in Fig. 1. The 2D XFEM models introduced in Gustafsson et al. (2018a) were used to simulate crack propagation around an osteon. An initial screening analysis including all 14 material parameters from Gustafsson et al. (2018a) was performed to identify the most influential material parameters through analysis of variance (ANOVA). A subset of 7 parameters were then analyzed using a response surface design, where nonlinear dependencies and interactions between parameters were evaluated. The final step consisted in mapping the effect of material toughness, critical interface strains and cement line stiffness on the crack pattern.

Cortical bone damage model
A damage model was recently developed for cortical bone to simulate crack propagation around an osteon at the microscale using XFEM (Gustafsson et al. 2018a). Three models were created representing the cortical bone microstructure with different orientations, where each model comprised one osteon embedded in an interstitial matrix and surrounded by a cement line interface. The Haversian canals were neglected in this study and filled with osteon material, as they were previously shown to induce unrealistic effects in a 2D model (Gustafsson et al. 2018a). All microstructural features, i.e., interstitial matrix, osteons and cement lines, were modeled as linear isotropic elastic materials. The models are referred to as longitudinal, radial and transversal models, according to the orientation of the osteon (Fig. 2). Plane stress 4-node bilinear elements with reduced integration (CPS4R) were used and longitudinal, radial and transversal models contained 27,116, 11,915 and 18,555 elements, respectively. The initial cracks were inserted in the left edge of all models. Tensile tests until failure were simulated using displacement-controlled loading in a quasi-static analysis. Boundary conditions and model dimensions are illustrated in Fig. 2. More detailed information about the models and the mesh design, including a convergence study, can be found in Gustafsson et al. (2018a).
Crack propagation was modeled using the XFEM framework implemented in Abaqus Standard (v2017, Dassault Systemes) with the cohesive segments approach. Damage initiation in matrix and osteon was modeled using the maximum principal strain criterion (MAXPE). The fracture criterion f is defined as where max is the maximum principal strain and 0 max is the critical damage initiation strain. Damage is initiated when f MAXPE > 1 and the crack normal is given by the maximum principal strain orientation.
For the cement line, the interface damage model proposed in Gustafsson et al. (2018a) was used. In this model, the MAXPE criterion is used to model crack propagation through the cement line and the quadratic nominal strain criterion (QUADE) is used to model interface damage that causes crack deflection along the interface. Both criteria were implemented in the user-defined damage initiation subroutine UDMGINI, and damage was initiated when  In all models, the osteon was embedded in a square matrix and surrounded by the cement line interface separating the osteon from the interstitial matrix. The initial crack length was 10% of the side of each model. In the longitudinal and transversal model, the osteon was modeled as a rectangle rotated 10° with two half circles as endpoints, and in the radial model the osteon was modeled as a circle. The model dimensions were L =1 mm, h = 650 µm, d = 150 µm, t = 5 µm, and the thickness of the models was 100 µm max f QUADE , f MAXPE > 1 . The QUADE fracture criterion was defined as where n is the normal strain, s is the shear strain and 0 n and 0 s denote the critical interface strains in the cement line. A vector oriented perpendicular to the cement line interface was used as crack normal for the QUADE criterion.
For all damage criteria, the damage evolution was assumed to be mode independent and the degradation of the cohesive crack was modeled with an energy-based evolution law. A linear softening behavior was assumed for the traction-separation response inside the cohesive crack, and the strain energy release rate G was determined as the area under the traction-separation curve. The strain energy release rate (material toughness) corresponds to the fracture energy needed to completely open the crack.

Literature study
A literature study was performed to estimate the range of values for each model parameter (Table 1). Values from experimental and numerical studies were separated to illustrate commonly used model parameters in relation to available experimental data. In cases where quantitative experimental data were not available, relative numbers have been reported. No experimental data exist comparing critical strains in the cement line to critical strains for damage initiation in matrix and osteons; instead, data were reported in terms of interfacial shear strength from osteon push-out tests.

Parameter studies using the design of experiments approach
The first two steps of the parameter study ( Fig. 1a, b) are based on a design of experiments approach using fractional factorial designs which enables the parameters that have the largest influence on the outcome variables to be identified. The method lends itself well to computational studies in biomechanics (Dar et al. 2002;Isaksson et al. 2008Isaksson et al. , 2009). The parameter set introduced in Gustafsson et al. (2018a) was used as a baseline in this study (Table 3). In total, 14 material parameters (factors) were identified and analyzed: Young's modulus, Poisson's ratio, the critical damage initiation strain and the strain energy release rate for matrix, osteon and cement line materials, and the critical normal and shear strains in the cement line interface. For the critical interface strains, the average value between crack penetration and deflection as reported in Gustafsson et al. (2018a) f QUADE = n 0 n 2 + s 0 s 2 for each osteon orientation was used as baseline value. Each factor was assigned a parameter space that was ± 20% of the corresponding baseline value (Table 3). The importance of each factor for different outcome parameters was determined using analysis of variance (ANOVA). The evaluated outcome parameters were the maximum force (peak value in load curve), fracture energy, crack length (from the coordinates of all crack segments) and crack score. The crack score was used to categorize the crack pattern by assigning a number between 1 and 5 to the crack pattern according to Table 2.
The importance of the factors was estimated based on the assumption that important factors give a large contribution to the variance. The analysis followed the same rational as presented in Isaksson et al. (2008) where the following parameters needed to be calculated for the analysis of variance. The total sum of squares of the deviation about the mean (SS T ) was calculated as where N was the number of treatment conditions, y i the outcome parameter for the ith treatment conditions and ȳ the mean of all y i . The sum of the squares of deviation about the mean for each factor (SS F ) was calculated as where n was the number of levels, N F,i was the number of treatment conditions at each level of each factor and ȳ F,i was the mean outcome parameter at each level of each factor. The percentage of the total sum of squares for each factor ( %TSS ) represented the contribution of each factor to the variance and was considered a measure of the importance of each factor.

Screening design
An initial two-level screening experiment was performed on longitudinal, radial and transversal models to determine the most important factors (X screen in Table 3, Fig. 1a). A Resolution IV Two-Level fractional factorial design was used, where the main effects were assumed to be clear of two-factor interactions (Montgomery 2005). The Resolution IV array contained 14 factors (material parameters) and 32 treatment conditions (simulations with different combinations of factors) and was generated with the software JMP (see Supplementary Table 1).

Response surface design
Next, a response surface design was used to evaluate potential interactions and to further analyze important parameters identified in the screening experiment for the longitudinal and radial models (Fig. 1b). A subset of factors was identified based on their importance for fracture energy, crack length and crack score in the screening study, as these were assumed to be of most importance for the global fracture resistance of bone tissue. A Box-Behnken design was used with 7 factors (X BB in Table 3), three levels and 62 treatment conditions. The three levels were created by adding the baseline value as a middle value between the low and high levels ( Table 3). The Box-Behnken surface design array was generated with the software JMP (see Supplementary Table 2). The critical interface normal and shear strains were in this case assumed to be equal. The mean value is given in case mean ± SD was reported in the cited article. Where no quantitative data were available, relative measures were reported. Values for damage initiation were considered in experimental data reporting stress intensity or fracture energy (strain energy release rate). The stiffness for the matrix was measured under wet conditions. Symbols used in the table: ~ approximate value interpreted from figures; K G calculated from a reported stress intensity factor K, where G = K 2 /E, assuming E = 20 GPa (Koester et al. 2008

Mapping of parameter effect
The final step was to map the effect of two (Fig. 1c) and three (Fig. 1d) important parameters on the crack pattern while keeping all other parameters at baseline values. The analyzed parameter values spanned a wider range compared to values used in the screening and response surface designs in order to capture the full model behavior as described in Table 2. First, the effect of the material toughness (strain energy release rate) and the critical interface strain was evaluated (Fig. 1c). To reduce the number of varied parameters, both the material toughness parameters and the critical interface parameters were grouped and varied simultaneously (i.e., G = G mat = G ost = G cl and 0 cl = 0 n,cl = 0 s,cl ). The material toughness was varied from 0.05 to 0.4 kJ/m 2 and the critical interface strains from 1e−5 to 0.0045. Simulations were run using all three osteon orientations. Next, the effect of the cement line stiffness was also included (Fig. 1d) and evaluated at two levels: equal to the matrix stiffness (E cl = E mat = 15 GPa) and 20% lower than the osteon stiffness (E cl = 0.8·E ost = 9.6 GPa). Simulations were run using the longitudinal model, and the material toughness was varied from 0.05 to 0.4 kJ/m 2 and the critical interface strains from 0.0005 to 0.003.

Screening design
ANOVA was used to identify the most important factors for four outcome criteria (maximum force, fracture energy, crack length and crack score) in three different model geometries representing longitudinal, radial and transversal osteons (Table 4, contributions higher than 5% are shown in bold font). E mat and 0 max,mat were important only for the maximum force outcome and therefore not considered for the response surface design. The fracture energy outcome was dominated by the strain energy release rate G of the different materials. Differences were observed when comparing the fracture energy outcome for the different osteon orientations, where G mat was most important for the longitudinal models and G cl for the transversal models. For the radial model, also G ost was found to be important. The crack length and crack score outcomes showed similar results where E cl , 0 max,cl , 0 s,cl and 0 n,cl were important factors. However, 0 s,cl scored high for longitudinal and radial models, whereas 0 n,cl was most important for the transversal model. E ost also influenced the crack length and crack score in the radial and transversal models. Poisson's ratios for the different materials were found to have no or very low importance for all evaluated outcomes. Factors included in Table 4 could explain 57-96% of the variation in the evaluated outcome criteria. In total, seven factors were identified as important for the longitudinal and radial models (E ost , E cl , 0 max,cl , 0 s,cl , G mat , G ost and G cl ) and included in the response surface design (X BB in Table 3).

Response surface design
Results from the response surface design experiment (Table 5) confirmed the trends observed in the screening experiment (Table 4). The interactions scored low in general, and significant interactions involved factors that were already identified as important. The exception was the G cl · G cl interaction that explained over 7% of the variation for crack length in the radial model, while G cl alone was less important. Factors and interactions included in Table 5 could explain 85-92% of the variation in all outcome criteria.
The prediction profiles of the seven factors used in the Box-Behnken surface design (Fig. 3) showed minor nonlinearity in several factors. Most importantly, it showed nonlinear behavior for factors associated with the cement line, where no effect was seen for high levels of the factors and sudden high effects were found when the factors were lower than a threshold value. Most obvious was the effect of the cement line on the crack length, where high damage initiation strain and low critical interface strains promoted crack deflection. The material toughness parameters showed modest effects at the evaluated levels.

Mapping of parameter effect
Both the material toughness and the critical interface strains affected the crack trajectory (Fig. 4). In general, low material toughness and high interface strength promoted crack penetration of the osteon, while high toughness and low interface strength promoted crack deflection in the cement line. The different crack patterns (crack scores) were found within continuous regions in Fig. 4, with nonlinear lines marking the borders between the different regions. Lower interface strengths were required for crack deflection in the longitudinal models compared to the transversal. Asymptotic trends reached high toughness values (G > 0.2 kJ/m 2 ) for all osteon orientations.
The effect of the cement line stiffness was seen as a shift in the crack trajectory (Fig. 5), where lower cement line stiffness promoted crack deflection at higher interface strengths. With low cement line stiffness, small deflections were more common (crack score 2, black region in Fig. 5), compared to models with higher cement line stiffness where the crack directly penetrated the interface (crack score 1, red region in Fig. 5).

Discussion
Local damage properties in cortical bone are difficult to measure experimentally, and hence, material parameters used in computer models are based on very general assumptions and not well established. The aim of this study is to identify important factors for crack propagation in cortical bone and to identify parameters that need further quantification. A comprehensive parameter study was performed using the interface damage model introduced in Gustafsson et al. (2018a), and the importance of 14 material parameters was evaluated in a screening design. This was followed by a response surface design including 7 factors and finally a more in-depth parameter study mapping the effect of a few important parameters. The results emphasized the importance of factors related to the cement line, which are all loosely defined in the literature.

Parameter studies using the design of experiments approach
In the initial screening experiment, the maximum force outcome (Table 4) showed that the importance of the osteon and cement line properties was low when the osteon was oriented parallel to the applied load (longitudinal model) and high when oriented perpendicular to the load (transversal model). Similarly, the damage initiation strain in the matrix was a dominant factor in the longitudinal and radial   Table 3 Material parameters used in the screening experiment with two levels (low and high) and 14 factors (X screen ) and the Box-Behnken surface design with three levels (low, baseline and high) and 7 factors (X BB ) Critical interface strains are specified for each osteon direction: L longitudinal, R radial, T transversal. Baseline values are based on values from the literature (Table 1) and introduced in a previous study (Gustafsson et   models, while instead all critical strains showed a modest importance (%TSS ~ 5%) in the transverse models. When analyzing the fracture energy and the crack trajectory outcomes, similar trends were seen in the screening study and the Box-Behnken surface design experiment (Tables 4, 5).
In the screening study, the material toughness parameters (G) were most important for the fracture energy, while in the three-level surface design the cement line parameters also scored high. The Young's moduli were not important for the fracture toughness (Table 4). This is in line with experiments, showing that the microstructure (osteon density) and tissue porosity are more important for the fracture toughness than the material heterogeneity (stiffness) in cortical bone (Granke et al. 2016). Tissue porosity has been shown to correlate both to crack initiation toughness and to total dissipated energy during fracture at the tissue scale (Granke et al. 2015(Granke et al. , 2016. These effects were captured by the critical damage initiation strains and the strain energy release rates for the different materials in the cohesive damage model even though the porosity was not explicitly modeled in this study. At lower length scales, periodic variation in elastic modulus has been shown to increase the fracture toughness   (Fratzl et al. 2007). However, as the lamellar bone structure was not included in our models, this effect was also incorporated into the material toughness parameters. Surprisingly, the factors describing the material toughness were not important for the crack length or crack score outcome criteria in the screening and response surface studies. Still, there seemed to be a correlation between crack length and fracture energy (Fig. 6) and this trend was confirmed when analyzing a wide range of material toughness values (Fig. 4). Again, all factors related to the cement line were important for the crack length and the crack score outcomes (Tables 4,  5) and the results were similar when comparing the two outcome criteria. The competition between the crack driving force driving the crack through the osteon and the interfaces trying to deflect the crack, as discussed in Zimmermann et al. (2009Zimmermann et al. ( , 2010, was illustrated by the prediction profiles (Fig. 3), where high damage initiation strain (high 0 max,cl ) and weak interfaces (low 0 cl ) were both beneficial for crack deflection and increased fracture energy.
As many model parameters lack clear experimental quantification, all factors were varied ± 20% around the baseline value. In this way, all factors were given similar importance, which is beneficial when the range of possible values is not well established. The drawback is that the chosen approach is sensitive to the choice of baseline values and it is possible that another combination of material parameters or another parameter space could display another behavior. As there is no consensus regarding the cement line stiffness (Burr et al. 1988;Milovanovic et al. 2018;Montalbano and Feng 2011;Skedros et al. 2005), line (E cl = 0.8·E ost = 9.6 GPa), the middle plot to a cement line with the same stiffness as the matrix (E cl = E mat = 15 GPa) and the lower plot to a stiff cement line (E cl = 1.2·E mat = 18 GPa, same as shown in Fig. 3). Each point in the graph corresponds to a simulation, and colored regions are drawn to visualize the region of each crack score, as shown in the right the effect of the baseline value for the cement line stiffness was evaluated by repeating the initial screening study for the longitudinal model with E cl = E mat = 15 GPa. This had no effect on what parameters that were identified as important (data not shown). Furthermore, as both screening and surface design analyses predicted a wide range of different crack paths (crack score 1-4), it indicates that a sufficiently broad and relevant parameter space for capturing the model behavior was chosen. In our previous study (Gustafsson et al. 2018a), the critical normal and shear strains in the interface were assumed to be equal based on that the weaker of the two critical strains would determine the threshold value for a certain crack trajectory. The effect of this assumption was tested in the screening experiment, which showed that the normal and shear strains were important under different loading conditions (Table 4). The critical normal strain 0 n,cl was important in the transversal model, when the load was applied perpendicular to the cement line interface. The critical shear strain 0 s,cl was important when the load was applied parallel to the interface, as in the longitudinal model. In the remaining simulations performed in this study, the critical normal and shear strains were varied simultaneously (i.e., 0 cl = 0 n,cl = 0 s,cl ), as only one parameter at a time had an effect when separating them for the evaluated model geometries (Table 4).

Mapping of parameter effect
When evaluating the effect of the critical interface strains and the material toughness on the crack trajectory, we found that the critical interface strains were highly important for the crack score, as crack deflections could always be achieved with a sufficiently low interface strength, independent of the material toughness (Fig. 4). On the other hand, material toughness alone could not be used to enforce crack deflection. This is in agreement with what has been shown using cohesive elements (Mischinski and Ural 2011;Parmigiani and Thouless 2006). However, the material toughness had an effect on the crack trajectory at low levels (G ≤ 0.2 kJ/ m 2 ) and simulations with higher toughness values (up to G = 1 kJ/m 2 ) were run to confirm that the crack trajectory was not changed (data not shown). The asymptotic behavior for high toughness agreed with the study by Parmigiani and Thouless (2006). Furthermore, the effect of the toughness in each microstructural feature was evaluated by varying the parameters G mat , G ost and G cl separately and keeping the others fixed at 0.2 kJ/m 2 (data not shown). Varying only the matrix toughness gave the same results as varying all toughness parameters simultaneously (as shown in Fig. 4), while no effect was seen on the crack score when varying the material toughness in the osteon or cement line. This illustrates the importance of the matrix toughness and the mechanical state inside the crack when encountering an interface, where cracks in materials with low toughness have a higher risk of penetrating an interface. This aspect was not evaluated by the parameter studies using cohesive elements. Instead, Mischinski and Ural (2011) kept the fracture toughness in the matrix fixed at high values (fracture toughness for normal mode was G nc = 1.16 kJ/m 2 and shear mode G sc = 2.97 kJ/ m 2 ) and Parmigiani and Thouless (2006) modeled a crack already impinging on an interface, not including the effect of the matrix fracture toughness at all. Interestingly, the matrix toughness had large influence in our model at levels corresponding to damage initiation in cortical bone where higher values (G ~ 0.2 kJ/m 2 ) correspond to young bone and lower values (G ~ 0.05 kJ/m 2 ) for old bone (Nalla et al. 2004). The reduced material toughness could be one explanation to the straighter cracks seen in old bone (Chan et al. 2009;Koester et al. 2011) as a reduction in toughness would increase the risk of crack penetration of osteons (Fig. 4). This could be highly relevant for understanding crack propagation in cortical bone and changes in fracture resistance due to aging.
The stiffness of the cement line also had an impact on the crack trajectory (Fig. 5), which confirmed the screening results (Tables 4, 5, Fig. 3). With reduced interface stiffness, crack deflections were seen at higher critical interface strains (seen as a shift to the right in Fig. 5). Furthermore, small kinks in the interface (crack score 2) were more common with a more compliant cement line under conditions where direct crack penetrations (crack score 1) occurred in stiffer  (Fig. 5). Nevertheless, experimental findings suggest that interface strength rather than stiffness might be more relevant for determining crack paths in cortical bone. A study at the sub-osteonal scale compared osteonal lamellae and interlamellar areas and reported that both areas had similar mineralization (Katsamenis et al. 2013). Still, cracks preferably propagated through interlamellar areas, even when it meant major deflections from the initial crack path (Katsamenis et al. 2013). Katsamenis et al. (2013) suggested that the orientation of mineralized collagen fibers, rather than the level of mineralization, could be a predictor for the crack path. Similarly, results from osteon push-out tests suggest that the strength of the cement line interface depends on the collagen fiber orientation (Bigley et al. 2006). This means that brittle rather than compliant interfaces could be responsible for deflecting cracks in cortical bone, which is illustrated in Fig. 4. Still, further experimental characterization of the cement line is needed, with focus on the damage properties. Quantification of normal and shear damage properties in the interface would improve the predictive capability of computer models and our understanding of the role of microstructure in cortical bone.

Model limitations and future studies
Crack propagation, both in experimental setups and in vivo, is a 3D phenomenon that was simplified and analyzed in 2D with three different osteon orientations in this study. The literature considering simulation of crack propagation in bone in full 3D is scarce. Available works focus on analyzing crack initiation as converge problems do not allow the full crack paths to be captured. Models in 3D have been used to simulate crack initiation in femurs (Ali et al. 2014;Marco et al. 2018b), teeth (Zhang et al. 2016) and trabecular bone (Hammond et al. 2019). Another approach was used by Demirtas et al. (2016) who evaluated the possibility to use both XFEM and cohesive elements to simulate crack propagation in cortical bone 3D-geometries. Crack propagation in matrix and osteons was modeled with XFEM but limited to propagate in 2D within a given plane perpendicular to the cortical bone microstructure, while cement lines were modeled with cohesive elements following the cylindrical osteons. An important limitation of that study was that XFEM and cohesive element cracks were not connected; hence, the matrix crack was confined to propagate in the same 2D plane even when cracks followed the cement line out of plane. In all, it is a great challenge to simulate crack propagation in 3D. Our future studies will explore the possibility to extend our modeling framework in 3D to model continuous crack paths in cortical bone at the microscale. The interface damage model can be generalized to 3D by adding the second shear t direction as but simulating crack propagation in cortical bone, including large crack deflections out of plane, is still beyond the state of the art. In the meantime, important aspects, such as crack deflections and interaction with cortical bone microstructure, can qualitatively be evaluated in 2D.
Another limitation is the simplified material descriptions used in this study, with isotropic linear elastic material models and no open pores. Our previous study showed that Haversian canals exaggeratedly affected the structures in longitudinal and transversal models in 2D, where critical interface strains 50 times lower than the critical matrix strains were required for crack deflection (Gustafsson et al. 2018a). To keep the same model design for all orientations, all Haversian canals were filled with osteon material in this study. The main outcome of the study is not expected to change from excluding the Haversian canals, as relative rather than absolute numbers were of interest when identifying important parameters focusing on the interaction between the crack and the cement line interface. Furthermore, cortical bone displays anisotropic toughening, where the difference in toughness originates from extrinsic toughening mechanisms that are active in different orientations (Koester et al. 2008;Nalla et al. 2004;Zimmermann et al. 2009Zimmermann et al. , 2010. However, less variation is seen when comparing the damage initiation toughness in different directions (Koester et al. 2008), and therefore, the same material toughness values were used in all models. Intrinsic toughening mechanisms (e.g., micro-cracking) appear as plastic mechanisms at higher length scales, and the intrinsic mechanisms occurring in front of the crack tip were included in the damage law describing the cohesive crack. However, new cracks cannot nucleate when there is an active crack in the region. As a result, elements that fulfill the damage criteria away from the propagating crack are not able to fracture, and hence, the stiffness of the structure is overestimated. In complex models with several osteons, this limitation can influence the results, as strain concentrations then can occur at multiple locations in the model. Plastic material formulations for matrix and osteons could then be used to capture more realistic deformation patterns when limited to one single crack and the effect of plasticity could be of interest when looking at the effect of aging. However, for the simple case focusing on crack propagation around one osteon, the use of a simple elastic material formulation is assumed to have a minor impact on the result.
Besides the effect of material parameters on crack propagation around an osteon, there are several geometrical f QUADE = n 0 n 2 + s 0 s 2 + t 0 t 2 parameters that could be important (e.g., osteon density, size, shape, porosity, cement line thickness etc.). To study the effect of geometrical parameters, more realistic models including multiple osteons are needed. Such complex models are not feasible to use for large screening studies. Instead, the presented parameter study will serve as a useful base when looking in to more realistic scenarios focusing at the effect of aging.
It has been shown repeatedly that cortical bone demonstrates a rising R-curve behavior, where the toughness increases with increased crack length (Koester et al. 2008(Koester et al. , 2011Nalla et al. 2004Nalla et al. , 2005Zimmermann et al. 2010). This is due to extrinsic toughening mechanisms, such as crack deflection and crack bridging, which are activated during crack propagation. The toughening effect is largest when the crack propagates perpendicular to the long axis of the osteons (Koester et al. 2008) and diminishes with age (Chan et al. 2009;Koester et al. 2011;Nalla et al. 2004Nalla et al. , 2006). An extension of the current models including multiple osteons and realistic osteon geometries will be used to evaluate if the XFEM model can predict a similar rising R-curve behavior as seen in experiments with anisotropic effects when comparing different osteon orientations.

Conclusion
To model crack propagation in cortical bone, it is necessary to use a model that can predict realistic crack paths and capture crack deflections in cement lines. For this, XFEM models are beneficial as they do not need predefined crack paths, and in combination with the proposed interface damage model, realistic crack paths can be simulated. Overall, the cement line properties are crucial for the crack trajectory and both the critical strains in the interface and the stiffness affect the ability to deflect propagating cracks. Further characterization of the interface is needed where also the effect of aging should be investigated. Critical damage strains in the interface have not yet been measured experimentally, and these should be evaluated in relation to critical damage initiation strains in matrix and osteons. The bulk fracture toughness of cortical bone is well documented in the literature (e.g., (Chan et al. 2009;Koester et al. 2008Koester et al. , 2011Nalla et al. 2004Nalla et al. , 2005, and this is assumed to correspond fairly well to the matrix toughness. However, the toughness value for damage initiation should be used for microstructural models instead of the total toughness values that include extrinsic toughening effects. The toughness values for osteons and cement lines have not been well determined experimentally. Our results indicate that these values are important for the total fracture energy dissipated during crack opening but have little importance in determining the crack trajectory. Finally, the results presented in this study illustrate how both the cement line interface and the cohesive state inside the crack when encountering an interface affect the crack trajectory and their interplay could be highly relevant for understanding how aging affects crack propagation and fracture resistance in cortical bone. Funding This study was funded by the Swedish Foundation for Strategic Research (Grant No. IB2013-0021).

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.