A Probabilistic Approach for Studying the Reliability of Cementless Hip Prostheses in the Presence of Mechanical Uncertainties

This paper presents a probabilistic approach for studying the reliability of cementless hip prostheses in the presence of mechanical uncertainties and its application to the investigation of the influence of bone-implant interface properties. The non-linear deterministic model of the bone-implant coupled system and its finite element implementation are described, and the proposed reliability analysis is exposed. It is demonstrated that the distribution (uniform, truncated Gaussian and truncated lognormal distribution) of the two chosen parameters and the truncation lengths have a minor influence on the Hasofer-Lind index. This index logically increases as the failure threshold increases. FORM and SORM approximations are compared with the results obtained using a crude Monte-Carlo method for the estimation of failure probability. The performance of three Monte-Carlo methods is studied in terms of the necessary number of FE calculations. The method based on the Directional Simulation (DS) technique is efficient and less time-consuming. The validity and operational capacity of the proposed approach would not be compromised by an increase in the number of uncertain parameters.


Introduction
The amount of relative micro-movement between bone and implant plays a crucial role in the success of Total Hip Replacement (THR) [1].Ideally, to promote the Primary Stability of the prosthesis, the amount of relative micro-movement between bone and implant induced early after the operation, before any biological process takes place, must remain lower than a prescribed threshold.If not, the lack of Primary Stability may lead to the aseptic loosening of the implant [1].The primary stability of a cementless hip stem is not only affected by the implant design, but also by biomechanical properties, among which the coefficient of friction at the bone-implant interface or the elastic properties of the cancellous bone are known to play a crucial role.The natural variability of these parameters has to be taken into account to predict the Primary Stability of the THR.This is practically unfeasible expe-rimentally, due to the great number of tests which would have to be carried out.The combination of numerical models and probabilistic methods offers a promising alternative.
Probabilistic methods enable the estimation of the effects of parameter variability in the resulting statistical variation of the system response to be determined.Each parameter is generally represented as a random variable, and the propagation of their randomness into the system response is predicted.By understanding the distribution of performance, evaluations of quality and risk assessment can be performed.Recent studies have proposed probabilistic approaches to assessing the structural integrity of orthopaedic implants [2][3][4][5][6][7][8].Previously, Browne et al. [9] applied reliability theory to aid in fracture mechanics-based life prediction procedures for a tibial tray component represented as a cantilever beam subjected to constant amplitude loading.Dar et al. [10] demonstrated how Taguchi and probabilistic methods could complement each other to account for uncertainties when predicting stresses with finite element analysis, in a study of a fixation plate represented as a cantilever beam.The femoral stem component of a total hip replacement has been the subject of several probabilistic structural integrity studies [11].Bah and Browne [12] used an idealized cylindrical finite element model to represent an implanted cemented hip stem in order to assess the most likely failure mode and to identify which parameters had the largest contribution, where geometry, material properties and loading were considered as random variables.The latest studies concerning the application of probabilistic approaches to cementless hip prostheses have focused on the effect of femur characteristics and implant design geometry [13] or implant positioning on Primary Stability [14].
This paper presents a probabilistic approach for studying the reliability of cementless hip prostheses in the presence of mechanical uncertainties and its application to the investigation of the influence of bone-implant interface properties.In Section 2, the nonlinear deterministic model of the bone-implant coupled system and its finite element implementation are concisely presented.The uncertain parameters of the problem, the control variable and the associated failure criterion defining Primary Stability are then detailed.Section 3 is dedicated to the probabilistic modelling of the problem.Section 4 exposes the proposed reliability analysis.In Section 5, numerical experiments are presented in order to quantify the influence of the statistical characteristics of the random parameters, the Young modulus of the cancellous bone and the friction coefficient at the interface between the cancellous bone and the stem on the reliability of the bone-implant coupled system.

Deterministic Modelling
In this section we first concisely present the deterministic model used in this study, from its mechanical formulation to its practical implementation, as well as its finite element formulation.Then we highlight the uncertain parameters of the model, we introduce a failure criterion for the latter and we express this criterion in terms of the selected uncertain parameters.

