On role of matrix behavior in compressive fracture of bovine cortical bone
aa r X i v : . [ c ond - m a t . s o f t ] J un On role of matrix behavior in compressive fracture of bovine cortical bone
Ashwij Mayya, ∗ Anuradha Banerjee, † and R. Rajesh
2, 3, ‡ Department of Applied Mechanics, Indian Institute of Technology-Madras, Chennai 600036, India The Institute of Mathematical Sciences, C.I.T. Campus, Taramani, Chennai 600113, India Homi Bhabha National Institute, Training School Complex, Anushakti Nagar, Mumbai 400094, India (Dated: May 16, 2018)In compressive fracture of dry plexiform bone, we examine the individual roles of overall meanporosity, the connectivity of the porosity network, and the elastic as well as the failure properties ofthe non-porous matrix, using a random spring network model. Porosity network structure is shownto reduce the compressive strength by upto 30%. However, the load bearing capacity increaseswith increase in either of the matrix properties – elastic modulus or failure strain threshold. Tovalidate the porosity-based RSNM model with available experimental data, bone-specific failurestrain thresholds for the ideal matrix of similar elastic properties were estimated to be within 60%of each other. Further, we observe the avalanche size exponents to be independent of the bone-dependent parameters as well as the structure of the porosity network.
PACS numbers: 87.85.G-, 62.20.mm, 87.10.Hk
I. INTRODUCTION
Cortical or compact bone, found in the mid-shaft ofload bearing bones like femur and tibia, is a brittle,porous biomaterial. Being a living tissue, the local mi-crostructure and porosity network of the bone evolvesin response to the mechanical stresses that the bone issubjected to, and this in turn modifies the local mechan-ical properties. Understanding the relationship betweenmicrostructure and mechanical properties is crucial forapplications such as extraction of bone grafts [1, 2], indesign of mechanically compatible implants [3–5] andporous scaffolds for bone tissue engineering [6, 7], to in-terpret loading history [8–10], to evaluate the effective-ness of chemical and physical therapeutical measures forbone healing [11, 12] etc. An important aspect of thisunderstanding is the development and testing of mod-els that incorporate microstructural features and predictmaterial properties such as failure strength, elastic mod-ulus, fracture paths, etc. Such models, if general enough,would also be of use in understanding failure behaviorof a wider class of brittle materials with a well definedporosity network, such as wood, rock [13, 14] etc.Compared to the detailed experimental characteriza-tion of the microstructure-property relationship of thedifferent microstructures [15–28], the number of predic-tive microscopic models that are able to reproduce char-acteristic features of the complex fracture processes in-volved are few in number. This is primarily because thefracture process in complex heterogeneous quasi-brittlematerials like bone involves multiple pores and micro-cracks that interact and evolve stochastically prior tofinal failure [29], which is difficult to capture using de-terministic models. Thus, modeling such materials us- ∗ [email protected] † [email protected] ‡ [email protected] ing classical fracture mechanics theory or standard fi-nite element method restricts the scope mostly to find-ing effective elastic behavior or to finding the resistanceto growth of a single macroscopic crack [30–32]. Sta-tistical models, like the random spring network model(RSNM), that approximate the continuum by a networkof springs with statistically distributed characteristics,are much better suited for studying fracture in such sys-tems. RSNM has been successful in providing an insightinto the role of disorder in the fracture behavior of hetero-geneous material systems with no preexisting crack [33],reproducing features like transition from brittle to non-brittle macroscopic response [34], avalanche size distri-butions [35] and qualitative [36, 37] and quantitative [38]features of fracture of composite materials. In the con-text of bone, simple one-dimensional models of parallelsprings with statistically distributed properties have beenused to simulate the characteristic quasi-brittle softeningseen in tension and bending [39] as well as in compres-sion of cortical bone [40]. Including spatial effects, three-dimensional RSNM has been used to model cancellousbone, spongier bone found near joints, to account forpercolation effects in the power-law variation of strengthwith porosity [41, 42]. However none of these studiestake into account the role of the structure of the porositynetwork. This is of particular importance since the struc-ture of the porosity network, and not merely the meanporosity, is known to be important. Microstructures withlarger mean porosity have sometimes higher compressivestrengths than those with lower porosity [26, 28, 43].Cortical bone has predominantly two distinct mi-crostructures – plexiform bone that is brick shaped andhas woven bone and vasculature sandwiched within regu-lar lamellae, and Haversian bone that has cylindrical sec-ondary osteons that run along the length of the bone [44].Under compression it has been found that mechanisms offailure are noticeably different between the microstruc-tures. In plexiform bone damage is localized on weak ra-dial planes resulting in prismatic fracture surfaces whilein Haversian bone the presence of osteons deflects thecrack paths and leads to meandering crack paths [28].To model the splitting fracture of plexiform bone undercompression, recently we developed a RSNM model thatincorporates the details of the porosity network of thebone [45]. In our porosity-based RSNM (see Sec. II formore details), the local spring constants were determinedbased on the experimentally obtained porosity networkof the plexiform bone. The effect of porosity network onthe macroscopic compressive response of the bones wasexamined under the simplifying assumption that the ma-terial properties like Young’s modulus, Poisson’s ratio,and compressive strength are independent of the bone.While the model reproduced the overall force deflectioncurves, qualitative fracture paths, and avalanche expo-nents reasonably well, the correlation between the pre-dicted and observed compressive strengths was not good.However, material properties are expected to be depen-dent on the bone, as well as the history of loading, andthe discrepancy between model predictions and experi-mental data was attributed to ignoring bone dependentmaterial properties. In this paper we introduce bone de-pendent material properties into the model by perform-ing a detailed parametric study of the model by varyingthe model parameters systematically. From experimentalenergy dispersive spectroscopy, we find that elastic mod-ulus is more or less independent of the bone. However, toobtain the experimental macroscopic response, we haveto use sample-dependent strain thresholds. We concludethat samples with similar mean porosity and mineraliza-tion perhaps have different material organization, leadingto variation in load bearing capacity. Further, the roleof porosity network is examined by comparing the pre-dictions with results from homogenized distribution ofequivalent porosities. While homogenization leads to al-teration of the failure paths and strengths, it leaves theavalanche exponents largely unchanged.The remainder of the paper is organized as follows.In Sec. II, we describe in detail the model and howthe porosity based RSNM is obtained from the experi-mental CT scan images. In Sec. III, results from para-metric studies of the model parameters are discussed interms of their effect on failure paths and load bearingcapacity. Bone-specific properties are iteratively esti-mated to match with experiments and predictions arecompared with simulations using bone-independent pa-rameters. The exponents of avalanche size distributionsare also presented for networked and homogenized poros-ity. II. MODEL
In this section, we describe the formulation of thetwo dimensional spring network model for plexiformbone [45], and the generalizations that will be studiedin this paper. The incorporation of the experimentallyobtained CT scan images into a porosity based two di- mensional RSNM involves several steps that are summa-rized in the flowchart shown in Fig. 1. We describe eachof the steps below.
Plexiform bone3-dCT scan images10 representativevolumes2 dimensionalporosity network 2-d RSNMDetermine E, as afunction of porosityDetermine spring constantsas a function of porosityE, determinedspring constantsPorosity based RSNM
FIG. 1. Flow chart detailing the steps involved in developingporosity-based RSNM. (c) (d)(a)
L : Longitudinal R : Radial T : Tangential
L RT (b)
FIG. 2. (a) Cubic sample of size 5mm × × We use the experimental data, reported earlier inMayya et al. [45] in which cubical samples from the an-terior section of the mid-diaphysis of bovine femur [seeFig. 2(a)] were tested in compression along the length ofthe bone. CT scan images of the samples were obtainedbefore and after compression failure [45]. Sample CTscan images of the pre- and post-compression samplesare reproduced in Fig. 2(c) and (d). The samples werechosen from regions that have plexiform microstructure,as evident from the brick-shaped layered structure seen in
L : LongitudinalR : RadialT : Tangential
RL T L R (b)(a) (d)(c)
FIG. 3. (a) Stack of representative domains, each spanningover 500 µm obtained from porosity analysis. (b) Spring stiff-ness are based on the %porosity at the node (pixel). (c) Atypical representative domain and (d) corresponding homog-enized configuration for sample III. the optical micrograph shown in Fig. 2(b). The CT scanproduces 500 slices, each slice corresponding to 10 µm .The grayscale image in each slice is converted to a binaryimage (0 or 1) by setting a threshold determined fromthe valleys of the probability distribution function of thegrayscale. Here 1 corresponds to a pore and 0 to a non-pore. The 500 slices are now divided into 10 representa-tive volumes, each consisting of 50 slices, correspondingto 500 µm thickness [see Fig. 3(a)]. The thickness waschosen such that it is much larger than the mean poresize but less than the inter pore distance. A representa-tive volume was mapped onto a two-dimensional squarenetwork where the porosity at a given location was ob-tained as the average of the binary data over the 50 slices.Thus, for each location, a porosity value varying from 0to 100% was obtained. A typical representative volume,shown in Fig. 3(c), is fairly dense, having mean poros-ity lower than 5%, and consists of fine pores that areevenly distributed and extend across the height of thebone sample while being slightly inclined to the longitu-dinal direction. The two-dimensional porosity network,thus obtained from the CT scans, is an input for theRSNM model constructed in the following manner.A representative volume is modeled using a 150 × θ ijk shown in Fig. 4(b)] from the initial value of π/ V , has contribu-tions from extension of springs, distortion of angle be-tween springs and also a repulsive contact force modeledas Hertzian contact between particles: V = X h ij i k ij δr ij + X h ijk i c ijk δθ ijk + α X mn ( d − r mn ) Θ( d − r mn ) , (1)where h ij i denotes those pairs of particles that are con-nected by linear springs, the spring constant being k ij and extension being δr ij . c ijk is the spring constant thatresists the distortion in angle, δθ ijk , between triad ofparticles h ijk i . The contact between any two particles m and n is initiated by a Heaviside function, Θ( x ), onlywhen the inter particle distance, r mn , is less than theparticle diameter, d . The elastic contact force parameter α is a material constant [46]. If the springs constantsof horizontal/vertical, diagonal and torsional springs are2 k , k and c respectively, then the system is isotropic withelastic modulus, E , and Poisson’s ratio, ν , given by [47] E = 8 k ( k + ca − )3 k + ca − ; ν = k − ca − k + ca − , (2)where a is the lattice spacing. We account for the elas-tic compliance of the testing machine by connecting thesites in the top and bottom rows of the spring networkto springs whose compliance match with that of the ma-chine used in experiments. In the simulations, the dis- FIG. 4. (a) Schematic diagram of the spring network model.(b) Unit cell of the network showing the linear springs at-tached to a site. Also shown is an example of an angle θ ijk whose distortion is resisted by a bending spring. placements are applied incrementally from the top. Forevery increment, the system is equilibrated to its mini-mum energy configuration by numerically integrating theequations of motion in terms of the position vectors, r i : d r i dt = −∇ r i V − γ d r i dt , (3)using velocity-verlet algorithm [48]. The damping coef-ficient γ dissipates energy and brings the system to itsminimum energy configuration [35]. For each incrementin the applied downward displacement, after equilibra-tion, if any spring is stretched beyond its correspondingthreshold strain ǫ f , it and the bending springs associatedwith it are broken. The system is re-equilibrated till nofurther breakage takes place within the increment.To incorporate the experimentally obtained porositynetwork into the model, the spring constants and break-ing strain thresholds are assigned in accordance with thelocal porosity. To find the porosity dependent materialbehavior, first the effective behavior of a porous net-work that has homogeneous matrix properties is eval-uated. The response is evaluated under tension as thebone samples experience tensile transverse strains thatlead to splitting under longitudinal compression. ARSNM of given porosity, P , is constructed by remov-ing voids [of size 50 µm for 0 ≤ P ≤
5% and 100 µm for 5% ≤ P ≤ E ( P ), ν ( P ) and ǫ f ( P ) are obtained. The scalingof these curves depends on the values of the spring con-stants and breaking thresholds of the homogeneous ma-trix. These values are chosen such that at 4% porosity,the known experimental average macroscopic response ofbone is reproduced [26]. Next, for simulation of frac-ture in bone samples, the individual spring characteris-tics of the RSNM are assigned on the basis of the localporosity data from experiments using the correspondingcalibrated E ( P ), ν ( P ) and ǫ f ( P ) in Eq.(2) to achieve aporosity dependent RSNM.To simulate the fracture process under compressionloads, in a previous work, we developed a random springnetwork model where the matrix properties were takento be independent of bone and the failure was pre-dicted accounting only the structure of the porosity net-work [45]. While the porosity based RSNM was shownto capture the characteristic features of the quantitativemacroscopic response as well as the qualitative failurepaths during the fracture process, quantitative correla-tions were poor. To improve the predictive capability ofthe model, here, we evaluate the use of matrix propertiesthat are specific to a bone. For validation, we use ex-perimental data [45] on a total of six samples that wereharvested from three different bovine femurs and referredto in the remainder of the paper as sample I and II frombovine-1, sample III and IV from bovine-2, and sampleV and VI from bovine-3. III. RESULTS AND DISCUSSION
Initially, a parametric study is performed to evaluatethe sensitivity of the predictions to model parameters– elastic modulus and failure strain. In the study, thematrix properties of the homogeneous porous networkare varied by ±
10% of the bone independent parametersused in Mayya et al. [45] that are considered here as basevalues. Based on each combination of matrix properties,first the porosity dependent elastic modulus, Poisson’sratio and failure strain strain, as shown in Fig. 5, areevaluated from a porous network, as described in Sec II,for a range of overall porosity. For the cases when there Porosity % +10% E+10% f E l a s t i c m odu l u s ( G P a ) -3 F a il u r e s t r a i n +10% E+10% f Porosity % P o i ss on ’ s r a t i o +10% E+10% f Porosity % (a) (b)(c) Base values
Elastic modulusFailure strain 18 G Pa0.0122
FIG. 5. Variation of (a) elastic modulus (b) Poisson’s ratioand (c) failure strain with porosity for different modulus andfracture strain of the matrix. is 10% increase in E and 10% ǫ f from the base value,as expected stiffer matrix results in stiffer macroscopicresponse and affects only the effective elastic behaviorwhile change in the failure threshold correspondingly haseffect only on the effective limiting strain and not onthe elastic behavior. We note that, for thermodynamicstability, the spring constants k and c must be positivein Eq. (2), which bounds the Poisson ratio from aboveby 1/3. A. Effect on failure paths
Using a porosity dependent RSNM, the effect of themodel parameters: elastic modulus and failure strain onthe macroscopic response is shown for a typical sample(sample III) in Fig. 6. For clarity only the predictionsbased on the base value and 10% independent increase ineach of the model parameters is shown. The increase in E as well as ǫ f is shown to increase the overall load bearingcapacity, obtained from the macroscopic response of thesample that has been averaged over the 10 representativedomains, by approximately 8%. Most of the individualsets exhibit a similar increase in load bearing capacity asevident in Fig. 6(b) that shows the macroscopic responseof a typical representative domain. The correspondingfailure paths are presented in Figs. 6(c)-(e). The firstbonds to fail for any combination of parameters are fromthe regions of relatively higher porosity (10 − E could introduce differences in load transfer displacement (mm) f o r c e ( k N ) +10% E+10% f displacement(mm) f o r c e ( k N ) +10% E+10% f (b) (a) (d)(c) (e) FIG. 6. (a) Macroscopic response of sample III and (b) fora typical representative domain. The corresponding failurepaths at the first load drop with (c) base values (d) 10% E and (e) 10% failure strain as input parameters. paths and thereby leads to some variations in the split-ting paths. On rare occasions (one representative domainout of the 10 considered here), the damage at the criticaldefect localizes in a manner leading to the formation of alarge single crack resulting in splitting fracture at lowerstrains. B. Effect on macroscopic response
To develop a comprehensive understanding of the effectof model parameters on the macroscopic response undercompression, simulations were performed for all the 6samples, using the structure of their respective porositynetworks. Figure 7(a) shows the predicted load bear-ing capacity for elastic modulus ranging between ± P ea k f o r c e ( k N ) P ea k f o r c e ( k N ) bovine 1bovine 2bovine 3bovine 1bovine 2bovine 3base value -10% E +10% Ebase value -10% E +10% E Porosity networkHomogenised dist. (a)(b)
FIG. 7. Variation of load bearing capacity with mean porositylevels and ±
10% variation in elastic modulus, E , accounting(a) for porosity network and (b) for homogenized distributionof porosity. P ea k f o r c e ( k N ) P ea k f o r c e ( k N ) bovine 1bovine 2bovine 3bovine 1bovine 2bovine 3base value -10% f +10% f base value -10% f +10% f Porosity networkHomogenised dist. (a)(b)
FIG. 8. Variation of load bearing capacity with mean porositylevels and ±
10% variation in failure strain, ǫ f accounting (a)for porosity network and (b) for homogenized distribution ofporosity. as a result of increase in elastic modulus and therebyresulting in differences in crack paths as discussed earlier.Homogenization of porosity network results in an overallincrease in predictions of peak force values as shown inFigs. 7(b) and 8(b). As expected, the increase in the loadbearing capacity of a homogeneously porous network ismore comparable for the increase in either parameters:22% for increase in E and 18% for increase in ǫ f . C. Use of bone-specific properties
Element counts (Ca) P ( E l e m en t c oun t s ) bovine 1bovine 2bovine 3 FIG. 9. Probability density of element counts from EDS linescans of bone samples.
To determine the appropriate bone-specific matrixproperties for each bone, we first examine the elastic be-havior experimentally. The elastic properties of corticalbone are known to be influenced by the mineral den-sity [49, 50] which are typically characterized by differ-ent techniques including electron micro-probe techniquessuch as energy dispersive spectroscopy (EDS) [51, 52].Here, we estimate the differences in the elastic behaviorof the bone matrix between bones by performing EDSline scans using Quanta-200 FEI scanning electron mi-croscope. The sample measurements are grouped on thebasis of the bovine femurs from which they are harvested.The probability density of calcium counts obtained foreach bone, shown in Fig. 9, are bell shaped curves withsignificant overlap. The expectation values for the countsfor the bones are with in 4% of the overall average countwhich is possibly because the bones are from healthy an-imals of similar age. Thus, for further analysis the elasticbehavior of all the bones is taken to be the base value.To incorporate the differences in the fracture behaviorof bone matrix, the threshold failure strain for each bonewas estimated iteratively till it compared well with theexperimental data and is presented in Table I.For a given sample, 5 random realizations for eachof the 10 representative networks were performed. The
TABLE I. Bone specific input parametersBone Sample no. Parameter for simulationsbovine - 1 I E and 1 . × ǫ f IIbovine - 2 III E and 0 . × ǫ f IVbovine - 3 V E and 1 . × ǫ f VI* Base values : E = 18 GP a and ǫ f = 0 . strain threshold for each bond was taken from a Gaus-sian distributions with mean as in Table I and 5% stan-dard deviation to account for heterogeneity at lengthscales much smaller than the lattice parameter. Themacroscopic response, averaged over all the 50 simula-tions per sample are compared with the experimentaldata in Fig. 10. The characteristic features of the macro-scopic response such as the initial linear elasticity, themaximum load bearing capacity, multiple smaller eventsprior to final failure etc. are well reproduced. It is re-markable to note that for each bone, of the two samplestested, same model parameters are effective in predictingresponse that compares closely with experimental datafor both the samples.A comparative summary of the maximum load takenby any sample as per simulation based on bone specificproperties, base values and as from experimental datais presented in Fig. 11(a). Maximum loads as predictedusing base values for all samples, as seen earlier, tendto have a monotonic decrease with increasing porositywhich is in contradiction to the experimental data wherethere is no such consistent trend with respect to porosity.While the samples from bovine 2, bovine 1 and bovine3 appear to be in order of increasing porosity, the loadtaken by them has no direct porosity dependence, in factthe highest load is taken by the most porous of the sam-ples. However, when the bone specific properties are in-corporated in the simulations a significant improvementin predictions are seen for both samples of each bone asillustrated in Fig. 11(b). The concordance correlation co-efficients from the simulations with bone-specific param-eters (0.79) are also significant as compared to base val-ues (-0.12). The differences in load bearing capacity be-tween samples within a bone that are extracted from thesame anatomical site can be attributed primarily to theinherent differences in the porosity network. However,between samples from different bones with apparentlysimilar elastic behavior, as per the model of the presentstudy, the fracture behavior prediction requires 60-70%difference in matrix failure strain threshold. These in-ferences require validation through detailed characteri-zation of failure properties at smaller length scales. f o r c e ( k N ) SimulationExperiment f o r c e ( k N ) SimulationExperiment f o r c e ( k N ) SimulationExperiment f o r c e ( k N ) SimulationExperiment f o r c e ( k N ) SimulationExperiment0 0.2 0.4 0.602468 f o r c e ( k N ) SimulationExperiment (a) (c) (e)(b) (d) (f)displacement (mm) displacement (mm) displacement (mm)displacement (mm) displacement (mm) displacement (mm)
FIG. 10. Comparison of macroscopic response from simulations and experiments of bovine samples, where (a) is for sample I,(b) is for sample II and so on. experiments base value bone specific0.6 0.7 0.8 0.9 1 1.1 1.2345678 Porosity % P ea k f o r c e ( k N ) S i m u l a t i on : f o r c e ( k N ) (b)(a) bovine 1bovine 2bovine 3 bovine 1bovine 2bovine 3 FIG. 11. Load bearing capacity from experiments and simula-tions: (a) variation with porosity and (b) correlation betweenexperiment and simulations.
D. Avalanche size distribution
One of the quantities that has been used in the liter-ature to characterize fracture in heterogeneous media isthe avalanche size distribution which measures the incre-mental response of the system to incremental increasesin the external loading. In experiments, the response ismeasured through energy of acoustic emissions, E , thatoccur during the fracture process. In simulations, it ismeasured by the number of events (broken springs), s ,that occur per increment of strain. The probability dis-tributions for these quantities are known to be powerlaws, implying the absence of a typical avalanche sizethat is independent of system size. The exponents char-acterizing the acoustic emission and avalanche size dis- tribution may be related to each other. Let P E ( E ) and P s ( s ) denote the respective distributions. Asymptoti-cally, they behave as P E ( E ) ∼ E − τ E and P s ( s ) ∼ s − τ s .However, it is known that E ∝ s [35]. Using therelation P E ( E ) dE = P s ( s ) ds , arising from the con-servation of probability, it is straightforward to obtain τ E = (1 + τ s ) /
2. Remarkably, these exponents are quiteuniversal and independent of details of material or mod-eling. For example, τ E is known to be 1 . − . τ E = 1 . − . τ s from our simulations, and ask whether and howthe exponent τ s depends on the choice of bone-dependentmaterial properties, or the porosity network.The avalanche size distribution for a typical sample isshown in Fig. 12. For each of the representative domains,for all 5 realizations, avalanche sizes are obtained for allincrements in displacement. As can be seen from Fig. 12,the distribution is insensitive to whether the bone-specificparameters or base values are used. This is understand-able as the bone specific parameters involve only a changein ǫ f , and this does not appreciably change the load dis-tribution in the network during compression. For thesample shown in Fig. 12, we find τ s ≈ .
06 correspondingto τ E = 1 .
53. For the other samples, the exponents thatwe obtain are shown in Table II. The exponent values for τ E lie in the range 1 . − .
5, consistent with the valuesof 1 . − . −8 −6 −4 −2 Avalanche size, s P ( s ) base value : networkbone specific : networkbone specific : homogenized porosity τ s = 2.06 FIG. 12. Avalanche distribution P ( s ) for sample III fromsimulations. values of the exponent with the porosity of the sample.In Fig. 12, we also present the avalanche size distribu-tion associated with an equivalent homogenized porousnetwork. Surprisingly, the avalanche exponent does notchange noticeably. However, the distribution for smallavalanche sizes is different, reflecting the contrast in thenetwork at small scales. TABLE II. Power law exponents, τ s and τ E for all samples.Bone Sample no. τ s τ E bovine - 1 I 1.93 1.47II 1.80 1.40bovine - 2 III 2.06 1.53IV 1.93 1.47bovine - 3 V 1.69 1.35VI 1.93 1.47 IV. SUMMARY AND CONCLUSIONS
Complex fracture processes in porous brittle materi-als are controlled by several factors. Using a porositydependent RSNM, in the present study, we examine theindividual roles of overall mean porosity, the networkingof the porosity, the elastic behavior of non-porous matrixand the failure behavior of the non-porous matrix on thecompressive strength of dry plexiform bone. As per themodel, increasing mean porosity results in reduced com-pressive strength. Porosity network structure, as seen inplexiform bone, reduces the compressive strength furtherby upto 30%. While the initiation of the fracture process is typically at regions of highest porosity, the crack path-ways are predominantly controlled by the connectivityof the porosity network. Among the multiple competingporosity pathways, the damage localizes in only a fewas driven by the load distribution which in turn is influ-enced by the elastic properties of the model matrix andnot by the failure strain threshold of the matrix. How-ever, the load bearing capacity increases with increase ineither elastic modulus or failure strain threshold of thematrix.To validate the porosity-based RSNM model withavailable experimental data, bone-specific propertieswere applied. Of the six samples (two each from three dif-ferent bones), the elastic properties were found to differminimally as reflected in the mineral composition deter-mined by EDS. Assuming elastic similarity, bone-specificfailure strain thresholds for the ideal matrix were esti-mated to be within 60% of each other. It is to be notedthat the compressive strengths differ between the twosamples of the same bone and yet, identical bone-specificmodel parameters achieve excellent concordance correla-tion between experiment and simulations for both sam-ples of each bone. To provide an accurate input for themodel, it would therefore be of importance to develop anexperimental scheme that characterizes the strain thresh-old of the matrix at small length scales.We find that the avalanche size exponents are inde-pendent of the bone-dependent parameters and in therange of the experimental values obtained for porcinebone [58]. Surprisingly, the exponent is also independentof the structure of the porosity network as the avalanchesize distribution exponent for the equivalent homoge-nized network has similar value. This may be the reasonwhy the avalanche exponent is universal and has simi-lar numerical value for many different kinds of material.At the same time, the exponents obtained in this paper(1 . − .
05) are significantly different from that for homo-geneous RSNM under tension ( ∼ .
5) [35]. Whether thisdifference is due to the difference in the nature of load-ing, or it is because of the strong correlation between thevalues for the spring constant and its strain threshold, asin the current study, remains to be answered.It would be interesting to study the predictions of themodel for Haversian bone which has a distinctly differ-ent microstructure from plexiform bone. Unlike the lay-ered plexiform bone, both porosity network character-istics and material constitution are different that mayinfluence the load bearing capacity. A similar analysisfor Haversian bone is part of ongoing studies.
ACKNOWLEDGMENTS
We thank Prof. Krishnan Balasubramanian, Centrefor Non-Destructive Evaluation and Prof. V. S. Sarma,Scanning Electron Microscopy Lab, IIT Madras for pro-viding access to the experimental facilities. [1] A. Pearce, R. Richards, S. Milz, E. Schneider, S. Pearce, et al. , Eur. Cell Mater. , 1 (2007).[2] M. Conward and J. Samuel, J. Mech. Behav. Biomed.Mater. , 525 (2016).[3] K. A. Hing, Int. J. Appl. Ceram. Tec. , 184 (2005).[4] A. Bansiddhi, T. Sargeant, S. I. Stupp, and D. Dunand,Acta Biomater. , 773 (2008).[5] F. Libonati and L. Vergani, Compos. Struct. , 188(2016).[6] Z. Chen, D. Li, B. Lu, Y. Tang, M. Sun, and S. Xu,Scripta Mater. , 157 (2005).[7] X. Liu, M. N. Rahaman, and Q. Fu, Acta Biomater. ,4889 (2013).[8] K. E. Keenan, C. S. Mears, and J. G. Skedros, Am. J.Phys. Anthropol. , 657 (2017).[9] M. Zedda, M. R. Palombo, D. Brits, M. Carcupino,V. Sath´e, A. Cacchioli, and V. Farina, Zoomorphology ,1.[10] J. Mitchell, L. J. Legendre, C. Lef`evre, and J. Cubo,Zoology (2017).[11] L. E. Claes, H.-J. Wilke, and H. Kiefer, J. Biomech. ,1377 (1995).[12] T. P. Gocha and A. M. Agnew, J. Anat. (2015).[13] M. A. Meyers and K. K. Chawla, Mechanical behaviorof materials , Vol. 2 (Cambridge Univ. Press Cambridge,2009).[14] A. Shojaei, A. D. Taleghani, and G. Li, Int. J. Plast. ,199 (2014).[15] D. T. Reilly and A. H. Burstein, J. Biomech. , 393(1975).[16] S. F. Lipson and J. L. Katz, J. Biomech. , 231 (1984).[17] J. L. Katz, H. S. Yoon, S. Lipson, R. Maharidge, A. Me-unier, and P. Christel, Calcif. Tissue Int. , S31 (1984).[18] R. Shahar, P. Zaslansky, M. Barak, A. Friesem, J. Currey,and S. Weiner, J. Biomech. , 252 (2007).[19] R. Martin and J. Ishida, J. Biomech. , 419 (1989).[20] R. Martin and D. Boardman, J. Biomech. , 1047(1993).[21] T. L. Norman, D. Vashishth, and D. B. Burr, J. Biomech. , 309 (1995).[22] G. C. Reilly and J. D. Currey, J. Exp. Biol. , 543(1999).[23] J. H. Kim, M. Niinomi, T. Akahori, J. Takeda, andH. Toda, JSME Int J., Ser. A , 472 (2005).[24] J. Kim, M. Niinomi, T. Akahori, and H. Toda, Int. J.Fatigue , 1039 (2007).[25] D. Carter, W. C. Hayes, and D. J. Schurman, J.Biomech. , 211 (1976).[26] S. Li, E. Demirci, and V. V. Silberschmidt, J. Mech.Behav. Biomed. Mater. , 109 (2013).[27] A. Mayya, A. Banerjee, and R. Rajesh, Sci. Rep. (2013).[28] A. Mayya, A. Banerjee, and R. Rajesh, Mater. Sci. Eng.,C , 454 (2016).[29] Z. P. Baˇzant, Proc. Natl. Acad. Sci. U.S.A. , 13400(2004).[30] E. Budyn and T. Hoc, Int. J. Numer. Meth. Eng. , 940 (2010).[31] A. A. Abdel-Wahab, A. R. Maligno, and V. V. Silber-schmidt, Comput. Mater. Sci. , 128 (2012).[32] K. R. Brown, S. Tarsuslugil, V. Wijayathunga, andR. Wilcox, J. R. Soc. Interface , 20140186 (2014).[33] M. J. Alava, P. K. Nukala, and S. Zapperi, Adv. Phys. , 349 (2006).[34] W. Curtin and H. Scher, J. Mater. Res. , 535 (1990).[35] P. Ray, Comput. Mater. Sci. , 141 (2006).[36] C. Moukarzel and P. Duxbury, J. Appl. Phys , 4086(1994).[37] C. Urabe and S. Takesue, Phys. Rev. E , 016106(2010).[38] D. Boyina, T. Kirubakaran, A. Banerjee, and R. Velmu-rugan, Mech. Mater. , 64 (2015).[39] D. Krajcinovic, J. Trafimow, and D. Sumarac, J.Biomech. , 779 (1987).[40] J. Schwiedrzik, R. Raghavan, A. B¨urki, V. LeNader,U. Wolfram, J. Michler, and P. Zysset, Nat. Mater. ,740 (2014).[41] G. H. Gunaratne, C. S. Rajapaksa, K. E. Bassler, K. K.Mohanty, and S. J. Wimalawansa, Phys. Rev. Lett. ,068101 (2002).[42] C. S. Rajapakse, J. S. Thomsen, J. S. E. Ortiz, S. J.Wimalawansa, E. N. Ebbesen, L. Mosekilde, and G. H.Gunaratne, J. Biomech. , 1241 (2004).[43] H. Tr¸ebacz, A. Zdunek, J. Cybulska, and P. Pieczywek,Australas. Phys. Eng. Sci. Med. , 43 (2013).[44] J. D. Currey, Bones: structure and mechanics (Princetonuniversity press, 2002).[45] A. Mayya, P. Praveen, A. Banerjee, and R. Rajesh, J.Roy. Soc. Interface , 20160809 (2016).[46] L. D. Landau, A. Kosevich, L. P. Pitaevskii, and E. M.Lifshitz, Theory of elasticity (Butterworth, 1986).[47] L. Monette and M. Anderson, Modell. Simul. Mater. Sci.Eng. , 53 (1994).[48] L. Verlet, Phys. Rev. , 98 (1967).[49] M. B. Schaffler and D. B. Burr, J. Biomech. , 13 (1988).[50] J. D. Currey, J. Biomech. , 131 (1988).[51] K. ˚Akesson, M. Grynpas, R. Hancock, R. Odselius, andK. Obrant, Calcif. Tissue Int. , 236 (1994).[52] R. D. Bloebaum, J. Skedros, E. Vajda, K. N. Bachus,and B. Constantz, Bone , 485 (1997).[53] A. Petri, G. Paparo, A. Vespignani, A. Alippi, andM. Costantini, Phys. Rev. Lett. , 3423 (1994).[54] A. Garcimartin, A. Guarino, L. Bellon, and S. Ciliberto,Phys. Rev. Lett. , 3202 (1997).[55] A. Guarino, S. Ciliberto, A. Garcimartın, M. Zei, andR. Scorretti, Eur. Phys. J. B , 141 (2002).[56] C. Maes, A. Van Moffaert, H. Frederix, and H. Strauven,Phys. Rev. B , 4987 (1998).[57] S. Pradhan, A. M. Stroisz, E. Fjær, J. F. Stenebr˚aten,H. K. Lund, and E. F. Sønstebø, Rock Mech. Rock Eng. , 2529 (2015).[58] J. Bar´o, P. Shyu, S. Pang, I. M. Jasiuk, E. Vives, E. K.Salje, and A. Planes, Phys. Rev. E93