Two distinct geometrical models of bone at the nanoscale (collagen fibril and mineral platelets) are analyzed computationally. In the first model (model I), minerals are periodically distributed in a staggered manner in a collagen matrix while in the second model (model II), minerals form continuous layers outside the collagen fibril. Elastic modulus and strength of bone at the nanoscale, represented by these two models under longitudinal tensile loading, are studied using a finite element (FE) software abaqus. The analysis employs a traction-separation law (cohesive surface modeling) at various interfaces in the models to account for interfacial delaminations. Plane stress, plane strain, and axisymmetric versions of the two models are considered. Model II is found to have a higher stiffness than model I for all cases. For strength, the two models alternate the superiority of performance depending on the inputs and assumptions used. For model II, the axisymmetric case gives higher results than the plane stress and plane strain cases while an opposite trend is observed for model I. For axisymmetric case, model II shows greater strength and stiffness compared to model I. The collagen–mineral arrangement of bone at nanoscale forms a basic building block of bone. Thus, knowledge of its mechanical properties is of high scientific and clinical interests.

## Introduction

Bone is a hierarchically structured biological composite material, and thus its properties depend on its architectures at different structural levels [1–3]. Bone is composed of organic matter, inorganic hydroxyapatite (HA) nanocrystals, and water [1,4]. The organic phase consists of collagen (around 90 wt %) and noncollagenous proteins (10 wt %), which interact with both minerals and collagen [4,5]. These constituents are arranged in a complex way giving bone its excellent mechanical properties; bone is strong, stiff, tough, and lightweight. Collagen, which is a soft and highly deformable protein made of nanosized tropocollagen molecules, provides bone with its capacity for energy dissipation when bone is under mechanical deformation, while HA platelet-shaped crystals, which are strong and stiff but brittle, contribute to bone's high stiffness and strength [6–8]. At the nanoscale, the tropocollagen molecules are crosslinked to form collagen fibrils which are biomineralized by the HA crystals [4]. The mineralized collagen fibril forms a basic building block of bone. Thus, the structure and resulting properties of bone at the nanoscale are of high scientific and clinical interests as they influence the overall properties of bone. Mechanical properties of bone at the nanoscale can serve as valuable inputs for multiscale studies of bone strength and adaptation, among other applications.

The main deformation mechanisms in collagen fibrils are the breaking of cross-links and intermolecular sliding between collagen molecules which, in turn, play a vital role in bone deformation [6,9]. Another mechanism for energy dissipation in bone is sliding between collagen and HA platelets. This mechanism leads to the formation of plastic zones around defects and cracks, and it is proposed to be the dominant toughening mechanism in bone, allowing local energy dissipation and leading to protection of whole bone structure [10].

Bone structure–property relations at the nanoscale have been widely studied [4,11]. However, there is still no full consensus on the collagen/HA arrangement in bone. This is due to the limitations of current imaging techniques in probing submicron scale needed to visualize spatial arrangements of the nanoconstituents of bone, their dimensions, the nature of interfaces between the constituents, and interface energies which all significantly affect the mechanical behavior of bone [12]. Different possible interfaces may exist at subnanoscale/nanoscale such as collagen–HA, HA–HA, collagen–noncollagenous proteins, and HA–noncollagenous protein interfaces, with water playing a role in bonding. Other factors that make understanding of bone properties challenging are bone dynamics (self-healing mechanisms and remodeling processes), spatial variations in bone properties, changes due to bone disuse, gender, age, prior history of broken bones, and other factors [2,10].

