The effect of boundary and loading conditions on patient classification using finite element predicted risk of fracture

Background: Osteoporotic proximal femoral fractures associated to falls are a major health burden in the ageing society. Recently, bone strength estimated by finite element models emerged as a feasible alternative to areal bone mineral density as a predictor of fracture risk. However, previous studies showed that the accuracy of patients' classification under their risk of fracture using finite element strength when simulating posterolateral falls is only marginally better than that of areal bone mineral density. Patients tend to fall in various directions: since the predicted strength is sensitive to the fall direction, a prediction based on certain fall directions might not be fully representative of the physical event. Hence, side fall boundary conditions may not be completely representing the physical event. Methods: The effect of different side fall boundary and loading conditions on a retrospective cohort of 98 postmenopausal women was evaluated to test models' ability to discriminate fracture and control cases. Three different boundary conditions (Linear, Multi-point constraints and Contact model) were investigated under various anterolateral and posterolateral falls. Findings: The stratification power estimated by the area under the receiver operating characteristic curve was highest for Contact model (0.82), followed by Multi-point constraints and Linear models with 0.80. Both Contact and MPC models predicted high strains in various locations of the proximal femur including the greater tro- chanter, which has rarely reported previously. Interpretation: A full range of fall directions and less restrictive displacement constraints can improve the finite element strength ability to classify patients under their risk of fracture.


Introduction
Fragility hip fractures are a major public health problem frequently causing permanent disability among elderly people with a substantial economic burden to society. In the UK, fragility hip fractures cost NHS approximately £1.1 billion (Leal et al., 2015). Consequently, there is a need for improved prognosis, i.e. to more accurately recognise patients at risk of hip fracture. Bone strength is usually considered a good proxy for the risk of fracture (Riggs et al., 1982). In the current standard of care, dual X-ray absorptiometry (DXA)-based areal bone mineral density (aBMD), or its T-score, are frequently used to measure bone strength (World Health Organization, 1994). The most investigated alternative method is the quantitative-computed-tomography (QCT)-based subject-specific finite element modelling (FE) (Viceconti et al., 2018).
FE models can predict the mechanical strength of a femur as measured experimentally on a cadaveric bone under simulated side-fall conditions with excellent accuracy (Grassi et al., 2012;Pottecher et al., 2016;Schileo et al., 2014). However, the stratification power (the ability to separate fractured cases from controls) appears only slightly better than the aBMD (Enns-bray et al., 2019;Keyak et al., 2011;Kopperdahl et al., 2014;Nishiyama et al., 2014;Orwoll et al., 2009). Similarly, in a retrospective study carried out at our Institute on the same cohort used in this study, the stratification power of femur strength under side-fall conditions was 79% compared with 75% for aBMD (Qasim et al., 2016). While even small increases of predictive https://doi.org/10.1016/j.clinbiomech.2019.06.004 Received 15 March 2019; Accepted 4 June 2019 accuracy may have a significant clinical impact, additional the cost and complexity of obtaining the FE strength predictor when compared to aBMD might not justify its wide scale adoption. Given the % standard error of the estimate affecting the strength prediction when compared to cadaver experimental test is 15-16% (Pottecher et al., 2016;Dragomir-Daescu et al., 2011;Schileo et al., 2014;Dall'ara et al., 2013), it is reasonable to expect some room for improvement, that justify further investigation.
The relatively reduced stratification accuracy of femur strength when compared to its predictive accuracy against cadaver experiments could be attributed to number of modelling factors: among the others CT imaging resolution, the assumption of isotropic or anisotropic material, failure criterion, modelling of the cortical bone, considering or neglecting muscle forces and soft tissue and, boundary and loading conditions (displacement constraints and loading directions). In this study we focus on the boundary and loading conditions. The displacement constraints for side-fall simulation varied among previous studies. The trochanteric aspect of the proximal femur was usually constrained directly against the loading direction (Bessho et al., 2004;Dall'ara et al., 2013;de Bakker et al., 2009;Haider et al., 2013;Hambli et al., 2013;Keyak and Rossi, 2000;Koivumäki et al., 2012;Nishiyama et al., 2013). Other FE models included a frictionless contact at the impact location, which showed a promising improvement in the FE prediction (Ariza et al., 2015;Rossman et al., 2015). On the other hand, the methods used to constrain the distal end vary from either directly fully constrained (Dall'ara et al., 2013;de Bakker et al., 2009;Haider et al., 2013;Hambli et al., 2013;Keyak and Rossi, 2000;Koivumäki et al., 2012) or allowing rotation using a hinge node (Ariza et al., 2015;Dragomir-Daescu et al., 2011;Nishiyama et al., 2013;Rossman et al., 2015). Rossman et al. (2015) investigated the effect of various displacement constraints on the prediction of the femur stiffness. They stated that displacement constraints have a strong effect on the predicted femur stiffness. Furthermore, they reported that constraining the trochanteric aspect eliminates the rotation about the axis transverse to the load, leading to an overestimation of femur stiffness. However, the effect of displacement constraints in term of the stratification power (fracture and non-fracture cases) has not been investigated.
The use of multiple loading conditions has been reported to be beneficial (Falcinelli et al., 2014) to identify bone weakness by providing a wider spectrum of possible loading cases, and therefore to more accurately identify the weakest orientation in the event of a fall.
Posterolateral falls reported to be associate with the weakest structural conditions of the femur (Ford et al., 1996;Majumder et al., 2009;Pinilla et al., 1996), more than it is with anterolateral falls. Thus, most previous studies considered posterolateral falls only while investigating the stratification power of FE model (Falcinelli et al., 2014;Qasim et al., 2016), while very few studies considered both posterolateral and anterolateral falls (Nishiyama et al., 2014). Therefore, it is important to understand whether considering anterolateral falls moreover to posterolateral falls can improve the ability of the FE model to stratify patient under their risk of fracture.
The aim of this study is therefore, to investigate if the ability of the FE predicted strength (improved from the model reported in Qasim et al., 2016) in classifying fracture and non-fractured cases could be improved by: (1) the introduction of a non-linear large-sliding contact boundary condition, or a kinematically equivalent linearized version of such boundary conditions, and (2) considering both anterolateral and posterolateral falls.

