Can the diverse elastic properties of trabecular and cortical bone be attributed to only a few tissue-independent phase properties and their interactions?

Abstract

As candidates for tissue-independent phases of cortical and trabecular bone we consider (i) hydroxyapatite, (ii) collagen, (iii) ultrastructural water and non-collagenous organic matter, and (iv) marrow (water) filling the Haversian canals and the intertrabecular space. From experiments reported in the literature, we assign stiffness properties to these phases (experimental set I).

On the basis of these phase definitions, we develop, within the framework of continuum micromechanics, a two-step homogenization procedure: (i) at a length scale of 100–200 nm, hydroxyapatite (HA) crystals build up a crystal foam (“polycrystal”), and water and non-collagenous organic matter fill the intercrystalline space (homogenization step I); (ii) at the ultrastructural scale of mineralized tissues (5–10 microns), collagen assemblies composed of collagen molecules are embedded into the crystal foam, acting mechanically as cylindrical templates. At an enlarged material scale of 5–10 mm, the second homogenization step also accommodates the micropore space as cylindrical pore inclusions (Haversian and Volkmann canals, inter-trabecular space) that are suitable for both trabecular and cortical bone. The inputs for this micromechanical model are the tissue-specific volume fractions of HA, collagen, and of the micropore space. The outputs are the tissue-specific ultrastructural and microstructural (=macroscopic=apparent) elasticity tensors.

A second independent experimental set (composition data and experimental stiffness values) is employed to validate the proposed model. We report a small mean prediction error for the macroscopic stiffness values. The validation suggests that hydroxyapatite, collagen, and water are tissue-independent phases, which define, through their mechanical interaction, the elasticity of all bones, whether cortical or trabecular.

Appendix: Continuum micromechanics relationships for the definition of the two-step homogenization procedure of Figure 1

Homogenization step I

As for the polycrystal, RVE $\hat{V}_{{\text{p}}}$, the average local strains in the two phases, ε _HA and ε _uw, are related to the homogeneous strains E _p imposed at the boundary of $\hat{V}_{{\text{p}}}$, by the average of the localization tensor over $\hat{V}_{{{\text{HA}}}}$ and $\hat{V}_{{{\text{uw}}}}$, A _HA and A _uw,

$$ \begin{array}{*{20}c} {{{\user2{\varepsilon }}_{{{\text{HA}}}} = {\mathbf{A}}_{{{\text{HA}}}} :{\mathbf{E}}_{{\text{p}}} ,}} & {{{\user2{\varepsilon }}_{{{\text{uw}}}} = {\mathbf{A}}_{{{\text{uw}}}} :{\mathbf{E}}_{{\text{p}}} }} \\ \end{array} $$

(A1)

An estimate for the homogenized stiffness tensor is given by the classical relation of micromechanics (Zaoui 1997),

$${\mathbf{C}}^{{{\text{est}}}}_{{\text{p}}} = {\left\langle {{\mathbf{c}}:{\mathbf{A}}^{{{\text{est}}}} } \right\rangle }_{{\hat{V}_{{\text{p}}} }} = {\sum\limits_r {\hat{f}_{r} {\mathbf{c}}_{r} :{\mathbf{A}}^{{{\text{est}}}}_{r} } }$$

(A2)

where r∈[HA,uw] is the phase index, ${\left\langle {{\left( . \right)}} \right\rangle }_{V} = 1/V{\int {_{V} {\left( . \right)}dV} }$ stands for the volume average, c _r =3J k _r +2K µ _r denotes the (isotropic) stiffness tensor of phase r∈[HA;uw]; J, J _ijkl =1/3δ _ij δ _kl , is the volumetric part of the fourth-order unity tensor I, I _ijkl =1/2(δ _ik δ _jl +δ _il δ _kj ). δ _ij is the Kronecker delta; in other words δ _ij =1 for i=j, and δ _ij =0 for ij. The deviatoric part of I is K=I–J. $\hat{f}_{r} = \frac{{\hat{V}_{r} }}{{\hat{V}_{{\text{p}}} }}$ is the volume fraction of phase r, $\hat{f}_{{{\text{HA}}}} + \hat{f}_{{{\text{uw}}}} = 1$, and ${\mathbf{A}}^{{{\text{est}}}}_{r}$ is an estimate for the localization tensor of phase r. For the polycrystal estimate, ${\mathbf{A}}^{{{\text{est}}}}_{r}$ is given by an implicit relationship,

