Extended mechanical loads for the analysis of acetabular cages

To analyse the strength and mechanical behaviour of hip implants, it is essential to employ an appropriate loading model. Generating computational models supplemented with muscle forces is a complicated task, especially in the initial phase of implant development. This research aims to expand the possibilities of the simpler acetabular cage model based on joint loads without significantly increasing the demand for computing resources. A Python script covered and grouped the loads from daily activities. The ten calculated major loads were compared with the maximum of the walking and stair climbing loads through the finite element analyses of a custom-made acetabular cage. Sensitivity analyses were performed for the surrounding bones’ elastic modulus and the pelvis boundary conditions. The major loads can geometrically cover the entire load spectrum of daily activities. The effect of many high-magnitude force vectors is uncertain in the approach that uses the most common maximum loads. Using these resultant major loads, a new stress concentration area could be detected on the acetabular cage, besides the stress concentration areas induced by the loads reported in the literature. The qualitative correctness of the results is also supported by a control computed tomography scan: a fracture occurred in an extensive, high-stress zone. The results are not sensitive to changes in the elastic modulus of the surrounding bone and the boundary conditions of the model. The presented load vectors and the algorithm make more extensive static analyses possible with little computational overhead. The proposed method can be used for checking the static strength of similar implants.


Introduction
The hip joint is a spherical joint. The joint loads pass through the centre of the spherical head of the femur or through an artificial ball, which represents the head of the thigh bone when hip arthroplasty is required. When a large acetabular bone defect is presented, patient-specific solutions, such as custom-made implants, are mandatory. In such cases, geometric reconstruction and mechanical analysis are challenging problems (Ahmad and Schwarzkopf 2015;Paprosky et al. 1994).
To develop these types of implants, finite element (FE) analyses are used, among other methods. The FE models around the pelvis have evolved significantly over the years in terms of the material and the boundary conditions-loading models. Bergmann et al. (2001) suggest that for biomechanical and FE analyses, the maximum load during normal walk and the maximum load while walking upstairs and downstairs should be considered.
Many publications have used joint forces for the loading model, or the suggested maximum load of the gait cycle, or sometimes the maximum load from ascending stairs (Costin et al. 2014;Plessers and Mau 2016;Totoribe et al. 2018;Vogel et al. 2020). In one publication, the data were supplemented with loads for stumbling (La Rosa et al. 2016). The maximum normal walking load in a different patientspecific model was investigated in a previous publication (Dóczi et al. 2020), and the results were consistent with the 1 3 observed implant failure. Another publication preferred to use a model standing on one leg (Kawanabe et al. 2011), which others used but with hemipelvic boundary conditions (Moussa et al. 2020). The use of vertical loads of different magnitudes was also discussed by several authors (Du et al 2020;Fu et al. 2018;Iqbal et al 2017;Maslov et al. 2019).
However, if an implant can withstand a large load from a certain direction, it is not evident whether it would be strong enough to resist a smaller load from a completely different direction. Separate loads sampled from certain phases of the walking load are presented in Ma et al. (2013). Typically, five loads are modelled, for which individual analyses were carried out (Hsu et al. 2007;Maslov et al. 2021), but a 32-load model also exists , which is highly detailed, and examines the contact conditions of the acetabulum. One publication examines the importance of different load directions by considering five load cases (the maxima of one-leg stance, normal walking, ascending and descending stairs, and stumbling) (Del-Valle-Mojica et al. 2019). This concept is important for topology-optimized implants, because if some loads were not intended in the optimization task, it could provide poorly optimized results. Iqbal et al. (2019) used multiple load cases for their topology optimizations. Still, these were only the maximum values of the following daily activities: normal walking, standing on a single leg, standing up, sitting down, ascending stairs, and descending stair.
The hypothesis of our study is that these most common loads from the literature (the maximum load of walking and the loads of walking upstairs and downstairs) can only provide a narrow segment of sample for analysing acetabular implants. This will be supported by comparing these load magnitudes and directions with the force vectors from daily activities and with a patient-specific finite element model of a custom-made acetabular cage, where the stress concentration areas will be investigated. Sensitivity analyses will be carried out for the material model and the boundary conditions. Moreover, the von Mises stress field will be compared with the observed implant failure from the patient's control computed tomography (CT) data. The other goal of this study is to propose a heuristic algorithm for calculating the recommended loads for implants from any given dataset, albeit with some restrictions.

