## Abstract

We examine the specific role of the structure of the network of pores in plexiform bone in its fracture behaviour under compression. Computed tomography scan images of the sample pre- and post-compressive failure show the existence of weak planes formed by aligned thin long pores extending through the length. We show that the physics of the fracture process is captured by a two-dimensional random spring network model that reproduces well the macroscopic response and qualitative features of fracture paths obtained experimentally, as well as avalanche statistics seen in recent experiments on porcine bone.

## 1. Introduction

Cortical or compact bone, like most structural biomaterials, is a natural composite of a biopolymer (collagen) and a biomineral (hydroxyapatite) along with other non-collagenous proteins and water, hierarchically organized at several length scales [1]. It is typically found in load-bearing bones like femur and tibia, both along the mid-shaft as well as a protective outer covering for the spongier (cancellous) bone near the joints. Being a living tissue whose primary function is to bear mechanical stresses, it dynamically adapts its microstructure and porosity network to local biomechanical forces, which in turn alters properties such as strength and resistance to fracture [2,3]. For assessment of its structural health, therefore, it is vital to understand the structure–property relation of bone material. This understanding is also central to selection of sites for extracting bone grafts [4], in design of mechanically compatible implants [5,6] and porous scaffolds for bone tissue engineering [7,8], to evaluate the effectiveness of chemical and physical therapeutical measures for bone healing [9], etc. While there has been significant effort towards experimental characterization of the structure–property relationship, development of predictive models that are able to reproduce characteristic features of the complex fracture processes involved has been comparatively quiescent. Such models, if robust enough, would also be applicable in simulation of fracture for a wider class of brittle porous materials, such as wood and rock [10,11].

Bulk cortical bone of bovids and cervids has predominantly two distinct microstructures: plexiform bone that is brick shaped and has woven bone and vasculature sandwiched within regular lamellae, and Haversian bone that has cylindrical secondary osteons that run along the length of the bone [12]. The mid-diaphysis of a human femur, except in those of very young bones, is known to be Haversian in the entire cross-sectional area. While plexiform microstructure results from the layer by layer laying of bone material, Haversian bone is formed by remodelling existing bone through a process of resorption and deposition of bone material. The mechanical properties of the bone depend on the microstructure. For instance, elastic modulus and Poisson ratio of Haversian bone are lower than that of plexiform bone [13–16]. In tension, Haversian bone has lower strength [17,18], fracture toughness [19,20] as well as resistance to fatigue crack growth [21–23]. In compression, however, Haversian bone has higher strength that is proportional to the extent of remodelling [24,25].

A significant factor affecting the mechanical properties within the same type of bone, whether cancellous or compact, is its porosity. From experimental characterization, it is evident that elastic modulus as well as strength decrease with mean porosity as a power law [17,18,26–28]. On the contrary, in cortical/compact bone, Haversian microstructure has higher porosity than plexiform bone [29,30], yet has higher compressive strength [25,30]. Furthermore, it has been recently argued that plexiform bone is remodelled into Haversian bone only in regions where the bone experiences high compressive stresses [25,31]. If indeed the plexiform bone is unsuitable for bearing compressive loads, the structure and arrangements of its pores, and not just the mean porosity, must be playing a crucial role in the failure process under compression.

The fracture process in complex heterogeneous quasi-brittle materials like bone involves multiple porosities and micro-cracks that interact and evolve stochastically prior to final failure [32]. Modelling such materials using classical fracture mechanics theory or standard finite-element method restricts the scope mostly to finding effective elastic behaviour or to finding the resistance to growth of a single macroscopic crack [33–35]. Statistical models, like the random spring network model (RSNM), approximate the continuum by a network of springs with statistically distributed characteristics. Thereby, it has been successful in providing an insight into the role of disorder in the fracture behaviour of heterogeneous material systems with no pre-existing crack [36], reproducing features like transition from brittle to non-brittle macroscopic response [37], avalanche size distributions [38] and qualitative [39,40] and quantitative [41] features of fracture of composite materials. In the context of bone, simple one-dimensional models of parallel springs with statistically distributed properties have been used to simulate the characteristic quasi-brittle softening seen in tension and bending [42] as well as in compression of cortical bone [43]. Including spatial effects, three-dimensional RSNM has been used to model cancellous bone to show that the variation in strength with porosity as predicted by a power law is an overestimate for higher porosities as beyond bond percolation threshold the network loses its load-bearing capacity [44,45]. However, by assuming the porosity of the bone material to be homogeneously distributed, these studies completely ignore the role of the structure of the porosity network.