Researchers have used scanning electron microscopy, transmission electron microscopy, atomic force microscopy, and other experimental methods, to investigate the arrangements of constituents of bone, dimensions of the constituents, and interfaces between them [4,13–20]. Most studies agree that minerals are in the shape of elongated nanoscale polycrystalline platelets with dimensions of 3–5 nm × 25 nm × 50–100 nm [17,21,22]. Regarding the collagen/HA arrangement, some studies claim that the HA minerals are located in the gap zones inside the collagen fibrils [18–20,23] while others believe that minerals are both outside and inside the collagen fibrils, while yet others concluded that they are mainly situated outside the collagen fibrils [16,17,24–29]. When modeling a mineralized collagen fibril, Nikolov and Raabe [30] included both intrafibrillar and extrafibrillar minerals in their model. They assumed that ∼25% of minerals are located outside of collagen fibrils and used homogenization techniques to estimate the elastic properties. Hellmich et al. [31] modeled the elastic properties of bone at the nanoscale assuming three representations of collagen–mineral interaction including collagen inclusions in a mineral foam matrix, interpenetrating network of HA and collagen, and collagen–HA network in a mineral foam matrix. Recent studies suggest that minerals and collagen form interpenetrating phases, and the pores are filled with noncollagenous proteins and water [32,33]. Finally, most studies of bone at the nanoscale focus on predictions of bone's stiffness [11], while little attention has been given to modeling of strength and fracture.

In this paper, we study computationally the stiffness and strength of bone at the nanoscale assuming two different geometrical arrangements of HA minerals, one with a staggered arrangement of HA inside collagen fibrils and second with the HA minerals residing outside collagen fibrils (see Fig. 1). Comparison of the mechanical properties computed using these two models can provide a further understanding of the behavior of bone at the nanoscale.

## Geometry of Models and Assumptions

We study the mechanical properties of bone at the nanoscale by considering two cases: one when minerals are isolated and located inside collagen fibrils (model I) and second when minerals are connected and form sheets outside collagen fibrils (model II). The applied loading is a longitudinal tensile loading aligned along the collagen fibril direction. Thus, the computed properties are the longitudinal tensile elastic modulus and strength of bone at the nanoscale (at collagen and mineral level). The analysis employs a finite element (FE) method and a cohesive surface law for interfaces between the constituents of bone to account for bone fracture. Both models are studied under the assumptions of plane stress, plane strain, and axisymmetry. We include plane stress/strain cases since most studies on bone at the nanoscale used planar geometries. Second, since a three-dimensional (3D) arrangement of crystals in a mineralized collagen fibril is still not well characterized, we assume an axisymmetric case for simplicity. To facilitate comparison between these models, the volume fraction of the HA minerals in all considered models is 45%.

### Model I.

Model I follows the geometric arrangement proposed by Jäger and Fratzl [34], which is the most widely used model in the literature for bone at the nanoscale [6,35,36]. In this model, the HA platelets are periodically arranged in a staggered fashion inside a collagen fibril as shown in Fig. 2. For plane stress and strain cases, the dimensions of the HA platelets are taken as 100 nm in length and 4 nm in thickness following [6]. For axisymmetric case, the HA platelets are assumed to have the length of 100 nm while the thickness increases with the distance from the center (1 nm, 3 nm, and 3.42 nm). These mineral thicknesses were chosen such that all cases have the same volume fraction (45%) of the HA platelets. Figure 2 also shows the boundary conditions applied on model I. In computations, model I is fixed at the bottom surface and is displaced at the top surface. In the plane stress/strain cases, the left side is fixed (zero normal displacement) in the horizontal direction and has a zero traction in the vertical direction (zero shear traction) while the right side has zero tractions in both directions (see Fig. 2(a)). In the axisymmetric case, symmetric boundary conditions (fixed in the horizontal direction and no rotation about the other two directions) are applied on the centerline, and a zero traction boundary condition is applied on the other vertical side (see Fig. 2(b)).

### Model II.