Covering vectors
Calculating the mechanical response of a structure when its load vectors have the same initial point in every different load case would be computationally very challenging. A more effective way would be to groups the loads and define the so-called major loads, which may be used for simulations, where there are no other larger loads in their small surrounding areas.
The small surrounding area can be defined as a circular sector in the 2D case; in 3D, it is a spherical sector (see Fig. 1).
The apex angle is user-defined. With a smaller angle, the calculation is more detailed, but this can increase the number of simulations required. If the apex angle is large, the coverage area is larger, and fewer major loads are required. However, this method can produce poor simulation results because the direction difference between the major vector and a common vector can be significant.
Calculating the major loads is similar to the minimal disc covering problem, which is non-polynomial (NP) hard (Das et al. 2011). In our case, the task is to cover 3D vectors with spherical sectors with a given half-apex angle, which the authors selected at 10°. The motivation for this decision will be explained later. In this study, a possible heuristic solution will also be provided.

Heuristic algorithm for calculating the major loads
The dataset of the loads was the HIP98 dataset by Bergmann et al (2001). All the given loads from daily activities (slow, normal, and fast walking, walking upstairs, walking downstairs, standing up, sitting down, standing, knee bending) for an average patient were considered. These data contain 9 × 201 = 1809 3D vectors with their X, Y, and Z components. Our study focused on loads of the acetabular joint in the pelvic coordinate system (Bergmann et al. 2001). The loads have the same initial point, which is their origin. Thus, the load vectors are represented as their terminal points on the figures.
The heuristic algorithm was written in Python and can be divided into two parts (see flowchart in Fig. 2).
The images in Fig. 3 also demonstrate the process, using colour markings similar to Fig. 2.
In the first part, to calculate the directions of the major loads, the data's magnitudes are irrelevant, so the force vectors should be transformed into unit vectors. The next step is to calculate the mean of this unit vector data; then, the algorithm calculates the farthest data point from this global mean vector. This edge vector can be applied to separate the local set of the unit vectors to reduce computing time.
In the following step, the search vector was defined. The search vector is the axis of a cone, which has a userdefined half-apex angle ( = 10°). If this covers a vector, the following equation is satisfied [Eq. (1)]: where u is the vector be covered, c is the rotation axis of the cone as a vector, and is the cone's half-apex angle.
The search vector can also move around the edge vector, but it always contains the edge vector. The increment of the movement is user controlled. The increments define the number of different search positions (the search resolution) and the details of the calculation, including running time.
The next step is to choose an appropriate position from among these various positions: this will be defined by the search vector, and it will be the direction of a major load. Fewer major loads are produced if the algorithm is set to prefer to cover the unit vectors closer to the edge vector. This can be implemented as described below.
The distances from the edge vector to every local unit vector are calculated. Next, all the distances are normalized by the longest distance (d norm ). After that, the following penalty function is used [Eq. (2)]: Hence, all unit vectors have a certain value, (v). Each unit vector has a minimum value (v min ), and if the exponent (p) is greater than 0, the closer the unit vector is to the edge vector, the higher the v value will be.
The sum of these values is calculated in every search vector position. Where it is the maximum, that will be the direction of a major load, and the local unit vectors covered by this major load's cone are removed.
The process is then repeated in a loop to calculate the new mean values of the remaining unit vectors until all the unit vectors have been covered.
In summary, the user-defined values were the following. It was possible to cover the dataset with ten major loads, choosing the exponent of the penalty p = 6, the half-apex angle = 10°, and the increment of the movement was 1°. Figure 4 shows the chosen parameters. By modifying the exponent of the penalty function, different numbers of major load vectors can be used to cover the dataset. In the case of a half-apex angle of 10°, this number can be reduced to 10 (by setting the exponent of the penalty function to 6), which is comparable to the number of types of loads from daily activities in the data set (10 vs. 9). It was not advisable to run more simulations than this to demonstrate the method.
Selecting a 1° increment of the movement, the running time for calculating the major loads was feasible. Smaller increments did not significantly improve the algorithm's performance in terms of fewer major loads but increased the running time significantly.
The second part of the algorithm calculates the magnitudes of major loads. First, all the major load direction magnitudes are changed to the maximum magnitude from the dataset of the force vector, and these are then reduced step by step as described in the process above ( Fig. 2).

