Computationally efﬁcient modelling of hip replacement separation due to small mismatches in component centres of rotation Journal of Biomechanics

Patient imaging and explant analysis has shown evidence of edge loading of hard-on-hard hip replace- ments in vivo. Experimental hip simulator testing under edge loading conditions has produced increased, clinically-relevant, wear rates for hard-on-hard bearings when compared to concentric conditions. Such testing, however, is time consuming and costly. A quick running computational edge loading model (Python Edge Loading (PyEL) - quasi-static, rigid, frictionless), capable of considering realistic bearing geometries, was developed. The aim of this study was to produce predictions of separation within the typical experimental measurement error of (cid:1) 0.5 mm. The model was veriﬁed and validated against com- parable ﬁnite element (FE) models (including inertia and friction) and pre-existing experimental test data for 56 cases, covering a variety of simulated cup orientations, positions, tissue tensions, and loading envi- ronments. The PyEL model agreed well with both the more complex computational modelling and experimental results. From comparison with the FE models, the assumption of no inertia had little effect on the maximum separation prediction. With high contact force cases, the assumption of no friction had a larger effect (up to (cid:1) 5% error). The PyEL model was able to predict the experimental maximum separations within (cid:1) 0.3 mm. It could therefore be used to optimise an experimental test plan and efﬁciently investigate a much wider range of scenarios and variables. It could also help explain trends and damage modes seen in experimental testing through identifying the contact locations on the liner that are not easily measured experimentally.


Introduction
Separation of the femoral head and acetabular cup in total hip replacements (THRs) has been associated with increased wear and device damage (Kovochich et al., 2018;Manaka et al., 2004;Sanders et al., 2012). Experimental THR simulator testing can include a transverse force to cause this separation, representing soft tissue tension from slightly misaligned (mismatched) components, as well as a superoinferior (SI) force (ISO, 2018). The intention of these standardised edge loading tests is to capture an important aspect of in vivo contact area and pressure conditions, providing a method for comparison of new devices with well performing predicate devices.
Different surgical approaches have been linked clinically to changes in separation (Glaser et al., 2008a). Post-operative joint positioning alters soft tissue tensions and moment arms (Asayama et al., 2005;Cantin et al., 2015), which have been linked to changes in separation through computational models (LaCour, 2017). Understanding the connection between loading environment, separation, and the mechanisms of material damage under edge loading conditions, is the subject of on-going research. For mismatch driven separation, the transverse force dominates when the components are under lower SI loads, creating head-liner separation and moving the contact area towards the liner edge. The degree of head-liner separation, and the subsequent damage modes, are dependent on various interacting factors, such as the cup orientation, translational mismatch between the components, tension force, and loading/activity. The large number of test cases required to evaluate these effects can currently be analysed using an experimental biomechanical procedure (0.5 h per case). Longer term, more time consuming, wear and durability testing associate the mechanical scenario with damage of the device . The time required for these studies, https://doi.org/10.1016/j.jbiomech.2019.07.040 0021-9290/Ó 2019 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). however, limit the range and resolution of the test cases that can be considered and these trends need to be re-evaluated for any significant design changes (for example lipped liners, medialised liners, altered bearing sizes). The increasing need for more realistic and stratified implant testing conditions (Yamamoto et al., 2005), reflected in new testing standards (ISO, 2018), will progressively escalate the pre-clinical experimental testing work load to impractical levels. The provision of efficient simple computational methods would enable many more situations to be assessed and complex computational and experimental techniques to be focused on high risk conditions.
Computational modelling is an important tool in understanding the source of different types of material damage seen in THR devices. Simple computational approaches can capture information on the relationship between the force environment and contact position. More sophisticated computational approaches, one example being finite element methods, can be used to predict more detailed contact pressures and areas (Askari and Andersen, 2018;Donaldson et al., 2015;Genda et al., 2001;Hast et al., 2019), and possibly material strain, stress, and wear based damage metrics. They are relatively time consuming and specialised (and hence costly) which limits the number of variable combinations which can be investigated.
To overcome these issues this study aimed to develop and evaluate a new computational implementation of a quasi-static analytical model (Python Edge Loading, PyEL) to predict THR separation behaviour, which quickly provides information on geometryspecific contact position. The model was developed for hard-onhard bearing designs where large deformations and liner geometry changes, due to wear and creep/plastic deformation, are minimised (Hart et al., 2013;Manaka et al., 2004;Williams et al., 2007), and all components were assumed to behave as rigid bodies. The model is designed to replicate the action of current ISO standard separation testing, to enable detailed validation, and to facilitate understanding of the mechanisms connecting device geometry, orientation, and load environment to damage modes.
Three studies were undertaken to: (1) verify the implementation of the PyEL model via discretisation studies and comparison with equivalent finite element (FE) models; (2) compare the PyEL predictions to component separations that occur experimentally under a range of orientation and force conditions; (3) investigate the time penalties and any benefit of more sophisticated FE models which include dynamic behaviour and friction.