Model II assumes that the elongated HA mineral platelets are arranged outside the collagen fibril. This model is motivated by experimental studies of Refs. [13] and [17] and earlier reports that minerals reside mainly outside of collagen fibrils and form circumferential lamellae [27–29,37]. Since there is limited information in the literature on the nature of collagen–HA and HA–HA interfaces, we consider two interfacial conditions: (a) a noncollagenous matrix (interphase) exists between minerals [38] and (b) no matrix (no interphase) exists between the minerals. Figure 3 depicts model II under different geometrical assumptions. Figures 3(a) and 3(b) show the cases when there is no matrix between the minerals while Figs. 3(c) and 3(d) represent the cases when the matrix exists between the minerals. In all the considered cases, the HA platelets are located outside the collagen phase. In this model, the length of the HA platelets is assumed to be 100 nm, as in model I, while the thicknesses vary from one case to another. The thicknesses of geometric models shown in Figs. 3(a)–3(d) are 8 nm, 2 nm, 9 nm, and 4 nm, respectively. Although the thickness of the HA crystals varies from one case to another, the volume fraction of the HA platelets is 45% for all cases. Boundary conditions applied on model II are shown in Fig. 3. Model II is fixed at the bottom surface and is displaced at the top surface. In the cases of plane stress/strain, the left side is fixed in the horizontal direction (zero normal displacement) and has a zero traction in the vertical direction (zero shear traction) while the right side is traction-free (see Figs. 3(a) and 3(c)). For the axisymmetric case, symmetric boundary conditions are applied on the centerline, and a traction-free boundary condition is applied on the other vertical side (see Figs. 3(b) and 3(d)).

### Model Properties.

For simplicity, mechanical properties of the minerals, collagen, and matrix (noncollagenous interphase) are assumed to be linear elastic and isotropic. There is still no full consensus in the literature about the properties of collagen and minerals [6,35,39–42]. The elastic properties used in this paper are listed in Table 1; they are within a range of those reported in the literature.

In this paper, we assume that the failure mechanism in mineralized collagen fibril is a delamination of minerals and collagen/matrix, i.e., fracture happens at the interfaces between different phases. The interfacial characteristics of bone at the nanoscale have been studied based on experimental data, finite element results, and molecular dynamics results [6,43–47]. In order to investigate the effect of the interfacial properties between the constituents of bone in different models, we vary the total fracture energy between 0.01 and 1.0 $J/m2$. This range of the fracture energy values represents different types of interfaces: strong (ionic interactions) and intermediate (thin water layer) or weak (thick water layers) following [6,43–47]. Cohesion of collagen and mineral interface is due to the existence of a layer of structural water following [35,46,48,49]. Total fracture energy is defined as the total energy released before fracture occurrence. In this study, the strength of the interface is set to be 30 and 64 MPa representing weak and strong interfaces, respectively, following [6,43]. The stiffness of the interface is assigned to be 80 GPa; the choice of the stiffness is mainly based on numerical studies [6,43]. Generally, failure in bone occurs due to debonding of minerals and collagen (interfacial fracture) rather than the fracture of a single phase [6,35]. We assume that failure due to debonding happens before the failure of collagen or minerals, and thus, for simplicity, debonding is assumed to be the only failure mechanism.

## Governing Equations and Finite Element Analysis

where $\sigma $ represents the second-order stress tensor.

Fracture is assumed to occur only at the interfaces. Delamination is studied using a cohesive surface model [50,51]. It must be noted that this surface-based cohesive model is different than a cohesive element (zone) model [50]. The cohesive element (zone) behavior is defined by a thickness, and this may lead to erroneous results if used with plane stress assumption. This is due to the fact that in the case of plane stress the body exhibits change in its thickness, which are usually not accounted for in the cohesive element (zone) model. On the other hand, the cohesive surface behavior is defined as a surface interaction property in abaqus, it does not add mass to the model (contact property rather than material property), and it is typically used for cases in which the thickness of the interface is relatively small; therefore, it can be properly used for plane stress case.

The interface constitutive model used in this paper has a cohesive surface contact behavior which is based on the traction-separation response that relates the opening displacements to resisting tractions [50]. Figure 4 portrays the traction-separation law, which is a law allowing the separation to occur tangentially and normally. The traction-separation assumes an initial linear elastic response followed by a damage initiation and evolution. The initial linear elastic response can be described according to the following expression:

