The poro-viscoelastic properties of trabecular bone: a micro computed tomography-based finite element study

journal of the mechanical behavior of biomedical materials 44 (2015) 1 –9

Available online at www.sciencedirect.com

www.elsevier.com/locate/jmbbm

Research Paper

The poro-viscoelastic properties of trabecular bone: a micro computed tomography-based finite element study Clara Sandinoa,b, David D. McErlaina,b, John Schipilowa,c, Steven K. Boyda,b,c,d,n a

McCaig Institute for Bone and Joint Health, University of Calgary, Canada Department of Radiology, Faculty of Medicine, University of Calgary, Canada c Roger Jackson Centre for Health and Wellness Research, University of Calgary, Canada d Schulich School of Engineering, University of Calgary, Canada b

art i cle i nfo

ab st rac t

Article history:

Bone is a porous structure with a solid phase that contains hydroxyapatite and collagen. Due to

Received 18 October 2014

its composition, bone is often represented either as a poroelastic or as a viscoelastic material;

Received in revised form

however, the poro-viscoelastic formulation that allows integrating the effect of both the fluid

15 December 2014

flow and the collagen on the mechanical response of the tissue, has not been applied yet. The

Accepted 18 December 2014

objective of this study was to develop a micro computed tomography (mCT)-based finite element

Available online 24 December 2014

(FE) model of trabecular bone that includes both the poroelastic and the viscoelastic nature of

Keywords:

the tissue. Cubes of trabecular bone (N¼25) from human distal tibia were scanned with mCT and

Trabecular bone

stress relaxation experiments were conducted. The mCT images were the basis for sample

Poro-viscoelasticity

specific FE models, and the stress relaxation experiments were simulated applying a poro-

Micro computed tomography

viscoelastic formulation. The model considers two scales of the tissue: the intertrabecular pore

Finite element models

and the lacunar-canalicular pore scales. Independent viscoelastic and poroelastic models were also developed to determine their contribution to the poro-viscoelastic model. All the experiments exhibited a similar relaxation trend. The average reaction force before relaxation was 9.28  102 N (SD75.11  102 N), and after relaxation was 4.69  102 N (SD72.88  102 N). The slope of the regression line between the force before and after relaxation was 1.92 (R2 ¼0.96). The poro-viscoelastic models captured 49% of the variability of the experimental data before relaxation and 33% after relaxation. The relaxation predicted with viscoelastic models was similar to the poro-viscoelastic ones; however, the poroelastic formulation underestimated the reaction force before relaxation. These data suggest that the contribution of viscoelasticity (fluid flow-independent mechanism) to the mechanical response of the tissue is significantly greater than the contribution of the poroelasticity (fluid flow-dependent mechanism). & 2015 Elsevier Ltd. All rights reserved.

n Corresponding author at: McCaig Institute for Bone and Joint Health, University of Calgary, 3280 Hospital Drive NW, Room HRIC 3AC64, Calgary, Alberta, T2N 4Z6. Tel.: þ1 403 220 3664; fax: þ1 403 210 9573. E-mail address: [email protected] (S.K. Boyd).

http://dx.doi.org/10.1016/j.jmbbm.2014.12.018 1751-6161/& 2015 Elsevier Ltd. All rights reserved.

2

1.

journal of the mechanical behavior of biomedical materials 44 (2015) 1 –9

Introduction

Trabecular bone is a biological tissue with three main hierarchical levels of porosity: At the intertrabecular pore level, the tissue consists of an irregular network of trabeculae with blood and marrow filling the pores; at the lacunar-canalicular pore level, the pores enclose the bone cells (osteocytes) and the interstitial fluid; at the collagen-apatite level, pores with water surround the organic collagen and the inorganic hydroxyapatite (Fritton and Weinbaum, 2009). The hydroxyapatite provides the hardness and the rigidity to the tissue while the collagen provides its toughness (Boskey and Coleman, 2010; Viguet-Carrin et al., 2006). The viscoelasticity of the tissue has been associated with both the layers of intercrystalline water within the extrafibrillar clusters of hydroxyapatite (Eberhardsteiner et al., 2014) and the arrangement of the molecules of collagen and water (Shen et al., 2011). The mechanical stimuli in the marrow and the interstitial fluid are associated with bone adaptation and regeneration (Fritton and Weinbaum, 2009; Gurkan and Akkus, 2008). Due to the porous nature of bone tissue and its intrinsic viscoelasticity, trabecular bone can be characterized as a poro-viscoelastic material. The biphasic poro-viscoelastic formulation was introduced by Mak (1986a) to study the contribution from the intrinsic viscoelasticity of the matrix and the interstitial fluid flow to the apparent viscoelastic behavior of cartilage. This formulation has been used to study cartilage (Mak, 1986a; Mak, 1986b; Setton et al., 1993) and hydrogels for bone tissue engineering (Kalyanam et al., 2009; Noailly et al., 2008; Strange et al., 2013); however, we are not aware of it being applied to bone tissue. At the intertrabecular pore scale, trabecular bone has been modeled as either as a poroelastic or as a viscoelastic material. The poroelastic theory (Biot, 1941) has been used to study the role of marrow and interstitial fluid flow (Cowin, 1999; Lim and Hong, 1998; Lim and Hong, 2000; Ochoa et al., 1991a, 1991b) and to explain the time-dependent material properties of the tissue (Carter and Hayes, 1977; Hong and Song, 1998; Hong et al., 2001). The poroelastic models assume elasticity in the solid phase, in spite of the inherent viscoelasticity of bone. On the other hand, the viscoelastic formulation has been used in studies of creep (Bowman et al., 1994; Moore et al., 2004) and stress relaxation (Deligianni et al., 1994; Quaglini et al., 2009; Sasaki et al., 1993). Several studies applying the microporomechanics principles have been conducted lately with the objective of including the role of the structure at the collagen-apatite pore scale in the mechanical behavior of the bone (Brynk et al., 2011; Eberhardsteiner et al., 2014; Fritsch et al., 2009; Hellmich and Ulm, 2005; Hellmich et al., 2009). Hellmich and Ulm (2005) were able to determine the poroelastic properties of trabecular tissue as a function of the micropore space and the collagen and hydroxyapatite volume fractions. Recently, in order to investigate the mechanical properties of the mineral and the protein phases, Chen and McKittrick (2011) conducted a study of demineralized and deproteinized trabecular bone and showed that removing either the organic or the inorganic material significantly reduced the strength of the tissue. The mechanical response of trabecular bone depends also