Finite element model of the fixation
To examine the applicability of the major loads, a patientspecific finite element model was developed.
The patient had a large Paprosky 3/b acetabular bone defect. The treatment was a custom-made sheet metal acetabular cage provided by Sződy et al. (2017). The CAD models were created using the patient's anonymized pre-and postoperative CT data. The hemipelvis, the acetabular cage and its screws were segmented with threshold-based and manual segmentation using Slicer 3D. Then, the STL files were exported and repaired with Autodesk Meshmixer.
To produce the solid CAD model of the pelvis, Solid-Works 2018 Scanto3D module was used. Near the bone defect, the cortical and the trabecular bone cannot be segmented easily. Due to the migration of the original acetabular component, a significant change had taken place in the bone's structure. Therefore, the regions close to the acetabular defect were separated from the other parts of the pelvis, and different materials were used for the models, which will be discussed later.
The mid-surface CAD geometry of the sheet metal acetabular cage was also made in SolidWorks 2018.
The screws had a simplified cylindrical geometry, with a nominal diameter of 4.5 mm. The heads of the screws had simplified spherical geometry with dimensions matching the relevant standard (ISO 5835:1991). The liner and the femoral ball were modelled as revolved bodies.
The pre-processing of the finite element model was carried out using HyperMesh 2017.2. Multiple mesh refinement Fig. 2 The flowchart for the calculation of the major loads direction and magnitude was applied to the acetabular cage because this part is the major focus of this study. The maximum von Mises stress was evaluated in a stress concentration area with different element sizes for the normal walking load case. Between 0.3 and 0.7 mm element sizes, the deviation between the results was under 5% compared to the last mesh refinement (Fig. 5). Fig. 3 The endpoints of the force vectors are coloured by their magnitudes, some of them are marked with black arrows (a); transformation to unit vectors (b); calculating the mean and the edge vector, separating the possible coverable unit vectors (c); penalizing the coverable unit vectors depending on their distance from the edge vector (d); the search vector in its initial position (grey) and local best coverage (gold) (e); the subprocess repeated, previously covered unit vectors were removed (f); the direction results (g); the length-adjusted direction results, the major loads (h)  This element size range was used in all stress concentration areas. The final FE mesh data are given in Table 1.
Homogenous, linear elastic, and isotropic material properties were used everywhere, including the bone (Anderson et al. 2005). Near the acetabular defect, a homogenous material was used with an averaged Young's modulus (Ravera et al. 2015). A sensitivity test was performed for this part due to the versatile bone structure.
In the healthy areas of the pelvis, a separate material model was applied. The solid elements represented the spongious part of the pelvic bone, while shell elements represented the cortical layer (Anderson et al. 2005;Plessers and Mau 2016), as shown in Table 2. Table 3 shows the connections between the parts of the model.
There were fixed boundary conditions at the sacroiliac joint and the pubic symphysis (Clarke et al. 2016;Plessers and Mau 2016). A sensitivity analysis was performed on the boundary condition change in the pubic symphysis.
The first load step is a bolt pretension, where every bolt has a 50 N pretension force, modelling the surgical procedure and closing the contacts. Then, in the next load step, as the main load, a force vector is applied to the centre of the femoral head. Each load was applied in a separate static load case after the bolt pretension, as the aim of the study is to investigate the response of the cage to different load directions. The applied loads based on the literature and simulated by the algorithm can be seen in Table 4.
The components of the main loads for FE analyses in the pelvic coordinate system (Bergmann et al. 2001) are given in Table 4. The first five rows are the peak loads from different daily activities from the HIP98 database (slow, normal, and fast walking; walking upstairs and downstairs)   in BW% (percent of body weight). The patient's weight was approximately 80 kg, and the force components were calculated with a 9.81 m/s 2 acceleration of gravity. The next ten loads are the major loads that were given by the algorithm. Small displacement quasi-static nonlinear analyses were performed because the deflections of the model were relatively small compared to its main dimensions, whereas the dynamic effect of the forces was neglected, and the frictional contacts required nonlinear calculation. The FE solver was Optistruct 2017.2.
The Python script used to determine the major loads and the FE calculations were run on a personal computer (Processor: Intel ® Core™ i7-9700 CPU @3.00 GHz; RAM: 64 GB; System type: 64-bits operative system, × 64-based processor).
The Python script for the major loads finished in 32 s using 199.6 MB of memory, while the 1 + 15 = 16 FE (15load case after the pretension step) analyses were completed in 11.75 h.
The geometry and the FE model are depicted in Fig. 6. The acetabular cage is the major focus of this study. The von Mises stress of the cage was examined during the sensitivity analyses in one load case. The preload was the bolt pretension, and subsequently, the maximum amplitude from the normal walking gait cycle was applied to the model as a static load. The sensitivity analyses were performed in two dimensions.
The homogenous bone part was in connection with the acetabular cage. Here, Young's modulus value was changed by ± 30%, so the elastic moduli were 9100 and 4900 MPa.
The boundary condition at the pubic symphysis was changed to a 1 DOF constraint in the X-axis direction.
The patient's postoperative and control CT scans were utilized for the qualitative verification of the results. The expected failure must be in an extensive, high-stress zone, which the FE model must be able to predict. Figure 7 compares the loads in the literature and the results presented by the algorithm. It clearly shows that a significant proportion of the loads from daily activities, including high-magnitude loads, are not taken into consideration. For example, there is a large vector with a 252 BW% magnitude (marked with a red circle), which is from the downstairs activity. In a quite different direction, there is a load with a large BW% magnitude (marked with a blue circle). An approach that uses the peak loads from the daily activities (i.e. nine loads) cannot cover all the force vectors. However, with these ten loads, our algorithm can cover all the force vectors from the daily activities with the spherical sectors. This method requires only one more FE analysis than when the maximum force approach was employed. FE analyses confirmed the fact that this extension influences the stress state of the cage.