where $tN\xaf$ and $tT\xaf$ are the traction components predicted by the linear elastic traction-separation law at the current separations assuming no damage exists. *D* = 0 represents no damage case while *D* = 1 represents fully damaged case. Moreover, traction-separation assumes that there is no interfacial delamination due to compressive traction. In other words, delamination occurs due to tensile or shear loading at the interface. Softening is assumed to occur here based on fracture energy which is the energy dissipated due to failure, and it can be quantified by the area under the traction-separation curve.

## Results and Discussion

For axisymmetric models, eight-node biquadratic axisymmetric quadrilateral elements (labeled CAX8R in abaqus element library) were used while for plane stress and plane strain models, eight-node biquadratic plane stress (labeled CPS8R in abaqus element library) and plane strain (labeled CPE8R in abaqus element library) quadrilateral elements were used, respectively. Meshes of the different models were refined until the results insensitive to the mesh were obtained. Model I shown in Fig. 2(a) has 44,000 elements while model I shown in Fig. 2(b) has 25,000 elements. Model II shown in Figs. 3(a)–3(d) has 150,000, 14,400, 369,000, and 126,000 elements, respectively. Finite element simulations were done for models I and II with varying interfacial fracture energy under tensile loading in the longitudinal direction (along the length of collagen fiber).

### Model I.

Figures 5(a) and 5(b) portray the effect of fracture energy on the stress–strain curves of model I under the plane stress and plane strain conditions, respectively. It can be observed that the stress–strain curves for both cases are almost the same. Figure 5(c) represents the effect of fracture energy on the stress–strain curves under axisymmetric condition. It is seen that the axisymmetric assumption leads to lower stresses compared to the plane stress and plane strain cases. The strength of the interface is set to be 64 MPa in the calculations. Table 2 reports the elastic moduli for plane strain, plane stress, and axisymmetric conditions. The elastic moduli are almost the same for plane stress and strain conditions, while the elastic modulus for the axisymmetric case is lower than the ones for plane stress and strain cases. Moreover, it can be inferred from Fig. 5 that the fracture energy does not affect the stiffness, but it significantly influences the strength as is replotted in Fig. 6. The stiffness of model I is not affected by the fracture energy because of the same initial linear behavior of the traction-separation law. Once the maximum strength at the interface is reached, the crack starts to propagate according to different prescribed damage evolution laws, and this introduces the nonlinear mechanical response shown in Fig. 5. Figure 5 indicates that strain at failure increases with the increase of the interfacial fracture energy for all (plane stress, plane strain, and axisymmetric) cases. In addition, Fig. 6 shows that a plateau is obtained at high fracture energy values. Figure 7 portrays the von Mises stress in model I at a strain of 2%. Figures 5 and 6 and Table 2 show similar results for both plane strain and plane stress conditions. Hence, similar stress fields are expected for both plane elasticity cases as shown in Fig. 7. However, in the case of plane stress condition, the stress field in the collagen phase is more uniform than in the case of plane strain, and this indicates that collagen carries less load in the case of plane stress. In the axisymmetric case, the stresses are concentrated at the HA located in the center and are not distributed over the HA platelets as in the cases of plane stress and plane strain conditions. Hence, higher stress levels are observed compared to the plane stress and plane strain conditions at the same applied strain, which leads to faster failure for the axisymmetric condition. Generally, tensile loading is mainly carried by HA minerals while the collagen contributes by transferring the load between adjacent mineral platelets. Furthermore, the mismatch between the elastic properties of collagen and minerals causes significant sliding between them as implied by the discontinuity of the contours (see Fig. 7). The delamination starts at the ends of the crystals and then propagates along their sides.

### Model II.

In model II, the mineral crystals are arranged outside the collagen phase as shown in Fig. 2. The effect of the presence/absence of a noncollagenous organic layer between mineral crystals is investigated and compared with the results of model I. Figure 8 represents the stress–strain curves for the different models for the fracture energy of 0.2 $J/m2$ and interfacial strength of 64 MPa. Initially, a linear response is observed and then nonlinearity starts due to the damage evolution accounted for in the traction-separation law applied at interfaces. The elastic modulus, strength, and strain at failure significantly vary for the different assumptions and conditions. Plane stress and plane strain cases give very similar results. Figure 9 demonstrates the longitudinal elastic moduli for the different models. In Fig. 9, we only show plane stress case and omit plane strain case, for clarity. The highest elastic modulus is attained in the case of plane stress/strain without a matrix between minerals. The other models have almost the same elastic modulus. All different conditions of model II have a higher stiffness than the stiffness computed for model I.

