Consideration of multiple load cases is critical in modelling orthotropic bone adaptation in the femur

Functional adaptation of the femur has been investigated in several studies by embedding bone remodelling algorithms in finite element (FE) models, with simplifications often made to the representation of bone’s material symmetry and mechanical environment. An orthotropic strain-driven adaptation algorithm is proposed in order to predict the femur’s volumetric material property distribution and directionality of its internal structures within a continuum. The algorithm was applied to a FE model of the femur, with muscles, ligaments and joints included explicitly. Multiple load cases representing distinct frames of two activities of daily living (walking and stair climbing) were considered. It is hypothesised that low shear moduli occur in areas of bone that are simply loaded and high shear moduli in areas subjected to complex loading conditions. In addition, it is investigated whether material properties of different femoral regions are stimulated by different activities. The loading and boundary conditions were considered to provide a physiological mechanical environment. The resulting volumetric material property distribution and directionalities agreed with ex vivo imaging data for the whole femur. Regions where non-orthogonal trabecular crossing has been documented coincided with higher values of predicted shear moduli. The topological influence of the different activities modelled was analysed. The influence of stair climbing on the properties of the femoral neck region is highlighted. It is recommended that multiple load cases should be considered when modelling bone adaptation. The orthotropic model of the complete femur is released with this study. Electronic supplementary material The online version of this article (doi:10.1007/s10237-015-0740-7) contains supplementary material, which is available to authorized users.


