## Abstract

We model trabecular bone as a nanocomposite material with hierarchical structure and predict its elastic properties at different structural scales. The analysis involves a bottom-up multi-scale approach, starting with nanoscale (mineralized collagen fibril) and moving up the scales to sub-microscale (single lamella), microscale (single trabecula) and mesoscale (trabecular bone) levels. Continuum micromechanics methods, composite materials laminate theory and finite-element methods are used in the analysis. Good agreement is found between theoretical and experimental results.

## 1. Introduction

Bone is a connective tissue made of a cortical (compact) bone, forming a hard outer layer, and a trabecular (spongy) bone, filling the interior spaces and ends of long bones. Bone is a multi-phase composite material consisting of organic phase (32–44% bone volume (BV)), inorganic phase (33–43% BV) and water (15–25% BV) [1]. The organic phase is composed of collagen type I (approx. 90%) and non-collagenous proteins (NCPs) (approx. 10%) [1]. The inorganic (mineral) phase is made of calcium phosphate, which is similar to hydroxyapatite (HA), Ca_{10}(PO_{4})_{6}(OH)_{2}. The mineral phase is stiff and strong but brittle, whereas the collagen phase is soft and highly deformable [2]. Water plays an important role in the bio-mineralization process and serves as a plasticizer, enhancing the toughness of bone. These components are arranged into a complex hierarchical structure, which makes bone stiff, strong, tough and yet lightweight.

We define five levels of hierarchical organization in bone, which are described below with a focus on trabecular bone (figure 1).

*Nanoscale* (mineralized collagen fibril), ranging from few to several hundred nanometres, includes collagen molecules and HA crystals. Cross-linked collagen molecules, which are triple helical protein chains about 1.5 nm in diameter and 300 nm in length with 40 nm gaps between ends, are staggered in parallel to form fibrils of about 50–100 nm in diameter. Neighbouring collagen molecules are shifted by 67 nm resulting in a banded structure with a 67 nm periodicity, as illustrated in the electron microscopy images (figure 2*a*). Collagen is formed first, followed by mineralization, which involves the filling up of gaps and spaces in-between collagen molecules as well as those outside with HA nanocrystals. HA minerals have a plate-like shape with an average size of 50 × 25 × 3 nm^{3} [3,4]. Such mineralized collagen fibrils are the building blocks of both cortical and trabecular bone types.

*Sub-microscale* (single lamella), spanning from one to few micrometres, contains mineralized collagen fibrils which are aligned preferentially to form a single lamella of thickness approximately 3–7 μm (figure 2*b*). Spaces between the fibrils are filled with randomly arranged minerals forming a porous foam-like structure. Each lamella contains ellipsoidal cavities, typically 5–15 μm in cross section and 25 μm in length, called lacunae, which house osteocytes. Again, this scale has similar features in cortical and trabecular bones.

*Microscale* (single trabecula), ranging from tens to hundreds micrometres, represents a trabecular bone tissue. It is made of trabecular bone packets [5] (consisting of layers of lamellae oriented in different directions) and cement lines, which form trabeculae, and an interstitial bone at interconnects of trabeculae. A typical trabecular packet, which is formed during bone remodelling, has a crescent shape and is about 50 μm thick and 1 mm long (figure 2*c*). Trabeculae can be in the form of rods or plates. In cortical bone, the microscale consists of osteonal, interstitial and circumferential lamellae.

*Mesoscale* (trabecular bone), ranging from hundred micrometres to several millimetres, or larger, depending on the bone size, consists of a porous network of trabeculae (figure 2*d*). The pores, typically in the order of 1 mm, are filled with bone marrow, fat and bone cells. In cortical bone, this level represents randomly arranged osteons embedded in an interstitial lamella, with some resorption cavities, all surrounded by a circumferential bone.

*Macroscale* (whole bone level), spanning from several millimetres to several centimetres or more, depending on species, consists of both cortical and trabecular bones.

Various analytical and computational models were proposed to predict elastic properties of trabecular bone at these different structural scales. At nanoscale, bone was mainly modelled by using continuum mechanics approaches and was represented as a composite material with collagen matrix and reinforcing HA inclusions [6–8]. More recent studies incorporated the effect of water and NCPs [9–12]. Computational models, using a finite-element method (FEM), were used in references Ji & Gao [13], Siegmund *et al*. [14] and Yuan *et al*. [15]. Molecular dynamics (MD) simulations were also used to study fibrillar collagen [2,16] and collagen–crystal interactions [17–20]. At sub-microscale, a single lamella was modelled computationally as a random network of preferentially oriented mineralized collagen fibrils [21] and analytically as a matrix-inclusion composite [9,10,12]. At microscale, properties of trabecular bone tissue were mainly studied experimentally using microtensile test [22–24], bending test [25,26], ultrasound [23,27,28] and nanoindentation [27,29–33]. Rice *et al*. [34] experimentally obtained Young's modulus of trabecular bone and used the data together with Christensen's model [35] for low-density materials (LDM) to back-calculate tissue properties of bone. Similarly, van Rietbergen *et al.* [36] used a three-dimensional FEM model along with experimental data for apparent modulus to back-calculate lower and upper bounds for trabecular bone tissue modulus. At mesoscale, elastic behaviour of trabecular bone was studied using several different approaches, involving analytical and computational techniques. Analytical studies represented trabecular bone as a cellular solid and expressed its Young's modulus by power law relations in terms of density [37–42]. Although density is a key parameter in determining properties of trabecular bone, it alone cannot fully capture the mechanical behaviour of bone. Other researchers defined a fabric tensor, which characterizes the textural or structural anisotropy of trabecular bone, and found the relationships between the elastic constants of trabecular bone and its fabric tensor and density [43–45]. Trabecular bone's architecture, characterized by thickness, number and separation distance of individual trabecula as well as their three-dimensional connectivity, plays an important role in its response. Thus, high-resolution imaging techniques, such as micro-computed tomography (μCT), accounting for actual trabecular bone architecture, were used in combination with the FEM to predict elastic moduli of trabecular bone [31,46–51]. Computational studies involving idealized periodic geometry included [52–57].