$${\mathbf{A}}^{{{\text{est}}}}_{r} = {\left[ {{\mathbf{I}} + {\mathbf{S}}^{{{\text{Esh}}}}_{r} :{\left( {{\mathbf{C}}^{{{\text{est,}} - 1}}_{{\text{p}}} :{\mathbf{c}}_{r} - {\mathbf{I}}} \right)}} \right]}^{{ - 1}} :{\left\langle {{\left[ {{\mathbf{I}} + {\mathbf{S}}^{{{\text{Esh}}}} :{\left( {{\mathbf{C}}^{{{\text{est,}} - 1}}_{{\text{p}}} :{\mathbf{c}} - {\mathbf{I}}} \right)}} \right]}^{{ - 1}} } \right\rangle }^{{ - 1}}_{{\hat{V}_{{\text{p}}} }} $$

(A3)

Eshelby’s tensor S ^esh for spherical inclusions in a matrix with ${\mathbf{C}}^{{{\text{est}}}}_{{\text{p}}}$ reads as (p. 300 of Zaoui 1997):

$${\mathbf{S}}^{{{\text{Esh}}}}_{{{\text{HA}}}} = {\mathbf{S}}^{{{\text{Esh}}}}_{{{\text{uw}}}} = {\mathbf{S}}^{{{\text{Esh}}}}_{{{\text{sph}}}} = \alpha ^{{{\text{est}}}} {\mathbf{J}} + \beta ^{{{\text{est}}}} {\mathbf{K}}$$

(A4)

with

$$\begin{array}{*{20}c} {{\alpha ^{{{\text{est}}}} = \frac{{3k^{{{\text{est}}}}_{{\text{p}}} }}{{3k^{{{\text{est}}}}_{{\text{p}}} + 4\mu ^{{{\text{est}}}}_{{\text{p}}} }},}} & {{\beta ^{{{\text{est}}}} = \frac{{6{\left( {k^{{{\text{est}}}}_{{\text{p}}} + 2\mu ^{{{\text{est}}}}_{{\text{p}}} } \right)}}}{{5{\left( {3k^{{{\text{est}}}}_{{\text{p}}} + 4\mu ^{{{\text{est}}}}_{{\text{p}}} } \right)}}}}} \\ \end{array} $$

(A5)

For the solution of this implicit problem we refer to Hellmich and Ulm (2002b).

Homogenization step II

The “global” microstructural strains E _m are related to the average local strains in the phases of the RVE V _m (Fig. 1c) by the average concentration tensors A _M, A _col, and A _por,

$$ \begin{array}{*{20}c} {{{\user2{\varepsilon }}_{{\text{M}}} = {\mathbf{A}}_{{\text{M}}} :{\text{E}}_{{\text{m}}} ,}} & {{{\user2{\varepsilon }}_{{{\text{col}}}} = {\mathbf{A}}_{{{\text{col}}}} :{\text{E}}_{{\text{m}}} ,}} & {{{\user2{\varepsilon }}_{{{\text{por}}}} = {\mathbf{A}}_{{{\text{por}}}} :{\mathbf{E}}_{{\text{m}}} }} \\ \end{array} $$

(A6)

The estimate for the homogenized stiffness tensor of the microstructure is given by

$${\mathbf{C}}^{{{\text{est}}}}_{{\text{m}}} = {\left\langle {{\mathbf{c}}:{\mathbf{A}}^{{{\text{est}}}} } \right\rangle }_{{V_{{\text{m}}} }} = {\sum\limits_r {f_{r} {\mathbf{c}}_{r} :{\mathbf{A}}^{{{\text{est}}}}_{r} } }$$

