Evaluating the macroscopic yield behaviour of trabecular bone using a nonlinear homogenisation approach

Computational homogenisation approaches using high resolution images and ﬁ nite element (FE) modelling have been extensively employed to evaluate the anisotropic elastic properties of trabecular bone. The aim of this study was to extend its application to characterise the macroscopic yield behaviour of trabecular bone. Twenty trabecular bone samples were scanned using a micro-computed tomography device, converted to voxelised FE meshes and subjected to 160 load cases each (to de ﬁ ne a homogenised multiaxial yield surface which represents several possible strain combinations). Simulations were carried out using a parallel code developed in-house. The nonlinear algorithms included both geometrical and material nonlinearities. The study found that for tension-tension and compression-compression regimes in normal strain space, the yield strains have an isotropic behaviour. However, in the tension-compression quadrants, pure shear and combined normal-shear planes, the macroscopic strain norms at yield have a relatively large variation. Also, our treatment of clockwise and counter-clockwise shears as separate loading cases showed that the differences in these two directions cannot be ignored. A quadric yield surface, used to evaluate the goodness of ﬁ t, showed that an isotropic criterion adequately represents yield in strain space though errors with orthotropic and anisotropic criteria are slightly smaller. although the isotropic yield surface presents itself as the most suitable assumption, it may not work well for all load cases. This work provides a comprehensive assessment of material symmetries of trabecular bone at the macroscale and describes in detail its macroscopic yield and its underlying microscopic mechanics.


1.
Introduction Exponential growth of older population implies that problems associated with deteriorated mechanical capabilities of bone need urgent attention. Computational modelling to examine the mechanical response of musculoskeletal systems requires the mechanical behaviour of bone to be defined satisfactorily (Pankaj, 2013). A continuum description of bone that can be related to its microstructure and includes its anisotropy and its yield behaviour will go a long way in predicting failure of bone and bone-implant systems.
The macroscopic elastic behaviour of bone has been mostly modelled using isotropic linear elasticity. Often, bone macroscopic properties are assumed to be homogeneous with separate elastic properties being assigned to cortical and trabecular bone (Completo et al., 2009;Conlisk et al., 2015). Sometimes, subject specific macroscopic elastic properties are assigned using computed tomography (CT) scans, which permit inhomogeneity in the material properties on the basis of CT attenuations (Helgason et al., 2008;Tassani et al., 2011). That said, since CT attenuations can only provide scalar values, assumption of isotropy needs to be made. However, it is well recognised that the macroscopic behaviour of bone is not isotropic. For trabecular bone, which resembles open cell foams, the anisotropy is largely a consequence of its anisotropic microarchitecture Turner et al., 1990). An ultrasonic approach proposed by van Buskirk et al. (1981) was shown to provide a good approximation of nine orthotropic elastic constants if a heterogeneity correction were included. In general, experimental mechanical techniques are unable to provide the complete stiffness tensor at the resolution required for modelling (Odgaard et al., 1989).
Image based computational approaches have been successfully applied for the evaluation of the macroscopic stiffness tensors (Donaldson et al., 2011;van Rietbergen et al., 1995). In these, micro-CT (or micro-magnetic resonance imaging) scans of bone are converted into high resolution 3D finite element (FE) meshes, with a detailed geometry of its microstructure. The solid phase (or bone tissue) is assigned isotropic elastic properties and the volume element (VE) is then computationally subjected to six strain/stress states (three normal and three shear). The response enables evaluation of the full macroscopic elastic stiffness tensor using the standard mechanics methodology (van Rietbergen et al., 1996). Previous studies have extensively employed these homogenisation approaches, and relationships between stiffness and microarchitectural indices (volume fraction and fabric tensor) have also been established (Cowin, 1986;Turner and Cowin, 1987;Turner et al., 1990;Zysset and Curnier, 1995).
While modelling bone as an elastic material may be adequate for a few applications, a significant proportion of applications requires evaluation of post-elastic response, e.g. to evaluate implant loosening resulting in its failure. Many studies still continue to employ elastic analyses to predict arbitrarily post-elastic behaviour (Falcinelli et al., 2014).
Both stress-and strain-based criteria have been used to describe the macroscopic yield surface of bone (Keaveny et al., 1994;Keller, 1994;Kopperdahl and Keaveny, 1998). In recent years a consensus appears to be emerging that strainbased criteria are easier to apply as trabecular bone behaviour in this space is "more isotropic" and density independent than in stress space (Bayraktar et al., 2004;Chang et al., 1999;Pankaj and Donaldson, 2013). There is also now some evidence to suggest that failure of bone is strain-controlled rather than stress-controlled (Nalla et al., 2003) (Bayraktar et al., 2004;Sanyal et al., 2015;Wolfram et al., 2012). All these studies used a simple bilinear criterion to represent the solid phase of bone. Wolfram et al. (2012) used a limited number of load cases which can lead to loss of information on physiologically possible complex load cases, while both Sanyal et al. (2015) and Wolfram et al. (2012) made a priori assumptions with regard to macroscopic yield surface symmetries; the former assumed it to be transverse isotropic and the latter orthotropic. Nanoindentation experiments on bone suggest that the solid phase of bone has a pressure-dependent yield surface (i.e. its yielding depends on hydrostatic stress), which arises because of bone's cohesive-frictional behaviour (Tai et al., 2006). Due to this reason, bone tissue (or the solid phase) can be modelled using classical criteria, such as Mohr-Coulomb or Drucker-Prager (Carnelli et al., 2010;Tai et al., 2006). On the macroscale, high density bone is prone to tissue yielding, while low density bone is likely to fail via a mixture of large deformation failure mechanisms and tissue yielding (Bevill et al., 2006;Morgan et al., 2004;Stolken and Kinney, 2003  2.