Can the specific role of the structure of the porosity network of plexiform bone in its vulnerability against compressive failure be identified? In this paper, we take a combined experimental and theoretical approach. Using computed tomography (CT) scan imaging, we obtain the detailed porosity network of plexiform bone, both before and after failure under compressive loads. On the basis of the observed periodicity in the features of the porosity network in the tangential direction as well as the failure planes being along the same direction, it is argued that a two-dimensional model in the longitudinal–radial plane is sufficient to capture the physics of the fracture process provided that the sample is oriented so that its axes are well aligned with the porosity. A detailed two-dimensional RSNM is developed, which is able to successfully reproduce both the qualitative fracture paths as well as the quantitative stress–strain curves data. Furthermore, the significance of the spatial correlations in the porosity network in lowering the compressive strength is shown by comparing response of plexiform bone with that of an equivalent homogeneously porous material.

## 2. Micro-computed tomography scans and compression testing

Cubical samples of size 5 × 5 × 5 mm (figure 1*a*), a total of six, were machined from the anterior cortices at the mid-diaphysis of dry bovine femurs that were extracted from four different healthy animals (samples I and II from one, samples III and IV from the second, sample V from the third and sample VI from the fourth animal). The cross-sectional dimensions of the specimens in this study are similar to those used in previous studies on bovine femur [46,47]. Loss of water content in bone is known to result in the decrease of strength and fracture toughness of bone material [48]. However, several studies have used dry bone properties as an indicator for comparative bone quality considering that the testing methodologies for dry bone are significantly simpler [49–52]. In the optical micrograph of the longitudinal face, shown in figure 1*b*, layers of woven bone were present between arrays of primary osteons that is typical of plexiform bone [14,25,31].

To obtain the three-dimensional porosity network, the samples were scanned using a GE Phoenix v/tome/x s micro-CT machine. One thousand radiographic projections of isotropic resolution, 8 µm, were acquired with fixed exposure time of 1000 ms and frame averaging time of 4 s for each projection. After scanning, 500 equally spaced slices were reconstructed using Feldkamp algorithm. The reconstruction of sample III before compression testing is shown in figure 1*c*. It shows the presence of an intricate porosity network composed of several long thin pores, many of which appear to stretch across the thickness of the sample and are at a slight inclination to the longitudinal axis.

The same sample was then tested under quasi-static compression using smooth parallel platens along the longitudinal direction at a rate of 0.1 mm min^{−1} until the load-bearing capacity dropped to 80–90% of the peak load. The compliance in the experimental set-up was 2 ± 0.2 × 10^{−5} mm N^{−1}, which is of similar order as reported in [53]. Additionally, the predictions of compressive strength from simulations, as described later in §4, were found to be largely insensitive to the minor variations observed in machine compliance. The sample develops very well-defined failure planes as evident from the CT scan image of the sample post-compressive failure shown in figure 1*d*. In spite of the initial wide spread porosity network (figure 1*c*), interestingly, the extensive damage is largely localized only in a few planes that extend through the thickness in the tangential direction. These planes are parallel to each other and, notably, oriented at the same angle as that of long pores extending across the thickness prior to testing.