Mechanical Formulation and Associated Numerical Model
The mechanical problem to solve consists of the contact with friction of two three-dimensional deformable bodies B 1 and B 2 occupying the domains V 1 and Denoting by c S the contact surface between B 1 and B 2 , and considering the prin-ciple of virtual work, this problem may be formulated in the following way (variational formulation): σ ε e δu e δu e δu (1) where dV and dS are the volume and surface measurements, respectively, σ is the Cauchy stress vector defined as ( ) is the associated virtual strain vector corresponding to the imposed virtual displacement vector ( ) is the vector of known externally applied force per unit volume, ( ) = − δu δu δu denotes the relative virtual displacement vector field of the two bodies in the contact zone.
The contact tractions acting on the contact surface c S may be decomposed into normal and tangential components as follows: 21 21 where λ and t are the normal and tangential tractions, and 21 n and 21 s are the normal and tangential unit vector fields associated with c S .Note that at this step all the parameters used are fields depending on the spatial coordinates ( ) , , ξ ξ ξ .
Let g be the gap function for the contact surface pair.It is a scalar function defined as follows: Let ( ) , , ξ ξ ξ ξ = be the coordinates vector of point ξ is a function of 1 ξ .Then ( ) g ξ is defined as: where 21  n is the unit normal outward vector to P .With this definition, the normal contact conditions can be written: where the last equation expresses the fact that the normal traction component is compressive, and exists only if the contact is effective i.e. if 0 g = .The tangential contact conditions are taken into account by assuming that Coulomb's law of friction holds locally on the contact surface.Let µ be the coefficient of friction.Coulomb's law of friction is usually formulated as a function of the local tangential relative velocity at the contact interface.However, in the static case it can be reformulated using the relative tangential displacement ( ) 1 , v ξ e due to a loading increment Δe of V e and S e : where 21 s is the unit tangential vector at point 2 P .
Let us also define the non-dimensional variable τ given by: With this definition, Coulomb's law of friction can be written: ( ) ( ) The constraint function method proposed by Bathe [15] is used to impose the normal and tangential contact conditions.
Let n w be a real-valued function of g and λ such that the solution of ( ) , 0 n w g λ = satisfies the normal contact conditions, and let s w be a real-valued function of v and t such that the solution of ( ) satisfies the tangential contact conditions.Thus the initial problem (1) takes the form: σ δε e δu e δu n s δu (9) The finite element formulation of this problem is achieved using the standard procedures exposed in reference [16].The displacement field u is approximated using an interpolation procedure based on nodal point displacements, leading to the classical expression: ) ( ) , where H is the dis- placement interpolation matrix and U is the vector of the nodal displacements.
Concerning the discretization of the contact conditions, the constraint function method proposed by Bathe [15] is adopted.This results in adding a nonlinear term to the initial problem, dependent on the unknown contact forces vector at the contact nodes of body 1  B .Denoting by c f this term, by l ∈  d the global unknown vector gathering U and c f , and by l ∈  f the prescribed conditions vector gathering the external nodal prescribed loads and the contact constraints, the displacement-based finite element formulation of the initial problem then takes the form: where the nonlinear function K from l  into l  is the mechanical operator, and the integer l is the dimension of the discretized problem; that is to say, of the finite element model.Equation ( 10) is solved using an incremental procedure coupled with the Newton-Raphson method.This leads to the calculation of the unknown vector 1 n+ d at loading increment 1 n + from the solution n d at increment n .At each step of this procedure, the incre- mental equation is obtained by linearizing Equation (10) about the last calculated state.Using the procedure described in [16], the matrix incremental equation obtained at iteration step ( ) where

Δd
is the increment of unknown d at iteration step ( ) i , and ( ) K is given by: where ∇ a denotes the gradient with respect to a .The solution to Equation (11) gives the increment ( ) i ∆d and therefore determines the value of the state vector 1 n+ d at iteration step ( ) The iterations are continued until a specified convergence criterion is satisfied.
Practically, this problem is solved using the FEMAP-NASTRAN software package [17].The stem and the bone are discretized using four-node tetrahedral elements leading to a model of 79067 elements and 18701 nodes (Figure 1).
A contact pair between the stem and the bone is defined in order to take into account the contact with friction between these two bodies.The host bone is subdi- vided into two regions, corresponding respectively to the cortical and the cancellous zones.The femoral head and femoral stem are tied together.The femur is fully constrained at its distal part.Frictional contact at the bone-implant interface is taken into account, leading to nonlinear behaviour.
The position of the femur reproduces the standing position in vivo.The applied force is vertical and is progressively increased to a maximum value corresponding to the maximum load during a walk cycle of a 75 kg person in the unilateral compression position.The numerical results obtained are validated by comparison with in vitro test results [18].