Finite element results
Investigating the von Mises stress results and the loads given in the literature (Bergmann et al. 2001;Costin et al. 2014;Plessers and Mau 2016;Totoribe et al. 2018;Vogel et al 2020;), they are found to be very similar. Four main stress concentration areas can be detected (see Fig. 8 Viewing the results induced by the major loads, a new stress concentration area, which the loads from the research of the literature cannot provide, can be detected, in addition to all the previously mentioned stress concentration areas (see Fig. 9).

Sensitivity analyses
Overall, the examined stresses showed low sensitivity to the changes mentioned in Table 5. The acetabular cage's nodal maximum von Mises stresses were compared to the original model. The change in the bone's elastic modulus had no significant effect on the von Mises stresses of the acetabular cage because it is much stiffer than the bone and the supporting areas around it.
Different boundary conditions did not give significantly different results either; in both cases, almost the same kind of displacements was produced.

Comparison of postoperative and control CT data
The comparison of the STL meshes from the postoperative and control CT data can be seen in Fig. 10. The primary failure of this implant was a fracture in one of the supports, where two stress concentration areas (zone 3 and zone 4) were in a wide section with a high von Mises stress. A fracture was detected in the control CT in this part of the cage. This failure was indicated by stress responses both to loads commonly found in the literature (the first five loads from Table 4) and stress responses to major loads modelled by our method. The authors want to emphasize that in this case, this type of qualitative verification is suitable for showing the consistency of the results from the loads from the literature with the results provided by the new method we propose.

Limitations of the results
The major vector results are presented only for the 10° halfapex angle, and it is not proven that there is no coverage with less major load vector with this angle. Another limitation is that only the HIP98 database was covered because this contained the forces in the pelvic coordinate system. The FE results only focus on the von Mises stresses of the cage. Moreover, the calculated von Mises stress results are greater than the yield stress of the unformed sheet metal. This is the consequence of using a linear elastic material model for the acetabular cage. The cage was a highly deformed sheet metal part. Due to the extensive cold work, the yield stress in the highly formed areas could be much higher than in the original. There was no annealing process after the forming so that residual stress could remain in the cage. However, modelling this phenomenon is challenging and has no other important effect on the result during the qualitative comparison of the stress concentration areas.

Limitations of the study
The study and the methodology used also have their limitations, and some criticisms may be made of them, which we tried to reduce with countermeasures but which cannot be eliminated entirely.
The muscle loads cannot be included because the algorithm only deals with vectors with the same initial point. The major loads provided do not belong to any measured data point, and the other joint load is unknown. However, for implant development, it can produce viable results, and for topology optimization, simplified models are recommended.
The 10° half-apex angle selected was the author's choice. A smaller angle could provide more detailed results, but it would mean more major loads and more simulations.
The heuristic algorithm cannot provide a globally optimal solution. It is worth running it multiple times with several parameters specified by the user. Nevertheless, it is guaranteed to provide a solution, and if the number of major loads is acceptable, it can be used immediately. Only static analyses can be performed with these forces.
Regarding the FE model, only a hemipelvis model without pelvic discontinuity is suitable for the analyses; the load on the other leg and the muscle forces and ligaments were replaced by boundary conditions. Accordingly, a sensitivity test was carried out for the boundary condition, but no significant change in the stress state of the cage was observed.
The starting point of the geometric reconstruction was a clinical CT scan with scattered X-ray photons. Segmentation, STL mesh formation, conversion into a CAD model, and finite element meshing again involve the possibility of geometric inaccuracies, which may arise due to the differences in the position and orientation of the acetabular cage in the model compared to in reality. To eliminate this potential error, we used screw pretensioning, but the pretension forces were also estimated, and the geometry of the screws was simplified. There may also be inaccuracies in the friction coefficients from the literature.
Based on the simplified calculation, a homogeneous, linear elastic, isotropic material model was used throughout. Where possible, separated material properties were used when specifying the material model of the pelvis, but the exact cortical thickness was only an estimate supported by the literature. Due to the limitations of the clinical CT Fig. 8 The stress concentration areas and the von Mises stress distributions for different peak loads produced by daily activities. Slow walking (a); normal walking (b); fast walking (c); ascending stairs (d); descending stairs (e) data, it was impossible to separate the cortical and the spongiosis parts in the areas around the cage. To compensate for this uncertainty, a sensitivity test was performed on the modulus of elasticity of the parts here, which had no significant effect on the stress state of the acetabular cage either.
Elsewhere, other authors have been able to use the inclusion of muscular forces successfully (Clarke et al. 2016;Fig. 9 The von Mises stress distributions for the ten major load cases (abbreviated ML#)  Ravera et al. 2015;Zaharie and Phillips 2018) and make sensitivity analyses to investigate the changes in the results. It is important to note that these publications primarily focused on the pelvis bone structure and the stresses around it. Considering all these muscular forces for simulations entails difficult and computationally expensive tasks, especially in the early phase of implant development. Overall, these loads were able to predict the failure's location, as was also seen from comparing the CT data. Hence, these results can be deemed to have been qualitatively verified. In this case, in terms of failure, the stress patterns given for the loads by the literature and those given by our proposed method match. Still, it is evident that further tests are needed to show the new stress concentration area, because the small sample size and the lack of a diverse range of implant designs and materials limit the study. Expanding the study to include a greater variety of implants would allow for a more comprehensive understanding of the stress concentration areas and potential failure locations. To achieve experimental quantitative verification, additional data are needed. Strain gauge stamp measurements could also confirm the stress concentration areas but that would require in vitro experiments.

Conclusion
For the stress analyses of implants, the most common loads given in the literature can only provide a tight inspection segment with the FE method. During daily activity, large magnitude force vectors occur in quite different directions than the maximum peak loads of walking and stair climbing. Stress concentration areas may remain hidden if these loads alone are used for the investigation of implant stresses and the development of implants based thereon. This paper proposes a new approach and method for the loading model of FE simulations of hip joint implants. This approach is mainly based on the statement that if many force vectors have the same initial point, the loads can be grouped by spherical sectors with small half-apex angles. The axes of these spherical sectors can be the direction of major loads. The magnitudes are the maximum magnitude of the covered vectors. Using the published heuristical method for an NPhard covering problem, the major loads can be calculated quickly with quite good efficiency.
The other useful result from using the major loads is that it made it possible to select a zone for development because all the major loads predicted this zone as a stress concentration area (zone 4).
The methodology proposed here can be applied in concept formation, using simple models, in preliminary static checks. Subsequently, detailed calculations can be carried out for the promising product variants with more advanced models.
The algorithm can also cover other vectors with the same initial point, for example, the force vectors of a femoral implant.
In the algorithm, adding step zero could remove the force vectors that do not have a larger magnitude than a certain value. Suppose there is any previously given load with a specified half-apex angle (such as the stumbling with its great magnitude). In that case, the inner force vectors can also be removed, and the remaining vectors could be covered. This method and its results can be used for further implant analyses, even with patient-specific joint loads from inverse dynamics muscle models.
The FE model could be improved, for example, by including the applicable ligaments around the hemipelvis, a CTbased material assignment for the bone, and the modelling of the cage's residual stresses during manufacturing, nonlinear material models and the detailed modelling of bone-screw connections, even with submodels. In addition, the authors will use the approach of major loads to study the topology optimization capabilities of acetabular cages.