Compression in the longitudinal direction leads to development of tensile transverse strains both in the radial and tangential directions. The failure planes observed are aligned in the tangential direction as per the known orthotropy in elastic as well as fracture properties. At the microscopic scale, the pores are significantly more aligned in the tangential direction such that micro-splitting cracks develop and propagate in planes that stretch across thickness in the tangential direction. To explore this further, we analyse the existing correlations within the porosity network. The output of the CT scan image is greyscale in nature; however, the mean intensity varies with depth due to non-uniform lighting. This bias is eliminated by converting the data to binary. In each of the 500 CT scan slices in the tangential direction, we observe that the greyscale values are distributed around two widely separated peaks. We introduce a cut-off that lies between the peaks and map the greyscale values less than the cut-off to 0 (vacancy) and the others to 1 (occupied site). We make a list of all pores that are clusters of vacancies with cubic connectivity. The largest 50 pores are shown in figure 2*a*. They run through the length of the sample in the longitudinal direction at a slight angle to the vertical as in the output of the CT scan. The two-dimensional longitudinal view of the clusters is straight lines aligned in the radial direction (figure 2*b*), implying that the clusters or pores are contained within the longitudinal–radial plane. In addition, the experimental fracture planes lie in the longitudinal–tangential planes. This allows us to construct a simplified two-dimensional model for the failure process from the three-dimensional data by coarse graining along the tangential direction.

## 3. Two-dimensional modelling