Uncertain Parameters of the Model
Until now, all the quantities appearing in Equation (10) (i.e. the operator K and the vectors d and f) have been supposed to be deterministic, that is to say not affected by uncertainties.
We now assume that the known prescribed conditions vector f is deterministic and that two mechanical parameters of the operator K are uncertain: the Young's modulus of the cancellous bone and the coefficient of friction at the interface between the cancellous bone and the femoral stem, respectively denoted 1 y and 2 y .Therefore K is uncertain through these two parameters.Let ( ) , y y = y be the vector of 2  gathering them.By hypothesis, K depends on y .Consequently, from Equation (10), the unknown vector d also depends on y and is therefore uncertain.

Control Variable and Associated Failure Criterion
Formally, the solution to Equation ( 10) can be written: (14) where S is a nonlinear operation from l  into l  such that: A control variable is a real variable, for example a scalar displacement (or a stress, or a strain) at a particular point of the studied mechanical system, whose value must be controlled within the framework of the reliability analysis of the system.
In this work, the chosen control variable is the Euclidean norm of maximum relative displacement between the femoral stem and the cancellous bone.
Let w be such a variable.It depends on the unknown d through a relationship of the form: where C is a known function from l  into +  , called observation operator.Inserting Equation ( 14) into Equation ( 16) yields: where D C S =  is also a known function from l  into +  .As a result, w can also be viewed as a function of f .Note that the functions C and D are known, in the sense that they are described by known numerical models.
As seen in Section 2.2, the mechanical operator K depends on the uncertain vector parameter ( ) , y y = y .As a result, its inverse 1 − = S K also depends on y , as does the function D , as can be seen from Equation (17).Indicating this fact by the new notation y S and D y in the place of S and D , Equations ( 14) and ( 17) be- come: and, since f is a constant vector of l  , can be re- written: where Let 0 w be a given admissible value of w , called the failure threshold, and let u be the real variable such that: Such a variable is called the safety margin associated with the control variable w , and its sign defines two fundamental states for the mechanical model: the safe state if > 0 u and the failure state if 0 u ≤ .The value X. S. HU ET AL. 16 0 u = characterizes a limit situation called the failure limit state.According to Equation (19), u can be re- written: where F is a function from 2  into  called the limit state function of the mechanical model.This function defines three specific subsets of 2  : which verify: and which are respectively called the safe domain, the failure domain and the limit state curve.The condition 0 u ≤ defining the failure state is the failure criterion associated with the control variable w and its limit value 0 w .The failure domain f D is the geometric representation of this criterion.
Let us note that, since the function F is only known numerically (through the finite element model described in Section 2.1), the sets s D , f D and C can only be determined numerically.

Probabilistic Modelling
In the following, all the considered random variables (RV) are assumed to be defined on the same probability space ( ) where Ω is a sample space,  is a σ -algebra of subsets of Ω and  is a probability on  .