models using micro computed tomography (mCT) images (Müller et al., 1994; van Rietbergen et al., 1998). However, in most of the mCT-based studies, trabecular bone has been characterized as linear elastic or elasto-plastic in order to determine bone strength (Bevill and Keaveny, 2009; Guillén et al., 2011; MacNeil and Boyd, 2008; Sanyal et al., 2012). Characterization of the mechanical properties of bone tissue is important to fully understand its mechanical behavior, and for the study of the effect of bone-related diseases on the biomechanical response of the tissue. The objective of this study was to develop a mCT-based FE model of trabecular bone at the intertrabecular pore scale that accounts for the poro-viscoelasticity of the tissue. This formulation provides the ability to include the intrinsic viscoelastic response of the bone matrix in addition to the poroelastic nature of the tissue. Validation will be achieved comparing sample specific mCT-based FE models with mechanical tests.

2.

Methods

2.1.

Mathematical description of the model

The constitutive equations for a linear isotropic poroelastic material can be written as follows (Biot, 1941; Rice and Cleary, 1976): The constitutive equation for the porous solid is:
 2G εkk  δij αp ð1Þ σ ij ¼ 2Gεij þ δij K 3 where σij and εij are the stress and strain tensors, p is the pore pressure, K and G are the bulk and the shear moduli, α is the Biot coefficient of effective stress, and δij is the Kronecker delta. The Biot coefficient can be expressed as: α ¼ 1

K Ks

ð2Þ

where Ks is the bulk modulus of the solid. The constitutive equation for the fluid is: ζ ¼ αεkk þ

α2 p Ku  K

ð3Þ

where ζ is the variation of fluid content per unit volume of porous media and Ku is the undrained bulk modulus. The undrained bulk modulus can be expressed as: Ku ¼ K þ

α2 Ks Kf φKs þ ðα φÞKf

ð4Þ

where φ is the pore volume fraction and Kf is the bulk modulus of the fluid. The Darcy's law is expressed as:
 k ∂p ð5Þ qi ¼  μf ∂xi where qi is the rate of fluid volume crossing a unit area of porous solid, k the intrinsic permeability having dimension of length squared and μf the fluid viscosity. The equilibrium equation is defined as:

on its micro architecture. Recently, bone studies have included

3 X ∂σ ij

the micro architecture of the tissue into finite element (FE)

i¼1

∂xi

¼0

ð6Þ

journal of the mechanical behavior of biomedical materials 44 (2015) 1 –9

and the continuity equation as: 3 ∂ζ X ∂qi þ ¼0 ∂t i ¼ 1 ∂xi

ð7Þ

The poro-viscoelastic formulation is derived from the poroelastic model assuming that the solid matrix is linear viscoelastic rather than linear elastic. Defining the bulk and shear relaxation functions as: ^ ¼ K1 þ KðtÞ

nK X

Ki e  t=τi

K

ð8Þ

Gi e  t=τi

ð9Þ

i¼1

^ ¼ G1 þ GðtÞ

nG X

G

i¼1

where K1 and G1 are the bulk and shear moduli at equilibrium, Ki and Gi are the bulk and shear relaxation moduli, and τKi and τGi are the relaxation times, the constitutive equations for poroviscoelasticity can be written as follows (Abousleiman et al., 1996): The constitutive equation for the porous solid: ! Z t Z t 0 ^ 0 ^ ^ t0 Þ 2Gðt t Þ dεkk ðt0 Þ Þdεij ðt0 Þ þ δij Gðtt Kðt σ ij ðtÞ ¼ 2 3 0 0 Z t α^ ðt t0 Þdpðt0 Þ ð10Þ δij 0

The constitutive equation for the fluid: Z t Z t α^ ðt t0 Þ2 ^ dpðt0 Þ αðt t0 Þdεkk ðt0 Þ þ ζðtÞ ¼ 0 0 ^ ^ Þ 0 0 Ku ðt t Þ  Kðtt The Darcy's law:
 k ∂pðtÞ qi ðtÞ ¼  μf ∂xi