In this paper, we present a multi-scale analytical and computational approach to predict the elastic properties of trabecular bone. At the mesoscale, we account for actual bone geometry. Results are compared with our own experimental data at mesoscale and experimental results reported in literature.

## 2. Methods

### 2.1. Experiments

#### 2.1.1. Specimen preparation

Bone tissue was obtained from the Emory Body Donor Program (Emory University, Atlanta, GA, USA). Samples were extracted from proximal tibia (near knee joint) of an 88-year-old male. The donor expired due to a myocardial infarction and was otherwise a relatively healthy individual with no discernable signs of bone disease. The top portion, approximately 10 cm, of tibia was obtained. The head of tibia was cut off, about 3 cm from the top, exposing its greatest cross section. Four samples were then cored out using a trephine 8 mm in diameter. The cylindrical samples were taken such that their axis of symmetry was along tibia's long axis. These samples were next cut to a 1 : 1 aspect ratio using a metal mould and an irrigated saw. Afterwards, the samples were kept in saline moistened gauze wraps and stored at −20°C. Samples were thawed in saline for 3 h before imaging and testing. One bone sample was used for μCT imaging, while all four samples were used for mechanical testing.

#### 2.1.2. Micro-computed tomography imaging

A trabecular bone sample was imaged at the Georgia Institute of Technology Orthopaedic Bioengineering Laboratory (Atlanta, GA, USA) using a SCANCO Medical Products (Basserdorf, Switzerland), MicroCT 40 unit. The cylindrical sample (8 mm in length and diameter) was scanned with a 37 μm resolution and a 20 μm voxel size. A total of 173 two-dimensional slices were taken of the sample. Contours were drawn around the filled voxels to determine an outer boundary of the sample and reconstructed into a three-dimensional image. A μCT image captures the actual trabecular bone architecture from which three-dimensional connectivity, trabecular thickness, trabecular number and trabecular spacing can be obtained. The μCT imaging process discretizes a sample into cubic elements called voxels. Thresholding distinguishes bone voxels from pore voxels. The threshold value was selected by choosing a middle point between the peak corresponding to trabecular bone and the peak corresponding to pores in the image spectrum, using the procedure described in Basillais *et al*. [58] and Stock [59]. The BV and total volume (TV) of the sample are calculated, respectively, from the number of voxels assigned to bone and the total number of voxels and, then, the volume fraction (VF) of the sample can be determined as BV/TV.

#### 2.1.3. Compression test

Uniaxial compression test was used to obtain apparent Young's moduli of trabecular bone samples. The term ‘apparent’ is used because the sample size is smaller than the representative volume element (RVE) [60]. Four cylindrical trabecular bone samples were placed between polished steel plates at room temperature and loaded in the direction of their axis of symmetry. Specimens were not constrained at the platen from an in-plane motion. An Instron Mini-Bionix testing machine with a 1000 N load cell was used to apply a compressive strain rate of 0.01 s^{−1}. The applied strain rate was within the range of strain rates that occur *in vivo*, which are between 0.01 and 0.08 s^{−1} [61]. After obtaining the stress–strain curves of trabecular bone samples, their Young's moduli were calculated by fitting a linear regression through the initial linear portions of the curves.

### 2.2. Multi-scale modelling of trabecular bone

The multi-scale modelling approach consists of successive homogenization steps. Elastic properties of trabecular bone are calculated at each structural level, from nanoscale to mesoscale (figure 1). In the analysis, results from a lower level are used as inputs for a higher level. Continuum micromechanics approaches, laminate composite materials theory (LCMT) and FEMs are used to account for structures of bone at these different scales. Figure 3 summarizes all the modelling steps and shows how they are linked together at different scales. In the following sections *E*, *ν*, **C** and denote, respectively, Young's modulus, Poisson's ratio, elastic stiffness tensor and VF of a pertinent phase. We use the term ‘Young's modulus’ when a material is isotropic and, more generally, the term ‘elastic stiffness tensor’ when a material may be anisotropic. is the fourth-order Eshelby tensor accounting for the shape of phase *r* in a matrix with stiffness tensor **C**_{0}, where 0 is a generic subscript.