Figure 10 illustrates the effect of the interfacial fracture energy on the longitudinal tensile strength of bone at the nanoscale. The maximum strength is obtained in the axisymmetric model without matrix between minerals, followed by the axisymmetric model with matrix between minerals. On the other hand, the models under the plane stress or plane strain conditions with a matrix between minerals have lowest strengths. The strength of model II ranges from 140 to 370 MPa, while the strength of model I varies from 220 to 360 MPa based on the applied conditions (see Figs. 5 and 8). In addition, Fig. 11 reports the strength of bone at the nanoscale for the different cases for two different interfacial strengths, 30 and 64 MPa, representing weak and strong interfaces, respectively [6,43]. The results shown in Fig. 11 are computed assuming interfacial fracture energy of 0.2 $J/m2$. Here, again we only show the plane stress cases since the corresponding plane strain cases are very similar.

Figure 12 portrays von Mises stress contours for the different cases of model II. In all cases, collagen has a uniform stress distribution which implies homogeneous stretching of collagen. The load is mainly carried by minerals, and the maximum stress occurs at crystals within central lamella for all different cases (see Fig. 12). Moreover, fracture starts at the interfaces between the outer and inner crystals, and this is implied by the stress relaxation (low stress values) that takes place at fracture lines and the discontinuity of the contours as represented in Fig. 12. Similar to model I, debonding starts at the short sides of the mineral platelets and then sliding between collagen and minerals propagates to the long sides of the mineral platelets.

### Comparison and Limitations.

Interestingly, model II shows higher values for stiffness and strength for the axisymmetric case than for the plane stress and plane strain cases while model I shows the opposite trend (see Figs. 6 and 10). The analyses that are based on plane stress and strain assumptions are performed for completeness and because most of the earlier studies considered planar geometries. The axisymmetric models provide idealized representations of the 3D HA–collagen arrangement. Yuan et al. [52] stated that the axisymmetric case gives similar results to a detailed 3D model for the elastic case. Thus, the axisymmetric assumption is a good approximation for the elastic case and may lead to reasonable approximations for strength. Second, since it is not well understood how the HA crystals are arranged in 3D manner in a mineralized collagen fibril, these idealized geometries provide a start for future studies accounting for more detailed representation of crystals in bone. Also, based on the literature, it is believed that crystals reside both inside and outside of collagen fibrils but percentages are still not well quantified. This study addresses the two limit cases of extrafibrillar and intrafibrillar crystals. Further studies of more detailed representations of HA crystal arrangements and accounting for HA crystals both within and outside collagen fibrils could be done in the future.

The elastic moduli of bone at the nanoscale for three cases out of the four cases considered for model II are found to be around 30 GPa (see Fig. 9). These results are consistent with statistical nano-indentation results [53] and ultrasonic data [54]. Furthermore, these results agree with studies accounting for extrafibrillar mineralization [30,31,37,55–57]. Moreover, properties predicted by the axisymmetric model II are close to the results obtained from micropillar mechanical tests, which give elastic moduli around 30 GPa and uniaxial strength of around 400 MPa [58] (see Figs. 9 and 13). All these comparisons indicate that model II is more realistic than model I.

Nair et al. [26], among others, discussed the role of extrafibrillar mineralization. They performed atomistic simulation of 3D molecular structure of collagen and minerals with a staggered arrangement of minerals in a collagen matrix under compressive loading. They computed elastic modulus as a function of mineral content (up to 40% by volume) and found that the elastic modulus was significantly lower than the macroscopic elastic modulus of bone reported in literature. Thus, they concluded that intrafibrillar minerals alone are not sufficient and extrafibrillar mineralization plays an important role in mechanical properties of bone. Our findings agree with the above-mentioned study as the model II represents the extrafibrillar mineralization case. The model II geometry was motivated by experimental observations of Schwarz [13,16,17] and others reporting on extrafibrillar mineralization [24–29].

