|Home | About | Journals | Submit | Contact Us | Français|
The relationship between fabric (a measure of structural anisotropy) and elastic properties of trabecular bone (TB) was examined by invoking morphology and homogenization theory on the basis of micro magnetic resonance images (μMRI) from the distal tibia in specimens (N=30) and human subjects (N=16) acquired at a 160×160×160 μm3 voxel size. The fabric tensor was mapped in 7.5×7.5×7.5 mm3 cubic subvolumes by a 3D mean-intercept-length method. Elastic constants (three Young’s and three shear moduli) were derived from linear micro finite-element (μFE) simulations of 3D grayscale bone volume fraction mapped images. In the specimen data, moduli fit power laws of bone volume fraction (BV/TV) for all three test directions and subvolumes (R2 0.92 – 0.98) with exponents ranging from 1.3 to 1.8. Weaker linear relationships were found for the in vivo data due to a narrower range in BV/TV. When pooling the data for all test directions and subvolumes, BV/TV predicted elastic moduli less well in the specimens (mean R2=0.74) and not at all in vivo. A model of BV/TV and fabric was highly predictive of μFE-derived Young’s moduli: mean R2s of 0.98 and 0.82 (in vivo). The results show that fabric, an important predictor of bone mechanical properties, can be assessed in the limited resolution and SNR regime of μMRI.
Osteoporotic fractures typically occur at sites rich in trabecular bone (TB). The contribution of TB to load bearing has been shown to exceed that of cortical bone in the distal radius and tibia (1). The mechanical competence of the TB network depends on a number of parameters, chief among which are bone volume fraction, trabecular architecture (topology, orientation), and degree of mineralization. The clinical measure of bone mineral density (BMD) expressed in terms of areal or volumetric density is a combination of both volume fraction and degree of mineralization. Experimental studies have revealed that approximately 60–80% of the variance in bone’s mechanical properties can be attributed to BMD (2–4). BMD, however, offers no information on trabecular architecture, and thus, provides an incomplete picture of osteoporotic fracture risk (5,6).
An important characteristic of the TB structure is its directional dependence (structural anisotropy), an attribute recognized early by Wolff (see, for example, (7)). Whitehouse (8) first quantified structural anisotropy by conceiving the line intercept method, which is based on a measurement of the number of line intercepts across marrow along a particular direction, typically obtained from optical micrographs of sections of trabecular bone. When plotted as a function of direction the mean intercept lengths trace an ellipse (also called rose plot). In three dimensions, angular sampling of mean intercept length yields an ellipsoid, which can be expressed as a second-rank tensor known as a fabric tensor (9). Fabric tensors provide compact descriptions of orthotropic structural anisotropy in the form of 3×3 matrices, where the eigenvectors w1, w2, w3 define the main directions of the ellipsoid’s axes and the eigenvalues λ1, λ2, λ3 the axes’ lengths.
Work by Goldstein et al (10) revealed, on the basis of 3D reconstructions of microcomputed tomography (μCT) images of specimens of trabecular bone, that in normal bone, more than 80% of the variance in mechanical behavior can be explained by measures of density and orientation. Cowin linked structural and mechanical anisotropy by establishing a theoretical framework under the assumption that trabecular bone possesses orthotropic symmetry and that bone can be modeled by assuming an isotropic tissue modulus (11). Invoking Cowin’s model, Turner et al subsequently showed that fabric measures could explain 72 to 94% of the variance in elastic constants (12) in TB samples from bovine femur and human proximal tibia. Van Rietbergen et al also demonstrated in 3D images obtained from serial sections of TB from a whale vertebral body that inclusion of fabric substantially improved bone volume fraction-based prediction of elastic constants (13). Additional evidence of the role of structural anisotropy as a contributor to bone mechanical properties was provided by Oden et al (14), showing that BMD was only moderately predictive of failure stress (R2=0.44, p=0.004) while the mean intercept length measured along the loading direction was much more strongly associated with this parameter (R2=0.85, p<0.001). Another study, conducted on the basis of μCT images of femoral neck surgical specimens, indicated that TB structural anisotropy was markedly greater in individuals with hip fracture than in those without fractures (15).
Advances in high-resolution peripheral quantitative computed tomography (HR-pQCT) (16,17) and micro magnetic resonance imaging (μMRI) (18,19) permit noninvasive assessment of bone architecture and measures of biomechanical competence by subjecting the resulting images to micro-finite element (μFE) analysis (20–22). However, the role of fabric and its relationship to bone mechanical properties has not previously been examined under conditions of limited signal-to-noise ratio (SNR) and resolution characteristic of in vivo imaging.
Micro-FE, performed on the basis of high-resolution μCT images, has previously demonstrated to yield excellent agreement with mechanical testing (see, for example, Kabel, et al, (23)). Further, some of the present work’s authors have shown recently that μFE derived mechanical constants based on in vivo resolution images, are highly correlated with those derived from high-resolution μCT images (21).
In this work we examined the contributions from fabric, derived from specimen and in vivo μMR images of the human distal tibia, to the elastic constants estimated from image-based μFE models of the same images. The relationship derived by Cowin linking fabric to elastic measures was then examined relative to the dependence of the elastic measures on bone volume fraction alone. We hypothesize that by including TB orientation in the model linking bone volume fraction and elastic properties, an improved prediction of the elastic properties can be achieved.
High-resolution MR images of thirty left (15) and right (15) human distal tibia specimens of 25 mm axial length previously obtained at 1.5T field strength (21) were processed as described below. In addition, 3T MR images of the distal left tibia from sixteen subjects (seven healthy young and nine post-menopausal) were analyzed for the purpose of evaluating the method in images acquired in vivo at the same image voxel size. Images of the younger subjects (age range: 22–42 years, mean: 34 years) were taken from a recently completed reproducibility study in which two female and five male volunteers had been imaged repeatedly (24). Images of the older subjects represent a subset of an ongoing drug intervention trial in postmenopausal women (age range: 60–84 years, mean: 68 years) who were either osteoporotic or osteopenic, acquired at baseline. The study was approved by the University of Pennsylvania’s IRB and all subjects enrolled provided informed consent. Images of both subject groups were obtained with the same scan protocol as detailed in (24) to achieve an isotropic voxel size of 160μm in 16 minutes 40 seconds.
Each specimen image (Figure 1a) was bone volume fraction (BVF) processed to yield a grayscale image with pure bone and pure marrow corresponding to intensities of 100 and 0 (Figure 1b). The BVF-mapping procedure involves scaling of the image intensities relative to locally-calculated marrow intensity levels to remove intensity variations due to surface receive coils, as described by Vasilic et al (25). Bone volume divided by total volume (i.e. BV/TV) represents the average BVF of all voxels within a TB subvolume. The vertical trabecular plates in the tibia are approximately parallel to the endosteal cortical boundary (26,27), leading to substantially different structural anisotropies in anterior, posterior and lateral locations as visually apparent (Figure 1). Therefore, three 7.5×7.5×7.5 mm3 subvolumes of TB were selected at posterior (P), lateral (L), and anterior-medial (AM) locations (Figures 1c–e). To sample the heterogeneous TB architecture at this skeletal site, the three subvolumes were automatically positioned within the cortical boundary by moving away from the center of each tibia along the x (L), y (P), and −y, −x (AM) directions. The structural heterogeneity at this skeletal site is captured well by the 3D renditions of the three subvolumes (Figures 1c – e) and the resulting mean-intercept-length tensors (Figures 1f – h). Details on the computation of the mean-intercept-length tensor are provided below. A single 7.5×7.5×7.5mm3 subvolume of TB was then extracted from the anterior portion (where SNR was highest) of the BVF-mapped distal tibia within each of the in vivo images (Figure 2).
BV/TV and fabric tensor were calculated for each 7.5×7.5×7.5mm3 BVF-mapped subvolume (Figure 3a). MIL was sampled at 5° angular increments of polar angles ϕ and θ (Figure 3b), resulting in 905 MIL samples distributed evenly over the unit sphere. A threshold value of BVF ≥15% was used to identify the bone phase during intercept sampling. Each 3D Rose plot of the MIL sample points was then fit to an ellipsoid whose matrix representation denotes the MIL structural tensor M. The associated fabric tensor, denoted H, was calculated as the inverse square-root of M (28) so that the largest eigenvalue of H corresponded to the direction of largest MIL. An ellipsoid representation of H (after normalization of the radii) is shown in Figure 3c. Eigenvalues and eigenvectors of H were determined using spectral decomposition. The normalized eigenvalues (λ1, λ2, λ3) and their eigenvectors w1, w2, w3 were used to describe the TB orientation relative to the image coordinate system through the three angles θx, θy, θz defined in Figure 4a. The angle between the principal material axis w3 and z, (the direction of the imaging gradient, approximately parallel to the tibial axis, i.e. the infero-superior direction) was computed as θz = cos−1(w3•z). Angles θx and θy were defined similarly, but restricted to the transverse plane (i.e. xy plane, x and y corresponding to antero-posterior and medio-lateral directions, respectively) such that θy (θx) is the angle between the projection of w2 (w1) onto the xy plane and image y (x) axis.
BVF-mapped cubic subvolumes were subjected to six stress/strain simulations via μFE analysis (Figure 3d). The μFE approach used in the present work (21) closely follows the algorithm described previously by van Rietbergen and others (see, for example, (29)). Briefly, image voxels were modeled as hexahedral elements with a Poisson’s ratio of 0.3 and tissue moduli (Etissue) proportional to BVF such that for pure bone voxels (BVF=100%), Etissue=15 GPa (30), values which have been used widely (see, for example, (31). Three compressive and three shear displacements were individually applied upon the faces of a subvolume while partially restricting the opposing face. Simulations resulted in a 6×6 stiffness matrix defined relative to the laboratory reference frame (i.e. image coordinate system):
Off-diagonal elements, denoted δij, are zero when the three loading directions (i.e. the measurement coordinate axes) parallel the normal vectors of the three planes of orthotropic symmetry. In practice, as shown in Figure 3e, values of δij are non-zero but small relative to the Sij terms. Through tensor rotation of S from the laboratory coordinate system to the coordinate system of approximate orthotropic symmetry, So= RTθ x,θ y,θ z S Rθ x,θ y,θ z, δijs of So will be minimized (zero if the structure were strictly orthotropic). Thus, an orthotropic stiffness matrix So (Figure 3f) can be estimated through transformation of S by minimizing the residual function (29),
The residual function was minimized relative to three rotation angles about x, y, and z axes of the image coordinate system, denoted θx, θy, θz, which were allowed to vary from −45° to 45° in 1° increments. Vectors u1, u2, u3 define the mechanical axes of approximate orthotropic symmetry, based on the rotation of the image coordinate system via Rθ x,θ y,θ z. The compliance matrix was then calculated by C = So−1 (Figure 3g).
The nine orthotropic elastic constants (three Young’s moduli, Eii (i=1, 2, 3), three shear moduli, Gij (i, j=1, 2, 3; i≠j), and three Poisson’s ratios, vij (i, j=1, 2, 3; i≠j)) were determined from the expressions in Figure 3h for the orthotropic compliance matrix Co (i.e. Co with all δij replaced by zeros). The error in the stress-strain relationship (29) caused by forcing orthotropic symmetry was calculated as:
where I is the 6×6 identity matrix and || ||is the matrix norm operator.
Cowin’s model (11) expresses the nine coefficients of the orthotropic compliance tensor Co in terms of fabric eigenvalues λ1, λ2, λ3, the second invariant of the fabric tensor Π= λ1λ2+ λ1λ3+ λ2λ3, and nine functions of the bone volume fraction k1-k9 (12):
where i, j=1,2,3; i≠j. Eq. [4d] expresses the power-law dependence of the elastic modulus on the material volume fraction ρ where kam, kbm are eighteen empirically determined constants and m=1, 2, …, 9. Prior studies found that bone stiffness can be described as a power of material volume fraction with an exponent α in the interval of 0.5 to 3.0 (32–34).
The alignment between the fabric eigenvectors w1, w2, w3 and orthotropic mechanical axes u1, u2, u3 (as determined from minimization of Eq. 2) was evaluated and expressed in terms of offset angles Ω1, Ω2, Ω3 (Figure 4b). For simplicity, the average offset angle Ωavg was computed from Ω1, Ω2, Ω3. Assuming Ω1, Ω2, Ω3 are small per Cowin’s assumption, the fabric eigenvalues, BV/TV (in place of ρ), and the nine μFE-calculated compliance entries were input into the matrix formulation of Eq. :
Here, A is a 9×9 matrix containing the fabric information and engineering strain conversion factors of Eq. . The 18 unknown constants kαam, kαbm were found using a least-squares approximation for each exponent α, which was varied between 0.8 and 3 in increments of 0.1. Goodness-of-fit expressed in terms of R2 was calculated for each correlation and then adjusted, when necessary, to account for the degrees of freedom in multiple parameter-fits according to R2adj=1− (1− R2) (N−1)/(N−K−1), where N is the number of measurements and K is the number of unknowns (12). The joint fit of Eq.  relative to the pooled specimen or in vivo data involved N=90×9=810 (16×9=144) and K=18+1=19. The 18 unknown constants kαam, kαbm and exponent α were calculated using the joint fit. When considering regressions of specific elastic constants, e.g. the three Young’s moduli, E11, E22 and E33, N is 90×3=270 for the specimen, and 16×3=48 for the in vivo subset, while K = 19. For comparison, the nine μFE-derived elastic constants were also correlated to linear and power-law models of BV/TV.
Finally, analysis of variance (ANOVA) was performed to assess the ability of the best-fit orthotropic model to distinguish young, healthy subjects from the post-menopausal sub-group. Significance levels for structural measures, μFE-derived elastic constants, and model predicted elastic constants were compared.
Means and standard deviations of structural measures for specimen and in vivo datasets are summarized in Table 1a. Fabric eigenvalues λ1, λ2, λ3 and orientation angles θx, θy, θz confirm the directional dependence visually apparent in the specimen subvolumes (Figure 1). Fabric ellipsoids indicate that the principal direction of trabeculae coincides with the longitudinal direction of the tibia (i.e. image z axis, parallel to the external field Bo). This is shown by λ3 being the largest eigenvalue and θz being small (i.e. w3 is closely aligned with z, Table 1a). The minor axes differ substantially between the three anatomical regions (labeled L, P, and AM in Figure 1). The average offset angle (Ωavg: mean of Ω1, Ω2, Ω3) indicates that the fabric axes and orthotropic mechanical axes are approximately 11.2° (specimen) and 9.2° (in vivo) offset from each another.
Means and standard deviations of the orthotropic elastic constants and orthotropic errors for specimen and in vivo datasets are summarized in Table 1b. Orthotropic errors ε based on S are 30±17% and 24±11% for specimen and in vivo data, respectively. After calculation of So via minimization of Eq. 2 values are 18±9% (specimen) and 16±7% (in vivo), corresponding to a 40% and 33% reduction in orthotropic error.
In Figure 5, elastic moduli are plotted as a function of BV/TV for the three sub-volumes in the specimen images. We note that the computed Young’s moduli fit power laws of bone volume fraction very well, with three distinct families of curves, each pertaining to a particular test direction (R2 = 0.92 – 0.98, exponent α = 1.3 – 1.8 where larger α values are found for transverse Young’s moduli). As expected, E33, the Young’s modulus approximately in the infero-superior loading direction, is largest, whereas medio-lateral and antero-posterior directions yield similar but lower stiffness values. The three distinct sets of curves are a direct manifestation of the structural arrangement of the trabecular structure, paralleling the results of the fabric calculations.
Table 2 lists coefficients of determination of all nine elastic constants (individually and pooled) from BV/TV in specimen and in vivo images. Linear models were used for fitting in vivo measures Eii and Gij to bone volume fraction because of the much smaller range in the in vivo TB sub-volumes, yielding generally weaker correlations than did the specimen data (Table 2). With the exception of shear moduli, pooled constants were not significantly predicted by BV/TV.
Goodness-of-fit (R2) for the models incorporating fabric is plotted in Figure 6 as a function of exponent α. An optimum R2 relative to the combined elastic constants was found for α=1.5. Table 3 summarizes the model predictions for the pooled μFE-computed elastic constants by (BV/TV)1.5 and fabric measures.
Figure 7 plots μFE-computed elastic moduli for both specimen and in vivo data versus those predicted on the basis of the model considering bone volume fraction and fabric (Eq. , using α=1.5). This model predicts the elastic moduli very well, regardless of loading direction for both the specimen and in vivo data sets (Radj2=0.975 and 0.928, respectively). In addition, in vivo μFE-derived elastic constants are well predicted (Radj2=0.884) by the model derived using the specimen data (see Figure 7c and far right column of Table 3).
Lastly, Table 4 lists compressive and shear moduli for the two in vivo sub-groups. Except for one (E11) all compressive and shear moduli were lower in the postmenopausal sub-group. These differences were significant for E22, G12 and G23 (p<0.05 to <0.0005) for both μFE and model-derived parameters. The difference in moduli was supported by lower BV/TV in the post-menopausal sub-group, 9.52±0.95 % versus 11.3±0.73 % (p<0.05), however, the mechanical parameters were stronger group differentiators.
This study investigated the relationship between morphology and finite-element derived elastic moduli in μMR images of TB subvolumes from the distal tibia both in specimen images acquired under in vivo scanning conditions and in vivo in a small cohort of human subjects. Elastic constants derived from μMRI are comparable in magnitude to values reported by other investigators, including data obtained by direct mechanical testing of samples from the proximal tibia (elastic moduli range: 10–1500 MPa) (35) and elastic constants determined from μFE analysis of HR-pQCT images of the distal tibia (mean elastic modulus: 720±297 MPa) (1). Our results show that variations in μFE-derived elastic constants associated with anatomical location and test direction were not adequately explained by bone volume fraction in either the specimen or in vivo data, exhibiting a strong regional and test-direction dependence. These findings agree with previous studies conducted on the basis of high-resolution μCT imaging evaluating the relationship between morphological measures and elasticity (2–4,36).
When including the MIL-derived fabric tensor information via Cowin’s model of orthotropic materials, the nine orthotropic μFE-calculated elastic constants could be predicted with high coefficients of determination (mean R2adj: 0.92 in specimen, 0.82 in vivo) independent of sub-volume location and test direction. Additionally, the model-predicted elastic properties were found to distinguish post-menopausal women with reduced bone mass from younger healthy subjects with significance levels similar to those achieved by μFE modeling.
The dependence of the orthotropic elastic constants on volume fraction and fabric has previously been shown in high-resolution reconstructions of TB specimen image data (R2adj ranging from 0.88 to 0.99) (13,35,36). However, until now, the validity of Cowin’s model had not been demonstrated on the basis of in vivo images.
The 18 empirically determined model parameters did not appear to change substantially with input images from the same anatomical location, as evidenced by their ability to predict elastic constants derived from in vivo images using the specimen-derived model parameters. This result suggests the possibility of developing a model on the basis of a large set of specimen datasets, potentially in concert with mechanical experiments, and then applying the model to structural measures assessed from images acquired in vivo to predict bone mechanical properties in patients. This finding is of particular relevance considering the spread in BV/TV values and moduli between specimen and in vivo images.
The difference in mean bone volume fraction (and thereby elastic moduli) found between specimen and in vivo data is somewhat surprising and demands some scrutiny. One possible reason is that spin dephasing at trabecular bone boundaries caused by the susceptibility difference between the slightly paramagnetic fixing solution and diamagnetic bone is likely to inflate BV/TV values in the specimens relative to values obtained in vivo (37) (since the bone-fixing solution susceptibility difference is greater than that between bone and marrow).
The validity of Cowin’s model for the current data is particularly noteworthy. The degree to which TB complies with the notion of orthotropic symmetry was evaluated in terms of the orthotropic error ε (see Eq. ). Forcing orthotropic symmetry was associated with a ~17% error in the calculated elastic constants. This error is considerably larger than the 6% error found in μCT images of (4mm)3 TB specimens from seven distinct human metaphyseal locations by Zysset et al (35). One explanation for the relatively large mean error in the present work is the proximity of the TB subvolumes to the cortical boundary where the structural, and thus the mechanical environment, is expectedly more heterogeneous. Trabecular bone subvolumes used in Zysset et al’s work were likely extracted from more centrally located regions of the skeleton where the architecture is more homogeneous. Deviations from orthotropic symmetry also impact the degree of alignment between fabric and orthotropic symmetry axes. The average offset angles Ωavg were 11.2° (specimen) and 9.2° (in vivo), data that are in excellent agreement with those found by Zysset et al (~11°) (35). Increasing angular offsets Ωavg were associated with increasing orthotropic error for specimen and in vivo data (see Figure 8). Therefore, the error in the orthotropic elastic constants and the model prediction were augmented by the degree to which the TB sub-volume adhered to orthotropic symmetry. This source of error was minimized by excluding TB sub-volumes with orthotropic errors greater than 4% (12). As mean errors in the measured Young’s (Shear) moduli of 9.5% (1.1%) were found for offset angles of 10° (38), it would have been necessary to exclude a substantial number of datasets from this study, as it was in Turner’s study, if utilizing such an extreme cut-off value. Instead, all datasets were retained at the expense of lower accuracy in the model prediction.
The spatial variations in architecture, notably in terms of structural orientation and trabecular density between the three subvolumes examined (Figure 1), is substantial. Thus, inclusion of a larger volume (> 1–2 cm in width) would likely violate the condition of orthotropy and the underlying theory would be inapplicable. For this reason, prediction of mechanical parameters in terms of structural orientation are limited to relatively small volumes (subvolumes with sides ranging from 0.4 to 1cm), which may not be representative of the bone’s overall mechanical properties. Moreover, larger volumes are now readily amenable to image-based FE analysis. However, the method provides insight into the regional structural and mechanical make-up of the trabecular network. By decomposing the mechanical behavior of TB into contributions from volume fraction and fabric, its underlying causes can be inferred. For example, it is understood that structural anisotropy increases with age and osteoporosis, thereby rendering the bone more likely to fail due to buckling (see, for example (15)). Further, it has previously been shown that Cowin’s model parameters do not substantially differ for normal and osteoporotic bone on the basis of micro-computed tomography (39). It would therefore be of interest to examine these structural alterations during disease progression and treatment.
This work is, to the authors’ knowledge, the first application of Cowin’s model to images of the distal tibia acquired either under conditions matching those in vivo or actual images obtained in human subjects by μMRI. The data emphasize the role of fabric in bone mechanical properties, and that the relationship between fabric and mechanical constants can be assessed in the regime of limited resolution and SNR of high-resolution magnetic resonance. The findings should be equally relevant to other imaging modalities, notably HR-pQCT, that have similar point-spread function limited resolution.
The authors would like to acknowledge support from the National Institute of Health through NIH R01-AR41443, NIH RO1-AR53156, NIH R01-AR55647, and NIH F31-EB74482. We would also like to recognize members of the Bone Bioengineering laboratory at Columbia University for providing the distal tibia specimens.