Mechanical properties and VF of bone's constituents are the key parameters for determining its overall behaviour. Our choices of properties and VFs of bone's constituents, collagen, HA, water and NCPs, are listed in table 1. A wide range of values for Young's moduli of collagen and HA are reported in literature (see table 1 in Hamed *et al*. [10]). The reported values for Young's modulus of collagen range from less than 1 GPa [62] up to 12 GPa [63], depending on the state of hydration, species and experimental technique used to measure them. In order to assess the effect of collagen properties on the overall properties of trabecular bone, we select two values of collagen Young's modulus: 1.5 [6,64–66] and 5.4 GPa [67], which are, respectively, on lower and higher sides of the range. The Poisson ratio of collagen is assumed to be 0.28, so that the overall Poisson's ratio of wet collagen composite becomes about 0.35, as estimated by Katz [7] and used by Nikolov & Raabe [11]. For HA crystals, Young's modulus is chosen to be 114 GPa [68,69], while Poisson's ratio is set to be 0.23 based on the results of *ab initio* calculations by Snyders *et al*. [70]. To our knowledge, no data are available in the literature for the mechanical properties of NCPs. Since they consist of flexible coiling macromolecules, their Young's modulus must be lower compared with modulus of collagen with its relatively stiff triple-helical structure. Here, we assume that the NCPs have isotropic properties with Young's modulus equal to 1 GPa [11] and Poisson's ratio of 0.45 which is a typical value for soft polymers with flexible molecules [11]. For water, the bulk modulus is chosen to be 2.3 GPa with a Poisson's ratio of 0.49 corresponding to a nearly incompressible material. For simplicity, all components are assumed to have linear elastic and isotropic behaviour. This is an idealization. Many researchers [68–70] reported the properties of a single HA crystal to be anisotropic. Similarly, for collagen, Cusack & Miller [63] reported a transversely isotropic behaviour. However, while modelling bone as a collagen–HA composite, almost all researchers used isotropic behaviour for both collagen and HA crystals. This was due to simplicity and, more importantly, owing to lack of enough experimental data for anisotropic properties of collagen and HA in bone. Most references report isotropic properties of HA and collagen (giving a wide range of values), while the data on anisotropic constants are very limited (also showing a wide range of values). Our model is set up in such a way that it could easily handle anisotropic properties of collagen and HA. However, owing to the uncertainty of what the anisotropic constants are, we choose isotropic properties for constituents of bone.

Furthermore, researchers used VFs for HA minerals ranging from 32 to 52 per cent [10,11]. Here, we selected an intermediate value between the lower and upper bounds, 42 per cent which was also used by Jager & Fratzl [71] for a fully mineralized bone. The collagen VF is assumed to be 41 per cent [72]. Given that collagen comprises approximately 90 per cent of the organic matrix [1], the VF of NCPs is about 5 per cent. The remaining phase is water with a VF of 12 per cent. Moreover, table 2 presents the local VFs of phases used as modelling inputs at each hierarchical level.

#### 2.2.1. Nanoscale

Collagen molecules are assembled in a staggered array with gap and overlap zones to form a collagen fibril [3]. Mineral crystals are located in gaps and spaces between collagen molecules, as well as outside the collagen fibril, while water and NCPs fill remaining open spaces [12]. This structure is called a mineralized collagen fibril. We model collagen in its wet state as a collagen–water composite and represent a mineral, which contains intercrystalline pores within, filled with water and NCPs, as a HA–water mixture. We consider two different cases, using two assumptions about collagen–HA interactions: interpenetrating phases versus matrix-inclusion phases. The first one is motivated by the recent experimental observations of Chen *et al*. [73] which showed that HA crystals form a continuous phase in bone, while the second one is included for comparison since most previous studies used this model.

*Assumption 2.1.* Both phases, namely collagen–water and interfibrillar HA–water composites, interpenetrate each other. In this case, we use a self-consistent (SC) method [74,75], with two types of inclusions and no matrix, to model bone at nanoscale. Cross-linked collagen molecules are assumed to be cylindrical with an aspect ratio of 1000 : 1 : 1 following the approximately 100 μm length [76] and approximately 100 nm diameter of collagen fibrils [1,77], while HA crystals are represented as ellipsoids with an aspect ratio of 50 : 25 : 3 following [4]. Given the stiffness tensors of collagen–water composite, **C**_{colw}, and interfibrillar HA–water mixture, **C**_{HAw}, the effective elastic stiffness tensor of a mineralized collagen fibril, **C**_{fib}, is obtained as
2.1

Subscripts ‘colw’, ‘HAw’ and ‘fib’ denote collagen–water mixture, interfibrillar HA–water composite and mineralized collagen fibril, respectively. Superscripts ‘cyl’ and ‘ellips’ denote, respectively, cylindrical and ellipsoidal shapes of two phases, namely collagen and HA minerals. Since effective properties of a mineralized collagen fibril, **C**_{fib}, are not isotropic, the components of Eshelby tensor are evaluated numerically for an ellipsoidal inclusion embedded in a general anisotropic matrix using a Fortran code developed by Gavazzi & Lagoudas [78]. Equation (2.1) is solved iteratively to obtain the implicit unknown **C**_{fib}. Note that the Eshelby tensors and , which are dependent on **C**_{fib}, must be updated in each iteration.