ð11Þ

ð12Þ

The equilibrium equation: 3 X

∂σ ij ðtÞ ¼0 ∂xi i¼1

ð13Þ

The continuity equation: 3 ∂ζðtÞ X ∂qi ðtÞ þ ¼0 ∂t ∂xi i¼1

3

The cubes were scanned with mCT with an isotropic voxel size of 20 mm (Scanco Medical mCT 35, Brüttisellen, Switzerland; 70 kVp, 114 mA, 400 ms integration time). The images were smoothed (Gaussian filter; sigma¼ 0.8, support¼ 1.0) and segmented (threshold of 13% of the maximum grey value) (Fig. 1). The subsequent 3D images were re-scaled to an isotropic voxel size of 80 mm in order to reduce the average number of voxels from 125 to 2 million per model. Disconnected voxels due to image noise were removed. The bone volume ratio (BV/TV) was determined and the error of the rescaling was computed as: e¼

BV=TV20 mm  BV=TV80 mm BV=TV20 mm

ð17Þ

where e is the error, and BV/TV20 mm and BV/TV80 mm are the BV/TV computed for the 20 mm and 80 mm images, respectively. The cubes with an error higher than 1% were excluded from the study. The final sample size was 25 cubes (obtained from 15 tibiae). After the scanning the cubes were wrapped in gauze soaked in saline solution, to keep them hydrated, and frozen at 30 1C.

2.3.

Mechanical experiments

The cubes were thawed for 3 h at room temperature (24 1C) before the experiment and were immersed in saline solution during the test. Stress relaxation experiments were carried out on an MTS testing machine (Mini Bionix 858, MTS Corp.) equipped with a 10 kN load cell. A compressive deformation of 0.1 mm was applied at a strain rate of 1 mm/s and then maintained for 600 s (Fig. 2). The strain was applied in the direction aligned with the long axis of the tibiae (approximately the principal orientation of the trabeculae). The test was repeated with two specimens after 1 h of recovery, in order to maximize the repeatability of the experiment. After the loading experiment, a second scan was performed and the images obtained before and after the mechanical test were registered in order to ensure that trabecular cracks were

ð14Þ

where the Biot coefficient is defined as: αðtÞ ^ ¼ 1

^ KðtÞ Ks

ð15Þ

and the undrained bulk modulus as: K^ u ðtÞ ¼ K^ ðtÞ þ

2.2.

a^ ðtÞ2 Ks Kf φKs þ ða^ ðtÞ  φÞKf

ð16Þ

Sample preparation and lCT imaging

Human trabecular bone cubes of 10 mm side-length (n¼ 33) were obtained from 17 distal tibiae of cadavers (49 to 85 years of age, average¼ 64.35 years) using a linear precision saw (IsoMet 5000, Buehler, Canada). Bone tissue samples were provided by the Gross Anatomy Laboratory and the Joint Transplantation Program at the University of Calgary. All experimental procedures were approved by the Conjoint Health Research Ethics Board at the University of Calgary.

Fig. 1 – 3D image reconstruction of a representative specimen of trabecular bone (BV/TV¼28.5%). Cubes of 10 mm side length were cut from human distal tibiae and were scanned with micro computed tomography (lCT; nominal voxel size 20 lm). The lCT images were used for FE models.

4

journal of the mechanical behavior of biomedical materials 44 (2015) 1 –9

not visible. Finally, the experimental reaction force during relaxation was smoothed applying a moving average filter of 100 points. The reaction force was normalized as: NFðtÞ ¼

FðtÞ2Fð600Þ Fð600Þ

ð18Þ

where NF(t) is the normalized force and F(t) is the total reaction force at time t. The ratio between the force before and after relaxation, F(0.1)/F(600), was determined and Pearson's correlations between this ratio, the donors' age and the morphometric variables were determined. Correlations were considered statistically significant at po0.05. The correlations and their statically significance were determined in Minitab 16 (Minitab Inc, USA).

2.4.

FE Modeling

The meshes for the FE models consisted of eight-noded, cubic hexahedral elements of 80 mm with two sets of material properties discretized: the bone and the marrow phase. These sets were defined directly from the 3D images using a voxelbased method. The bone and the marrow phase were simulated with the poroelastic formulation (Biot, 1941) adding a viscoelastic component to the bone phase. Poroelastic material properties from the mechano-regulation approach of Lacroix and Prendergast (2002) used to simulate fracture healing were used (Table 1). For the viscoelastic component, a generalized Maxwell model was applied: An exponential decay equation of two terms was adjusted to the average curve of the normalized reaction force determined experimentally by using the iterative method. Finally the corresponding coefficients of a Prony series were determined. The simulated stress relaxation test applied a compressive displacement of 0.1 mm in the z direction to the nodes on the top surface of the cubes. This displacement was applied using a strain rate of 1 mm/s and then maintained for 600 s. The nodes on the bottom surface were fixed in the z direction and the pore pressure in the lateral walls of the cubes was set as zero. The total reaction force and the total stress were computed every time increment. Additionally, to determine the drained Young's modulus of the bone phase, linear elastic models of the bone phase were conducted. A total compressive strain equivalent of 1% was applied and the computed reaction force was fitted to the