Retrospective cohort
A retrospective cohort of 98 postmenopausal women was used in this study. The cohort was divided into two groups: a fracture group and a control group. The fracture group consisted of 49 women who had been diagnosed with low energy trauma fractures in the proximal femur. The control group consisted of 49 women who were pair-matched for age, height, and weight ( Table 1). Details of the cohort are reported previously (Yang et al., 2014).

Finite element modelling
The proximal femur of each patient was segmented using ITK-Snap 2.0.0 (University of Pennsylvania). The segmented bone was meshed with 10-node tetrahedral elements with an average element size of 3 mm (ICEM CFD 14.0, ANSYS Inc., PA, USA) (Helgason et al., 2008). However, for Contact model, mesh refinement was checked based on post-hoc verification of the strain field smoothness. This was monitored in term of strain energy error (ANSYS, Inc. Theory Reference), which in no case exceeded 8% within the volume of interest. For the model with the highest post-hoc error, a mesh convergence study was conducted, which confirmed that the values of first and third principal strains changed by only 4% and 1%, respectively, within the volume of interest between each mesh refinement (3 mm) and the finest mesh (1.5 mm). Material properties were mapped from the CT scans using Bonemat v3.0 (Schileo et al., 2008a). The anatomical coordinate system was defined as described in Qasim et al. (2016). Briefly, the origin is located at the centre of the femoral head, positive x pointed towards the distal end, positive y pointed laterally (towards the greater trochanter), and positive z pointed posteriorly.
Three different boundary conditions were used to represent the sideways fall configuration (Fig. 1): a) Linear model: the most lateral node on the greater trochanter was fixed (in the direction of loading), while being able to translate in the other two directions, simulating a non-friction slider. The distal end of the proximal femur was fully constrained in all directions ( Fig. 1a). This boundary condition has been used frequently in previous studies (Falcinelli et al., 2014;Qasim et al., 2016). b) Multi-point constraints model (MPC): this model employed the same boundary condition (non-friction slider) applied to the greater trochanter in the Linear model, with a relaxed constraint at the distal end. To achieve this, a pilot node was added to the finite element model representing the centre of the knee joint of each patient. The coordinates of the knee joint centre were obtained from the full femur anatomy, estimated using rigid-body registration and statistical shape model-guided fit (Qasim et al., 2016). A multi-point constraints (MPC) method was used to connect the distal end of the proximal femur with the centre of the knee joint. The linear constrained equation related all the degrees of freedom (DOFs; translation and rotation in x, y, and z directions) of nodes at the distal end to the pilot node to simulate a rotational hinge around the knee joint. The pilot node could rotate around the axis transverse to the applied load while all other DOFs were fixed (Fig. 1b). This was intended to mimic the knee joint motion during the fall. c) Contact model: this model employed the same boundary condition at the distal end as in the MPC model, but the proximal constraint was approximated with contact. To achieve a minimal constraint at the greater trochanter, a full non-linear surface-to-surface contact using augmented-Lagrangian algorithm was employed, where the surface area of the greater trochanter was in contact with an  (7) Z. Altai, et al. Clinical Biomechanics 68 (2019) 137-143 Fig. 1. Schematic of the FE model with three different displacement constraints: (a) Linear model with a non-friction slider on the greater trochanter which was free to translate in x and y directions; the distal end was completely fixed. (b) MPC model with non-friction slider on the greater trochanter which was free to translate in x and y directions while the distal end connected to a pilot node by a multi-point constraints (MPC) method; the pilot node was free to rotate around z axis. (c) Contact model with nonlinear surface-to-surface contact between the greater trochanter and an infinitely rigid plate representing the ground. MPC elements were used to connect the distal end of the femur and the pilot node, which was free to rotate around the z axis.