Properties of the muscles and ligaments
shows the values for the properties of the muscles and ligaments included in the model: peak contractile forces ( ), tendon slack length ( ) and reference stiffness values ( ). The stiffness values were distributed by the number of connectors defined for each group, .
26 muscles and 7 ligamentous structures ( Figure 1, in green) were represented as groups of spring elements, in number proportional to their insertion area.
Musculotendon stiffness was calculated as in Phillips (2009) based on the dimensionless force-strain relationship proposed by Zajac (1989) for the tendon using values of maximum isometric force and tendon slack length taken from the literature (Delp 1990).
The stiffness values for the ligaments included in the model were extracted from publically available experimental measurements and validated models for the cruciate and collateral knee ligaments (Butler 1989, Li et al. 1999, Mesfar and Shirazi-Adl 2006, the patellar ligament (Stäubli et al. 1996) and the iliotibial tract (Merican and Amis 2009).

Imaging Data
The results produced by the orthotropic adaptation process were compared with ex vivo imaging data of two femur specimens.

Bone density calculation
The predicted density distributions for the orthotropic models based on a single load case and multiple load cases were plotted and compared with a CT scan slice of the complete femur. From the FE simulation results, a mean Young's modulus, mean , was calculated for each element as the average of the three orthotropic Young's moduli (Equation S1).
Simultaneously, the mean shear modulus, , was calculated as the average of the three orthotropic shear moduli (Equation S2).
The density for each element was calculated using a modulus-density relationship measured for trabecular bone in the femoral neck (Morgan et al. 2003 where rep was taken as the maximum value between mean and the equivalent value associated with mean for each element, with ν taken as 0.3 (Equation S4), in order to account for areas where the shear modulus is predominant over the Young's modulus.
The conversion between orthotropic properties and density was adapted from an experimentally derived curve relating mass and elastic modulus for the femoral neck (Morgan et al. 2003). The choice of this equation is discussed in Geraldes (2013).
Other studies have used similar relationships to convert the orthotropic Young's moduli to density (Miller et al. 2002).
Since shear modulus adaptation was considered, it is possible for certain elements to converge with higher values of shear moduli than those that might be expected when basing their adaptation on the directional Young's moduli. This was taken into account in the attempt to calculate a representative value of Young's modulus, E rep , from which the equivalent density was calculated. The criterion presented was selected as the best representation of the element's Young's modulus after a sensitivity study involving different criteria was performed.
The use of a modulus-density relationship to produce the density plots is limiting.
However, since this relationship is used for both the single and multiple load case model results, it is reasonable to believe that the relative comparison between them holds. More rigorous micromechanics approaches have been put forward by Hellmich et al. (2008) and Blanchard et al. (2013) and would allow for a more accurate calculation of the density distribution in the femur. Alternatively, an extension of the method proposed by Fernandes et al. (1999) for calculating relative bone density from a unit cell with material volume in the three orthotropic directions proportional to the directional Young's moduli could be a more physiological approach and is presented in the next section.

Alternative approach for density calculation
An alternative approach to the use of an empirical relationship was also considered in deriving bone density from the orthotropic elastic constants, using a similar method to that adopted by Fernandes et al. (1999). Considering a cube of unit volume as shown in Figure S2, the relative density (ratio of solid volume to total volume) can be expressed based on the dimensions of each of the internal cuboids representing the solid phase, assumed to be isotropic. where p min is the minimum bone porosity, taken as 0.02 (Cooper et al. 2007) giving a t max of 0.7286. Each value of t ij can be calculated as: where E max and G max are the maximum permissible values of Young's moduli and shear moduli, respectively. The relative density based on the elastic constants can be calculated as: where ρ ̅ is the relative density. In order to compare the density values against those given by the empirical formula a tissue density of 2.11g/cm 3 was used (Morgan et al. 2003). Figure S3 shows the densities calculated for each element using both approaches. Figure S3: Comparison of the predicted densities using the cuboid and empirical approaches.
The R 2 value for a slope of 1 passing through the origin was 0.9512 with a root mean squared error (RMSE) of 0.1286. It can be seen that the cuboid approach tends to predicted lower densities in trabecular bone and higher densities in cortical bone than the adopted empirical formula. It is reasonable to conclude that while specific values of predicted density are sensitive to the chosen approach, differences in density distributions between the two loading scenarios will be similar.

Bone density comparison for the femur
A previous study by authors compared isotropic vs. orthotropic adaptation under a single load case (Geraldes and Phillips 2014) showing that the orthotropic assumption improved the prediction of bone density distribution when compared with the more commonly used isotropic approach. Inclusion of multiple load cases for a variety of frequent daily activities was suggested in order to improve the distribution of the material properties in the trabecular bone regions of the distal part of the femur and across the femoral head. Such quantitative comparison was beyond the scope of this manuscript but, for clarity purposes, we decided to include a comparison of grayscale profiles across 10 mm thick 2D slices of the proximal femur, cortical shaft and distal femur regions ( Figure S4 and Table S2).  The multiple load case model produced grayscale distributions with larger Pearson's product moment correlation coefficients than the single load case model, particularly in the distal region of the femur, showing better predictions in the spatial distribution of bone density, as expected.

The importance of a physiological loading environment of the femur
It is accepted that in order to model the behaviour of bone remodelling with FE models, the driving stimulus of the adaptation process needs to be physiologically meaningful (Bitsakos et al. 2005). This involves careful selection of the loading and boundary conditions applied and has been thoroughly explored in literature. Initial models of the femur were fixed at the distal end and with the hip joint reaction force and muscle forces included as point loads (Taylor et al. 1996). Lengsfeld et al. (1996) showed that the femoral strain pattern and principal stress orientation is sensitive to the resultant joint contact forces and muscle forces at the hip joint, particularly within the coronal plane. Duda et al. (1998) further demonstrated the importance of considering a fully balanced loading configuration in order to reproduce a physiological strain distribution, rather than other loading configurations. Polgar et al. (2003) found that when muscles with large attachment areas were It is clear from literature that boundary and loading conditions play an important role in the resulting predicted stress and strain distributions in the femur, with an increase in their physiological significance associated with more physiological results.
Therefore, the focus of this work was directed to the effect of multiple load cases in the predicted bone structure. The improvement to the adaptation results of physiological loading and boundary conditions when compared to simplified models is discussed in detail in Geraldes (2013). The authors therefore recommend the use of a balanced model to allow for the application of the adaptation process for the complete femur, without artefacts induced by non-physiological boundary conditions. This allows for reporting results for all regions of the femur in different anatomical planes, unlike other bone adaptation models. A similar approach that combined forces obtained from musculoskeletal analysis with finite element models has been used to expand a homogenization method to report density values for the complete scapula with reasonable agreement with CT scan data (Quental et al. 2014), study the accuracy of FE modelling of the scapula (Campoli et al. 2014), and a structural mesoscale model of the femur (Phillips et al. 2015). These studies further confirm the importance of careful physiological representation of joint loading and muscle forces in order to be able to report results for the complete bone geometry.