Stochastic Modelling of the Uncertain Parameters
In order to take into account the random variability of the uncertain parameters 1 y (the Young's modulus of the cancellous bone) and 2 y (the coefficient of friction between the cancellous bone and the femoral stem), these two parameters are modelled as continuous RVs, denoted Y respectively.The continuous 2   -valued random vector ( ) gathering these two scalar RVs is the probabilistic model of the uncertain vector parameter ( ) . In view of the intended numerical applications, this random vector is assumed to satisfy the following hypotheses: (H1) Its components 1 Y and 2 Y are independent.(H2) They both follow the same type of probability distribution.
(H3) Three types of probability distribution are ad-missible for these components: uniform, truncated Gaussian and truncated lognormal.
Concerning these hypotheses, we can make the following remarks: (R1) Let ( ) respectively.Then, from (H1): (R3) In connection with (H3), let j Y be a scalar continuous RV with pdf j Y p .Then j Y is said to follow: a) an uniform distribution with support , where j S 1 is the indicator function of j S (i.e. ( ) ( ) ≤ +∞ , and j m , j σ are the solutions of the nonlinear system: , .
where φ and Φ are respectively the one-dimensional standard Gaussian pdf and the one-dimensional standard Gaussian cumulative distribution function, such that: 1 e , d , 2π where j S is the interval ( ) ≤ +∞ , and j m , j σ are the solutions of the nonlinear system: j j j j j j j j j j j j j j j j j j

Safe and Failure Events
Equations ( 19) and (22) show that when the uncertain vector parameter ( ) is modelled as a random vector ( ) , the control variable w and the safety margin u , which depend on y , become scalar RVs.Let W and U be respectively these RVs.They are such that: and are both defined on the probability space ( ) Naturally, it is assumed here that the definition domains of H and F (which are coincident according to Equation ( 22)) contain the support of the probability distribution of Y .
In this random context, the safe state and the failure state are defined by two events, the safe event s E ∈  and the failure event f E ∈  , such that: and which verify: E s is the event associated with the safe domain D s and E f is the event associated with the failure domain D f .

Fundamental Objective
The fundamental objective of reliability analysis [19] is to calculate the probabilities that the events s E and f E will occur, that is ( ) . These probabilities are given by: and, according to Equation (37), satisfy: ( ) ( ) where ( ) From Equation ( 39), ( ) is known as soon as ( ) is known and vice versa.That is why in this work our attention is only focused on ( ) , called the failure probability and denoted f P in the following.This probability can be rewritten: and its calculation requires the use of a numerical procedure.For the numerical applications treated in this study, we have chosen a Monte Carlo method ( [19][20][21]).

Standard Formulation
In the classical reliability processes [19], it is customary to transform the initial formulation into a standard formulation in which the vector of random parameters follows a standard Gaussian distribution.This leads in our case to the construction of a regular transformation T , with inverse 1 − T , such that the random vector ( ) is a 2  -valued standard Gaus- sian random vector defined on ( ) , , Ω   .Once this transformation has been constructed, the failure event can be expressed in terms of X as fol- lows: where Γ is a mapping from 2  into  such that, 2 ∀ ∈ x : ( ) ( ) ( ) From Equation (32), the failure probability is then given by: ( ) where ( ) : 0 And 2 φ is the bidimensional standard Gaussian pdf, given by: ( ) ( ) Note that Equation (44) can also be obtained by carrying out the variable change ( ) = y T x in Equation (40).The set f ∆ is the failure domain in the standardized x -space, that is in the space (identified with 2  ) of the X random vector's realizations.Its complementary in 2  , that is the subset : > 0 , and the common boundary between f ∆ and s ∆ , that is the curve : 0 , are respectively the safe domain and the limit state curve in this standardized space (Figure 2).Equation ( 44) represents the standard formulation of the reliability problem.To estimate this integral, three Monte Carlo methods were used: the crude direct method and two more refined methods, based respectively on the importance sampling technique and the directional simulation technique ( [19][20][21]).
Such a formulation is completely defined as soon as transformation T is known, and the latter only depends on the probability distribution of the random vector ( ) . As the components 1 Y and 2 Y of Y are assumed independent, this transformation is of the form: where j T , , is the j -th coordinate of T and Y are both assumed to follow the same type of probability distribution): a) , follows a uniform distribution with support , , follows a truncated Gaussian distribution with mean where ( ) is the solution of the system (30).
where ( ) is the solution of the system (33).Note that the fact that 1 Y and 2 Y are both assumed to follow the same type of probability distribution is not a loss of generality.It is a simple working hypothesis, which could be changed without any incidence on the proposed methodology.

Hasofer-Lind Index
The Hasofer-Lind index HL β ( [19,22]) is a reliability indicator whose calculation is easier than that of failure probability f P .It is defined in the standardized x -space as follows: ( ) where ( ) denotes the usual Euclidean distance between the origin point O and the failure domain f ∆ .Such an index is thus given by: is the canonical Euclidean inner product on 2   .

FORM Approximation [19,22,23]
This approximation (FORM means First-Order Reliability Method) consists of replacing the failure domain f ∆ by the half-space: where FOTA Γ is the First-Order Taylor Approximation of Γ at point * M , such that: whose graph FORM Σ is a straight line tangent to Σ at * M (Figure 2).The failure probability (44) is then approximated by: and a simple calculation gives: ( ) where HL β is given by Equation (54).