*Assumption 2.2.* Cross-linked wet collagen molecules play the role of a continuous matrix, while interfibrillar HA–water minerals act as reinforcing inclusions. Here, we use a Mori–Tanaka (MT) scheme [79,80] to obtain **C**_{fib} as follows:
2.2

where subscripts and superscripts are as defined above.

Effective stiffness tensor of collagen–water composite, **C**_{colw}, is obtained by using the MT method, with the cross-linked collagen molecules being a matrix and the voids filled with water and NCPs being inclusions, as
2.3

Moreover, a mineral interacting with water can be represented as a porous solid with some intercrystalline voids within (internal defects), filled with water and NCPs. We use the MT method to predict the elastic stiffness tensor of HA–water mixture, **C**_{HAw}, as follows:
2.4

In equations (2.3) and (2.4), ‘col’, ‘w’ and ‘HA’ denote, respectively, dry collagen, water and NCPs and HA crystals. Also, superscript ‘sph’ denotes the spherical shape of pores. We assume that water VFs in collagen–water and HA–water composites are equal.

Note that HA minerals in bone are distributed non-uniformly owing to stochastic nature of mineralization process and owing to remodelling [81], leading to spatial heterogeneity at different structural scales. We use effective medium theories which allow us to specify the alignment and shapes of phases but not their specific arrangement. Our models assume that different phases are distributed randomly but in a statistically uniform way which enables us to obtain homogenized properties of a mineralized collagen fibril. Alternatively, one could assume a non-uniform distribution of minerals. However, since the exact distribution of minerals in collagen fibrils is still not well understood (experimental techniques are not yet available to make such measurements), we choose the above described approach for simplicity. Secondly, local mineral contents vary spatially owing to remodelling. In this paper, we assume a uniform spatial mineral distribution, again for simplicity.

#### 2.2.2. Sub-microscale

At sub-microscale, we use two different modelling steps [10]: (i) combining mineralized collagen fibrils with an extrafibrillar HA foam, and (ii) combining the composite obtained in step (i) with lacunar cavities to form a single lamella.

##### 2.2.2.1. Mineralized collagen fibrils combined with an extrafibrillar hydroxyapatite foam (coated fibrils)

X-ray diffraction, atomic force microscopy and transmission electron microscopy (TEM) confirmed the existence of randomly dispersed extrafibrillar minerals outside the fibrils [82–85]. Such HA crystals with intercrystalline pores in-between, filled with water and NCPs, can be thought of as an extrafibrillar HA foam [86,87]. This structure motivates the use of SC scheme with two interpenetrating phases, HA–water minerals and pores filled with water and NCPs, to determine the overall elastic stiffness of extrafibrillar foam, **C**_{Efoam}, as
2.5where subscript ‘Efoam’ refers to extrafibrillar HA foam. Random orientation of extrafibrillar HA crystals leads to isotropy of the homogenized material. Therefore, for simplicity, both phases are assumed to be spherical in shape, following [88].

We propose the following two methods to model the interaction between mineralized collagen fibrils and extrafibrillar HA foam:

*Assumption 2.3*. Mineralized collagen fibrils, having elastic properties obtained in §2.1, and extrafibrillar HA foam, with effective properties obtained from equation (2.5), interpenetrate each other. Fibrils are assumed to be cylindrical in shape and unidirectionally aligned, while extrafibrillar HA foam is assumed to be spherical. Subscript ‘cfib*’* denotes mineralized collagen fibrils combined with extrafibrillar HA (coated fibrils). The SC scheme is applied to obtain the effective elastic tensor of coated fibrils, **C**_{cfib}, as
2.6

In equation (2.6), superscripts ‘cyl’ and ‘sph’ represent cylindrical and spherical shapes of mineralized fibrils and extrafibrillar HA foam, respectively, while and denote Eshelby tensors for cylindrical (fibrils) and spherical (HA foam) inclusions.

*Assumption 2.4*. Extrafibrillar HA foam and mineralized collagen fibrils are represented, respectively, as a matrix and inclusions. This assumption motivates us to use the MT method to predict **C**_{cfib} as
2.7where subscripts ‘cfib’, ‘fib’ and ‘Efoam’ are as defined above, and denotes the Eshelby tensor corresponding to cylindrical inclusions (fibrils) embedded in the extrafibrillar HA foam.

##### 2.2.2.2. Single lamella

The matrix obtained in equations (2.6) or (2.7) is perforated by ellipsoidal cavities (lacunae) to form a single lamella. Effective stiffness tensor of a single lamella, **C**_{lamella}, is obtained by using the MT scheme as
2.8

Subscripts ‘lamella’ and ‘lac’ denote, respectively, a single lamella and lacunae. Lacunae are represented by ellipsoids with a 5 : 2 : 1 aspect ratio, following dimensions 25 × 10 × 5 µm^{3} [12,89], occupying 3 per cent VF.

