Jarunan Panyasantisuk Institute for Surgical Technology and Biomechanics, University of Bern, Stauffacherstr. 78, Bern CH-3014, Switzerland e-mail: [email protected]
Dieter H. Pahr Institute for Lightweight Design and Structural Biomechanics, Vienna University of Technology, Gusshausstr. 27-29/317, Vienna A -1040, Austria e-mail: [email protected]
Thomas Gross Institute for Lightweight Design and Structural Biomechanics, Vienna University of Technology, Gusshausstr. 27-29/317, Vienna A-1040, Austria e-mail: [email protected]
Philippe K. Zysset Institute for Surgical Technology and Biomechanics, University of Bern, Stauffacherstr. 78, Bern CH-3014, Switzerland e-mail: [email protected]
1
Comparison of Mixed and Kinematic Uniform Boundary Conditions in Homogenized Elasticity of Femoral Trabecular Bone Using Microfinite Element Analyses Mechanical properties o f human trabecular bone play an important role in age-related bone fragility and implant stability. Microfinite element (pFE) analysis allows computing the apparent elastic properties o f trabecular bone fo r use in homogenized FE (hFE) anal ysis, but the results depend unfortunately on the type o f applied boundary conditions (BCs). In this study, 167 human femoral trabecular cubic regions with a side length o f 5.3 mm were extracted from three proximal femora and analyzed using pFE analysis to compare systematically their stiffness with kinematic uniform BCs (KUBCs) and periodicity-compatible mixed uniform BCs (PMUBCs). The obtained elastic constants were then used in the volume fraction and fabric-based orthotropic Zysset-Curnier model to identify their respective model parameters. As expected, PMUBCs lead to more com pliant apparent elastic properties than KUBCs, especially in shear. The differences in stiffness decreased with bone volume fraction and mean intercept length (MIL). Unlike KUBCs, PMUBCs were sensitive to heterogeneity o f the biopsies. The Zysset-Curnier model fitted the apparent elastic constants successfully in both cases with adjusted coeffi cients o f determination (r2^ ) o f 0.986 fo r KUBCs and 0.975 fo r PMUBCs. The proper use o f these BCs fo r hFE analysis o f whole bones will need to be investigated in future work. [DOI: 10.1115/1.4028968]
Introduction
Mechanical behavior of trabecular bone plays an essential role in osteoporosis-related fractures and in stability of orthopedic or dental implants. FE analysis is becoming a standard technique to predict the biomechanical behavior of whole bones and boneimplant systems [1-3]. For trabecular biopsies, elastic and yield properties can be computed using /(FE analysis. However, highresolution computed tomography images of the trabecular micro structure and the required computing resources for /(FE may not always be available. On the other hand, CT-based hFE analysis becomes an accepted standard to predict stiffness and strength of bones at risk of fracture [4], The latter method relies on accurate homogenized material properties for trabecular bone that proved to be difficult to assess experimentally due to various artifacts [5,6]. Accordingly, several investigators applied homogenization analysis to compute the full elasticity tensors of trabecular bone with /(FE models subjected to different types of BCs [7-9]. Unfortunately, the used BCs such as kinematic uniform or peri odic BCs do not seem to deliver apparent elastic properties that are compatible with experimental validation tests [10]. The ques tion of proper apparent elastic properties of human trabecular bone for hFE analysis requires therefore further research. Homogenization [11,12] is a powerful theory to quantify the apparent elastic properties of trabecular bone. The homogeniza tion approach decomposes the total strain and stress in an average (apparent) plus a fluctuating component [13] and delivers apparent
Manuscript received May 30, 2014; final manuscript received October 23, 2014; accepted manuscript posted November 3, 2014; published online December 10, 2014. Assoc. Editor; Blaine Christiansen.
Journal of Biomechanical Engineering
elastic properties that relate the apparent strain and stress with Hooke’s linear constitutive law. These apparent properties can be computed by applying BCs on volume elements that should be representative but that are not necessarily periodic [14]. As an important energy equivalence principle, the Hill condition states that the average of the internal product of the stress and strain ten sors at the microlevel equals the internal product of their averages at the macrolevel [11,15]. For nonperiodic random media, the three BCs that fulfill Hill’s condition are as follows [16,17]: (1) kinematic uniform boundary conditions (KUBCs) in which the points of the boundaries are constrained to displace uniformly, (2) static uniform boundary conditions (SUBCs) in which the points are constrained with uniform traction, and (3) mixed uniform boundary conditions (MUBCs) which are com binations of uniform displacement and traction constraints. In general, for volume elements that are not representative, these BCs provide distinct apparent elastic properties, but for represen tative volume elements (RVEs), the results coalesce into effective properties. In particular, the exact effective properties can be obtained by applying periodic BCs [7] on periodic volume elements [18,19]. Hazanov and Huet [20] showed that KUBCs and SUBCs repre sent, respectively, the upper and lower bounds of all apparent stiffness tensors and, in particular, the ones obtained by MUBCs and periodic BCs. Nonetheless, a stiffness tensor based on MUBCs is generally not equal to the effective stiffness tensor. For trabecular bone, a characteristic length of a few millimeters was claimed to fulfill the continuum assumption for a RVE [21], but the microstructure is obviously not periodic. Apparent properties of trabecular bone were originally computed by pFE with KUBCs but
Copyright © 2015 by ASME
JANUARY 2015 , Vol. 137 / 011002-1
not acknowledged as an upper bound of the sought in situ elastic properties [8]. Applying periodic BCs on trabecular bone cubes may connect bone tissue to marrow on opposed surfaces and lead to biased elastic properties of the real connected bone [9], Pahr and Zysset [15] showed that a selected set of MUBCs, namely periodicity-compatible MUBCs (PMUBCs), are equiva lent to periodic BCs when the considered cubic volume element exhibits orthotropic symmetry and the planes of symmetry coin cide with the faces of the cube. Their investigation showed that cubic trabecular bone biopsies with a side length of 5 mm were an excellent choice, since the apparent stiffness resulting from PMUBCs matched with that from periodic BCs. However, only six samples were investigated in this study. Homogenized elastic properties of trabecular bone have a close relationship to the underlying morphology [22]. Zysset [23] reviewed the theoretical models relating volume fraction and fabric to trabecular bone apparent stiffness. These models include the orthotropic Zysset-Curnier model, which predicts a positive-definite apparent stiffness with a minimal set of five constants [24], In a recent study, computational homogenization was con ducted using /iFE with KUBCs on 701 trabecular bone samples from three clinically relevant anatomical locations, the proximal femur, thoracolumbar vertebra, and distal radius [25]. Prediction of apparent stiffness with the Zysset-Curnier model was outstand ing and highly similar for the three locations, but it remains unclear how the results are affected by BCs. Accordingly, in the present study, linear /iFE homogenization analysis of femoral trabecular bone was performed to compute apparent stiffness by applying KUBCs and PMUBCs on a same large set of trabecular bone biopsies. The objectives were (1) to compare apparent stiffness obtained from KUBCs versus PMUBCs and (2) to determine the model parameters of the Zysset-Curnier model for PMUBCs.
2
Theoretical Model
Briefly, the Zysset-Curnier model [24] relates bone volume fraction (BV/TV) and bone fabric or texture to the apparent stiff ness. The fabric tensor M is a positive-definite second-order tensor with eigenvalues m, and normalized eigenvectors m, [22,26]. M is made independent of bone volume fraction with the normalization tr(M) = 3. 3
M =
(1) i=i
where M, = m, ® m, are the dyadic products of the eigenvectors. The apparent stiffness tensor S is calculated using bone volume fraction, complete fabric information, elasticity parameters (Ao, Aq, and fi0), and the exponents (k and /) 3
S (p, M) =
(A0 + 2pi))pkm fM i M,1=1 3
+ ^ 2'0pkmlim'jMj ® My '■>/=1
¥J 3
+^2
2 p 0p km \m l]M i ® M j .
(2)
'7=1
¥i
where p is the compact mathematical expression of bone volume over total volume (BV/TV).
3
Materials and Methods
Bone sections from the head of three human proximal femora (3 female, age 62-75 years, mean age 66 ± 8 years) were scanned 0 1 1 0 0 2 -2 / Vol. 137, JANUARY 2015
E ast
N o rth
n
Top
ex u i # 0 tx = 0 tx = 0 e2 t2 = 0 u 2 = 0 t2 = o - - se3 t3 = 0 t3 = 0 u3= 0
o
n
o
- 3u - J ; : ex l *- [ ; 1 ! d o ci o u ■ Bottom; Top / , ■ ............ .............‘ ________
E a st
N o rth
Top
e i ti = 0 ti = 0 Ux 0 e 2 u 2 = 0 u 2 = 0 U2= 0 e3 U3i* o t3 = o f c - °
|
/
*"3*/
West |/ -
'T
!
-
J
/!
!
/Bottom Fig. 1 Top: a /iFE mesh converted directly from a /iCT image with a resolution of 3 6 /im; middle: the uni-axial loading case in 1-direction applying PMUBCs; bottom: the shear loading case in 1- and 3-direction applying PMUBCs; and e, are displacem ents, tractions, and unit directions, respectively. N orth -S ou th plane is normal to W es t-E a s t and B ottom -Top planes.
with microcomputed tomography (/iCT 40, SCANCO Medical AG, Briittisellen, Switzerland) with a resolution of 18 pm. The cortical shell was removed from the images by segmentation and a total of 167 intact trabecular regions of interest with a side length of 9.2 mm were extracted. PMUBCs require that the volumes of interest are aligned with the orthotropic material symmetry axes. Therefore, trabecular regions of interest were rotated iteratively so that the principal axes of the fabric tensor of a cubic subregion with 5.3-mm side length were aligned with its edges. The MIL method [26,27] was used to evaluate the fabric tensor. The tolerance error of the prin cipal axes to the cube edges was 1 deg. A linear interpolation scheme was applied in the image rotation. After rotation, the cubic subregions with a side length of 5.3 mm were cropped, coarsened to a resolution of 36 pm (see Fig. 1 top), segmented using a single-level threshold [28], and subsequently all of the uncon nected bone regions were removed. To generate /iFE models, image voxels were converted directly to mesh elements. Image processing and FE model generation were done by MEDTOOL.r Linear homogenization analysis of these trabecular bone cubes was performed using /iFE analysis with KUBCs and PMUBCs. A Young’s modulus of lOGPa and a Poisson’s ratio of 0.3 were assigned to each element [29]. Six independent load cases were applied on the cubes, including three uni-axial and three shear cases. Uni-axial and shear loading cases of PMUBCs are illus trated in Fig. 1 middle and bottom, respectively. From homoge nized /iFE analysis, stress and strain fields were obtained, from which the overall elastic stiffness tensor S fe^ was evaluated using stress and strain volume averages [15]. Subsequently, an orthotropic representation SFEortlo of the stiffness tensor was deter mined by neglecting the nonorthotropic tensor components in the coordinate system of the fabric tensor [15,25]. The norm error cor responding to this assumption can be calculated with 1www.dr-pahr.at
Transactions of the ASME
Fig. 2 Exam ples of highly heterogeneous bone biopsies which were not included in the filtered data set: BV/TV = 11.21% (left) and 24.07% (right)
NE aniso—ortho
®
The orthotropic qFE stiffness components were then fitted to the Z ysset-C um ier model using multiple linear regressions, which yield the m odel parameters. In this study, multiple linear regres sions were solved by a linear solver routine (numpy.linalg.solve(), p y t h o n 2.6.6). In general, the (global) linear system includes 12 nonzero components of the orthotropic stiffness tensor of each sample, w hich were aligned in a vector y y = Xc + e
(4)
where X is a 12 x p matrix containing p and m,-, c is a vector containing p constants of the model, and e is the vector of the residuals. The global linear system is shown in Appendix in matrix form. In an attempt to generalize the fabric-elasticity model without losing positive-definiteness, the three shear com po nents were split from the global linear system
yx = X aca +
(5)
y„ = X ,,^ + e„
(6)
where the subscription A denotes the linear system of normal components o f a stiffness tensor and the subscription p denotes the linear system o f shear components. The split linear systems require two additional exponents k and / as seen in Appendix. A fter computing the model param eters, model stiffness could be predicted by the Z ysset-C um ier model using the obtained model param eters, bone volume fraction, and the complete fabric information. The difference between the model stiffness S m0dei and orthotropic pFE stiffness SpEortho is expressed by the model norm error NEortho_model NE ortho—model
II^FEotat,
§ m o d e l||
ll§FEortJ |
calculated for each subcube. Mean and standard deviation of Psubcube were computed and the coefficient o f variation o f p Subcube calculated with S t d ( /L u b c u h e
(1) The complete data set: a data set w hich includes all the trabecular bone biopsies (2) The filtered data set: a data set which excludes highly heter ogeneous bone biopsies and includes only the trabecular bone cubes with CV below the computed threshold For the complete data set, the m odel parameters were calculated using the global linear system. For the filtered data set, a set of model parameters were calculated using the global linear system and another set using the split linear systems. Linear regressions in log scale were performed between the model stiffness S modei and the orthotropic pFE stiffness S FE„„h„• KUBC-based stiffness of the current methodology was then compared to that o f Gross et al. in linear regressions using the model parameters from Ref. [25]. For the filtered data set, pFE stiffness tensor components based on PMUBCs and KUBCs were compared in a one-to-one relation ship. The norm error N E KUBC- p m u b c describes the difference between pFE stiffness of the two BCs
(7)
||sL UBCs ™rC.0rth0 UBC|||| I rC .ortho_____
The total norm error NE^so. moliei describes the difference of the m odel stiffness S m0dei with respect to the full anisotropic stiffness tensor § fe^ „ IIS f
^>model 11
— ‘
(8)
IIS f b _ . I I
Unlike KUBCs, PMUBCs showed sensitivity to heterogeneity, which required further investigation. Among other tested m eth ods, the coefficient o f variation (CV) o f bone volume fraction was retained to assess heterogeneity. Each bone cube was divided into eight identical subcubes and bone volume fraction (psubcube) was
Journal of Biomechanical Engineering
(9)
The bone cubes were sorted by CV in descending order. Trabecu lar bone cubes with the highest CV were excluded one-by-one. At each step, coefficients of determination (r2) was obtained from the linear regression o f model stiffness to pFE stiffness based on PMUBCs. In the beginning of the cube exclusion, r2 rose steeply. Then, it reached a peak value after which it stabilized and started to fall. The coefficient of variation at the peak value was taken as the threshold to exclude highly heterogeneous bone samples (illustrated in Fig. 2). Hence, two data sets were analyzed in this study.
N E k u b c - pm u b c
NE,aniso—model
)
m ean(psubcube)
PMUBC
||S FEortho
(10)
Table 1 Mean values (±s ta n d a rd deviation) of bone volum e fraction (BV/TV) and degree of anisotropy (DA), and num ber of fem oral trabecular bone samples (n) of the exam ined data sets
Gross et al. [25] Complete data set Filtered data set
BV/TV
DA
n
0.19 (±0.10) 0.25 (±0.08) 0.27 (±0.08)
1.67 (±0.34) 1.54 (±0.20) 1.57 (±0.18)
264 167 126
JANUARY 2 0 1 5 , Vol. 137 / 011002-3
• ■
g -3 .0
Gross etal. (2012) Complete data set Filtered data set
a
O «
‘c 2.5
while the lowest difference is in the diagonal terms of the normal components Xu. Moreover, the rela tive difference decreases with MIL and is therefore minimal along the major fabric direction A33 (see Fig. 5). The exponents of volume fraction k range from 1.87 to 2.01, which agree well with the power coefficient of 2 for the apparent density resulting from a statistical review by Rice et al. [30]. The exponents of fabric l range from 1.00 to 1.68, which are bounded by those from Zysset [23]. Interestingly, the exponents k and l based on PMUBCs are higher than the ones based on KUBCs. Using the split linear systems, the exponents k and / are around 1.3 times higher, except the exponent / of the shear components which is around 1.6 times higher. This indicates that, based on PMUBCs, bone volume fraction and fabric information play a stronger role in the apparent properties, in particular, fabric eigen values have a strong effect on the shear components. Linear regressions of the Zysset-Cumier model to /dFE stiffness yield r^dj from 0.945 to 0.975 for PMUBCs, and from 0.984 to 0.986 for KUBCs. Exclusion of highly heterogeneous bone samples, which violate the notion of RVE, improved significantly the norm errors of both KUBCs and PMUBCs. Splitting linear systems between the normal and shear components improved sig nificantly the norm errors only for PMUBCs, whereas for KUBCs not significant benefit could be established. There were some limitations in this study. First, the PMUBCs are based on orthotropic symmetry of the RVE but trabecular bone cubes are only approximately orthotropic as witnessed by the approximately 5% error on the tensor norm. However, PMUBCs do not require periodicity, which was shown to induce artificial connections between bone and marrow, especially at low volume fraction [9]. Second, trabecular bone images had to be rotated because PMUBCs require orthotropic microstructures aligned with the planes of orthotropic symmetry. The morphology of the rotated bone images deviated slightly from that of the original images. Bone volume fraction and degree of anisotropy were on average 0 1 1 0 0 2 -6
/ Vol. 137, JANUARY 2015
higher after rotation with the mean deviation of 3.46% (±4.29) and 1.09% (±6.57), respectively. Among other morphology infor mation, the highest mean deviation was found in connectivity with an average decrease of 17.37% (±7.38). Fortunately, this unexpected finding did not affect the morphology-elasticity relationships. Third, PMUBCs are sensitive to heterogeneity, which lead to the analysis of two data sets. The filtered data set that excludes highly heterogeneous bone cubes also excludes outliers of stiff ness tensor components (see Fig. 6). The Zysset-Cumier model tends to overestimate elastic properties for highly heterogeneous samples, especially for the ones with low bone volume fraction. In future studies, these two BCs will need to be compared in FE analysis of whole bones. KUBCs represent an upper bound of the apparent stiffness [15,20], They may be appropriate for a trabecu lar bone volume element with a rigid neighborhood, e.g., regions surrounded by cortex or high BV/TV. On the other hand, if a vol ume element is surrounded by similar morphology, PMUBCs may yield more appropriate apparent elastic properties. A solution will need to be found for homogenization of trabecular bone material properties in the presence of marked heterogeneity, i.e., of high gradients of volume fraction or fabric. In particular, the influence of BCs on the apparent yield properties of trabecular bone remain to be explored. Finally, proper homogenization of trabecular bone is expected to improve the accuracy of both linear and nonlinear hFE analysis of bone or bone-implant systems.
Appendix: Multiple Regression Models The global linear system has the form of
M A 22 )
0 0 ln(p) 1 0 0 ln(p)
ln(A33)
1
0 0
ln(p)
ln(mf)
ln(A32)
0
1
0
ln(p)
ln(m3m2)
ln (2i3)
0
1
0 ln(p)
ln(mi/w3)
ln(A2i)
0
1
0 ln(p)
ln(m2m i)
K /% )
0
0
1
ln(p)
ln(m2m3)
M P i 3)
0 0
1
ln(p)
ln(m im 3)
ln(Pi2)
0
0
1
ln(p)
ln(/wi/??2)
ln(/l23)
0
1
0 ln(p)
ln(/??2 /H3)
ln(A3i)
0
1
0 ln(p)
ln(m3/wi)
1
0 ln(p)
ln(/Wi7M2) /
/ In (^ n ) \
Vln(Ai2) )
ln(m f) \
(1
u
In (m\)
(ei \
Dieter H. Pahr Institute for Lightweight Design and Structural Biomechanics, Vienna University of Technology, Gusshausstr. 27-29/317, Vienna A -1040, Austria e-mail: [email protected]
Thomas Gross Institute for Lightweight Design and Structural Biomechanics, Vienna University of Technology, Gusshausstr. 27-29/317, Vienna A-1040, Austria e-mail: [email protected]
Philippe K. Zysset Institute for Surgical Technology and Biomechanics, University of Bern, Stauffacherstr. 78, Bern CH-3014, Switzerland e-mail: [email protected]
1
Comparison of Mixed and Kinematic Uniform Boundary Conditions in Homogenized Elasticity of Femoral Trabecular Bone Using Microfinite Element Analyses Mechanical properties o f human trabecular bone play an important role in age-related bone fragility and implant stability. Microfinite element (pFE) analysis allows computing the apparent elastic properties o f trabecular bone fo r use in homogenized FE (hFE) anal ysis, but the results depend unfortunately on the type o f applied boundary conditions (BCs). In this study, 167 human femoral trabecular cubic regions with a side length o f 5.3 mm were extracted from three proximal femora and analyzed using pFE analysis to compare systematically their stiffness with kinematic uniform BCs (KUBCs) and periodicity-compatible mixed uniform BCs (PMUBCs). The obtained elastic constants were then used in the volume fraction and fabric-based orthotropic Zysset-Curnier model to identify their respective model parameters. As expected, PMUBCs lead to more com pliant apparent elastic properties than KUBCs, especially in shear. The differences in stiffness decreased with bone volume fraction and mean intercept length (MIL). Unlike KUBCs, PMUBCs were sensitive to heterogeneity o f the biopsies. The Zysset-Curnier model fitted the apparent elastic constants successfully in both cases with adjusted coeffi cients o f determination (r2^ ) o f 0.986 fo r KUBCs and 0.975 fo r PMUBCs. The proper use o f these BCs fo r hFE analysis o f whole bones will need to be investigated in future work. [DOI: 10.1115/1.4028968]
Introduction
Mechanical behavior of trabecular bone plays an essential role in osteoporosis-related fractures and in stability of orthopedic or dental implants. FE analysis is becoming a standard technique to predict the biomechanical behavior of whole bones and boneimplant systems [1-3]. For trabecular biopsies, elastic and yield properties can be computed using /(FE analysis. However, highresolution computed tomography images of the trabecular micro structure and the required computing resources for /(FE may not always be available. On the other hand, CT-based hFE analysis becomes an accepted standard to predict stiffness and strength of bones at risk of fracture [4], The latter method relies on accurate homogenized material properties for trabecular bone that proved to be difficult to assess experimentally due to various artifacts [5,6]. Accordingly, several investigators applied homogenization analysis to compute the full elasticity tensors of trabecular bone with /(FE models subjected to different types of BCs [7-9]. Unfortunately, the used BCs such as kinematic uniform or peri odic BCs do not seem to deliver apparent elastic properties that are compatible with experimental validation tests [10]. The ques tion of proper apparent elastic properties of human trabecular bone for hFE analysis requires therefore further research. Homogenization [11,12] is a powerful theory to quantify the apparent elastic properties of trabecular bone. The homogeniza tion approach decomposes the total strain and stress in an average (apparent) plus a fluctuating component [13] and delivers apparent
Manuscript received May 30, 2014; final manuscript received October 23, 2014; accepted manuscript posted November 3, 2014; published online December 10, 2014. Assoc. Editor; Blaine Christiansen.
Journal of Biomechanical Engineering
elastic properties that relate the apparent strain and stress with Hooke’s linear constitutive law. These apparent properties can be computed by applying BCs on volume elements that should be representative but that are not necessarily periodic [14]. As an important energy equivalence principle, the Hill condition states that the average of the internal product of the stress and strain ten sors at the microlevel equals the internal product of their averages at the macrolevel [11,15]. For nonperiodic random media, the three BCs that fulfill Hill’s condition are as follows [16,17]: (1) kinematic uniform boundary conditions (KUBCs) in which the points of the boundaries are constrained to displace uniformly, (2) static uniform boundary conditions (SUBCs) in which the points are constrained with uniform traction, and (3) mixed uniform boundary conditions (MUBCs) which are com binations of uniform displacement and traction constraints. In general, for volume elements that are not representative, these BCs provide distinct apparent elastic properties, but for represen tative volume elements (RVEs), the results coalesce into effective properties. In particular, the exact effective properties can be obtained by applying periodic BCs [7] on periodic volume elements [18,19]. Hazanov and Huet [20] showed that KUBCs and SUBCs repre sent, respectively, the upper and lower bounds of all apparent stiffness tensors and, in particular, the ones obtained by MUBCs and periodic BCs. Nonetheless, a stiffness tensor based on MUBCs is generally not equal to the effective stiffness tensor. For trabecular bone, a characteristic length of a few millimeters was claimed to fulfill the continuum assumption for a RVE [21], but the microstructure is obviously not periodic. Apparent properties of trabecular bone were originally computed by pFE with KUBCs but
Copyright © 2015 by ASME
JANUARY 2015 , Vol. 137 / 011002-1
not acknowledged as an upper bound of the sought in situ elastic properties [8]. Applying periodic BCs on trabecular bone cubes may connect bone tissue to marrow on opposed surfaces and lead to biased elastic properties of the real connected bone [9], Pahr and Zysset [15] showed that a selected set of MUBCs, namely periodicity-compatible MUBCs (PMUBCs), are equiva lent to periodic BCs when the considered cubic volume element exhibits orthotropic symmetry and the planes of symmetry coin cide with the faces of the cube. Their investigation showed that cubic trabecular bone biopsies with a side length of 5 mm were an excellent choice, since the apparent stiffness resulting from PMUBCs matched with that from periodic BCs. However, only six samples were investigated in this study. Homogenized elastic properties of trabecular bone have a close relationship to the underlying morphology [22]. Zysset [23] reviewed the theoretical models relating volume fraction and fabric to trabecular bone apparent stiffness. These models include the orthotropic Zysset-Curnier model, which predicts a positive-definite apparent stiffness with a minimal set of five constants [24], In a recent study, computational homogenization was con ducted using /iFE with KUBCs on 701 trabecular bone samples from three clinically relevant anatomical locations, the proximal femur, thoracolumbar vertebra, and distal radius [25]. Prediction of apparent stiffness with the Zysset-Curnier model was outstand ing and highly similar for the three locations, but it remains unclear how the results are affected by BCs. Accordingly, in the present study, linear /iFE homogenization analysis of femoral trabecular bone was performed to compute apparent stiffness by applying KUBCs and PMUBCs on a same large set of trabecular bone biopsies. The objectives were (1) to compare apparent stiffness obtained from KUBCs versus PMUBCs and (2) to determine the model parameters of the Zysset-Curnier model for PMUBCs.
2
Theoretical Model
Briefly, the Zysset-Curnier model [24] relates bone volume fraction (BV/TV) and bone fabric or texture to the apparent stiff ness. The fabric tensor M is a positive-definite second-order tensor with eigenvalues m, and normalized eigenvectors m, [22,26]. M is made independent of bone volume fraction with the normalization tr(M) = 3. 3
M =
(1) i=i
where M, = m, ® m, are the dyadic products of the eigenvectors. The apparent stiffness tensor S is calculated using bone volume fraction, complete fabric information, elasticity parameters (Ao, Aq, and fi0), and the exponents (k and /) 3
S (p, M) =
(A0 + 2pi))pkm fM i M,1=1 3
+ ^ 2'0pkmlim'jMj ® My '■>/=1
¥J 3
+^2
2 p 0p km \m l]M i ® M j .
(2)
'7=1
¥i
where p is the compact mathematical expression of bone volume over total volume (BV/TV).
3
Materials and Methods
Bone sections from the head of three human proximal femora (3 female, age 62-75 years, mean age 66 ± 8 years) were scanned 0 1 1 0 0 2 -2 / Vol. 137, JANUARY 2015
E ast
N o rth
n
Top
ex u i # 0 tx = 0 tx = 0 e2 t2 = 0 u 2 = 0 t2 = o - - se3 t3 = 0 t3 = 0 u3= 0
o
n
o
- 3u - J ; : ex l *- [ ; 1 ! d o ci o u ■ Bottom; Top / , ■ ............ .............‘ ________
E a st
N o rth
Top
e i ti = 0 ti = 0 Ux 0 e 2 u 2 = 0 u 2 = 0 U2= 0 e3 U3i* o t3 = o f c - °
|
/
*"3*/
West |/ -
'T
!
-
J
/!
!
/Bottom Fig. 1 Top: a /iFE mesh converted directly from a /iCT image with a resolution of 3 6 /im; middle: the uni-axial loading case in 1-direction applying PMUBCs; bottom: the shear loading case in 1- and 3-direction applying PMUBCs; and e, are displacem ents, tractions, and unit directions, respectively. N orth -S ou th plane is normal to W es t-E a s t and B ottom -Top planes.
with microcomputed tomography (/iCT 40, SCANCO Medical AG, Briittisellen, Switzerland) with a resolution of 18 pm. The cortical shell was removed from the images by segmentation and a total of 167 intact trabecular regions of interest with a side length of 9.2 mm were extracted. PMUBCs require that the volumes of interest are aligned with the orthotropic material symmetry axes. Therefore, trabecular regions of interest were rotated iteratively so that the principal axes of the fabric tensor of a cubic subregion with 5.3-mm side length were aligned with its edges. The MIL method [26,27] was used to evaluate the fabric tensor. The tolerance error of the prin cipal axes to the cube edges was 1 deg. A linear interpolation scheme was applied in the image rotation. After rotation, the cubic subregions with a side length of 5.3 mm were cropped, coarsened to a resolution of 36 pm (see Fig. 1 top), segmented using a single-level threshold [28], and subsequently all of the uncon nected bone regions were removed. To generate /iFE models, image voxels were converted directly to mesh elements. Image processing and FE model generation were done by MEDTOOL.r Linear homogenization analysis of these trabecular bone cubes was performed using /iFE analysis with KUBCs and PMUBCs. A Young’s modulus of lOGPa and a Poisson’s ratio of 0.3 were assigned to each element [29]. Six independent load cases were applied on the cubes, including three uni-axial and three shear cases. Uni-axial and shear loading cases of PMUBCs are illus trated in Fig. 1 middle and bottom, respectively. From homoge nized /iFE analysis, stress and strain fields were obtained, from which the overall elastic stiffness tensor S fe^ was evaluated using stress and strain volume averages [15]. Subsequently, an orthotropic representation SFEortlo of the stiffness tensor was deter mined by neglecting the nonorthotropic tensor components in the coordinate system of the fabric tensor [15,25]. The norm error cor responding to this assumption can be calculated with 1www.dr-pahr.at
Transactions of the ASME
Fig. 2 Exam ples of highly heterogeneous bone biopsies which were not included in the filtered data set: BV/TV = 11.21% (left) and 24.07% (right)
NE aniso—ortho
®
The orthotropic qFE stiffness components were then fitted to the Z ysset-C um ier model using multiple linear regressions, which yield the m odel parameters. In this study, multiple linear regres sions were solved by a linear solver routine (numpy.linalg.solve(), p y t h o n 2.6.6). In general, the (global) linear system includes 12 nonzero components of the orthotropic stiffness tensor of each sample, w hich were aligned in a vector y y = Xc + e
(4)
where X is a 12 x p matrix containing p and m,-, c is a vector containing p constants of the model, and e is the vector of the residuals. The global linear system is shown in Appendix in matrix form. In an attempt to generalize the fabric-elasticity model without losing positive-definiteness, the three shear com po nents were split from the global linear system
yx = X aca +
(5)
y„ = X ,,^ + e„
(6)
where the subscription A denotes the linear system of normal components o f a stiffness tensor and the subscription p denotes the linear system o f shear components. The split linear systems require two additional exponents k and / as seen in Appendix. A fter computing the model param eters, model stiffness could be predicted by the Z ysset-C um ier model using the obtained model param eters, bone volume fraction, and the complete fabric information. The difference between the model stiffness S m0dei and orthotropic pFE stiffness SpEortho is expressed by the model norm error NEortho_model NE ortho—model
II^FEotat,
§ m o d e l||
ll§FEortJ |
calculated for each subcube. Mean and standard deviation of Psubcube were computed and the coefficient o f variation o f p Subcube calculated with S t d ( /L u b c u h e
(1) The complete data set: a data set w hich includes all the trabecular bone biopsies (2) The filtered data set: a data set which excludes highly heter ogeneous bone biopsies and includes only the trabecular bone cubes with CV below the computed threshold For the complete data set, the m odel parameters were calculated using the global linear system. For the filtered data set, a set of model parameters were calculated using the global linear system and another set using the split linear systems. Linear regressions in log scale were performed between the model stiffness S modei and the orthotropic pFE stiffness S FE„„h„• KUBC-based stiffness of the current methodology was then compared to that o f Gross et al. in linear regressions using the model parameters from Ref. [25]. For the filtered data set, pFE stiffness tensor components based on PMUBCs and KUBCs were compared in a one-to-one relation ship. The norm error N E KUBC- p m u b c describes the difference between pFE stiffness of the two BCs
(7)
||sL UBCs ™rC.0rth0 UBC|||| I rC .ortho_____
The total norm error NE^so. moliei describes the difference of the m odel stiffness S m0dei with respect to the full anisotropic stiffness tensor § fe^ „ IIS f
^>model 11
— ‘
(8)
IIS f b _ . I I
Unlike KUBCs, PMUBCs showed sensitivity to heterogeneity, which required further investigation. Among other tested m eth ods, the coefficient o f variation (CV) o f bone volume fraction was retained to assess heterogeneity. Each bone cube was divided into eight identical subcubes and bone volume fraction (psubcube) was
Journal of Biomechanical Engineering
(9)
The bone cubes were sorted by CV in descending order. Trabecu lar bone cubes with the highest CV were excluded one-by-one. At each step, coefficients of determination (r2) was obtained from the linear regression o f model stiffness to pFE stiffness based on PMUBCs. In the beginning of the cube exclusion, r2 rose steeply. Then, it reached a peak value after which it stabilized and started to fall. The coefficient of variation at the peak value was taken as the threshold to exclude highly heterogeneous bone samples (illustrated in Fig. 2). Hence, two data sets were analyzed in this study.
N E k u b c - pm u b c
NE,aniso—model
)
m ean(psubcube)
PMUBC
||S FEortho
(10)
Table 1 Mean values (±s ta n d a rd deviation) of bone volum e fraction (BV/TV) and degree of anisotropy (DA), and num ber of fem oral trabecular bone samples (n) of the exam ined data sets
Gross et al. [25] Complete data set Filtered data set
BV/TV
DA
n
0.19 (±0.10) 0.25 (±0.08) 0.27 (±0.08)
1.67 (±0.34) 1.54 (±0.20) 1.57 (±0.18)
264 167 126
JANUARY 2 0 1 5 , Vol. 137 / 011002-3
• ■
g -3 .0
Gross etal. (2012) Complete data set Filtered data set
a
O «
‘c 2.5
while the lowest difference is in the diagonal terms of the normal components Xu. Moreover, the rela tive difference decreases with MIL and is therefore minimal along the major fabric direction A33 (see Fig. 5). The exponents of volume fraction k range from 1.87 to 2.01, which agree well with the power coefficient of 2 for the apparent density resulting from a statistical review by Rice et al. [30]. The exponents of fabric l range from 1.00 to 1.68, which are bounded by those from Zysset [23]. Interestingly, the exponents k and l based on PMUBCs are higher than the ones based on KUBCs. Using the split linear systems, the exponents k and / are around 1.3 times higher, except the exponent / of the shear components which is around 1.6 times higher. This indicates that, based on PMUBCs, bone volume fraction and fabric information play a stronger role in the apparent properties, in particular, fabric eigen values have a strong effect on the shear components. Linear regressions of the Zysset-Cumier model to /dFE stiffness yield r^dj from 0.945 to 0.975 for PMUBCs, and from 0.984 to 0.986 for KUBCs. Exclusion of highly heterogeneous bone samples, which violate the notion of RVE, improved significantly the norm errors of both KUBCs and PMUBCs. Splitting linear systems between the normal and shear components improved sig nificantly the norm errors only for PMUBCs, whereas for KUBCs not significant benefit could be established. There were some limitations in this study. First, the PMUBCs are based on orthotropic symmetry of the RVE but trabecular bone cubes are only approximately orthotropic as witnessed by the approximately 5% error on the tensor norm. However, PMUBCs do not require periodicity, which was shown to induce artificial connections between bone and marrow, especially at low volume fraction [9]. Second, trabecular bone images had to be rotated because PMUBCs require orthotropic microstructures aligned with the planes of orthotropic symmetry. The morphology of the rotated bone images deviated slightly from that of the original images. Bone volume fraction and degree of anisotropy were on average 0 1 1 0 0 2 -6
/ Vol. 137, JANUARY 2015
higher after rotation with the mean deviation of 3.46% (±4.29) and 1.09% (±6.57), respectively. Among other morphology infor mation, the highest mean deviation was found in connectivity with an average decrease of 17.37% (±7.38). Fortunately, this unexpected finding did not affect the morphology-elasticity relationships. Third, PMUBCs are sensitive to heterogeneity, which lead to the analysis of two data sets. The filtered data set that excludes highly heterogeneous bone cubes also excludes outliers of stiff ness tensor components (see Fig. 6). The Zysset-Cumier model tends to overestimate elastic properties for highly heterogeneous samples, especially for the ones with low bone volume fraction. In future studies, these two BCs will need to be compared in FE analysis of whole bones. KUBCs represent an upper bound of the apparent stiffness [15,20], They may be appropriate for a trabecu lar bone volume element with a rigid neighborhood, e.g., regions surrounded by cortex or high BV/TV. On the other hand, if a vol ume element is surrounded by similar morphology, PMUBCs may yield more appropriate apparent elastic properties. A solution will need to be found for homogenization of trabecular bone material properties in the presence of marked heterogeneity, i.e., of high gradients of volume fraction or fabric. In particular, the influence of BCs on the apparent yield properties of trabecular bone remain to be explored. Finally, proper homogenization of trabecular bone is expected to improve the accuracy of both linear and nonlinear hFE analysis of bone or bone-implant systems.
Appendix: Multiple Regression Models The global linear system has the form of
M A 22 )
0 0 ln(p) 1 0 0 ln(p)
ln(A33)
1
0 0
ln(p)
ln(mf)
ln(A32)
0
1
0
ln(p)
ln(m3m2)
ln (2i3)
0
1
0 ln(p)
ln(mi/w3)
ln(A2i)
0
1
0 ln(p)
ln(m2m i)
K /% )
0
0
1
ln(p)
ln(m2m3)
M P i 3)
0 0
1
ln(p)
ln(m im 3)
ln(Pi2)
0
0
1
ln(p)
ln(/wi/??2)
ln(/l23)
0
1
0 ln(p)
ln(/??2 /H3)
ln(A3i)
0
1
0 ln(p)
ln(m3/wi)
1
0 ln(p)
ln(/Wi7M2) /
/ In (^ n ) \
Vln(Ai2) )
ln(m f) \
(1
u
In (m\)
(ei \