Fig. 2 – Schematic of the stress relaxation experiment. 1% of compression (0.1 mm) was applied at a strain rate of 1 mm/s and maintained for 600 s.

force after relaxation of the sample specific mechanical experiments. In order to understand the contribution of the poroelastic and the viscoelastic components to the poro-viscoelastic model, the stress relaxation curves were simulated using: (1) only the poroelastic formulation applied to the bone and the marrow phases, and (2) only the viscoelastic formulation applied to the bone phase. The FE simulations of stress relaxation were conducted in reduced meshes of 4 mm side-length taken from the center of the cubes. This model size reduction (from 1.95  106 to 1.25  105 elements) was made in order to reduce the amount of memory and computation time required to simulate the relaxation curves. In order to correct the error caused by the size reduction and to be able to compare the FE results with the experimental curves, a geometry correction factor for each cube was determined as:   RFLE model;10 mm =100   ð19Þ CF ¼ RFLE model;4 mm =16 were CF is the correction factor, and RFLE model; 10 mm and RFLE model; 4 mm are the reaction forces computed using the linear elastic model described above with the 10 mm and the 4 mm side length cubes, respectively. Then, the total stress of the FE models conducted with 4 mm side length meshes was adjusted multiplying by CF. FE models were solved in Abaqus (Dassault Systèmes, USA) on a supercomputing system (four 64 bit 2.4 GHz AMD Istanbul processors; 15 GB of memory per 4 processor node; average of 2 h per simulation; Westgrid, Canada). A summary of all simulations conducted is shown in Table 2.

3.

Results

3.1.

Sample preparation and lCT imaging

No significant differences between the mCT derived morphometric variables measured before and after the stress relaxation experiment were detected (Table 3). Additionally, no trabecular fractures were visible when inspecting the registered images before and after the experiment (resolution of the images¼ 20 mm). These findings suggest damage to the cubes at the micro-scale was minimal. The average BV/TV of the sample was 29% with a standard deviation of 8%.

3.2.

Mechanical experiments

All the cubes exhibited a similar mechanical trend under relaxation. The range in reaction force values before relaxation (at 0.1 s) for the 25 cubes was 1.06  102-2.43  103 N, with an average of 9.28  102 N (SD75.11  102 N). The range after relaxation (at 600 s) was 5.38  101-1.34  103 N, with an average of 4.69  102 N (SD72.88  102 N). In general, 70% of the relaxation occurred in the first 20 s, 80% in the first 80 s, and 90% in the first 200 s. At 600 s almost 100% of the relaxation was achieved (Fig. 3a). The exponential decay curve adjusted from the average normalized reaction force (Fig. 3b) was: NFðtÞ ¼ 1 þ 0:60e  0:25t þ 0:31e  0:0054t

ð20Þ

5

journal of the mechanical behavior of biomedical materials 44 (2015) 1 –9

Table 1 – Material properties. Property

Young's modulus (MPa) Poisson's ratio Porosity Permeability (m/s) Grain bulk modulus (MPa) Fluid bulk modulus (MPa)

Linear elastic

Viscoelastic

Poroelastic

Poro-viscoelastic

Bone phase

Bone phase

Bone phase

Marrow phase

Bone phase

Marrow phase

Determined from experiments 0.325

Determined from experiments 0.325

Determined from experiments 0.325b 0.04c 9.79  10  14b

2a

2a

0.167a 0.8a 9.79  10  11a

Determined from experiments 0.325b 0.04c 9.79  10  14b

0.167a 0.8a 9.79  10  11a

17660d

2300e

17660d

2300e

2300e

2300e

2300e

2300e

a

Lacroix and Prendergast (2002). Cowin (1999). c Schaffler and Burr (1988). d Smit et al. (2002) e Anderson (1967). b

Table 2 – Summary of FE models. Model

Phases

Formulation

Model size

Applied strain and boundary conditions

1

Bone phase Marrow phase Bone phase

Poro-viscoelastic Poroelastic Linear elastic

4 mm

1% of compression during 600 s, pressure zero in the lateral sides 1% of compression

Bone phase Bone phase Marrow phase

Viscoelastic Poroelastic Poroelastic

10 mm 4 mm 4 mm 4 mm

2 3 4

1% of compression during 600 s 1% of compression during 600 s, pressure zero in the lateral sides

Table 3 – Measured morphometric variables before and after the stress relaxation experiment. Morphometric Variable

Bone volume ratio Connectivity density Structure model index Trabecular number Trabecular thickness Trabecular separation

BV/TV [%] ConnD [1/mm3] SMI [1] TbN [1/mm] TbTh [mm] TbSp [mm]

Before the mechanical experiment

After the mechanical experiment

average

stand dev

average

stand dev

28.71 6.52 0.47 1.60 0.22 0.61

8.03 2.63 0.54 0.27 0.05 0.09

28.14 6.23 0.41 1.60 0.22 0.60

7.98 2.22 0.63 0.23 0.05 0.07

