Preventing femoral fractures is an important goal in osteoporosis research. In order to evaluate a person's fracture risk and to quantify response to treatment, bone competence is best assessed by bone strength. Finite-element (FE) modelling based on medical imaging is considered a very promising technique for the assessment of in vivo femoral bone strength. Over the past decades, a number of different FE models have been presented focusing on the effect of several methodological aspects, such as mesh type, material properties and loading conditions, on the precision and accuracy of these models. In this paper, a review of this work is presented. We conclude that moderate to good predictions can be made, especially when the models are tuned to specific loading scenarios. However, there is room for improvement when multiple loading conditions need to be evaluated. We hypothesize that including anisotropic material properties is the first target. As a proof of the concept, we demonstrate that the main orientation of the femoral bone structure can be calculated from clinical computed tomography scans. We hypothesize that this structural information can be used to estimate the anisotropic bone material properties, and that in the future this could potentially lead to a greater predictive value of FE models for femoral bone strength.
Osteoporosis is a very frequent skeletal disorder that is characterized by compromised bone strength predisposing to an increased fracture risk (NIH 2000), mainly for the proximal femur, vertebrae and distal radius. Lifetime fracture risk for women and men older than 50 rises to 40–50 and 13–22 per cent, respectively (Johnell & Kanis 2005) and the resulting osteoporotic fractures cause increased disability, morbidity and mortality (Cummings & Melton 2002). Owing to ageing of the world population, osteoporosis will become not only an increasing social but also an increasing economic problem. For example, European medical care costs associated with osteoporosis were estimated at €31.7 billion in 2000 and are expected to rise to €76.7 billion in 2050 (Kanis & Johnell 2005). This prognosis points out the importance of preventing and reducing osteoporotic fractures.
Today, most treatment strategies for osteoporosis focus on medication, be it antiresorptive such as bisphosphonates (Black et al. 2000), anabolic such as parathyroid hormone (Bauer et al. 2006) or combination therapy (Bone et al. 2000). Recently, weight-bearing and high-impact exercises have received some attention because of their potential effect on preserving bone density and preventing osteoporosis (Heinonen et al. 1996). However, a major disadvantage of these exercises is that they may increase the risk of fracture because of their high impact. Therefore, other training methods with less fracture risk are being sought. Whole-body vibration (WBV) is such an alternative method that combines different aspects of impact exercises and power training. It has been shown that high-frequency WBV can stop bone loss or even slightly increase bone mass (Verschueren et al. 2004).
Regardless of what treatment strategy is used, its effect on bone quality should ideally be evaluated in terms of bone strength, the best bone parameter to quantify fracture risk. At present, bone mineral density (BMD) measurements with dual-energy X-ray absorptiometry (DXA) are being used as a surrogate measure of bone strength. However, the value of BMD measurements at a single time point is limited, because it has been shown that bone density explains only 70 per cent of the variance in the mechanical properties of bone (Goulet et al. 1994). This result is in line with the finding that several clinical trials, investigating the effect of antiresorptive drugs, have shown that BMD cannot predict bone fracture risk accurately (Cummings et al. 2002; Graeff et al. 2007). Furthermore, the difference between two consecutive BMD measurements provides only poor indications of changes in femoral strength due to osteoporosis treatment. A 40–50 per cent reduction in the risk of fracture through the use of calcium and vitamin D supplementation has been reported, whereas BMD increased by only 1 per cent (Chapuy et al. 1992; Dawson-Hughes et al. 1997). Hence, although for large populations DXA measurement of BMD may provide an indication of the changes in femoral strength induced by osteoporosis treatments, for individual patients the predictive capacity of bone density is of more limited value. BMD explained only approximately 4–28 per cent of the 35–50 per cent reduction in vertebral fracture risk after antiresorptive treatment (Ettinger et al. 1999; Black et al. 2000; Reginster et al. 2000). Sarkar et al. (2002) even found that BMD can explain only 4 per cent of the reduction in the risk of vertebral fracture after a treatment with raloxifene; 96 per cent remains unexplained.
Hence, other aspects such as bone geometry and its internal architecture need to be considered too. Computed tomography (CT) is a methodology to measure both density and structure in a single measurement. Especially, the combination with finite-element (FE) analysis, the most widely used computational technique for structural analysis in engineering, seems promising. CT-based FE models can provide insights into load transfer through the bone architecture and with that help our understanding of how differences in bone microarchitecture influence bone strength. They integrate bone density measurements with bone geometry in a three-dimensional fashion, thereby effectively compensating for the limitations of DXA. As a consequence, CT-based FE models have been proposed as an alternative to BMD in the assessment of bone status. The potential of CT-based FE models has been clearly recognized, and over the past years a number of papers have been published on several methodological aspects of FE models, such as mesh type, material properties, failure mechanisms and loading conditions.
In this paper, we will present a review of the different aspects of CT-based FE modelling of the femur. Based on this review, we hypothesize that it could be possible to improve the predictive value of these FE analyses for in vivo femoral strength in individual patients. One of the aspects that has been rather unexplored until now is the use of anisotropic material properties in a multi-level FE model. We hypothesize that the anisotropic nature of bone is directly related to trabecular orientation and that it is possible to derive it from clinical CT data. The methodology to derive this trabecular orientation based on local grey values will be described in this paper. In the future, it might be possible to use this microstructural information to create patient-specific FE models that incorporate anisotropic material properties in a hierarchical fashion.
2. Creating subject-specific models of the human femur
To build an FE model, a three-dimensional representation of the patient's femur is needed. Such a representation can be acquired from CT scan data. In short, the different attenuation of X-rays by bone and soft tissue leads to images in which different tissues have different grey values. After segmentation, the femoral bone can be digitally extracted from the surrounding tissues. These segmented femur slices are then joined together in a three-dimensional representation of the femur, from which an FE model can be built. In the creation of accurate models, several methodological aspects need to be addressed. These aspects are discussed below.
(a) FE mesh
A basic concept in FE modelling is to divide a complicated object artificially into a finite number of small and manageable pieces, so-called elements; the behaviour of each element is described in mathematical equations, which are then solved using computers. To generate an FE mesh from CT scan data, there are basically two possible meshing techniques: geometry-based meshing and voxel-based meshing. An example of both techniques is shown in figure 1.
(i) Geometry-based meshing
Geometry-based meshing encompasses a two-step procedure. First, the contours of the femur are extracted from the CT scan. Based on these contours, the surface of the bone is reconstructed, from which, in a second step, an FE mesh with tetrahedral (Taddei et al. 2006) or hexahedral (Lotz et al. 1991a,b) elements is built. Comparison of the two element types showed little difference between the respective models (Ramos & Simoes 2006; Yosibash et al. 2007), as long as the meshes are sufficiently well refined.
(ii) Voxel-based meshing
Voxel-based meshing (Keyak et al. 1990) comprises only one step: the segmented voxel data are directly converted into an FE model with brick elements, hence meshes can be generated very quickly. The main advantage of this approach is that it is a simple automated technique that can deal with any shape of arbitrary complexity. However, curved surfaces cannot be represented properly by brick elements. Hence, where voxel meshes can provide accurate internal stresses and strains, the jagged-edged surface causes peak stresses and strains, which is a disadvantage when accurate surface stresses and strains are needed. Post hoc smoothing algorithms can alleviate these surface errors (Charras & Guldberg 2000). A more fundamental approach consists in smoothing the surface before performing the actual FE calculations. However, conflicting results have been reported. Where Boyd & Muller (2006) showed that peak stresses on the surface of a sphere model were lowered by a factor of 4 after smoothing, Arbenz & Flaig (2007) found only minimally lower stresses for a similar model. Hence, further investigation on smoothing of voxel-based meshes is needed.
Only a limited number of voxel-based models for the prediction of femoral bone strength have been presented. Hence, it seems that the full potential of voxel-based meshes has not been explored. However, the ease in generating voxel models, combined with the recent massive increase in speed and size of very large FE models, appear to provide good conditions for further exploration.
Unlike macroscopic-scale CT-based FE modelling, voxel meshes are already generally accepted as a meshing technique for μCT-based modelling of trabecular bone. For these complex structures, it has been shown that tetrahedral and voxel-based meshes yield almost identical results at higher resolutions (Ulrich et al. 1998) and that there is a very good agreement between the predictions of a voxel-based μFE model and experimental measurements (Ladd et al. 1998; Kabel et al. 1999; Homminga et al. 2003).
(b) Material properties
The use of appropriate material properties in FE analyses has received a lot of attention. Today, it has become the state of the art to take into account the inhomogeneity of bone, based on CT numbers, the attenuation coefficient of the material present in each voxel. Depending on the research question, a number of relationships need to be defined, which include: the conversion of CT number to elastic modulus; the use of a linear or a nonlinear stress–strain relationship; the assignment of isotropic or anisotropic material properties; and the implemented failure theory. When generating such heterogeneous models, it has to be noted that in some FE codes only a limited number of materials can be defined for each model. Therefore, the number of materials often gets reduced by dividing the range of elastic moduli into a limited number of intervals (Zannoni et al. 1998).
(i) From CT number to Young's modulus
To estimate the elastic modulus of the elements based on CT numbers, a two-step procedure is normally taken: (i) convert the CT numbers into bone density and (ii) estimate Young's modulus from bone density.
Bone density can be directly assessed from the CT scan data by its linear correlation with the CT number (McBroom et al. 1985). This relationship can be obtained by scanning a phantom, with known material densities, together with the patient. Bone density is often expressed as apparent density, which is calculated by dividing the bone's wet weight by the total volume of the bone under investigation, including its pores. Also wet density, calculated as the wet weight divided by bone tissue volume, and ash density, calculated as ash weight over bone tissue volume, are being used.
In contrast to the first step in the conversion of CT data into material properties, the relationship between bone density and bone elasticity is less clear. In a recent literature review, Helgason et al. (2008) have investigated density–elasticity relationships derived from experimental tests. After selecting and normalizing 22 relationships (figure 2), they found considerable differences between the investigated studies that could only be partially explained by methodological inter-study differences such as boundary conditions and specimen geometry. Morgan et al. (2003) investigated whether site specificity exists in the density–elasticity relationships. They tested on-axis bone samples from different anatomical sites with a protocol that minimized the end artefacts. They found that, when the site-specific trabecular structure was taken into account, there were no inter-site differences in tissue moduli. Hence, the differences in the Young's modulus–density relationship are due to site-specific variations in trabecular architecture. This implies that, in the FE modelling of bone, site-specific density–elasticity relationships should be used.
(ii) Isotropic or anisotropic properties
It is generally known that bone material, cortical as well as cancellous, is orthotropic rather than isotropic (Yang et al. 1998). Cortical or compact bone is the dense material that forms the outer shell of the bone. Cancellous or trabecular bone is the porous structure enclosed by the cortical shell. The main orientations of this structure appear adapted to the main loading directions, a feature known as Wolff's Law. A meta-analysis showed that, for trabecular bone samples, approximately 64 per cent of the variabilities in bone stiffness and strength are determined by its apparent density (Wehrli et al. 2002). Bone microstructure explains a significant part of the remaining variation: by including architectural parameters in the analysis, the predictive power increases to over 90 per cent for bone stiffness (Turner et al. 1990; van Rietbergen et al. 1998; Yang et al. 1998). Nevertheless, in most femoral FE models of the human femur, isotropic material properties are assigned to the elements (Keyak et al. 1990; Lotz et al. 1991a; Taddei et al. 2006). This approach is driven by the fact that CT scans provide only information about bone density in every voxel, and do not contain information on local anisotropy. Two recent studies have indicated that this is an appropriate methodology (Peng et al. 2006; Baca et al. 2008). By comparing models with isotropic and orthotropic material properties, Peng et al. and Baca et al. concluded that they gave very similar results. However, it has to be noted that Peng et al. did not vary the orthotropic orientation for each element; instead, they assumed that the principal directions of the bone are identical throughout the entire bone. Baca et al. introduced more diversity by defining 30 orthotropic orientations. However, for each element, the orientations were based on the organization of the central vascular canals of the Haversian lamellar systems in the cortical bone. Therefore, they gave an indication of the anisotropic structure of cortical bone but not of the internal trabecular architecture. Hence, none of these models represent local anisotropy of the more complex internal architecture and it is questionable whether their conclusions hold. In addition, the small differences between the isotropic and anisotropic models that were observed could also be related to the applied loading. Both Peng et al. (2006) and Baca et al. (2008) were simulating stance phase, which is a physiological loading condition. Verhulp et al. (2006) suggested that the transverse modulus, introduced by assigning orthotropic material properties, does not play a very important role when physiological loading conditions such as single-legged stance are considered, since bone architecture is adapted to those conditions, but that differences between the isotropic and orthotropic models can be substantially larger for non-physiological loading such as falls to the side.
(iii) Nonlinear models and failure theories
At yielding and subsequent failure, bone loading is no longer linearly related to displacement and may show post-yielding behaviour (Keaveny et al. 2001) or immediate failure (Schileo et al. 2008). In most studies, however, linear material properties are assigned to the elements of FE models to predict the load at which a femur will break and the location where a failure will initiate. The use of these linear models is in principle limited, because they are valid only up to the onset of local bone failure, whereas the use of nonlinear models is also justified past the onset of bone failure until total failure.
Lotz et al. investigated the predictive value for fracture load of linear (Lotz et al. 1991a) and nonlinear (Lotz et al. 1991b) models. In addition, they studied the predictive value for femoral fracture load of stress- and strain-based theories in linear models (Lotz et al. 1991a). For both linear and nonlinear models, they found a very good agreement between measured and calculated fracture loads. Additionally, the nonlinear model could predict bone yielding well. The results also showed that application of a strain-based failure criterion improved the prediction of the failure load. Similar to Lotz et al., Schileo et al. (2008) found that strain-based criteria resulted in a better prediction of the fracture onset and the level of failure risk. These results seem obvious, as it is generally accepted that strain parameters describe bone failure. However, in another comparative study of failure theories, Keyak & Rossi (2000) found just the opposite after investigating eight different failure theories; when using linear FE models, stress-based theories predicted fracture load better than strain-based theories. Hence, no consensus seems to exist on which failure criterion predicts femoral fracture load best.
Keyak and co-workers also investigated the effect of linear and nonlinear material properties on the accuracy of their FE models. For the linear models, they found significant relationships between the measured and predicted fracture loads (r2=0.75 for single-legged stance; r2=0.90 for fall; Keyak et al. 1998). For a single-legged stance, precision for identifying subjects with below-average fracture loads was found to be insufficient. FE models with nonlinear material properties overcame this problem (Keyak 2001). These models were, in contrast to the linear models, precise enough to identify 20 per cent of the subjects with below-average fracture loads with a reliability of at least 97.5 per cent. Hence, these data show that, by using nonlinear material properties, more subjects at risk of breaking a hip could be identified, especially those with low fracture loads, who are at the greatest risk.
(c) Loading conditions
Fractures occur when the loads acting on a bone are higher than the strength of that bone. Hence, besides an accurate assessment of bone strength, information on the loading conditions that encompass the greatest risk of femoral fracture is also needed.
Most femoral fractures are caused by a fall to the side, but spontaneous fractures also account for a significant amount of hip fractures. Consequently, one or both of these situations are simulated in FE analyses, by modelling the impact of a fall to the side and/or by modelling single-legged stance (Lotz et al. 1991a,b; Keyak et al. 1998). Keyak et al. confirmed the importance of these two loading conditions by studying the effect of loading direction on the femoral fracture load. They found that fracture load (i.e. bone strength) depends on the direction of applied loading; furthermore, they found that fracture loads were the smallest for three loading conditions in particular: a fall to the posterolateral aspect of the femur; single-legged stance; and stair climbing (Keyak et al. 2001). They also showed that the fracture load for one of these loading conditions accounts for only 81 per cent of the variance in the fracture load for the other loading condition, leaving 19 per cent unexplained (Keyak 2000). This is probably related to the complex geometry and internal trabecular architecture of the bone, and indicates the limitations of isotropic material properties.
Although abnormal contraction of major rotator muscles could induce hip fracture (Yang et al. 1996), it has been shown that excluding the abductor muscle forces in FE models does not have a notable effect on the predicted femoral strength and fracture location (Keyak et al. 2005; Cristofolini et al. 2007). Hence, muscle forces can be omitted in FE models of femoral bone strength.
3. Experimental validation
The gold standard of measuring bone stiffness and strength is through biomechanical testing. Hence, FE models are best validated by mechanical testing of the bone under investigation. Of course, the boundary conditions as applied in these FE models should mimic the loads as applied in the experimental set-up. Depending on which research question is to be answered, different set-ups simulating different loading conditions have been investigated. It is beyond the scope of this paper to discuss all these different set-ups in detail, and to comment on their precision and accuracy. Taking into account all inter-study differences, it is difficult to compare the accuracy of one FE model relative to another and to determine the effect of the number of femurs used on the accuracy of the model. To the authors' knowledge, the most rigorous validation in the prediction of femoral fracture risk was performed by Keyak et al. (2005). They used a group of 18 femora to fine-tune their parameter settings. Afterwards, a group of 26 other femora were used to evaluate the accuracy and the precision of the model. They found that the second group predicted fracture load for a single-legged stance with nearly the same accuracy (r2=0.88 and 0.83, respectively) as the development group.
4. Improved capacity for anisotropic FE models?
In most FE models heterogeneous isotropic material properties are assumed, accounting for the bone density distribution but not for internal architecture. However, it has been shown that the accuracy and the precision of these isotropic models for the prediction of femoral strength are limited (Cody et al. 1999; Keyak et al. 2005). This is most probably due to the fact that bone is orthotropic rather than isotropic.
The orthotropic nature of the internal trabecular architecture is a very important structural parameter that contributes to the global femoral strength. However, from the preceding literature review, it can be deduced that only very little research has been undertaken on anisotropic FE modelling of the proximal femur. Studies that did account for anisotropy either used methods to derive internal architectural parameters that cannot be used in vivo for patient-specific models (Wirtz et al. 2003) or assigned orthotropic material properties to all elements in a global coordinate system not accounting for local differences in orthotropy (Taylor et al. 2002; Peng et al. 2006). Hence, no valid method exists to derive local anisotropy parameters in vivo. As a consequence, it has not been possible to analyse patient-specific FE models with anisotropic material properties. However, based on our review, we hypothesize that inclusion of local anisotropy could help to improve the prediction of femoral strength in vivo.
5. Determination of trabecular orientation from clinical CT scans
The anisotropic nature of trabecular bone is directly related to the anisotropic trabecular architecture. Unfortunately, state-of-the-art technology does not allow detailed images of the trabecular architecture. This begs for the development of multi-level technology to incorporate anisotropic properties. In this section, the methodology to derive trabecular orientation from clinical CT data and the concept of applying these data to build patient-specific FE models are described. An overview of this procedure is represented in figure 3.
From clinical CT scans, it is not possible to derive trabecular architecture in detail. However, the distribution of the grey values in the images gives an indication of the main orientations of the internal structure, containing a measure of local anisotropy. To derive these orientations from two-dimensional CT slices, the method of line fraction deviation (LFD) developed by Geraets (1998) for small square areas was extended in order to be applicable to CT scans of whole femora. This new method scans the image with a squared region of interest, for which the LFD is calculated in every pixel. As the region of interest covers several pixels in every position, orientation will be calculated more than once for every pixel. The resulting orientation for a pixel is calculated by averaging all orientations found for that pixel.
As a proof of principle, we applied our method to a frontal slice of an in vivo CT scan of the right femur of a 63-year-old woman (Siemens Sensation 16; 120 kV; 140 mA s; pixel size 1 mm), participating in a study on whole-body vibration training (Verschueren et al. 2004). The participant gave her informed written consent before enrolment. The CT data were segmented (figure 4a,b) and the cortical bone was separated from the trabecular bone. Subsequently, the trabecular and cortical parts of the bone were analysed separately. During the analysis of the trabecular bone, a 10×10 pixel region of interest was moved over the slice, calculating the local orientation for each pixel with the LFD method. For the cortical bone, the region of interest was 12×12 pixels. As the cortical bone of the proximal femur is only a shell, this region of interest would often include soft tissue and cortical bone. Therefore, the grey values of the surrounding tissues were set to zero during segmentation. Subsequently, the LFD was calculated only for those regions of interest that contained pixels representing cortical bone. The results of these analyses (figure 4c) show good correspondence with the structural orientations in the regions where structure is most pronounced, which is the entire cortical bone, and the trabecular bone in the femoral head. In the other parts of the femur, structural orientations were less definable. This is probably due to the limited resolution of the CT scan, as the pixel size (1×1 mm) is too large to capture the thin trabecular bone structures. At the borders of the cortical bone, some calculated orientations differ from the expected direction, probably due to partial volume effects.
6. Towards improved prediction of bone strength
The newly developed method is now being expanded to three dimensions in order to calculate the spatial orientation of bone tissue for every voxel of a femur. These orientations can then be incorporated into multi-level FE models of the femur, as they account for local anisotropy based on structural parameters at the microscopic scale. These might have the potential to predict femoral strength more accurately than isotropic FE models because they mimic the structural and material properties of the bone better.
With the recent introduction of a new generation of clinical CT scanners that are able to image the proximal femur with a resolution of 0.3–0.4 mm isotropic resolution, there is potential to refine the assessment of trabecular orientation, although it will still not be good enough to visualize individual trabeculae. A single trabecula can be imaged with high-resolution peripheral quantitative CT (HR-pQCT, XtremeCT, Scanco Medical), providing an isotropic resolution of 80 μm. Unfortunately, this scanning technique is applicable only to ‘peripheral’ bones.
As osteoporosis is becoming an increasing social and economic problem (Kanis & Johnell 2005), a diagnostic tool that can detect osteoporosis in an early stage and predict bone strength accurately in individual patients is indispensable. Also, in view of studying the effects of new osteoporosis treatments, a method is needed that is more sensitive to small changes in bone strength than DXA, the current standard diagnostic tool. Rather than using a large population of patients for a long time until osteoporotic fracture appears, such a technique should allow for a much earlier detection of changes in fracture risk. FE modelling is such a technique. CT-based FE models provide the means for a mechanistic understanding of bone strength, which goes beyond the statistical correlations based on BMD; they have the potential to become a standard medical tool in the prediction of femoral strength, the best parameter to quantify fracture risk (Cody et al. 1999). The aim of this study was to review the current standing in FE modelling of the human femur and to identify which aspects of the modelling procedure should be improved to enhance the predictive value of FE models in the diagnosis of osteoporosis.
It was found that geometry-based as well as voxel-based FE models can provide moderate to good predictions of femoral bone strength, especially when the models are tuned to specific loading scenarios. However, there is room for improvement when multiple loading conditions need to be evaluated. Improvement may be achieved by a better description of bone failure, as no standardized algorithms exist for assessing fracture from FE analyses. It should be noted that the precise failure mechanisms and associated material properties are still not well understood.
This review also showed that very few studies addressed the problem of incorporating local anisotropy into their models and those that did based the anisotropic parameters on a priori assumptions (Taylor et al. 2002; Wirtz et al. 2003; Peng et al. 2006; Baca et al. 2008). For that purpose, we developed a new method that analyses two-dimensional slices of a CT scan using the LFD. Although the results of the analysis are promising, this method has some limitations. First, the results can vary through changes in the size of the region of interest. To find an optimal size for this region of interest, more CT images need to be analysed. Second, for now, only two-dimensional images can be analysed. At the moment, the method is being extended to three-dimensional images. Third, the method could not be validated, as no detailed information on the main orientations of the bone was available. However, visual inspection could be considered as a first validation. In a region similar to the femoral head, where the orientation of the bone is clearly visible, the new method is able to detect this direction. Fourth, owing to the limited resolution of the CT scan, it is not possible to distinguish trabecular structures in some parts of the femur. Hence, it is also difficult to calculate the main direction in those parts. However, with modern CT scanners, patients can be scanned with a higher resolution, as a result of which the orientation of the internal structures will be better visible and can be calculated more accurately. Hence, we expect that, when extended to three-dimensional images, the newly developed method will be able to calculate the spatial orientation of the bone tissue in every voxel of a high-resolution CT scan. These orientations can then be integrated as anisotropic material properties into a heterogeneous FE model. As local orientations are structural parameters on the microscopic level, which will be incorporated into a macroscopic femur model, the development of these types of model might be the beginning of the use of multi-level technology in the in vivo prediction of femoral strength. We hypothesize that the predictive value for fracture risk of the femur will be improved by including the calculated anisotropic material properties.
One contribution of 15 to a Theme Issue ‘The virtual physiological human: tools and applications I’.
- © 2009 The Royal Society