The total volume of a given sample is divided in the tangential direction into a set of 10 independent representative domains. The thickness of each is 500 µm (50 pixels) and is chosen in accordance with the two relevant length scales of the physical problem: larger than the largest pore size and of the order of the mean distance in tangential direction between pores. To obtain these distances, in a radial–tangential section as shown in figure 2*c*, we characterize the spatial distribution of vacancies by evaluating the probability that a site at a distance (*r*, *t*) from a given vacancy is also a vacancy, *g*(*r*, *t*). The probability in tangential direction, *g*(0, *t*), as presented in figure 2*d*, decreases exponentially from unity and has well-defined peaks at larger *t*. The first peak in *g*(0, *t*), indicative of the mean distance between pores, is at ≈500 µm. For smaller *t*, *g*(0, *t*) is fitted with an exponential exp(−*t*/*ξ*) to obtain *ξ* ≈ 10 µm (inset of figure 2*d*), which represents the dimension of a typical pore in the tangential direction.

The spatial distribution of the porosity within each domain is next reduced to a two-dimensional representation in R–T plane by averaging the binary data through the thickness that each of domain represents. The two-dimensional porosity distribution for 5 out of the 10 sets that were analysed are presented in figures 3*a* and 4*a* for samples II and III, respectively. Note that now the value associated with each pixel is between 0 and 100%. Most of the domain is fairly dense, having porosity lower than 5% including the porosity structure associated with the layering inherent in the plexiform bone. Porosities larger than 20% appear as few short pores with larger width, aligned vertically along the longitudinal axis, and do not extend across the height of the specimen (figure 3*b*) or are absent (figure 4*b*). By contrast, porosities in the 5–20% range (figures 3*c* and 4*c*) are finer, significantly larger in number, present in all the sections of both bones, and most of them extend across the height while being slightly inclined to the longitudinal axis. These inclined pores are perhaps part of the primary vasculature including blood vessels, etc. A few of them are vertically aligned and shorter in length (250 µm maximum).

To identify the features of the porosity network that contribute to the fracture mechanism, we examine the micro-CT scans of the fractured specimen for patterns in the preferred planes for splitting in terms of their orientations with respect to the structure of the initial porosity network. The two-dimensional representations of the porosities and failure planes are shown in figures 3*d* and 4*d*. It is observed that failure planes tend to follow those longest pores that have similarly aligned pores in neighbouring domains. The pores having higher porosity (more than 20%), even though they are weaker themselves and therefore more likely to act as stress concentrators, tend to have only a local effect on the fracture paths as they do not have any obvious correlation with the rest of the porosity network. Apart from the primary failure paths which follow the porosity network, there are also auxiliary failure paths which form as the load transfers differently post-formation of primary failure paths.

## 4. Random spring network model

We model the complex fracture processes of the two-dimensional representations using RSNM of size 150 × 150 on a square lattice with the nearest and next-nearest neighbour interactions, as presented in figure 5. The potential energy of the network, *V*, has contributions from extension of springs, distortion of angle between springs and the Hertzian contact between particles:
4.1where the sums *i*, *j*, are over linear springs, torsional springs and pairs of particles; *k _{i}*,

*c*are spring constants;

_{j}*δr*and

_{i}*δθ*are the extensions and distortions, respectively;

_{j}*α*is a material constant;

*d*is the diameter;

*r*is the distance between particles

_{ij}*i*,

*j*; and

*Θ*(

*x*) is the Heaviside theta function. If the spring constants of horizontal/vertical, diagonal and torsional springs are 2

*k*,

*k*and

*c*, respectively, then the system is isotropic with elastic modulus,

*E*, and the Poisson ratio,

*v*, given by [54] 4.2where

*a*is the lattice spacing. The sites at the top and bottom rows are connected to springs whose compliance matches with that of the machine used in experiments. After a displacement, the system is equilibrated to its minimum energy configuration by numerically integrating the equations of motion 4.3using velocity-Verlet algorithm [55], where the friction

*γ*dissipates energy [38]. For each increment in the applied downward displacement, after equilibration, if any spring

*i*is stretched beyond its threshold strain , it and the torsional springs associated with it are broken. The system is re-equilibrated until no further breaking takes place.

The experimentally obtained porosity network is incorporated into the RSNM by assigning the spring constants and breaking thresholds in accordance with the local porosity. Assuming that at small length scales the material behaviour is a function only of its local porosity, the porosity-dependent material constants are calibrated by evaluating them separately from the tensile macroscopic response of a porous network that has homogeneous matrix properties. The spring constant and failure threshold of the matrix are chosen such that the known properties of cortical bone in tension for approximately 4% porosity are well reproduced [29]. Using two sizes of porosities (50 µm for 0 ≤ *p* ≤ 5% and 100 µm for 5% ≤ *p* ≤ 20%), the elastic modulus, the Poisson ratio and failure strain are estimated from 10 realizations for each porosity configuration. The mean values from simulations are used to curve fit a power-law behaviour over a large range of porosities as shown in figure 6. With increase in porosity, the elastic modulus decreases while the Poisson ratio increases. The coupled effect is overall reduction in bulk modulus (approx. 12–13% reduction at 10% porosity), which implies increased compressibility as expected for a porous material. For *p* > 50%, failure strain is considered as zero as the bond percolation threshold for square lattice is 0.5. For thermodynamic stability, both the extensional spring constant *k* and the bending constant *c* must be positive which as per equation (4.2) puts an upper bound of one-third on the Poisson ratio. The above calibration assigns elastic modulus, Poisson ratio and strain threshold for a given porosity in a deterministic manner. To account for uncertainties in the calibration, stochasticity was introduced in the fracture behaviour by taking the strain threshold from a Gaussian distribution centred about the calibrated deterministic value and 5% standard deviation. The macroscopic response of each two-dimensional representative network was then obtained by averaging over five independent realizations of the strain thresholds.

The model predictions of the macroscopic response, generated by taking the mean of the averaged response of all 10 representative networks for each sample, shown in figure 7, reproduce well the elastic regime as well as the compressive strength. The experimental curves do not start at the origin as the toe regions were removed to account for the error in stiffness during the initial contact. The quantified data comparison for compressive strength and corresponding strain is presented in table 1. The model predictions for compressive strength fall within 20% of the experimental values except for sample V which is under predicted by 29%. The predictions for the corresponding strains fall within 25%.

We generate the corresponding failure paths by representing sites with higher (lower) probability of multiple bond breakage by white (blue), as shown in figures 3*e* and 4*e*. The predicted crack paths capture the characteristic features of the experimental data: damage growth is localized and follows the slightly inclined fine porosity network in most of the sections. Interestingly, in this study the samples are assumed to have similar bone matrix behaviour such that they differ only in terms of their porosity structure, and yet the predictions of the macroscopic response (concordance correlation coefficient for the compressive strength, *ρ*_{c} = 0.21) as well as failure path are very realistic. The predictions can be further improved by acquiring and incorporating bone-specific matrix behaviour.

Furthermore, to determine the influence of the structure of the porosity network, a corresponding model was developed that had the same porosity but distributed homogeneously, as shown for sample II at *t*/*T* = 0.05 in figure 8*a,b*. Comparison of the model predictions, shown in figure 8*e*, indicates the behaviour in the elastic regime to be similar in both cases, however, the failure behaviour is distinctly different. While the failure is localized in the original model along the porosity network, failure in homogenized model is diffused resulting in higher load-bearing capacity by as much as 25%. The porosity network in plexiform bone, perhaps arising due to physiological reasons of growth requirements, also results in localization of damage in lattice model and lower load-bearing capacity in compression implying the need for smaller and shorter porosities, similar to the remodelled Haversian bone.

We also measure the avalanche size distribution *P*(*s*), where the avalanche size *s* is defined as the number of bonds broken in an increment of the applied displacement. The avalanche size distribution follows a power law, , with *τ _{s}* varying between 1.8 and 2.06, as shown in the insets of figure 7. The energy emitted during an avalanche is known to be proportional to the square of the avalanche size (e.g. [38]). Let

*P*(

*E*) be the distribution of the energy emitted during an increment. From the probability conservation relation, , it follows that

*P*(

*E*) is also a power law with

*τ*= (1 +

_{E}*τ*)/2.

_{s}*τ*thus varies between 1.4 and 1.5. This is in excellent agreement with acoustic emission experimental data, recently reported for porcine cortical bone, where

_{E}*τ*varies from 1.3 to 1.7 [56].

_{E}## 5. Conclusion

We examined, both experimentally and computationally, the role of the inherent porosity network of bovine plexiform bone in its susceptibility to compressive failure. CT scans of samples show that plexiform bone is particularly vulnerable to splitting fracture primarily as a result of its well-defined porosity network. We developed a two-dimensional square RSNM incorporating the effects of the local porosity network that was able to reproduce both the quantitative macroscopic response as well as the qualitative features of the fracture paths observed in the experiments. Homogenization of the porosity showed a considerable increase in load-bearing capacity as the damage growth was now not localized suggesting advantages in remodelling of the plexiform bone to Haversian bone for increase in its load-bearing capacity in regions of high compression. It is to be noted however, weak planes in tangential direction are not necessarily detrimental in all failure modes; for instance, in tension weaker planes are known to deflect growth of cracks leading to higher toughness and strength.

The random spring network based predictive modelling, in this study, has been shown to be able to capture quantitative features of splitting fracture of plexiform bone with a given intricate pore structure. The damage initiation and growth in the bone material involves interaction between multiple defects and their evolution with increasing applied strain which result in a large number of dissipative micro-fracture events leading to progressive loss in load-bearing capacity prior to final failure. Such quantifiable details are outside the scope of classical approaches based on basic material science, stress analysis and fracture mechanics and therefore the model predictions make significant additions to the present understanding of bone fracture. However, the study is limited in its scope as it investigates dry plexiform bone. Extending its applicability to wet bone and Haversian bone, more relevant for human bone, would require significant modifications that are interesting areas for future studies. For further refining the results of this study, other effects need to be included such as friction at contacts, local orthotropy of bone material, bone-specific material behaviour, etc., which are part of ongoing projects.

## Authors' contributions

A.M., A.B. and R.R. contributed equally to the paper. P.P. participated in the initial data analysis.

## Competing interests

The authors have no competing interests.

## Funding

No funding has been received for this article.

## Acknowledgements

We thank Prof. Krishnan Balasubramanian for providing access to the experimental facilities of the Centre for Non-Destructive Evaluation, IIT Madras.

- Received October 6, 2016.
- Accepted November 8, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.