In situ parameter identification of optimal density-elastic modulus relationships in subject-specific finite element models of the proximal femur

Abstract

Quantitative computed tomography based finite element analysis of the femur is currently being investigated as a method for non-invasive stiffness and strength predictions of the proximal femur. The specific objective of this study was to determine better conversion relationships from QCT-derived bone density to elastic modulus, in order to achieve accurate predictions of the overall femoral stiffness in a fall-on-the-hip loading configuration. Twenty-two femurs were scanned, segmented and meshed for finite element analysis. The elastic moduli of the elements were assigned according to the average density in the element. The femurs were then tested to fracture and force-displacement data were collected to calculate femoral stiffness. Using a training set of nine femurs, finite element analyses were performed and the parameters of the density-elastic modulus relationship were iteratively adjusted to obtain optimal stiffness predictions in a least-squares sense. The results were then validated on the remaining 13 femurs. Our novel procedure resulted in parameter identification of both power and sigmoid functions for density-elastic modulus conversion for this specific loading scenario. Our in situ estimated power law achieved improved predictions compared to published power laws, and the sigmoid function yielded even smaller prediction errors. In the future, these results will be used to further improve the femoral strength predictions of our finite element models.

Figure 1

Photograph of the experimental setup. The distal end of the femur was potted in a block of dental cement and clamped in a fixture. The fixture was initially placed at an angle of 10° with respect to the y-axis, and could rotate about the x-axis. The trochanter was fixed to a cup with dental cement, which was connected to a load cell fixed to a metal frame. The femoral head was placed underneath an aluminum cup with curved surface. This cup was connected to a load cell which could move in the x and y directions using low-friction linear bearings.

Figure 2

Representative force-displacement curve of a femur with normal aBMD. Trochanter load is plotted versus femoral head displacement. The slope of the linear elastic part of the curve is measured as the overall femoral stiffness.

Figure 3

Surface and cut planes of finite element meshes used in the convergence analysis. Three uniform meshes with different maximum element edge lengths were used, and a smart mesh was developed with finer elements in the head/neck/trochanter region and coarser elements in the shaft. Moreover, the smart mesh contained smaller elements at the cortical shell and larger elements in the trabecular region.

Figure 4

Boundary conditions were applied to the finite element model to mimic a fall on the side. A vertical force was applied to the femoral head and vertical displacements were set to zero at the greater trochanter. The distal end of the femur was connected with beams to a rotation point. At the rotation point, only rotation about the x-axis was allowed; the other five degrees of freedom were set to zero.

Figure 5

Results of the mesh convergence study. Mean ± SD of the errors relative to the 1.5 mm meshes are plotted for the nine femurs in the training group for the smart mesh, 2.5 mm mesh, and the 5.0 mm mesh.

Figure 6

Measured vs. FEA predicted overall femoral stiffness using density-elastic modulus relationships previously published by (A) Morgan et al. (2003), (B) Keller (1994) and (C) Carter and Hayes (1977). For consistency, we differentiated between the estimation (solid circles) and validation (open squares) groups used in the parameter estimation procedure developed in this paper.

Figure 7

Convergence of the error between FEA prediction and experimental measurements (Eq. 6) during the parameter estimation procedure for the power law. After 62 iterations, the change in the objective function was below 10^-9.

Figure 8

Measured vs. predicted stiffness for the estimation and the validation groups, based on the optimized power law (A), and sigmoid functions with E_max = 15 GPa (B), 20 GPa (C) and 25 GPa (D).

Figure 9

Elastic (Young’s) modulus-density relationships for (A) three previously published power laws, (B) the estimated power law and sigmoid function with E_max = 25 GPa, and (C) all three estimated sigmoid functions.