#### 2.2.3. Microscale

Elastic properties of a single trabecula are obtained by following the homogenization scheme of Sun & Li [90] developed using LCMT. In their model, a thick laminate was composed of several laminas with different fibre orientations. In our problem, a single trabecula plays the role of a thick laminate and it consists of *k* lamellae, each having a preferential collagen fibril orientation. In trabecular bone, the lamellae are arranged into orthogonal, rotated and twisted motifs [91,92], similarly as in cortical bone [93]. Since we do not have enough information on the actual arrangement of lamellae in each trabecula, we assume for simplicity that the lamellae are oriented randomly, spanning whole set of directions, which gives rise to an isotropic response. In reality, trabeculae may be anisotropic and their properties may change from one trabecula to another.

Three vectors *V*_{1}, *V*_{2} and *V*_{3} are used to define a preferential orientation of the *k*th lamella. The vector *V*_{1} is obtained through the following equations:
2.9where *R* is a random number, with values 0 ≤ *R* ≤ 1, generating a randomly oriented pattern of lamellae. Since vectors *V*_{2} and *V*_{3} are contingent on vector *V*_{1}, they can be obtained as follows:
2.10

A transformation matrix, *T _{ij}*, is used to account for different fibril orientation in each lamella [94]

2.11

where *m _{i}, n_{i}* and

*p*are direction cosines of axis

_{i}*i*(

*i*= 1

*,*2

*,*3), that is, 2.12

After transformation, the stiffness tensor of the *k*th lamella, **C**^{(k)}, is obtained as
2.13where **C** is the stiffness tensor of a single lamella as obtained in equation (2.8). These *k* lamellae are then stacked together, according to Sun & Li's formulation [90], to build a single strut. A nearly isotropic response is obtained using a large number of randomly oriented lamellae.

#### 2.2.4. Mesoscale

We use the μCT-based FEM to study the elastic behaviour of trabecular bone. To create the FEM model from the digitized geometry, first the surfaces are exported from the μCT data into a file (.stl format) consisting of nodes and two-dimensional tetrahedral surfaces. These two-dimensional surface elements are then meshed into three-dimensional tetrahedral elements using HyperMesh v. 6.0 (Altair Engineering, Inc.). The μCT data are originally in the form of voxels giving rise to pixellated edges in the structure. However, after meshing with tetrahedral elements, the edges and boundaries become smooth. HyperMesh uses two types of smoothing algorithms (which are specially designed to smooth jagged edges of structures created from topology optimization): a modified Laplacian over-relaxation algorithm for size-correcting smoothing and a modified isoparametric-centroidal over-relaxation algorithm for shape-correcting smoothing. The FEM mesh, illustrated in figure 4*a*, is made of 6 063 971 elements and 1 558 400 nodes. After meshing, the quality of elements is checked, and for the elements that fail the quality criteria, node locations or element edges are swapped in order to improve quality. Detailed section of the tetrahedral mesh is shown in figure 4*b*.

After the mesh is constructed, linear elastic FEM analysis is conducted using OptiStruct v.11.0 (Altair Engineering, Inc.) to compute the apparent Young's modulus of trabecular bone in the direction of loading. The bone tissue is assumed to be linear elastic and isotropic with properties obtained in §2.2.3. Boundary conditions for the FEM model idealize those of a uniaxial compression test: a uniaxial displacement (uniform strain) is applied to the top surface of the cylindrical bone sample, the bottom surface is kept fixed, while the sides are taken to be traction-free. Rough surfaces at the top and bottom of the cylindrical sample are trimmed to obtain straight surfaces to facilitate the application of displacement boundary conditions. The FEM run for a single case takes about 20 min on a computer with CPU AMD Athlon 64 X2 5200 at 2.6 MHz, 8 GB RAM, Linux Fedora 14. The run requires 3.5 GB RAM and 20 GB disk space.

An energy approach is used to calculate the apparent Young's modulus of the bone sample, *E*_{bone}. Elastic strain energy density is defined as
2.14where *W* is the total elastic strain energy of the system and *V* its volume. and are, respectively, average stress and strain in the direction of loading, which are related to each other through equation . According to the average strain theorem, the volume average of strain is equal to the applied strain, , i.e. . Thus, the bone apparent Young's modulus is
2.15

Alternatively, one could use a direct method, instead of the energy approach, to obtain the apparent modulus as .

For comparison, we use the analytical model proposed by Gibson [40] to predict Young's modulus of trabecular bone, *E*_{bone}, as a function of its relative density:
2.16where *E*_{strut} is Young's modulus of a single trabecular strut as obtained in the previous level, *ρ*_{bone} and *ρ*_{strut} are, respectively, the density of trabecular bone and the density of solid bone (struts), and *A* is a constant of proportionality. Gibson [40] developed two types of structures for trabecular bone: a network of rod-like elements at low relative densities (less than 0.2) and a network of plate-like elements at higher relative densities (greater than 0.2). The first structure forms an open cell, while the later one forms a closed cell. Here, *n* was determined to be equal to 2 for an open cell and 3 for a closed cell [40]. The relative density, *ρ*_{bone}/*ρ*_{strut}, is equal to bone volume fraction determined by μCT.