SORM Approximation [19,24]
This approximation is a refinement of the previous method (SORM means Second-Order Reliability Method), which is based is based on the replacement of the failure domain f ∆ by the set: where SOTA Γ is the Second-Order Taylor Approximation of Γ at point * M , that is the quadratic function given by: whose graph SORM Σ is a second-order curve tangent to Σ at * M (Figure 2).In Equation (60), 2  ∇ Γ is the Hessian operator associated with Γ and is assumed to exist at point * M .Using this substitution, the failure probability (44) is then approximated by: and this integral can be evaluated by using approximate formulas, such as the Breitung asymptotic formula [25]: or the Hohenbichler asymptotic formula [26]: is an improvement of the previous one.In these formulas, χ is the curvature of the limit state curve Σ at

Numerical Experiments
The numerical applications presented in this section concern the model of cementless hip prosthesis described in Section 2.1 and aim at estimating its reliability in various situations, either through the Hasofer-Lind index HL β (calculated from the Rackwitz-Fiessler algorithm)   or by means of the failure probability f P (estimated from Monte Carlo simulations or by using the FORM and SORM approximations).All these methods have been implemented by the authors.We recall that the boneprosthesis system is considered as reliable if the Euclidean norm of the maximum relative displacement between the femoral stem and the cancellous bone do not exceed a given admissible value 0 w (the failure threshold).

Influence of Probability Distributions and Scatterings
The purpose of this application is to highlight the influence of distributions and scatterings of random parameters on bone-prosthesis system reliability, expressed in terms of the Hasofer-Lind index HL β .Y and 2 Y are truncated lognormal with sup- , , a b respectively, where in each case, for 1, j = : The chosen values for The results obtained are summarized in Table 1.In each case, we can see the significant effect of the coefficients of variation on the index HL β .We can also ob- serve that the effect of 2 Y cv (coefficient of variation of the coefficient of friction) is greater than that of 1 Y cv (coefficient of variation of the Young's modulus of the cancellous bone).Finally, we can point out that the probability distributions of 1 Y and 2 Y have a minor influence on HL β .

Influence of Truncation Lengths
In this application, the random parameters 1 Y and 2 Y are assumed to be independent and distributed according to truncated distributions with variable supports and constant first two moments.Under this assumption, the objective is to evaluate the effect of the support length of each distribution (so-called truncation length) on the reliability of the hip prosthesis expressed in terms of the Hasofer-Lind index HL β .
The chosen values for the statistics ( ) and the failure threshold 0 w are the following: The results obtained are listed in Table 2. First we can see that the behaviour of HL β with respect to the trun- cation lengths depends on the nature of the considered probability distributions.Indeed, in the truncated Gaussian case this behaviour is non-monotonous ( HL β de- creases then increases) in each considered situation (TG1 and TG2), while in the truncated lognormal case, it is either increasing monotonous (situation TLN1) or decreasing monotonous (situation TLN2).In each case, however, the influence of the truncation lengths on HL β is not very significant.

Influence of the Failure Threshold
In this section, we study the influence of the failure threshold 0 w on the reliability of the prosthesis, again using the Hasofer-Lind index HL β as a reliability indi- cator.
To this end, the random parameters 1 Y and 2 Y are assumed to be independent and uniformly distributed with means and standard deviations given by Equation (67); five values are successively considered for the failure threshold w (in m µ ): 270, 275, 28 0, 285, 290.The obtained values of HL β are given in Table 3 which also contains, for information, the corresponding values for the FORM approximation FORM f P of failure probability, given by Equation (58).
We can first observe that the values of HL β are quite low for all the considered values of 0 w , which are nevertheless common values for this threshold.This comes quite obviously from the chosen values for the coefficients of variation of 1 Y and 2 Y , which are relatively high ( ) . We can also underline that these results are logical in the sense that HL β increases as 0 w increases.It is interesting, moreover, to note that this increase is quasi-linear over the range of the considered 0 w values.

Comparison of Various Approximations of P f
In this application we consider the same example as in the previous application (see Section 5.3), and for each value of the failure threshold 0 w we compare the FORM and SORM approximations FORM f P and SORM f P of f P (given by Equations ( 58) and (63) respectively) with the target value of this probability, obtained by using the crude Monte Carlo method taken as the reference me-thod.
Estimating f P from the latter requires here the simulation of 4 10 realizations , of the standard Gaussian vector X , then, for each simulated realization j x , the calculation of the corresponding value ( ) j Γ x of the limit state function Γ .Such an estimation requires thus 4 10 finite element calculations.This is why two other less expensive Monte Carlo methods were used: one based on the importance sampling (IS) technique, the other on the directional simulation (DS) technique ( [19][20][21]).
The values of f P given by these two improved Monte Carlo methods and by the crude Monte Carlo method are respectively denoted  2) are not accurate enough to adequately represent Σ in a sufficiently large neighbourhood of the design point * M .The results in Table 5 confirm that the Monte Carlo methods based on the IS and DS techniques are more efficient than the crude Monte Carlo method.We can also observe that the method based on the DS technique is less time-consuming, and therefore more efficient, than the one based on the IS technique.