Material and methods
The tensorial notation used in this study largely follows the notation used by Schwiedrzik et al. (2013). A first-order tensor (or vector) is denoted by a lowercase bold letter (e.g. m), a second-order tensor is denoted by an uppercase bold letter (e.g. A) and a fourth-order tensor is denoted by a doublebarred uppercase letter (e.g. A). The tensor operations which will appear throughout the text are: single contraction of two second-order tensors, e.g. AB (or A ik B kj in indicial notation); double contraction of two second-order tensors, e.g. A : B ðA ij B ij Þ; double contraction of a fourth-order tensor and a second-order tensor, e.g. A : B ðA ijkl B kl Þ; double contraction of a second-order tensor and a fourth-order tensor, e.g. B : A ðB ij A ijkl Þ; tensor product of two first-order tensors, e.g. m m ðm i m j Þ; tensor product of two second-order tensors, e.g. A B A ij B kl À Á ; and the symmetric tensor product of two second-order tensors, e.g.

Imaging and finite element meshing
Ten trabecular bone specimens were extracted from bovine trochanters and femoral heads (young cattle, o2.5 years old). Coarse micro-CT images of the central femoral head and trochanter regions were taken for three specimens prior to coring. All subsequent samples were then cored with respect to the visually ascertained trabecular directions. The extracted cylindrical specimens had a diameter of 10.7 mm and a length of 29.9 mm. Diamond-tipped cores (Starlite Industries, Rosemont PA, USA) were used in the extraction of the specimens and the top and bottom edges of the cores were cut with a slow speed saw (Isomet 1000, Buehler, Düsseldorf, Germany) by using a diamond wafering blade designed for bone. All these procedures were performed under constant irrigation to avoid excessive abrasion and overheating. The specimens were submerged into phosphate buffered saline and scanned using micro-CT (Skyscan 1172, Bruker, Zaventern, Belgium) with a resolution of 17.22 mm. The scanning parameters were 94 kV, 136 mA and 200 ms integration time; 4 scans in 720 equiangular radial positions which were averaged. The grey scale images were binarised with an automatic thresholding script (Gomez et al., 2013).
Twenty virtual cubes of 5 mm length were extracted from the scanned and segmented cylinders. The Mean Intercept Length (MIL) fabric tensor (Harrigan and Mann, 1984) was evaluated using BoneJ (Doube et al., 2010) and then used to align the coordinate axes of the images with the eigenvectors of the fabric. This approach has been employed in a recent study (Wolfram et al., 2012). After the 5 mm cubes were cropped, the alignment was rechecked to ensure that there was no misalignment larger than 81 (Sanyal et al., 2015;Wolfram et al., 2012). MIL is known to approximate the elastic orthotropic directions of trabecular bone (Odgaard, 1997). By undertaking such alignment there was an expectation that these axes may also represent the orthotropic directions of the yield criterion if the criterion was orthotropic. The bone tissue volume (solid phase) over total volume ratios (BV/TV) had a range from 13.7% to 30.3%. The degree of anisotropy of these cubes, which is the ratio between the largest and smallest eigenvalues of the fabric tensor, ranged from 1.52 to 3.86. A 5 mm length Volume Element (VE) which has been previously considered appropriate to capture the features of trabecular bone (Harrigan et al., 1988;Sanyal et al., 2015;van Rietbergen et al., 1995) was employed for all simulations.