Fig. 2.
Multiple loading conditions (through point force applied to the centre of the femoral head) were used to represent various sideways fall configurations (28LCs), by varying the force direction from (0°t o +30) on the medial-lateral plane and (−30°to +30°) on the anterior-posterior plane. The lateral angles were not considered as this fall direction was not plausible because within this orientation the knee will the ground but not the hip. Region of interest used for analysing the FE strength of the femur is highlighted in blue, while the areas of the BCs that was excluded from the analysis is highlighted in grey. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) infinitely rigid plate (Fig. 1c). This allowed the femur to rotate and slide freely without any imposed displacement constraint on the greater trochanter.
A fitted sphere in the femoral head was used to estimate its centre. A concentrated point load of 1 kN was applied at the centre of the femoral head (an equivalent distributed load was tested with a small difference of only 1.3% found in minimum third principal strain). The 1 kN was an arbitrary load as a load sensitivity study has been performed for the three models and a linear relationship has been found between the applied load and the predicted peak strain; hence the resulting strains can be scaled linearly. All finite element models were solved with ANSYS Mechanical APDL (Ansys Inc., PA, USA).
To cover all possible fall directions, two sets of loading conditions (LCs) were investigated in this study. The first set consists of 28 LCs in which the load angle was incremented by 10 o from 0°to 30°in the frontal plane (adduction) and -30°to +30°in transverse plane (internal-external rotation). This range covers falls on both anterolateral and posterolateral directions. The second set consisted of 16 LCs, which was a subset of the first one to represent posterolateral falls only, with loading angles varied from 0°to 30°in both frontal and transverse planes (see Fig. 2).
Maximum principal strain criterion, which has been validated in vitro and successfully applied in vivo, was used to estimate the femur strength under each loading direction (Falcinelli et al., 2014;Schileo et al., 2008aSchileo et al., , 2008bSchileo et al., , 2014 within the region of interest, excluding areas where the BCs were applied (Fig. 2). The threshold values were 0.73% and 1.04% for tensile and compressive strains, respectively (Bayraktar et al., 2004). The principal strains were averaged at the surface nodes on a circle of 3 mm radius, to follow the continuum hypothesis avoiding local effects of the load (Qasim et al., 2016). Minimum fall strength (MFS) was calculated for each patient by identifying the minimum value of the strength across all simulated loading directions.
Contact model was solved with large deflection option on and the automatic time stepping for the Newton-Raphson scheme available in Ansys for Augmented-Lagrangian contact problems (ANSYS, Inc. Help Guide). All simulations were solved using the preconditioned conjugate gradient-iterative solver (PCG) with a tolerance value of 1E-8. The FE models were run on a high-performance computing cluster machines at the University of Sheffield (Iceberg). The average running time (for one loading direction) for Linear, MPC and Contact models was 3 min, 8 min, and 1 h, respectively.

Statistical analyses
Femur strength predicted by Linear, MPC, and Contact models was compared to investigate the effect of boundary conditions on strength estimation. Using paired t-test, logistic regression models were used to determine the ability of MFS to classify fractures and controls. Area under the receiver operating characteristic (ROC) curve (AUC) was calculated to compare the ability of MFS predicted using different FE models to differentiate between the fracture and control groups using (1) 28 LCs, considering both anterolateral and posterolateral falls, and (2) 16 LCs, considering posterolateral falls only.

Results
MFS predicted by the MPC and Contact models was smaller than the Linear model, with averaged values of 25% and 21%, respectively. MFS was significantly different between control and fracture groups for all FE models (Table 2). Logistic regression showed MFS calculated using Linear, MPC, and Contact models to be significantly associated with the fracture status (p < 0.0001). However, all AUC have been compared statistically (Delong and Carolina, 1988) and no statistically significant difference was found, which indicates that it is not a substantial improvement. Therefore, a larger cohort is needed to confirm our finding or to prove significance.
The highest stratification power was achieved using Contact model considering both anterolateral and posterolateral falls. When the largest loading spectrum of 28 LCs (first set) was considered, AUC was 0.80, 0.80 and 0.82 for Linear, MPC and Contact models compared to 0.75 using aBMD. The stratification power is slightly reduced when anterolateral falls were removed from the simulation, with reported AUC of 0.79 for each of Linear, MPC and Contact models (Fig. 3).
Mean and standard deviation (SD) of femur strength were reported in Table 3 for the three models. The weakest orientation was found when the applied loading direction was at 30°in both frontal and transverse plane. Posterolateral falls were the most critical loading directions for most of the cases under which the MFS was predicted, however, for a considerable number of cases, MFS was predicted under anterolateral falls (Table 4).
The regions of highest strains were consistently predicted by MPC and Contact models across various regions of the proximal femur (Fig. 4). These strains were most frequently located at the femoral neck; although 10 MPC model cases and 13 Contact model cases predicted high strain regions at the greater trochanter. Whereas, Linear model did not predict any high strains at the greater trochanter.

Discussion
The aim of this study was to investigate if the ability of the FE predicted strength in classifying fracture and non-fractured cases could be improved by refining the boundary and loading conditions.
The results show that, in a retrospective cohort of elderly women, considering the full spectrum of fall produced a higher stratification power of the FE model (classification of fracture and non-fracture cases) than considering a range of specific fall directions (posterolateral falls). Previous studies have addressed the importance of considering multiple loading conditions while classifying patients under their risk of fracture based on the estimated FE femur strength (Falcinelli et al., 2014;Qasim et al., 2016). However, for those studies only posterolateral falls were considered. This is consistent with other published results (Ford et al., 1996;Majumder et al., 2009;Pinilla et al., 1996). For most of the cases reported in the current study, posterolateral falls were indeed the critical fall conditions producing minimum load to fail. However, for considerable amount of cases (39 by Linear model, 27 by MPC model, and 23 by Contact model, out of a total of 98 cases), MFS was associated with a fall to the front (anterolateral falls). It should be acknowledged that anterolateral falls may result in the knee or the upper body hitting the ground first, especially in larger angles. However, there may still be cases where an elderly fail to react in time, resulting in the hip falling directly to the ground in this direction (Yang et al., 2018). This is plausible for elderlies with limited mobility. This suggested that it would be desirable to also consider anterolateral falls in FE simulations for a more complete analysis, in order to provide an accurate prediction of the minimum femur strength. In fact, anterolateral falls have been considered in at least one previous study (Nishiyama et al., 2014), which has reported a high predictive accuracy (AUC = 0.94). However, it should be noted that this paper used multivariate regression method (instead of univariate regression), which could lead to a higher AUC.
Regardless of the loading conditions, boundary conditions showed a slight variance on the stratification power for the FE models. The stratification power of the Contact model (AUC = 0.82) was higher than MPC and Linear model (AUC = 0.80). However, that was with the trade-off of computational time. MPC model is considerably less computationally expensive by comparison (8 min versus 1 h for a typical simulation). Although the 2.5% increase in stratification power is subtle, this could be substantial depending on the application, especially in term of cost-effectiveness. A recent review showed that aggregating the results of over 600 tests of human cadaver femurs, the accuracy of CT-based finite element models in predicting the peak force in side fall loading conditions was reliably between 84% and 85% (Viceconti et al., 2018). It is unexpected that the same method might have an accuracy larger than this, when used to stratify fractured and non-fractured cases. The systematic review done by the International Society for Clinical Densitometry in 2015 suggested that the best clinical stratification accuracy that could be achieved was not significantly different from that provided by aBMD (Zysset et al. 2015), that for the current cohort was 75%. Qasim et al. improved this to 79% by addressing some methodological issues (Qasim et al., 2016). In this study we make another important step: using more representative boundary conditions, we achieve a stratification accuracy of 82%. Even if the improvement is small in absolute terms, it brings us much closer to the upper limit of 84-85% which is the accuracy ex vivo.
Femur strength predicted by the MPC and Contact model was   Z. Altai, et al. Clinical Biomechanics 68 (2019) 137-143 substantially (26%, 23% respectively) lower than that predicted by the Linear model, although the three models have the same stratification power (AUC = 0.79) considering posterolateral falls only. This suggests that the inclusion of the MPC method at the distal end of the femur had a similar effect (although not identical) to the contact interaction, and both models were less constrained compared with the Linear model (also used in Qasim et al., 2016). Other studies (Dragomir-Daescu et al., 2011;Keyak, 2001;Rossman et al., 2015) also reported an overestimation of femur strength when direct displacement constraints were applied to the greater trochanter. Among these studies, Rossman et al. (2015) investigated the effect of using contact and MPC methods. However, because they used slightly different boundary criteria (modelled two contact regions at both the femoral head and greater trochanter), their results were not directly comparable with the current study. This highlights the importance of using appropriate boundary conditions (to avoid over-constraining the bone) at the trochanteric and distal aspect of the proximal femur in order to represent a more realistic sideways fall. Maximum principal strains of the MFS loading direction were predicted at different regions of the proximal femur, the locations of which indicate the potential fracture initiation. The majority of patients (72 in the Contact model and 65 in the MPC model) were predicted with maximum compressive strains in the superior region of the neck (Verhulp et al., 2008). However, 13 cases using Contact model and 10 cases using MPC model, the highest strain regions were at the greater trochanter, either the superior intertrochanteric region (compressive strains), or the inferior region (tensile strains) (Fig. 4). The results showed that an improved boundary conditions (either non-linear contact or a kinematically equivalent linearized version using MPC) can successfully predict high strain regions at the greater trochanter, whereas the Linear model failed to predict any. It is well known that the location of the maximum principle strains may provide the approximate site of the crack initiation. Thus, this could explain the lack of per trochanteric fracture predictions in previous FE models. A good agreement between the predicted strain locations by FE using those boundary conditions (Contact and MPC) and the experiments has been reported previously (Ariza et al., 2015;Dragomir-Daescu et al., 2011). Unfortunately, we are unable to comment whether per-trochanteric fracture has been predicted by FE models in these studies. Therefore, to the authors' knowledge, this is the first FE study where maximum strain was predicted in the greater trochanter region, indicating a potential of producing pertrochanteric fracture, which has been observed in the clinic, but not usually predicted by FE models. However, it should be noted that these strains may propagate towards different direction. This could be investigated by considering the crack propagation in the FE model, however, this is beyond the scope of the current study.
One limitation of the Contact model is the assumption of a direct contact between the bone and ground. This does not consider the soft tissue surrounding the greater trochanter, which is likely to have some damping effects and absorb some energy upon impact (Enns-bray et al., 2019;Kerrigan et al., 2003).
Muscle and joint contact forces could be an important factor that may affect the FE model prediction. However, to consider this factor, a personalized musculoskeletal model is needed to calculate these forces and then include them in the FE model. Since the data needed to generate the musculoskeletal model of the subjects of the current study is not available, these forces have not been considered in this study. Future work could investigate the effect of the muscle and joint reaction forces on the FE model prediction.
The size of the cohort of the current study could be one of the limitations. The major influence of falls in the hip may requires much larger groups to generate significant results. However, the cohort of this study is larger than the cohorts used in the previous studies (Falcinelli et al., 2014;Nishiyama et al., 2014).
In conclusion, this study showed that it is important to account for both posterolateral and anterolateral falls for an accurate classification, with posterolateral falls remained as the critical fall conditions for most cases. Furthermore, MPC and Contact models predicted more variations than the most commonly used Linear model in terms of both stratification power and the location of the predicted high strains. That is in agreement with what was reported clinically in terms of fracture types resulted from falls Mokawem et al., 2012). Contact model has higher classification power than MPC model and is more computationally expensive. Therefore, caution is needed in selecting the best method depending on the study objective and resources available. Finally, both Contact and MPC models were able to predict high strain region in the greater trochanter where a fracture might initiate to produce per-trochanteric fracture that has been rarely predicted by the previous FE models. These findings would help to improve hip fracture risk prediction using the FE approach.

Data policy statement
The full set of results can be freely downloaded DOI: 10.15131/shef. data.5901364. Please contact the corresponding author for further information on data access policies.

Declaration of Competing Interest
The authors declare that there is no conflict of interest.