Conclusions
The aim of this paper was to present a methodology to enable the prediction in a probabilistic context of the primary stability of cementless hip prostheses in the presence of uncertainties.Based on a reliability approach, this methodology supposes the numerical model describing the mechanical behaviour of the prosthesis-cancellous bone couple, the probabilistic models of uncertain parameters and the failure criterion expressing the loss of stability of the prosthesis to be known.The underlying mechanical problem is a classical problem of contact with friction between two three-dimensional deformable bodies: a Depuy Corail® prosthesis is implanted without cement in a femur.This problem is solved using a finite element model based on Bathe's formulation.Due to their variability in the literature, we chose to consider the Young's modulus of the cancellous bone and the coefficient of friction between the cancellous bone and the femoral stem as uncertain parameters.They were modelled as random variables, assumed successively to be uniform, truncated Gaussian and truncated lognormal.Finally, failure was defined as the exceeding of a given limit value by the maximum relative displacement between the femoral stem and the cancellous bone.The reliability of the coupled system was estimated using two indicators: the Hasofer-Lind index and the failure proba-bility.The first was calculated using the Rackwitz-Fiessler optimization algorithm.The second was estimated first by using approximate formulas derived from FORM and SORM approaches, then by means of three Monte Carlo methods: the crude Monte Carlo method, the method based on the importance sampling technique and that based on the directional simulation technique.
The presented numerical applications have shown the relevance of the proposed methodology and given an indication of its great potential in the field of prosthesis reliability.A different failure criterion could be envisaged, such as the acceptable maximum stress or the optimization of the contact area, for instance.
This methodology has been presented in the case where only two parameters are uncertain.It is clear, however, that its validity and its operational capacity would not be compromised by an increase in the number of uncertain parameters.This personal probabilistic model enables clinicians to compare the primary stability of cementless prostheses, and to choose the optimal implant design while taking into account uncertainties in the model input data.

Figure 1 .
Figure 1.Mesh of the FE model used for the study of thefemur-cementless prosthesis coupled system: DePuy Corail®.

2 Yp
and p Y be the probability den- sity functions (pdf) of 1 Y , 2 Y and uniform, or truncated Gaussian or truncated lognormal.

Figure 2 .
Figure 2. Geometric representation of the HL β -point

2 ⋅(Figure 2
denotes the canonical Euclidean norm on 2  .The Rackwitz-Fiessler algorithm ([23,24]) is well suited to solving the constrained nonlinear optimization problem (52) and it was chosen to treat the numerical applica-tions presented in Section 5. Note that for these applications the problem (52) has a unique solution.The solution point: called the design point, or HL β -point, and is located on the boundary Σ of f ∆ where Γ is the limit state function defined by Equation (43), ∇Γ is its gradient (assumed to exist at point *

φ
and Φ are the functions given by Equations (31), and the following conditions are assumed to be verified:

1 Y 2 Y 2 Y 1 Y
quantify the scattering of 1 Y and 2 Y respectively, three values are considered in each j C case according to following strategy: 1) cv is fixed to its reference value 0.10 and cv successively takes the values 0.03 , 0.10 , 0.30 ; 2) cv is fixed to its reference value 0.10 and cv successively takes the values 0.03 , 0.10 , 0.30 .