We also evaluate effective Young's modulus of trabecular bone, *E*_{bone}, by using Christensen's result [35] for isotropic low-density materials (LDMs), which is given as
2.17where *ν*_{strut} is the Poisson ratio of a single trabecular strut, obtained in §2.2.3, and *Φ*_{bone} is the VF of bone material, which is the same as VF parameter obtained by μCT.

Finally, as an alternative approach, we use the MT method with trabecular struts being a matrix and trabecular pores being inclusions. The effective stiffness tensor of trabecular bone, **C**_{bone}, is then found as
2.18where subscripts ‘strut’ and ‘pore’ denote, respectively, trabecular struts and pores. Properties of a single trabecula, **C**_{strut}, are taken from §2.2.3 and the pores are assumed to be spherical with VF as obtained from μCT imaging.

#### 2.2.5. Summary of modelling approaches

In §2.2.1, we introduced two different methods to model collagen–HA interactions: one using SC scheme for modelling interpenetrating phases and the other using MT scheme based on matrix-inclusion geometry. Similarly, two modelling methods were discussed in §2.2.2 accounting for different types of interactions between fibrils and extrafibrillar HA foam. These, altogether, give rise to four sets of approaches to model bone at nano- and sub-microscale and, subsequently, at higher scales. Such approaches, called as approaches I, II, III, and IV, are summarized in table 3.

## 3. Results

### 3.1. Experiments

A trabecular bone specimen was imaged using μCT to provide an input for FEM analysis. Bone volume fraction was determined to be approximately 8.1 per cent. Other material characteristics calculated from μCT imaging are listed in table 4.

Compression test was used to obtain the stress–strain curves of four trabecular bone samples and, then, their Young's moduli were calculated. Results are tabulated in table 5. The first sample was digitized and used to create the μCT-based FEM model.

### 3.2. Modelling

#### 3.2.1. Nanoscale

We assume that 75 per cent of total HA crystals are interfibrillar and the remaining crystals (25%) form an extrafibrillar HA foam. Cylindrical collagen molecules as well as ellipsoidal mineral crystals are aligned along the axis 1 of the Cartesian coordinate system. Effective stiffness tensor of a mineralized collagen fibril, obtained using approaches I–IV, are listed in tables 6–9. The results are reported for both values of collagen modulus: 1.5 and 5.4 GPa (shown inside the parentheses). The corresponding engineering elastic constants, namely Young's moduli, shear moduli and Poisson's ratios, are also given in tables 10–13.

#### 3.2.2. Sub-microscale

##### 3.2.2.1. Mineralized collagen fibrils combined with an extrafibrillar hydroxyapatite foam

Cylindrical mineralized collagen fibrils, with the elastic constants given in tables 6–9, are combined with extrafibrillar HA foam to form coated fibrils. We assume that extrafibrillar mineral crystals comprise 25 per cent of the total HA. Again, four different approaches are used to find effective elastic constants of coated fibrils and the results are tabulated in tables 6–9 for both values of collagen Young's modulus. Tables 10–13 list the corresponding engineering moduli.

##### 3.2.2.2. Single lamella

Components of stiffness tensor of a single lamella are given in tables 6–9 for different modelling approaches I–IV. The results corresponding to collagen Young's modulus of 1.5 and 5.4 GPa are listed, respectively, outside and inside parentheses. Equivalently, the engineering elastic parameters of a single lamella are given in tables 10–13.

#### 3.2.3. Microscale

A large number of lamellae, with the elastic properties given in tables 6–9, are randomly oriented to form a single trabecula. Isotropic elastic stiffness components of a single strut (trabecula) are found for such lamellar arrangement, following the modelling procedure mentioned in §2.2.3, and are tabulated in tables 6–9 for different modelling approaches and two different values of collagen Young's modulus. Also, tables 10–13 list Young's modulus and Poisson's ratio of a trabecula.

#### 3.2.4. Mesoscale

The μCT-based finite-element model is used to determine the apparent Young's modulus of trabecular bone by using as inputs the material properties obtained at the previous level and the mesostructural geometry obtained from the μCT images. Computed elastic constants of trabecular bone are given in tables 6–9 for different modelling approaches and collagen Young's moduli. Young's modulus and Poisson's ratio of trabecular bone for different modelling approaches are also listed in tables 10–13. Figures 5 and 6 illustrate, respectively, FEM results for a displacement magnitude, , and a displacement component in the direction of loading, *u _{z}*. The displacement magnitude, which captures both axial and lateral displacement components, tends to increase spatially towards the centre of cylinder and goes to zero at the fixed bottom edge. The axial displacement is non-uniform, as seen in figure 6, which means that the bone structure behaves like a loose collection of weakly connected pieces of bone. Also, figure 7 shows an elastic strain energy density distribution. In this figure, localized hot spots indicate bending as opposed to uniform tension/compression. Furthermore, most of the elastic strain energy is concentrated at thin strut connections. These imply that bending is a primary deformation mode for low-density trabecular bone and show that trabecular connectivity has a large effect on the elastic properties of trabecular bone [95–98].