(A7)

where r∈[col;M;por] is the phase index, c _r denotes the stiffness tensor of phase r; f _r =V _r /V _m is the volume fraction of phase r. The matrix stiffness is known from the previous homogenization step, ${\mathbf{c}}_{{\text{M}}} = {\mathbf{C}}^{{{\text{est}}}}_{{\text{p}}} {\left( {\hat{f}_{{{\text{HA}}}} } \right)}$ (Eq. A2). Using the Mori-Tanaka (1973) scheme, the localization tensors for two families of inclusions (collagen and pores) in the HA polycrystal matrix can be estimated as

$$\begin{array}{*{20}l} {{{\mathbf{A}}^{{{\text{est}}}}_{{{\text{col}}}} } \hfill} & { = \hfill} & {{{\left[ {{\mathbf{I}} + {\mathbf{S}}^{{{\text{Esh}}}}_{{{\text{col}}}} :{\left( {{\mathbf{c}}^{{ - 1}}_{{\text{M}}} :{\mathbf{c}}_{{{\text{col}}}} - {\mathbf{I}}} \right)}} \right]}^{{ - 1}} :} \hfill} \\ {{} \hfill} & {{} \hfill} & {{{\left\langle {{\left[ {{\mathbf{I}} + {\mathbf{S}}^{{{\text{Esh}}}} :{\left( {{\mathbf{c}}^{{ - 1}}_{{\text{M}}} :{\mathbf{c}} - {\mathbf{I}}} \right)}} \right]}^{{ - 1}} } \right\rangle }^{{ - 1}}_{{V_{{\text{m}}} }} } \hfill} \\ \end{array} $$

(A8)

$$\begin{array}{*{20}l} {{{\mathbf{A}}^{{{\text{est}}}}_{{{\text{por}}}} } \hfill} & { = \hfill} & {{{\left[ {{\mathbf{I}} + {\mathbf{S}}^{{{\text{Esh}}}}_{{{\text{por}}}} :{\left( {{\mathbf{c}}^{{ - 1}}_{{\text{M}}} :{\mathbf{c}}_{{{\text{por}}}} - {\mathbf{I}}} \right)}} \right]}^{{ - 1}} :} \hfill} \\ {{} \hfill} & {{} \hfill} & {{{\left\langle {{\left[ {{\mathbf{I}} + {\mathbf{S}}^{{{\text{Esh}}}} :{\left( {{\mathbf{c}}^{{ - 1}}_{{\text{M}}} :{\mathbf{c}} - {\mathbf{I}}} \right)}} \right]}^{{ - 1}} } \right\rangle }^{{ - 1}}_{{V_{{\text{m}}} }} } \hfill} \\ \end{array} $$

(A9)

$${\mathbf{A}}^{{{\text{est}}}}_{{\text{M}}} = {\left\langle {{\left[ {{\mathbf{I}} + {\mathbf{S}}^{{{\text{Esh}}}} :{\left( {{\mathbf{c}}^{{ - 1}}_{{\text{M}}} :{\mathbf{c}} - {\mathbf{I}}} \right)}} \right]}^{{ - 1}} } \right\rangle }^{{ - 1}}_{{V_{{\text{m}}} }} $$

(A10)

For cylindrical (collagen or micropore) inclusions in the polycrystal matrix, the non-zero components of Eshelby’s tensor ${\mathbf{S}}^{{{\text{Esh}}}}_{{{\text{col}}}} = {\mathbf{S}}^{{{\text{Esh}}}}_{{{\text{por}}}} = {\mathbf{S}}^{{{\text{Esh}}}}_{{{\text{cyl}}}}$ read as