The ratio between the reaction force before and after relaxation, which may be related to the relative content of collagen and mineral phase, was linearly distributed (slope¼ 1.92, R2 ¼0.96) varying between 1.59 and 2.59 and with an average of 2.05 (SD70.24) (Fig. 4). A statistically significant correlation (po0.05) was found between this ratio and the trabecular separation (r¼ 0.41). Slightly lower correlations were found between the relaxation ratio and the BV/TV (r¼  0.39, p-value¼ 0.06), the trabecular number (r¼  0.39, pvalue¼ 0.06), and the structure model index (r¼0.36, pvalue¼ 0.09). The correlation between the ratio and the age was not statistically significant (r¼ 0.22) (Table 4).

3.3.

p-value

0.80 0.68 0.71 0.96 0.97 0.87

FE modeling

Fitting the linear elastic models to the total reaction force after relaxation, the drained Young's modulus was determined to be 3.50 GPa (R2 ¼0.41) (Fig. 5). Therefore, the drained bulk modulus, K ¼E/3(1 2ν), was 3.33 GPa and the drained shear modulus, G¼ E/2(1þν), was 1.32 GPa. From the exponential equation described above (Eq. (20)), the relative shear and bulk modulus coefficients, g and k, and the relaxation time constant, t, for each Prony series term were determined as: g1 ¼ k1 ¼ 3.15  10  1, t1 ¼ 4 s and g2 ¼k2 ¼ 1.66  10  1, t2 ¼ 1.85  102 s.

6

journal of the mechanical behavior of biomedical materials 44 (2015) 1 –9

Table 4 – Correlation between decay ratio and morphometric variables.

Fig. 3 – Reaction force measured during the stress relaxation experiments. Average, minimum and maximum smoothed (a) and normalized (b) forces. All the experiments had a similar trend, with 70% of the relaxation obtained in 20 s and almost 100% in 600 s. The average normalized force was used to adjust an exponential decay equation of two terms.

Variable

Correlation

p-value

Bone volume ratio Connectivity density Structure model index Trabecular number Trabecular thickness Trabecular separation Age

 0.39  0.32 0.36  0.39  0.30 0.41 0.22

0.06 0.13 0.09 0.06 0.15 0.04 0.30

Fig. 5 – A regression plot comparing the force predicted with the linear elastic models with the force experimentally measured after relaxation. These FE models were used to estimate the drained Young's modulus used for the poroviscoelastic and the poroelastic models (E¼ 3.501 GPa). The line indicates the unity showing a linear correlation (R2 ¼ 0.41). respectively), while the relaxation predicted with the poroelastic formulation was significantly lower (average relaxation ratio of 1.03) (Fig. 7).

4.