2.2.
Constitutive model and computational procedure The twenty specimens were meshed using a voxelised mesh, where every voxel corresponds to a trilinear hexahedron, with an in-house developed script that meshes in parallel using the Message Passing Interface (MPI) (Forum, 1994). The meshing procedure was performed using 30 cores on a cluster at The University of Edinburgh, which is called Eddie (Edinburgh Compute and Data Facilities, ECDF). The largest mesh had around 9 million nodes. FE models for two of the twenty specimens are shown in Fig. 1.
The solid phase was modelled as an isotropic elastoplastic material. It is important to mention that trabecular bone at the tissue level is actually transverse isotropic or orthotropic (Hellmich et al., 2004;Malandrino et al., 2012;Wolfram et al., 2010). However, as pointed out by Cowin (1997), there is little to no error in assuming tissue isotropy. This is because the trabecula are composed of laminated material about their axes, which implies transverse isotropy or orthotropy; since the axis of the trabecula is the same as the loading axis, a beam made of orthotropic material can be reduced to a beam made of isotropic material.
The elastic regime was modelled using Hencky hyperelasticity, which restricts this material model to isotropic, with a Poisson's ratio of 0.3 and a Young's Modulus of 12700 MPa (Wolfram et al., 2012). A quadric yield surface (Schwiedrzik et al., 2013) given by was employed, where τ is the Kirchhoff stress and G and G are, respectively, a fourth-order tensor and a second-order tensor defined by and and I is the second-order unit tensor, s þ 0 and s À 0 are the tensile and compressive yield stresses, respectively, and ζ 0 is an interaction parameter. Eq.
(1) approximates a Drucker-Prager criterion when ζ 0 ¼ 0.49. Recent studies using indentation tests on bone tissue have suggested that the solid phase of bone can be represented using a Drucker-Prager type criterion (Carnelli et al., 2010;Tai et al., 2006). Uniaxial yield strains of 0.41% in tension and 0.83% in compression (Bayraktar and Keaveny, 2004) were converted to yield stresses by simply multiplying them by the Young's modulus of the solid phase (Schwiedrzik et al., 2015). Although there have been some experimental studies that have evaluated hardening of the extracellular matrix (Luczynski et al., 2015;Schwiedrzik et al., 2014), there is no agreement on the solid phase hardening behaviour of trabecular bone. Some studies have assumed linear hardening (Bayraktar and Keaveny, 2004;Bevill et al., 2006). In this study perfect plasticity was assumed (Carnelli et al., 2010), though small hardening (0.02% of the Young's modulus) was included to aid prevention of loss of ellipticity. In order to ensure global convergence of the Newton-Closest Point Projection Method (Newton-CPPM) scheme, a line search procedure was implemented as in the primal-CPPM algorithm proposed by Perez-Foguet and Armero (2002).
Each cubic specimen was subjected to 160 strain-controlled load cases as described in Table 1. The boundary conditions used to constrain the VE were kinematic uniform boundary conditions applied as described by Wang et al. (2009). It is recognised that these boundary conditions provide an upper bound for trabecular bone stiffness and also for yield (Panyasantisuk et al., 2015;Wang et al., 2009). The simulations were run on a Cray XC30 supercomputer hosted by ARCHER, the UK National Supercomputing Service. The analyses were carried out with an in-house parallel implicit finite strain solver, developed within the context of ParaFEM (Margetts, 2002;Smith et al., 2014), which uses an Updated Lagrangian formulation. This code uses MPI to perform the parallelisation Margetts, 2003, 2006). The high scalability of the code has already been demonstrated in previous work Margetts et al., 2015). Each of the 160 simulations per sample took approximately 12 min using 1920 cores; therefore the total number of core hours employed in this study is approximately 1.2 million.
A Newton-Raphson scheme was used as the solution tracking technique and a preconditioned conjugate gradient solver was used to solve the resulting linear algebraic systems. They are fast and if there are any convergence problems, they arise from the same origin (e.g. due to loss of positive definiteness of the stiffness matrix). However, convergence problems were only encountered in few of the porous samples (in 20 out of 3200 simulations) and can be related to a limit point or large-deformation related failure mechanisms (Bevill et al., 2006;de Souza Neto et al., 2008). Table 1 -Description of the load cases undertaken. Clockwise and counter-clockwise shear are differentiated by the sign of the off-diagonal terms of the homogeneous strain: clockwise corresponds to positive sign and counter-clockwise corresponds to negative sign. If T and C represent tension and compression respectively in normal strain space and clockwise and counter-clockwise shear in shear strain space, then the quadrants comprise of C-C, T-T, C-T and T-C; and octants comprise of C-C-C, C-C-T, C-T-C, C-T-T, T-C-C, T-C-T, T-T-C and T-T-T.  These points were marked differently in the figures and included in the fitting procedure as yield points. The initial load increment size corresponded to 0.1% macroscopic strain norm and could decrease to a minimum of 0.001% if global convergence was not achieved in larger load increments.