Moreover, to understand how trabecular bone architecture affects its elastic behaviour at mesoscale, we compare the results of our µCT-based FEM with some simpler models such as the Gibson model, the Christensen's LDM method and the MT scheme. The results are listed in table 14 for different modelling approaches I–IV and both values of collagen Young's modulus.

## 4. Discussion

We defined four distinct hierarchical levels in trabecular bone, namely the mineralized collagen fibril, single lamella, single trabecula and trabecular bone levels, and used several micromechanics methods, LCMT and FEM to obtain the elastic constants of trabecular bone at these different structural scales. Alternative methods could be used at each scale.

Table 15 lists Young's modulus of a single trabecula obtained by using our different modelling approaches and compares the results with those obtained by experiments. As mentioned before, different experimental techniques were used to measure the mechanical properties of individual trabeculae including the microtensile test [22–24], bending test [25,26], ultrasound [23,27,28] and nanoindentation [27,29–32,99]. Because of the complex geometry of a single trabecula, which has a curved shape and a varying cross section, machining bone samples for tensile and bending tests may cause significant surface defects [42]. Furthermore, the deformations are typically so small that any artefactual displacement, such as slipping at loading points, results in an underestimation of the measured value [42]. Therefore, the results of microtensile and bending tests of trabecular struts are lower than the true values [42]. That is the reason why the experimental data for a single trabecula listed in table 15 exclude the tensile and bending results, while includes the more reliable results of ultrasound and nanoindentation techniques. Note that most of the nanoindentation data on trabecular bone reported in literature are obtained using embedded samples in dry condition. Our analytical results for Young's modulus of a single strut are in a good agreement with experimental data. Gibson *et al*. [42] discussed different methods for measuring Young's modulus of wet human trabeculae in the longitudinal direction and concluded that 18 GPa is the best estimation for Young's modulus of solid trabeculae, which is the same as our modelling prediction for approach I and collagen Young's modulus of 1.5 GPa.

Also, the second part of table 15 compares the results of our µCT-based FEM analysis for different modelling approaches with experimental data. Since trabecular bone volume fraction significantly affects its Young's modulus, the values of bone volume fraction for experimental data are also included. The experimental data clearly show an increase in bone modulus with bone volume fraction. Our experimental results compare well with the experimental values reported in literature for trabecular bone of comparable porosity [102]. However, both our modelling and experimental results for Young's modulus of trabecular bone are on the lower side of the values reported in literature. This might be due to the fact that the specific bone samples used to create the FEM model and to do compression testing were from an old (88-year-old) donor and, consequently, were very porous. Another reason for the poor results of our experimental data may be the systematic errors in the platen compression test of trabecular bone, which occur owing to end artefacts [107]. Such errors lead to an underestimation of compressive Young's modulus of trabecular bone [107]. Also, the comparison between our modelling results and experimental data is performed at micro- and mesoscales but not at nano and sub-microscales, owing to lack of experimental data at lower scales.

We proposed different modelling approaches (approaches I–IV). The results, shown in table 15, suggest that the assumption of interpenetration for collagen and HA at nanoscale (in approaches I and II) gives rise to higher Young's moduli for trabecula and trabecular bone compared with the matrix-inclusion assumption (in approaches III and IV). When we use the SC method to obtain the effective elastic properties of mineralized collagen fibrils, the elastic properties of stiff HA are more dominant than the properties of soft collagen, considering the fact that their VFs are almost equal. However, the collagen properties contribute more to the overall elastic properties when the MT method is used rather than the SC method. Such a trend can be verified by the results obtained for different values of collagen Young's modulus: increasing collagen Young's modulus from 1.5 to 5.4 GPa increases Young's modulus of single trabecula, and accordingly of trabecular bone, more significantly in approaches III and IV (by 128 and 70%, respectively) than approaches I and II (by 40 and 27%, respectively).

Challenges and issues in creating FEM model might give rise to some errors in estimating Young's modulus of trabecular bone. The first challenge is how to set a threshold value for µCT images to accurately capture bone architecture and porosity. Different grey levels in the scan represent fully filled, partially filled or empty voxels. Thresholding determines whether a partially filled voxel is considered as bone or void. This might cause some errors in calculating a bone volume fraction which could be carried over into a finite-element model. Moreover, since side surfaces of cylindrical sample were rough, the volume of cylinder was approximated, which could cause inaccuracy in predicting bone's apparent Young's modulus by using equation (2.15). Also, at ends, we applied displacement boundary conditions for simplicity. However, more realistic boundary conditions should include friction. All the elements used in FEM were selected to be linear tetrahedral elements, for simplicity. Trabecular struts have evident bending behaviour giving rise to linear stresses, while tetrahedral elements are constant stress elements. However, we used a very fine mesh consisting of the largest number of elements manageable with our computer power that could minimize a possible error caused by the use of tetrahedral elements. We did a mesh convergence study in our previous work on modelling of trabecular bone as an idealized three-dimensional periodic cellular network [57]. The result showed that using five linear elements per width of struts improved the results by 13.17 per cent compared with the case of using two linear elements per width of struts. Finally, the resolution of our μCT image was selected to be 37 µm, based on sample dimensions, degree of X-ray penetration and limitations of the equipment used for imaging. Finer resolution would certainly capture better the trabecular bone architecture and lead to more accurate FEM predictions.