In summary, for the axisymmetric case model II has higher tensile longitudinal elastic modulus and strength than model I (see Fig. 13). This indicates that in model II collagen/crystal arrangement is more optimal for mechanical performance. Such geometry should also give enhanced bending and torsion responses (not considered here) since stresses under such loadings will be maximum at the outer lamellae. Thus, minerals which reside at the outer shell and are stiffer and stronger would carry such loads. Such geometry is also present at the whole bone level, which is designed by nature as a shaft with stiff and strong outer core formed by dense cortical bone to effectively resist complex mechanical loads.

This paper has several limitations. Models I and II represent idealized limit cases, while actual structure of bone at the nanoscale may involve a combination of both, with minerals being outside and within the collagen fibril. Models I and II were analyzed using axisymmetric, plane stress and plane strain assumptions for simplicity. More realistic 3D models should be explored in the future. Unfortunately, there is no full consensus in literature about the distribution of collagen and minerals in the 3D space [41,52,59–64]. This leads to difficulties in developing a realistic 3D model representing bone at the nanoscale.

Another limitation of the current work is that collagen and nanocollagenous proteins (interphase) were assumed to have linear elastic responses while more realistic models would involve nonlinear behaviors. For analysis, we selected representative values for collagen and HA elastic properties which are within a range of those reported in literature (see for example Table 1 in Ref. [65]). Parametric study could be done using other values for elastic constants. Our interest was to compare models I and II, representing intrafibrillar and extrafibrillar cases, respectively. Furthermore, certain values were assumed for the interfacial strength and fracture energy, based on values reported in literature [6,43], obtained from experiments or molecular level simulations. However, further experimental and atomistic simulation studies could be done in the future to obtain further insights on interfacial strengths. Also, we did not account for the presence of water at interfaces which would lead to a viscoelastic interfacial response [46,49]. It has been shown that water has a stabilizing effect of collagen triple helix as tropocollagen molecules that are hydrated require more energy to untie from surface of HA [1,2]. As a polar solvent with capability of making strong hydrogen bonds, water has an affinity for both collagen and HA; thus, it acts as a bridge between collagen and HA and strengthens the interface [3]. In addition, water delays failure of the HA–collagen as it acts as a glue between tropocollagen molecules [4,5,66–70]. In the analysis, we used mixed boundary conditions at the outer surfaces of the mineralized collagen fibril. More realistic boundary conditions would account for constraints from neighboring fibers [63]. Besides, our model does not consider effect of cross-linking of collagen molecules. All these could be subjects of future studies. Thus, modeling of bone at the nanoscale is a rich problem with numerous extensions and topics open for future explorations.

## Conclusions

In this study, bone at the nanoscale is analyzed, under a uniaxial tension in the longitudinal direction assuming two different geometrical arrangements of HA minerals and collagen. In the first arrangement (model I), minerals are assumed to lie inside collagen fibrils, while the second arrangement (model II) assumes that minerals are arranged outside the collagen fibrils. Models are analyzed using a finite element method and a traction-separation law employed between phases of bone. The failure mechanism is assumed to be debonding as it is the most dominant mode of failure. Under the axisymmetric assumption, model II gives higher elastic modulus and strength than model I. The results obtained from model II are more realistic than the ones obtained from model I, and they are in good agreement with experimental results reported in the literature. This study provides new insights and contributes to a deeper understanding of the mechanics of bone at the nanoscale.

## Acknowledgment

We thank Professor Henry Schwarcz for sharing his findings and stimulating discussions on model II. This research was partially supported by the National Science Foundation DMR Program Grant No. 15-07169. The findings, conclusions, and recommendations expressed in this manuscript are those of the authors and do not necessarily reflect the views of the NSF.