Theory and calculation
3.1.

Definition of the macroscopic yield points
The yield points were described in the plane where the abscissa is the Frobenius norm of the applied macroscopic strain, which corresponds to the Green-Lagrange strain, and the ordinate is the Frobenius norm of the homogenised Second Piola-Kirchhoff stress (Fig. 2), described by where there is no summation implied on repeated indices, V 0 is the initial volume of the VE, nel is the number of elements in the FE simulation, nip is the number of Gauss integration points in a trilinear hexahedra, J is the Jacobian, S is the Second-Piola Kirchhoff stress and w are the weights corresponding to the specific Gauss integration point. The 0.2% criterion was used to define the yield points (Wolfram et al., 2012), as shown in Fig. 2, and the elastic slope was obtained from the first two load increments, which were always fully elastic.

Formulation for macroscopic yield surface
Although a key aim of this study was to assess how different bone samples yield when subjected to the wide array of 160 load cases, we also examined the macroscopic yield surface fit using a quadric surface. The choice was based on its simplicity, because it has been previously related to the fabric tensor of bone (Cowin, 1986;Wolfram et al., 2012), and because it is a smooth surface. The quadric yield surface is described in strain space as where E is the Green-Lagrange strain tensor, F and F are used to define the shape, directionality and eccentricity of the yield surface. F is a fourth-order tensor, which has major and minor symmetries ( allowing it to be defined on a symmetric matrix space (Sym 6 ), by 21 coefficients (Mehrabadi and Cowin, 1990). As stated in Schwiedrzik et al. (2013), convexity of Eq. (6) is ensured with positive semi-definiteness of F. F is a symmetric second-order tensor, thus being described by 6 coefficients.

Different symmetries of the yield surface
Three different cases were investigated: isotropy, orthotropy and full anisotropy. Details about isotropic and orthotropic formulations can be found in (Schwiedrzik et al., 2013); only anisotropy is discussed in the following.
In the case of full anisotropy, or triclinic symmetry, the material can have different shear yield strains clockwise and counter-clockwise. For an anisotropic quadric, normal strains can interact with shear strains and shear strains can interact amongst themselves (Theocaris, 1992;Tsai and Wu, 1971). This means that in the triclinic case, F and F have 21 and 6 independent coefficients respectively.
By performing uniaxial strain load cases, several coefficients of F and all the coefficients of F can be determined. For F, the coefficients are In the case of F, the six diagonal coefficients of the projection of F onto Sym 6 are The 15 remaining parameters to be determined correspond to three normal strain interaction parameters, three shear strain interaction parameters and nine normal-shear strain interaction parameters. These parameters have expressions in the coefficients of F which are related to the previously stated diagonal coefficients, as shown in Table 2.
These, together with the six uniaxial normal strains and six uniaxial shear strains, add up to a total of 27 parameters. Calculating the determinant of 1 Â 1 and 2 Â 2 principal minors of the projection of F onto Sym 6 allows establishment Fig. 2 -Determination of the yield points by using the 0.2% criterion for the tensile and compressive uniaxial load cases of one sample. As it can be seen, the tensile and compressive uniaxial cases have the same elastic slope (dashed red line), as expected.
of basic restrictions on some of the coefficients to ensure that F is positive semi-definite, which are ε 7 ij Z 0;jζ kl j r1; i; j ¼ 1; 2; 3; k; l ¼ 1; 2; …; 6 The remaining restrictions on the coefficients are not expressed analytically but checked after the minimisation procedure to ensure positive semi-definiteness of F. For every symmetry case and for every sample, the eigenvalues of the projection of F onto Sym 6 were checked to ensure they were non-negative.

Evaluating goodness of fit
The macroscopic yield envelope was fitted by using a minimisation procedure in MATLAB (Mathworks, Natick MA, USA).
To evaluate the goodness of fit, the error was evaluated as where ‖Á‖ is the Frobenius norm of the corresponding macroscopic strain and N is the cardinality of a specific set of load cases. This error was evaluated for four different sets: for all the load cases; for the load cases which entirely lie on the normal strain space; for the load cases which entirely lie on the shear strain space; and for the strain cases that have one component in the normal strain space and one component in the shear strain space (which from now onwards will be referred to as combined normal and shear strain space).

Macroscopic yield strains
Macroscopic yield points in strain space for all twenty considered samples in the normal-normal and shear-shear planes are shown in Fig. 3. These represent 36 of the 160 load cases analysed for each sample, and no projections have been made, i.e. the yield strains only contain out-of-plane components equal to zero. Results show that the macroscopic yield surface of bone has a higher yield strain in compression than in tension in the normal strain space. This is expected due to the characteristics of its solid phase. It can be seen that the tensile quadrant displays quasiuniform macroscopic yield strains across samples (upper right quadrant of Fig. 3a, b, c). The compressive yield strains have some variability across samples, as can be seen from the spread of yield points in the lower left quadrant of Fig. 3a, b, c. The largest variation in the normal-normal planes is in the tensile-compressive quadrants (upper left and lower right quadrants, Fig. 3a, b, c).
Macroscopic yield strains in the shear-shear planes show a large variation of yield strains for different specimens (Fig. 3d, e, f). It can also be observed that the shear yield strains of bone are different in clockwise and counterclockwise directions, with these absolute differences ranging from 0.0034% to 0.4463%. A statistical comparison between these yield strains was performed for all pure shear cases with a paired t-test. This test suggests that paired clockwise and counter-clockwise shear yield strains are statistically different (po0.01).
Macroscopic uniaxial (tensile, compressive and shear) yield strains were related to BV/TV and fabric through multilinear regressions performed in log space. No relationship between yield strains and BV/TV and fabric was found. Only compressive uniaxial yield strains were mildly related to BV/ TV and fabric (R 2 ¼ 0.44, p-0).
In order to examine macroscopic yield strains in tensiontension, compression-compression and tension-compression regimes, we evaluated the mean of the macroscopic Green Lagrange strain norm for each of the above three regimes, as shown in Fig. 4. As expected, the mean of the norms is the lowest for tension-tension, highest for compressioncompression and in between for tension-compression. Fig. 4 also shows the standard deviation in the evaluated norms. It can be observed that the deviation is relatively small for the tension-tension regime, higher for compression-compression regime and the highest for tension-compression regime.

Solid phase strains
We examined strains at the microscale (solid phase strains). Under uniaxial macroscopic tension, more localised strains were found to occur at the solid phase level, and there were mostly no compressive solid phase strains anywhere in the specimens. However, under uniaxial macroscopic compression, the compressive solid phase strains were more diffused and found to occur throughout the geometry. Further, under macroscopic compression, large tensile solid phase strains were found to arise, due to bending and buckling of trabeculae. Fig. 5 and Fig. 6 show this for a slice from one typical porous sample.    Table 3. Plots for the two samples with the highest and lowest densities are shown in Fig. 7. The lower density sample shows a higher level of anisotropy in comparison to the higher density sample. Further, if consecutive macroscopic yield points of the porous sample are joined up, then the homogenised yield envelope does not always remain entirely convex in some of the planes (e.g. In general, the fitting errors were not found to correlate with bone density (Fig. 9). In normal strain space (Fig. 9a), the assumption of an isotropic quadric led to consistently higher errors, while the orthotropic and anisotropic assumptions resulted in smaller fitting errors. In the shear strain space (Fig. 9b), in the combined normal and shear strain space (Fig. 9c), and in the general strain space (Fig. 9d), the assumption of anisotropic quadric had the smallest errors. A mild trend of errors decreasing with increasing density was observed for the isotropic and orthotropic assumptions in shear strain space. In the general strain space, the errors and the error differences between assumptions tend to reduce with increasing density.

Discussion
Our study shows that the macroscopic yield surface of bone in normal strain space is fairly uniform across a wide range of samples; this confirms findings of previous research (Bayraktar et al., 2004;Lambers et al., 2014;Pankaj and Donaldson, 2013).
Our results also demonstrate that the full threedimensional macroscopic yield behaviour of trabecular bone can be reasonably well described using the isotropic quadric yield surface, though orthotropic and anisotropic surfaces lead to smaller errors. This is in agreement with previous studies in which the strain space yield surface was reported to be isotropic (Bayraktar et al., 2004) and more recently transversely isotropic (Sanyal et al., 2015) and orthotropic (Wolfram et al., 2012). However, unlike the above studies, our   treatment of clockwise and counter-clockwise shears as separate loading cases showed that the differences in the two directions are significant. This is probably because trabeculae are not symmetrically aligned with respect to the axes of the material, which in this case were assessed through the eigenvectors of the MIL fabric tensor. It is important to note that the assumption of identical macroscopic yield points in clockwise and counterclockwise directions restricts the system to orthotropy at best (Theocaris, 1992;Tsai and Wu, 1971). We also observed predominance of tensile solid phase strains in pure macroscopic shear, which is consistent with Sanyal et al. (2012). Multilinear regressions suggest that uniaxial yield strains are not correlated with BV/TV and fabric (Matsuura et al., 2008; Morgan and Keaveny, 2001;Panyasantisuk et al., 2015). Only a mild dependence was found for the uniaxial compressive yield strains (R 2 ¼0.44, p-0), with a positive slope for density and a negative slope for fabric, which suggests that long trabeculae, i.e. associated with a high fabric eigenvalue, have lower macroscopic yield strain, as suggested by Matsuura et al. (2008). Since no clear relationship between a fabric tensor and the yield strains was found, use of an isotropic yield surface formulation in strain space is more practical for real applications.
In uniaxial macroscopic tension, solid phase strains are almost exclusively tensile and independent of density (Bayraktar and Keaveny, 2004;Lambers et al., 2014). The highly oriented structure of trabecular bone results in the yield strains at the solid phase and at the macroscale being very similar in tension. When trabecular bone is loaded in macroscopic compression, yield mechanisms are different: in this case, yielding at the solid phase was found to occur both due to tension (arising from bending and buckling of trabeculae) and compression. As expected, we found considerable tensile strains in trabeculae for low density samples as has been previously reported (Bevill et al., 2006;Morgan et al., 2004;Stolken and Kinney, 2003). This density dependence results in macroscopic yield variation being displayed via a small spread of yield points in the compression-compression quadrants (lower left quadrants of Fig. 3a-c), as shown by the mild relationship between compressive uniaxial yield strains and density and fabric. This also implies that solid phase uniaxial yield strain asymmetry is not fully maintained at the macroscale and generally reduces with increased porosity and increased fabric eigenvalues. Our results are consistent with the experimental results of Lambers et al. (2014) in the sense that the number of microscopic yielded sites in macroscopic compression and in macroscopic tension are similar in number, but in macroscopic tension, the microscopic yield zones have more localised strains, which could be related to microcrack propagation.  The few previous studies that have evaluated the macroscale yield surface of bone from its microstructure have all used a bilinear criterion, with yield strain asymmetry, to define the solid phase of bone with a reduced stiffness beyond defined tissue yield values (Bayraktar et al., 2004;Sanyal et al., 2015;Wolfram et al., 2012). In our study, the solid phase of bone was modelled using a Drucker-Prager type criterion which has been validated via experimental studies (Carnelli et al., 2010;Tai et al., 2006). Our study considered yield points arising from 160 different load cases. Some recent studies have been limited to 17 load cases (Panyasantisuk et al., 2015;Wolfram et al., 2012). In order to compare our results we considered 17 strain cases similar to those in the cited studies. We found the errors to be 11.4% for the isotropic case and 10.3% for the orthotropic case for the 17 cases. The errors, taking into account all 160 cases, are 11.2% for the isotropic case, 10.4% for the orthotropic case and 7.8% for the anisotropic case. In other words, the 17 load cases lead to errors of a similar magnitude to those obtained using all 160 load cases. To further examine the effect of the considered load cases on the fitting error, we considered a single sample with different load cases. We evaluated the fitting errors considering all normal load cases proposed by Wolfram et al. (2012) (14 cases) and all our shear cases (26 cases). The errors were 18.0%, 16.3% and 7.6% for isotropic, orthotropic and anisotropic assumptions respectively. With the 17 load cases mentioned previously, the errors reduce to 9.8% and 6.5% for isotropic and orthotropic assumptions respectively for this sample. With all the 160 load cases, the errors are 13.3%, 12.9% and 9.2% for isotropic, orthotropic and anisotropic assumptions respectively. This illustrates that shear cases contribute importantly to anisotropy, and that the fitting errors clearly depend on the considered load cases, which illustrates the importance of examining a large range of complex load cases. With respect to the combined normal and shear strain spaces, the shear component of the macroscopic yield strain is often found to increase when there is a compressive normal component, indicating typical cellular solid behaviour (Fenech and Keaveny, 1999;Gibson and Ashby, 1997).
Our study has a number of limitations. We chose to use a quadric due to its simplicity, because it has been used in previous studies, because it requires fewer parameters than higher order criteria, and because it is a smooth surface, not requiring multiple plastic multipliers. Although our primary aim was to examine the effect of material symmetry assumptions on the macroscopic yield surface, the fitting errors clearly depend on the shape of the chosen surface and on the considered load cases. While we examined full anisotropy with a quadric surface, a previous study employed a higher order polynomial surface, a quartic, but restricted it to transverse isotropy (Sanyal et al., 2015); a restriction that is not fully supported by our results.
Although the solid phase constitutive law has been validated (Carnelli et al., 2010;Tai et al., 2006), there is no experimental validation for the macroscopic yield surfaces. In fact, such a validation is impossible as samples tested once cannot be retested and it is not possible to obtain numerous (or even two) identical samples. There have, however, been some attempts to on trabecular bone wherein the loading cases have been limited to triaxial compression Rincon-Kohli and Zysset, 2009). Thus, the range of complex strain cases we have tested can only be performed numerically.
We considered homogeneous tissue elastic properties while some previous studies have assessed the effect of heterogeneous mineral density on the macroscopic stiffness of bone (Blanchard et al., 2013;Renders et al., 2008). Renders et al. (2008) found a decrease of 21% in apparent stiffness when considering mineral heterogeneity. However the effect of heterogeneities at the solid phase on finite element models with geometrical nonlinearities is still unclear. Since we wanted to be able to compare our macroscopic yield strains with previously published results, we kept our solid phase elastic properties as homogeneous.
The study assumes that the solid phase of bone can be modelled as a plastic material, which is not entirely true as high localised strains can cause microcracks and eventual fracture; effects that full plasticity based models may not be able to capture. Furthermore, we did not consider hardening (Carnelli et al., 2010) because there is no agreement on the hardening law at the scale of the solid phase we are considering. However, previous experimental and theoretical studies such as Schwiedrzik et al. (2014), Luczynski et al. (2015) and Fritsch et al. (2009) showed that the extracellular matrix of bone has a hardening behaviour after yield.

Conclusions
Trabecular bone has fairly uniform macroscopic yield behaviour across samples in normal strain space. Thus, modelling it by using strain-based plasticity makes sense. For tensiontension and compression-compression quadrants in normal strain space, the strain norm at yield shows little variation, indicating an isotropic behaviour in these regimes.
In the tension-compression quadrants, pure shear and combined normal-shear planes, the macroscopic strain norms at yield have a relatively large variation, indicating a possible absence of isotropy. Further, differences in yield strain values in clockwise and counter-clockwise shear may indicate a possible anisotropy for the macroscopic behaviour of trabecular bone. However, due to the difficulties of formulating a non-isotropic closed-form yield surface in strain space due to the weak relationships between fabric and yield strains, and due to the small difference in fitting errors between isotropic and orthotropic or anisotropic considerations, an isotropic criterion presents itself as the most suitable approximation. However, for some load cases, considerable differences between the closed-form yield criterion and the actual yield strain may arise.
With respect to the yield surface, an eccentric-ellipsoid may adequately represent the macroscopic yield surface of bone as the fitting errors for all the considered symmetries are reasonably small. However, it is important to be mindful of the asymmetry in shear yield strains and that in the normal-shear load cases, the quadric may not be able to represent the macroscopic yield behaviour of trabecular bone.
This work provides a comprehensive assessment of material symmetries of trabecular bone at the macroscale and describes in detail its macroscopic yield and its underlying microscopic mechanics.