We made several simplifying assumptions and selections at different stages of modelling. First of all, the classification of scales in the way done here is not unique since the transitions between different hierarchies from nanoscale to macroscale in real bone are continuous rather than discrete. We accounted for four length scales while other choices could be made. Also, for simplicity, we assumed that an RVE exists at nano, sub-micro and microscales and that the features at the previous scale are much smaller than those at the next scale. In reality, bone features at previous scale are not necessarily infinitesimal relative to a larger scale.

An ambiguity exists in selecting elastic properties and VFs of collagen and HA from a wide range of values reported in literature [10], while little data are available for NCPs. Different choices may give rise to very different results. To illustrate this point, we used two different values for collagen Young's modulus, 1.5 and 5.4 GPa.

Furthermore, at first three scales, we made simplifying assumptions on bone's actual microstructures. At nanoscale, we did not model collagen cross-linking explicitly. Also, we assumed perfect bonding between collagen and HA and ignored the presence of an interphase layer in between them. Non-uniform distribution of HA crystals was also neglected in our model. At sub-microscale, fibrils were aligned unidirectionally, while in bone they are misaligned, but oriented in a preferential direction. We assumed that extrafibrillar HA foam fully filled spaces between mineralized fibrils. However, the foam might not fill all empty spaces between the fibrils leaving some voids. Also, we assumed that parallel layering of HA crystals in one collagen fibril is aligned with crystal layers in adjacent fibrils, while some TEM images suggest that crystals have a random arrangement rather than an orderly alignment in neighbouring collagen fibrils [108,109]. At microscale, we assumed bone tissue properties to be homogeneous and isotropic. Nevertheless, several studies showed that elastic properties of trabecular bone at this level are heterogeneous and change across thickness of a single trabecula (strut) [31,48,110–113]. Harrison *et al.* [31] performed nanoindentation along the strut length and across the strut width of trabecular samples and concluded that Young's modulus across all specimens varied from 3.14 GPa at strut exteriors to 19.75 GPa at strut centres. Such results verify the trend of gradual reduction in tissue modulus from centre to surface of trabeculae, where older bone at strut interior is more mineralized compared with the newly formed bone at strut surface. Also, for simplicity we assumed that bone tissue is isotropic at microscale, while some experimental data indicate anisotropy [114,115].

Since the main focus of this paper is on a multi-scale modelling approach, for simplicity, we used a single realization of trabecular bone architecture in the analysis. Then, we compared our computational results with the experimental data obtained using the same specimen. This was motivated by the fact that properties of trabecular bone are dependent on its architecture. More comprehensive analyses, considering a number of realizations, could be done in future. Also, nanoindentation technique could be employed to measure local elastic properties of trabecular bone samples used in the experiment to verify our theoretical predictions at microscale. Chevalier *et al*. [30] used such approach, a combination of theory (µCT-based FEM) and experiments (nanoindentation technique and macroscopic mechanical testing) to characterize trabecular bone at tissue and bone levels.

All the uncertainties and open issues discussed above along with some other factors ignored in this work make the multi-scale modelling of trabecular bone a challenging problem with much potential for future studies. Also, the modelling process discussed here is specifically applicable to a healthy mature trabecular bone. However, the proposed multi-scale model can capture elastic behaviour of aged or diseased trabecular bone by changing the input parameters and bone microstructures at different scales. Ideally, a comprehensive multi-scale experimental characterization of bone should be performed to provide accurate inputs for modelling.

## 5. Conclusions

We modelled trabecular bone as a hierarchical material and predicted its effective Young's moduli. Our analysis involved a multi-scale approach, starting with nanoscale (mineralized fibril level) and moving up the scales to sub-microscale (single lamella level), microscale (single trabecula level) and finally mesoscale (trabecular bone level). The selection of scales is not unique and other choices could be made. At nanoscale and sub-microscale, we used two alternative approaches, one assuming interpenetrating phases and the other involving matrix-inclusion geometry. At mesoscale, we used µCT-based finite-element model and compared the results with those obtained using micromechanics methods. We also considered two different values of Young's modulus of collagen, owing to a wide range of data available in literature. In the analysis, we used the models of micromechanics, LCMT and FEM. Good agreement was found between theoretical and experimental results. This research sets a framework for multi-scale modelling of materials with hierarchical structures.

## Acknowledgements

We acknowledge support of the National Science Foundation (CMMI 09–27909 ARRA, to Dr Ken Chong). We also thank Prof. William Hutton, Emory University Spine Center, for valuable discussions and sharing his laboratory for compression testing.

- Received November 23, 2011.
- Accepted January 4, 2012.

- This journal is © 2012 The Royal Society