Methods
In all studies a 36 mm ceramic-on-ceramic (CoC) Pinnacle BIO-LOX Ò Delta bearing (DePuy Synthes, Leeds, UK) (radial clearance 0.05 mm) was used. The primary output was the maximum separation in the mediolateral (ML) direction between the head and liner centres (linked to wear rates and which geometric regions of the liner are contacted Sariali et al., 2012)).
The experimental methodology, upon which the computational models were based, was reported in a separate study (Ali et al., 2017). Briefly, the experimental testing was performed on a six station ProSim Electromechanical hip simulator (EM13, Simulation Solutions, Stockport, UK) with three stations used per case. For practicality in terms of simulator design, the joint load was applied in the SI direction instead of 10°medially from the vertical (as is seen clinically). To maintain the correct load vector orientation relative to the liner, an equivalent 'clinical' inclination angle of 45°w as therefore represented by an applied inclination angle of 35°. Cases are described in this paper using 'clinical' liner inclinations and the full case list is provided in Table 1.
A spring was used in the ML axis to generate a transverse separation force (ISO, 2018;Nevelos et al., 2000). The component mismatch was represented by the level of spring compression with a concentric head and liner. The SI loading followed the ISO14242 gait profile (ISO, 2014) at 1 Hz. At peak SI loads, the head and liner come close to alignment, the contact area approaches the liner apex ( Fig. 1, a), and the lateral spring force is at its highest. At lower SI loads, the spring force causes separation and movement of the contact area toward the liner edge ( Fig. 1, b).
A quasi-static, rigid, frictionless, model was developed using Python version 3.5 (Python Software Foundation) to replicate the experimental methodology (PyEL model). This model used static mechanics (Fig. 2) to approximate the relationship between load and displacement (separation). The liner was represented as a three-dimensional point cloud surface created from a tetrahedral finite element mesh generated using Abaqus (v2017, Dassault Systémes, France). This point cloud was provided to the PyEL model as a text file and could be generated by any convenient software, including open source options. The other inputs, which were fixed for a given case, were the cup inclination angle, the spring constant, the magnitude of the spring compression (ML mismatch) and the magnitude of the axial load at each point in the gait cycle. Ultimately, the tool generated a head-liner contact position and a magnitude of ML separation for each point in the gait cycle. These values were calculated in two main stages (Fig. 2). In the first stage, the head-liner contact position and axial force required for static equilibrium were calculated at a series of separation positions. In the second stage, the contact positions through the gait cycle were re-constructed from the axial force value and the relationships found in the first stage.
Three additional rigid-body FE models were generated to compare to the PyEL model; the different modelling approaches allowed investigation of the effects of inertia and friction on separation behaviour.
Parametric tests were undertaken in both the experimental and computational models to understand the separation behaviour under different conditions. The variables were: liner inclination angle; translational mismatch (generating the ML force); and the SI load during swing phase of gait. Cases are listed in Table 1. For modelling methods with a greater analysis time, a smaller subset of cases were analysed (Table 1).
Study 1: The quasi-static rigid Python model (PyEL) was verified in terms of discretisation error and by comparison with an equivalent quasi-static rigid FE model (runtime: 1 h). The FE model was developed using the software package Abaqus (v2017) and meshed using quadratic tetrahedral elements (target element size 1 mm). Penalty-based hard normal head-liner contact was defined. Two cases were used (inclination angles of 45°and 65°), both with a 4 mm mismatch, 100 N/mm spring stiffness, and 70N swing phase load (SwPL). There were two main sources of potential discretisation error in the PyEL model: the resolution of the point cloud used to represent the geometry of the liner and the number of separation locations evaluated. Convergence testing was performed for both of these parameters on a generalised liner design, with the same inputs as above except with anatomical versions of 10°and 30°. This generated four test cases inducing contact across two different liner regions. Point clouds were created from element sizes of 0.27, 0.18, and 0.13 mm. The through-cycle separation behaviour was generated using 500, 1000, and 2000 linearly spaced separation positions. Changes in minimum and maximum predicted separation, and the integration through time of the SI and ML forces during edge loading (representing behaviour throughout the cycle) (O'Dwyer Lancaster-Jones et al., 2017), were compared.
Study 2: PyEL results were compared with experimental results for all parametric cases (two inclination angles, four mismatches, and seven SwPLs) described in Table 1 (56 cases total, average experimental testing time per case: 0.5 h). The experimental testing was performed as a separate study (Ali et al., 2017). Results presented are the mean from three repeats under each test condition (maximum separation standard deviation (SD) range 0.02-0.31 mm).
The aim of Study 2 was to compare identical, force-matched, scenarios and characterise the predictive ability of PyEL. PyEL cases were generated to match the mean measured experimental SwPL for each case. Due to simulator bending affecting the rate of spring compression (O'Dwyer Lancaster-Jones, 2017), the maximum experimental ML forces were used to calculate a modified spring constant for the PyEL cases (mean 63 N/mm, range 41-117 N/mm, SD 14 N/mm). This enabled the PyEL model to experience the same range of ML forces as the experimental case, across the same range of bearing positions.
Study 3: Separation was simulated using rigid-body explicit FE models which included (1) dynamic behaviour (inertia) (critically damped, runtime: 10 h) and (2) dynamic behaviour with bearing surface friction (critically damped, l = 0.05 (Brockett et al., 2007), runtime: 10 h). Results were compared to PyEL to assess the value of including additional physics. The femoral head mass was increased to 2.5 kg to include the mass of the attached experimental fixture. The mesh matched that of the quasi-static FE model. Mesh sensitivity testing gave a difference in maximum separation of <0.025 mm when comparing target element sizes of 0.75 mm and 1.25 mm. Two load cycles were simulated and the second cycle, where the response had settled, was used for data collection.

Results
The associated data are openly available from the University of Leeds data Repository (Etchels and Jones, 2019).

Study 1: Verification of PyEL model via discretisation sensitivity and comparison with quasi-static rigid FE
Across four different liner orientation cases (inclination 45°/65°, anatomic version 10°/30°) reducing the PyEL point cloud element size from 0.18 to 0.13 mm changed the calculated minimum and maximum separation results by less than 0.04 mm. Typical maximum separation results were in the range 0.5 -4 mm. Compared to an element size of 0.13 mm, 0.18 mm provided reasonable convergence for the integral of the ML and SI forces during edge loading with differences <3.5%.
Changing the number of contact position points evaluated by the PyEL model for the four cases from 500 to 1000 points caused changes to the minimum and maximum separation of less than 0.004 mm. For the integration of the ML and SI forces during edge loading, differences of <1.0% were seen.
Cases using an element size of 0.18 mm and 500 position points required a typical run time of <10 s.
A comparison between the PyEL model and a quasi-static, rigid, frictionless FE model is given in Fig. 3. At both inclination angles, the PyEL model predicted the FE maximum separation well (underestimate of 0.02 mm). There were two additional differences. The  1. Demand load profile applied to all models and the experimental test specimens. Describes how (a) a large axial load is able to force the components into alignment despite the mediolateral spring but that (b) a small axial load will allow for separation. first occurred at time t 1 , Fig. 3A. The PyEL model produced multiple separation positions for a single point in time when transitioning between different geometric liner regions. This was due to the relatively large change in surface normal from one region to the next. Positions with slightly less separation required the same axial force to achieve equilibrium as a position with slightly more separation. Points can therefore be seen with a separation 0 mm occurring at the same time as points with a separation of 1 mm. The second occurred at time t 2 , Fig. 3A. The PyEL model underestimated the lateral separation. This was due to the FE rigid contact implementation allowing for small nodal penetrations. PyEL used perfectly rigid point contact. Similar results were found at the higher inclination angle (Fig. 3B).

Study 2: Characterisation of PyEL predictive ability via comparison with experimental data
Dependent on the applied mismatch, the lowest separation value that the experimental simulator could measure was limited (0.3, 0.7, 1.1, and 1.5 mm for 1, 2, 3, and 4 mm mismatches, respectively). Measurements from cases expected to produce separations below these values were instead insensitive to all of the tested variables. For example, the experimental results between cases 2 and 4 ( Fig. 4) were insensitive to an increase in SwPL of approximately 500%. PyEL, however, was able to clearly differentiate these cases, based on the resultant force direction. As cases with a relatively low maximum separation for the applied mismatch were most affected by this issue, the greatest effects were seen under high SwPLs.
From both the experimental and PyEL results, the highest maximum separation occurred for cases with the lowest SwPL and highest peak ML force (corresponding to the highest mismatch) (Fig. 4A, case 1). Conversely, the lowest maximum separations were seen with a high SwPL and low peak ML force (Fig. 4A,  case 4). A steeper inclination angle caused a higher maximum separation.
When excluding cases where the simulator had lost measurement sensitivity, the PyEL model provided a strong prediction of the experimental result, with a maximum error of 0.3 mm and a slight overestimation. The overestimation was more pronounced with a lower inclination angle (e.g. case 1 in Fig. 4B compared to C), which may be related to increased frictional forces in the experiment (lower inclination angles reduced the separation, thereby increasing the spring force and contact force).

Study 3: Effects of inertia and friction on the computational results using the PyEL and dynamic FE models
The PyEL model and dynamic rigid, with friction, FE model were compared with a 4 mm mismatch, 100 N/mm spring constant, and SwPL from 70 to 300 N for clinical inclination angles of 45°and 65°.  Mean PyEL maximum separation error was a 0.08 mm overestimation (SD ± 0.05 mm) across these cases, which is smaller than typical experimental measurement error (Blumenfeld et al., 2011).
Separations through the gait cycle, for the PyEL and dynamic rigid FE (with and without friction) models in the more extreme 65°case, are given in Fig. 5. Inertia added a delay to the head movement into and out of the liner, up to approximately 0.025 s (2.5% of the total cycle time) at heel-strike and 0.01 s (1% of the cycle time) at toe-off, but had a negligible effect on the maximum separation (<0.01 mm difference across all inertia but no friction cases). Delaying the head movement at heel-strike extended the length of edge loading and therefore the peak axial force experienced (from 1400N to 2400N for the case in Fig. 5A).
Friction had a variable effect that was dependent on the SwPL, at 300N the maximum delay at heel-strike was 0.025 s. Friction reduced the maximum separation by 0.03 mm (0.8%) and 0.17 mm (5.4%) with 70 N and 300N SwPLs, respectively. Friction had a negligible effect on the time point at which the head left concentricity and moved into edge loading (toe-off).

Discussion
This study aimed to demonstrate the ability of a rapid, quasistatic, rigid analysis model, to predict separation between the head and liner, under various orientations and load conditions relevant to pre-clinical edge loading evaluation. PyEL replicated the experimental trends regarding the effect of inclination, mismatch, and swing phase load. As differences between the PyEL and experimental testing data were <0.3 mm, and typical experimental error is approximately 0.5 mm, PyEL pre-screening data could be used to select appropriate cases for wear testing when the practical limitations of the experimental testing process are taken into account.
In addition to a reduction in cost and time, by either replacing or refining experimental or FE parameter sweep testing, an additional advantage is the provision of context to the results in terms of the contacted liner features. The relationship between headliner separation (which is measureable experimentally) and liner contact location (which is not), is heavily dependent on the liner geometry and orientation. For the bearing used in this study, maximum separations >0.01 mm induced contact between nonconforming surfaces, and separations >1.6-1.9 mm (45-65°inclination) induced translational sliding of the contact area over the rim. This richer data can be used to answer research questions about the source of particular device damage modes.
The PyEL model compared well with FE models which included inertia, in terms of maximum separation, and therefore in terms of the prediction of the extremes of edge loading position. However, inclusion of inertia caused a delay to head movement at heel-strike and toe-off, which could have a large effect on the peak forces experienced under edge loading conditions (70% increase). In consequence, although it would be trivial to calculate resultant contact force data from PyEL, care should be taken in interpretation of that information from those static, rigid models alone.

Model roles, efficiency and robustness
The more complex FE models, and the experimental testing, offer advantages for specific scenarios in terms of wear/damage results or detailed contact area/pressure information. However PyEL is valuable for design screening, characterising trends, and explaining underlying mechanisms. FE analysis is the current standard method for numerical prediction of mechanical behaviour in this field. The quickest running FE solutions (static) took 1 h for a single case, compared to 10 s for the PyEL code. In addition, the FE approach was less robust, with occasional solver convergence issues requiring user intervention, making it less suitable for large parametric tests. The longest run time FE solutions (including inertia and friction) took 10 h to solve one case. For context, an experimental study by Partridge et al. (2019) compared two different bearing designs with a wide, but relatively low resolution, parameter sweep that resulted in a total of 300 tests. Liu et al. (2018) used a multi-body dynamics model to investigate a single implant design using 50 cases. Further optimisation of the FE models to minimise run times, via mesh optimisation and time incrementation scaling, could have reduced the computational time, with an associated compromise on results quality, but not to anywhere near equivalent to the PyEL model.
The robustness and relevance of pre-clinical testing is dependent upon the ability to confidently identify high-risk situations. Increasing the number of situations assessed in a reasonable time by a factor of 360-3600 increases the likelihood of identifying those which require further investigation using either computer based or experimental techniques. The PyEL model is also not dependent on the specialised skill sets, equipment, or commercial software needed for FE modelling and would therefore be more directly accessible to experimental design and test engineers.

Limitations and future work
The effect of inertia will be sensitive to loading variations. The rate of change of load magnitude tested here is relevant to walking gait. Higher rates are certainly possible, in events such as stumbling (Bergmann et al., 1993), and the use of a quasi-static model may be more limited.
The analysis of the effect of friction used a relatively arbitrary friction coefficient of 0.05 which was generated for standard, concentric, loading conditions (Brockett et al., 2007). It will not represent the complexity of the experimental lubrication conditions when the head slides into and out of the liner (Jin et al., 1997;Myant et al., 2012) and may alter the magnitude of the changes, if not the trends.
PyEL has only been investigated for a single CoC bearing design and further validation would be required for substantially different geometries and hard-on-soft bearing combinations (which currently represent 80% of UK THR surgeries (National Joint Registry, 2018)).
The ISO edge loading test described in this study is valuable as a standardised comparison point, but is a simplification of a complex in vivo mechanism. Assumptions made during the model development allowed it to solve quickly and extremely robustly across very large case lists. These assumptions restrict the flexibility of the model, however, meaning that it would not be possible to solve for predictions of hip separation based on in vivo measurements of joint and tissue forces and motions. Separations generated in this testing were in line with in vivo measurements of hip separation though, which to date have recorded maximum values between 0 and 5.6 mm, with means across multiple patients <4 mm (Blumenfeld et al., 2011;Dennis et al., 2001;Glaser et al., 2008a;Glaser et al., 2010;Komistek et al., 2002;Lombardi et al., 2000;Tsai et al., 2014).
Ongoing clinical (Boese et al., 2016;Glaser et al., 2008b) and in vivo modelling (LaCour, 2017) findings can be used for continuous improvement of the test cases and force environments used in pre-clinical testing. For example, subject-specific change in pelvic tilt angles during an activity could also alter the separation behaviour (Vasiljeva et al., 2018).

Conclusion
The results from this study suggest that a rapid, quasi-static mechanics based model can provide a good estimate of the maximum separation and contact location generated under a total hip replacement ISO standard edge loading protocol. Critically, models provide a direct link from the orientation and force environment to the position of the contact at the bearing surface, for a given device design. Additionally, these predictions were not substantially improved by the use of more time consuming FE methods or the inclusion of inertial effects. The use of the PyEL methodology allows for more extensive evaluation of edge loading kinematics under ISO 14242:4, thereby increasing the likelihood of identifying high risk situations within practical timescales.