Fig. 4 – A regression plot comparing the force experimentally measured before and after relaxation. A linear correlation between these forces with a ratio of 1.92 (R2 ¼ 0.96) was obtained. The FE poro-viscoelastic models explained 49% of the variability of the experimentally measured stress before the relaxation (at 0.1 s) (Fig. 6a) and 33% of the variability of the stress after the relaxation (at 600 s) (Fig. 6b). The trend of the stress relaxation curves predicted with this FE models was consistent with the curves obtained experimentally (Fig. 6c and d) presenting an offset over time. At the end of the relaxation curves, this offset varied between 4% (Fig. 6c) and 92% (Fig. 6d). The differences in the relaxation curves predicted with the poro-viscoelastic and the viscoelastic formulations were almost negligible (average relaxation ratios of 1.88 and 1.91,

Discussion

We proposed a mCT-based FE model of trabecular bone that accounts for both the micro architecture of the tissue and the combined poroelastic and viscoelastic nature of the bone, which has been applied in bone tissue for the first time. Previous biphasic poro-viscoelastic models have been applied to cartilage (Mak, 1986a, 1986b; Setton et al., 1993) and hydrogels (Kalyanam et al., 2009; Noailly et al., 2008; Strange et al., 2013), but did not consider the sample specific micro architecture that we have from our bone tissue. With our models we were able to simulate the exponential decay trend observed in the stress relaxation experiments and we demonstrated that at the intertrabecular pore and the lacunarcanalicular pore scale the effect of the viscoelasticity was considerably higher than the effect of the poroelasticity on the mechanical response of the tissue. Throughout the stress relaxation experiment, all the specimens presented a relaxation curve that was well characterized by an exponential equation of two terms, with relaxation time constants of 4 and 185 s. These results were similar to the experimental results of Deligianni et al. (1994) and Quaglini et al. (2009) who also adjusted exponential decay curves of two terms. However, our time relaxation

journal of the mechanical behavior of biomedical materials 44 (2015) 1 –9

7

Fig. 6 – Total stress predicted with the poro-viscoelastic models compared with experimental stress. Regression plots before (a) and after (b) relaxation indicating the line of unity, and relaxation curves for specimen A (c) and specimen B (d). The poroviscoelastic models captured 49% of the variability of the experimental stress before relaxation and 33% after relaxation. The error ranged from 4% to 92%.

Fig. 7 – Representative relaxation curves predicted with poro-viscoelastic, poroelastic, and viscoelastic FE models. The relaxation curves obtained with the viscoelastic formulation were similar to the curves obtained with the poro-viscoelastic formulation while the relaxation predicted with the poroelastic formulation was significantly lower. The average ratios of relaxation predicted were 1.88, 1.91 and 1.03 for the poro-viscoelastic, viscoelastic, and poroelastic models, respectively. constants were higher than the constants reported by others. Deligianni et al. (1994) reported values of 0.1 and 100 s, and

Quaglini et al. (2009) determined that the first constant varied between 2.18 and 4.25 s according to the applied initial load and the second constant was 95 s. No statistically significant differences between the morphometric variables taken before and after the experiments were found, and no differences in the registered images were detected on visual inspection; therefore, we are confident that we were measuring the actual tissue relaxation and not appreciable damage accumulation. This aspect is critical since Nagaraja et al. (2007) demonstrated that strain relaxation experiments in bone can cause damage accumulation. The ratio between the stress before and after relaxation, F(0.1)/F(600), varied between 1.59 and 2.59. The difference in the relaxation between specimens may be related to differences in the relative content between the collagen and the mineral phase. We did not measure the collagen and the mineral content of the specimens, but it is established that the mineral content and the collagen cross-links vary with age, affecting mechanical properties of the tissue (Boskey and Coleman, 2010). In this study the correlation between the ratio of relaxation and the age of the donors was not statistically significant. However, this ratio presented a positive correlation with the trabecular separation. In other words, even though we did not find evidence that the viscoelastic response of the bone vary directly with age, we found that the viscoelasticity of the tissue decreased with variations in the micro architecture that are typical of age. The statistically significant correlations between the

8

journal of the mechanical behavior of biomedical materials 44 (2015) 1 –9

relaxation ratio and the trabecular separation suggest that the viscoelasticity of the tissue may be possible to predict from the micro architecture of the tissue. Further mCT based studies that include sample specific measures of the viscoelastic properties of trabecular bone at both the intertrabecular and the lacunar-canalicular pore scale are necessary to clarify the relationship between the tissue micro architecture and its viscoelasticity. Additionally a higher sample size is needed. The computational approach proposed in this study did not include subject specific material properties. The average relaxation curve from the experimental results was used to adjust the viscoelastic parameters for the models, which in addition to the poroelastic properties were kept constant for all the cubes. Therefore we believe that the predicted variability (49% before relaxation and 33% after relaxation) is explained by the micro architecture of the tissue. In order to improve the accuracy of the models additional studies of the subject specific material properties are needed. The objective of implementing the poro-viscoelastic formulation was to include both the effect of the poroelasticity (fluid flow-dependent mechanism) and the effect of viscoelasticity (fluid flow-independent mechanism). Even though the time dependency of the material properties of bone has been previously explained with the poroelastic theory, the ratio between the force before and after relaxation measured in our experiments (41.59) cannot be modeled using only the poroelastic formulation, as the maximum ratio that can be obtained analytically using a positive poisson ratio is 1.5 (Armstrong et al., 1984). Our computational results suggest that at the intertrabecular pore and the lacunar-canalicular pore scale the viscoelastic component had a higher contribution to the relaxation than the poroelastic component, and therefore the mechanical response of the tissue can be simulated using a viscoelastic model. The advantage of including the poroelasticity in the model is that it allows determination of the pore pressure and fluid flow, which will provide useful information for the study of bone mechanotransduction. Additionally, studies at the molecular scale, have suggested that viscoelasticity may be governed by either the rearrangement of collagen molecules and water (Shen et al., 2011) or the intercrystalline layers of water (Eberhardsteiner et al., 2014). Therefore the poroelasticity could become the dominant mechanism at the molecular scale. The mechanical response of a poroelastic material depends on the boundary conditions. Under unconfined compression, the fluid can flow across the boundary and consequently the pore pressure relaxes; however, under confined compression no relaxation due to the pore pressure takes place. In contrast, the response of a viscoelastic material is independent of the pore pressure boundary conditions. In our study, the unconfined boundary conditions had an effect on the stiffness measured before and during the relaxation; however, these boundary conditions do not contribute to the relative effect of the poroelastic and the viscoelastic components on the poro-viscoelastic response. In our model, the viscoelastic parameters were adjusted according to the experimental results and the parameters in shear and in bulk were assumed to be equal. This assumption reduced the number of material coefficients in the formulation and decreased the complexity of the problem. The poroelastic

properties were taken from Lacroix and Prendergast (2002) approach which has been successfully used to predict the differentiation of tissue in fracture healing, but some of the properties they use (e.g. the Young's modulus and the poisson ratio of marrow) have not been validated. The use of these properties is a limitation of our study; however, due to the low effect of the poroelasticity compared to the effect of viscoelasticity, the sensitivity of these poroelastic parameters on the reaction force measured in the relaxation was low. The reduction of the model size from 20 to 80 mm and from 10 to 4 mm was conducted in order to reduce the total number of elements and to solve the complete relaxation curves. The error caused by applying a size reduction was previously studied (Sandino and Boyd, 2013). We found that changing the resolution of the images causes an error in the total stress of about 1%, and reducing the size of the cubes causes an error of about 3%. Our study also presented some experimental limitations, for instance, the stress relaxation tests were conducted with a 10 kN load cell (accuracy of 0.5%); therefore the error, which can be of about 50 N, is close to the magnitude of the inferior range of the measured loads after relaxation (5.38  101– 1.34  103 N). Moreover, the machine compliance was not taken into account. Additionally, we did not measure the marrow and the fluid flow across the boundaries of the specimens or the pressure inside the pores. In spite of these computational and experimental limitations we were able to explain up to 49% of the variability of the mechanical behavior of trabecular tissue under stress relaxation using mCT images. This demonstrates that the micro architecture of the tissue contributes to almost half of the variability in the mechanical response of the tissue and exhibits the potential of using CT images in bone studies. The other half of variability may be explained by the subject specific material properties. Models including the micro architecture of the bone tissue as well as using appropriate constitutive equations are necessary to improve the understanding of the behavior of this tissue.

Acknowledgments This study was supported by the Natural Science and Engineering Research Council of Canada (NSERC) and Alberta Innovates Health Solutions (AIHS). Dr. Dennis Coombe, from the Computational Modeling Group at the University of Calgary, Dr. Jérôme Noailly and Dr. Andrea Malandrino from the institute for Bioengineering of Catalonia, are acknowledged for their valuable discussions about bone permeability and the poro-viscoelastic formulation.

r e f e r e n c e s

Abousleiman, Y., Cheng, A.H.D., Jiang, C., Roegiers, J.C., 1996. Poroviscoelastic analysis of borehole and cylinder problems. Acta Mech. 119, 199–219. Anderson, C.B., 1967. Mechanics of fluids. In: Baumeister, T. (Ed.), Marks’ Standard Handbook for Mechanical Engineers. McGraw-Hill, New York, pp. 3.48–3.76.

journal of the mechanical behavior of biomedical materials 44 (2015) 1 –9

Armstrong, C.G., Lai, W.M., Mow, V.C., 1984. An analysis of the unconfined compression of articular cartilage. J. Biomed. Eng. 106, 165–173. Bevill, G., Keaveny, T.M., 2009. Trabecular bone strength predictions using finite element analysis of micro-scale images at limited spatial resolution. Bone 44, 579–584. Biot, M.A., 1941. General theory of three-dimensional consolidation. J. Appl. Phys. 12, 155–164. Boskey, A.L., Coleman, R., 2010. Aging and bone. J. Dent. Res. 89, 1333–1348. Bowman, S.M., Keaveny, T.M., Gibson, L.J., Hayes, W.C., McMahon, T.A., 1994. Compressive creep behavior of bovine trabecular bone. J. Biomech. 27, 301–310. Brynk, T., Hellmich, C., Fritsch, A., Zysset, P.K., Eberhardsteiner, J., 2011. Experimental poromechanics of trabecular bone strength: role of Terzaghi’s effective stress and of tissue level stress fluctuations. J. Biomech. 44, 501–508. Carter, D.R., Hayes, W.C., 1977. The compressive behavior of bone as a two-phase porous structure. J. Bone Jt. Surg. 59A, 954–962. Chen, P.Y., McKittrick, J., 2011. Compressive mechanical properties of demineralized and deproteinized cancellous bone. J. Mech. Behav. Biomed. Mater. 4, 961–973.. Cowin, S.C., 1999. Bone poroelasticity. J. Biomech. 32, 217–238. Deligianni, D.D., Maris, A., Missirlis, Y.F., 1994. Stress relaxation behaviour of trabecular bone specimens. J. Biomech. 27, 1469–1474.. Eberhardsteiner, L., Hellmich, C., Scheiner, S., 2014. Layered water in crystal interfaces as source for bone viscoelasticity: arguments from a multiscale approach. Comput. Methods Biomech. Biomed. Eng. 17, 48–63. Fritsch, A., Hellmich, C., Dormieux, L., 2009. Ductile sliding between mineral crystals followed by rupture of collagen crosslinks: experimentally supported micromechanical explanation of bone strength. J. Theor. Biol. 21, 230–252. Fritton, S.P., Weinbaum, S., 2009. Fluid and solute transport in bone: flow-induced mechanotransduction. Annu. Rev. Fluid Mech. 41, 347–374. Guille´n, T., Zhang, Q.H., Tozzi, G., Ohrndorf, A., Christ, H.J., Tong, J., 2011. Compressive behaviour of bovine cancellous bone and bone analogous materials, microCT characterisation and FE analysis. J. Mech. Behav. Biomed. Mater. 4, 1452–1461. Gurkan, U.A., Akkus, O., 2008. The mechanical environment of bone marrow: a review. Ann. Biomed. Eng. 36, 1978–1991. Hellmich, C., Celundova, D., Ulm, F.J., 2009. Multiporoelasticity of hierarchically structured materials: micromechanical foundations and application to bone. J. Eng. Mech. (ASCE) 135, 382–394. Hellmich, C., Ulm, F.J., 2005. Drained and undrained poroelastic properties of healthy and pathological bone: a poromicromechanical investigation. Transp. Porous Media 58, 243–268. Hong, J.H., Song, S.H., 1998. Poroelastic bahavior of trabecular bone: the effect of strain rate. KSME Int. J. 12, 421–428. Hong, J.H., Mun, M.S., Lim, T.H., 2001. Strain rate dependent poroelastic behavior of bovine vertebral trabecular bone. KSME Int. J. 15, 1032–1040. Kalyanam, S., Yapp, R.D., Isana, M.F., 2009. Poro-viscoelastic behavior of gelatin hydrogels under compression-Implications for bioelasticity imaging. J. Biomech. Eng. 131, 081005. Lacroix, D., Prendergast, P.J., 2002. A mechano-regulation model for tissue differentiation during fracture healing: analysis of gap size and loading. J. Biomech. 35, 1163–1171. Lim, T.H., Hong, J.H., 1998. Poroelastic model of trabecular bone in uniaxial strain conditions. J. Musculoskelet. Res. 2, 167–180. Lim, T.H., Hong, J.H., 2000. Poroelastic properties of bovine vertebral trabecular bone. J. Orthop. Res. 18, 671–677. MacNeil, J.A., Boyd, S.K., 2008. Bone strength at the distal radius can be estimated from high resolution peripheral quantitative

9

computed tomography and the finite element method. Bone 42, 1203–1213. Mak, A.F., 1986a. The apparent viscoelastic behavior of articular cartilage - The contributions from the intrinsic matrix viscoelasticity and interstitial fluid flows. J. Biomech. Eng. 108, 123–130. Mak, A.F., 1986b. Unconfined compression of hydrated viscoelastic tissues: a biphasic poroviscoelastic analysis. Biorheology 23, 371–383. Moore, T.L.A., O’Brien, F.J., Gibson, L.J., 2004. Creep does not contribute to fatigue in bovine trabecular bone. J. Biomech. Eng. 126, 321–329. Mu¨ller, R., Hildebrand, T., Ruegsegger, P., 1994. Non-invasive bone biopsy: a new method to analyze and display the threedimensional structure of trabecular bone. Phys. Med. Biol. 39, 145–164. Nagaraja, S., Ball, M.D., Guldberg, R.E., 2007. Time-dependent damage accumulation under stress relaxation testing of bovine trabecular bone. Int. J. Fatigue 29, 1034–1038. Noailly, J., Van Oosterwyck, H., Wilson, W., Quinn, T.M., Ito, K., 2008. A poroviscoelastic description of fibrin gels. J. Biomech. 41, 3265–3269. Ochoa, J.A., Heck, D.A., Hillberry, B.M., 1991a. The effect of intertrabecular fluid on femoral head mechanics. J. Rheumatol. 18, 580–584. Ochoa, J.A., Sanders, A.P., Heck, D.A., Hilberry, B.M., 1991b. Stiffening of the femoral head due to intertrabecular fluid and intraosseous pressure. J. Biomech. Eng. 113, 259–261. Quaglini, V., La Russa, V., Corneo, S., 2009. Nonlinear stress relaxation of trabecular bone. Mech. Res. Commun. 36, 275–283. Rice, J.R., Cleary, M.P., 1976. Some basic stress-diffusion solutions for fluid saturated elastic porous media with compressible constituents. Rev. Geophys. Space Phys. 14, 227–241. Sandino, C., Boyd, S.K., 2013. Trabecular bone poroelasticity for microCT-based FE models. In: Wittek, A., Miller, K., Nielsen, P. M.F. (Eds.), Computational Biomechanics for Medicine – Models, Algorithms and Implementation. Springer, New York, pp. 145–155. Sanyal, A., Gupta, A., Bayraktar, H.H., Kwon, R.Y., Keaveny, T.M., 2012. Shear strength behavior of human trabecular bone. J. Biomech. 45, 2513–2519. Sasaki, N., Nakayama, Y., Yoshikawa, M., Enyo, A., 1993. Stress relaxation function of bone and bone collagen. J. Biomech. 26, 1369–1376. Schaffler, M.B., Burr, D.B., 1988. Stiffness of compact bone: effects of porosity and density. J. Biomech. 21, 13–16. Setton, L.A., Zhu, W., Mow, V.C., 1993. The biphasic poroviscoelastic behavior of articular cartilage: role of the surface zone in governing the compressive behavior. J. Biomech. 26, 581–592. Shen, Z.L., Kahn, H., Ballarini, R., Eppell, S.J., 2011. Viscoelastic properties of isolated collagen fibrils. Biophys. J. 100, 3008–3015. Smit, T.H., Huyghe, J.M., Cowin, S.C., 2002. Estimation of the poroelastic parameters of cortical bone. J. Biomech. 35, 829–835. Strange, D.G.T., Fletcher, T.L., Tonsomboon, K., Brawn, H., Zhao, X., Oyen, M.L., 2013. Separating poroviscoelastic deformation mechanisms in hydrogels. Appl. Phys. Lett. 102, 031913. van Rietbergen, B., Majumdar, S., Pistoia, W., Newitt, D.C., Kothati, M., Laib, A., Ru¨egsegger, P., 1998. Assessment of cancellous bone mechanical properties from micro-FE models based on micro-CT, pQCT and MR images. Technol. Health Care 6, 413–420. Viguet-Carrin, S., Garnero, P., Delmas, P.D., 2006. The role of collagen in bone strength. Osteoporos. Int. 17, 319–336.