$$\begin{array}{*{20}l} {{S^{{{\text{Esh}}}}_{{{\text{1111}}}} } \hfill} & { = \hfill} & {{S^{{{\text{Esh}}}}_{{{\text{2222}}}} } \hfill} & { = \hfill} & {{\frac{9}{4}\frac{{k_{{\text{M}}} + \mu _{{\text{M}}} }}{{3k_{{\text{M}}} + 4\mu _{{\text{M}}} }}} \hfill} \\ {{S^{{{\text{Esh}}}}_{{{\text{1122}}}} } \hfill} & { = \hfill} & {{S^{{{\text{Esh}}}}_{{{\text{2211}}}} } \hfill} & { = \hfill} & {{\frac{1}{4}\frac{{3k_{{\text{M}}} + 5\mu _{{\text{M}}} }}{{3k_{{\text{M}}} + 4\mu _{{\text{M}}} }}} \hfill} \\ {{S^{{{\text{Esh}}}}_{{{\text{1133}}}} } \hfill} & { = \hfill} & {{S^{{{\text{Esh}}}}_{{{\text{2233}}}} } \hfill} & { = \hfill} & {{\frac{1}{2}\frac{{3k_{{\text{M}}} + 2\mu _{{\text{M}}} }}{{3k_{{\text{M}}} + 4\mu _{{\text{M}}} }}} \hfill} \\ {{S^{{Esh}}_{{1212}} } \hfill} & { = \hfill} & {{S^{{{\text{Esh}}}}_{{{\text{2121}}}} } \hfill} & { = \hfill} & {{\frac{1}{4}\frac{{3k_{{\text{M}}} + 7\mu _{{\text{M}}} }}{{3k_{{\text{M}}} + 4\mu _{{\text{M}}} }}} \hfill} \\ {{S^{{{\text{Esh}}}}_{{{\text{1313}}}} } \hfill} & { = \hfill} & {{S^{{{\text{Esh}}}}_{{{\text{3131}}}} } \hfill} & { = \hfill} & {{S^{{{\text{Esh}}}}_{{{\text{2323}}}} \; = \;S^{{{\text{Esh}}}}_{{{\text{3232}}}} \; = \;\frac{1}{4}} \hfill} \\ \end{array} $$

(A11)

Cylindrical inclusions are to be understood in the sense of Eshelby (1957), p. 386, as inclusions of prolate spheroid type with one axis of the ellipsoid being very much longer than the other two, which are of the same length. Homogenization of the ultrastructure (Fig. 1b) is performed by analogy to Eqs. A6–A11, replacing subscript “m” by “u”, f _r by $\bar{f}_{r}$, and r∈[col,M,por] by r∈[col,M].

We have also tested the possible significance of a deviation from the cylindrical pore shape, by considering inclusions of the prolate spheroid type with length a, oriented in the x ₃-direction (Fig. 1), and width b=c<a; or, in terms of a morphology parameter, ω=b/a=c/a<1. The limit case ω→0 refers to cylindrical inclusions. For prolate spheroid inclusions, the non-zero components of Eshelby’s tensor ${\mathbf{S}}^{{{\text{Esh}}}}_{{{\text{col}}}} = {\mathbf{S}}^{{{\text{Esh}}}}_{{{\text{por}}}} = {\mathbf{S}}^{{{\text{Esh}}}}_{{{\text{prsph}}}}$ read as (Eshelby 1957):

$$\begin{array}{*{20}l} {{S^{{{\text{Esh}}}}_{{{\text{1111}}}} } \hfill} & { = \hfill} & {{S^{{{\text{Esh}}}}_{{{\text{2222}}}} } \hfill} & { = \hfill} & {{Qb^{2} I_{{{\text{bb}}}} + RI_{{\text{b}}} } \hfill} \\ {{S^{{{\text{Esh}}}}_{{{\text{3333}}}} } \hfill} & { = \hfill} & {{Qa^{2} I_{{{\text{aa}}}} + RI_{{\text{a}}} } \hfill} & {{} \hfill} & {{} \hfill} \\ {{S^{{{\text{Esh}}}}_{{{\text{1122}}}} } \hfill} & { = \hfill} & {{S^{{{\text{Esh}}}}_{{{\text{2211}}}} } \hfill} & { = \hfill} & {{Qb^{2} I_{{{\text{bc}}}} - RI_{{\text{b}}} } \hfill} \\ {{S^{{{\text{Esh}}}}_{{{\text{1133}}}} } \hfill} & { = \hfill} & {{S^{{{\text{Esh}}}}_{{{\text{2233}}}} } \hfill} & { = \hfill} & {{Qa^{2} I_{{{\text{ab}}}} - RI_{{\text{b}}} } \hfill} \\ {{S^{{{\text{Esh}}}}_{{{\text{3322}}}} } \hfill} & { = \hfill} & {{S^{{{\text{Esh}}}}_{{{\text{3311}}}} } \hfill} & { = \hfill} & {{Qb^{2} I_{{{\text{ab}}}} - RI_{{\text{a}}} } \hfill} \\ {{S^{{{\text{Esh}}}}_{{{\text{1212}}}} } \hfill} & { = \hfill} & {{S^{{{\text{Esh}}}}_{{{\text{2121}}}} } \hfill} & { = \hfill} & {{Qb^{2} I_{{{\text{bc}}}} - RI_{{\text{b}}} } \hfill} \\ {{S^{{{\text{Esh}}}}_{{{\text{1313}}}} } \hfill} & { = \hfill} & {{S^{{{\text{Esh}}}}_{{{\text{3131}}}} } \hfill} & { = \hfill} & {{S^{{{\text{Esh}}}}_{{{\text{2323}}}} \; = \;S^{{{\text{Esh}}}}_{{{\text{3232}}}} \; = \;Q\frac{1}{2}{\left( {a^{2} + b^{2} } \right)}I_{{{\text{ab}}}} + R\frac{1}{2}{\left( {I_{{\text{a}}} + I_{{\text{b}}} } \right)}} \hfill} \\ \end{array} $$

(A12)

with an arbitrarily chosen value for a, and

$$\begin{array}{*{20}l} {b \hfill} & { = \hfill} & {{\omega \times a} \hfill} \\ {{I_{{\text{b}}} } \hfill} & { = \hfill} & {{\frac{{2\pi ac^{2} }}{{{\left( {a^{2} - c^{2} } \right)}^{{\frac{3}{2}}} }}{\left\{ {\frac{a}{c}{\left( {\frac{{a^{2} }}{{c^{2} }} - 1} \right)}^{{\frac{1}{2}}} - \cosh ^{{ - 1}} \frac{a}{c}} \right\}},\;\cosh ^{{ - 1}} {\left( {\cosh x} \right)} = x} \hfill} \\ {{I_{{\text{a}}} } \hfill} & { = \hfill} & {{4\pi - 2I_{{\text{b}}} } \hfill} \\ {{I_{{{\text{ab}}}} } \hfill} & { = \hfill} & {{\frac{{I_{{\text{b}}} - I_{{\text{a}}} }}{{3{\left( {a^{2} - b^{2} } \right)}}}} \hfill} \\ {{I_{{{\text{bc}}}} } \hfill} & { = \hfill} & {{\frac{\pi }{{3b^{2} }} - \frac{1}{4}I_{{{\text{ab}}}} } \hfill} \\ {{I_{{{\text{aa}}}} } \hfill} & { = \hfill} & {{\frac{{4\pi }}{{3a^{2} }} - 2I_{{{\text{ab}}}} } \hfill} \\ {{I_{{{\text{bb}}}} } \hfill} & { = \hfill} & {{\frac{{4\pi }}{{3b^{2} }} - I_{{{\text{ab}}}} - I_{{{\text{bc}}}} } \hfill} \\ {Q \hfill} & { = \hfill} & {{\frac{3}{{8\pi {\left( {1 - \nu _{{\text{M}}} } \right)}}}} \hfill} \\ {R \hfill} & { = \hfill} & {{\frac{{1 - 2\nu _{{\text{M}}} }}{{8\pi {\left( {1 - \nu _{{\text{M}}} } \right)}}}} \hfill} \\ \end{array} $$

(A13)

where ν _M=(3k _M–2µ _M)/(6k _M+2µ _M), Poisson’s ratio of the matrix.