A Genetic Algorithm Based Multi-Objective Shape Optimization Scheme for Cementless Femoral Implant Souptick Chanda Department of Mechanical Engineering, Indian Institute of Technology Kharagpur, Kharagpur, West Bengal 721 302, India
Sanjay Gupta1 Department of Mechanical Engineering, Indian Institute of Technology Kharagpur, Kharagpur, West Bengal 721 302, India e-mail: [email protected]
Dilip Kumar Pratihar Department of Mechanical Engineering, Indian Institute of Technology Kharagpur, Kharagpur, West Bengal 721 302, India
The shape and geometry o f femoral implant influence implantinduced periprosthetic bone resorption and implant-bone inter face stresses, which are potential causes o f aseptic loosening in cementless total hip arthroplasty (THA). Development o f a shape optimization scheme is necessary to achieve a trade-off between these two conflicting objectives. The objective o f this study was to develop a novel multi-objective custom-based shape optimization scheme for cementless femoral implant by integrating finite ele ment (FE) analysis and a multi-objective genetic algorithm (GA). The FE model o f a proximal femur was based on a subjectspecific CT-scan dataset. Eighteen parameters describing the na ture o f four key sections o f the implant were identified as design variables. Two objective functions, one based on implant-bone interface failure criterion, and the other based on resorbed proximal bone mass fraction (BMF), were formulated. The results predicted by the two objective functions were found to be contra dictory; a reduction in the proximal bone resorption vras accom panied by a greater chance o f inteiface failure. The resorbed proximal BMF was found to be between 23% and 27% fo r the trade-off geometries as compared to ~39% fo r a generic implant. Moreover, the overall chances o f interface failure have been mini mized fo r the optimal designs, compared to the generic implant. The adaptive bone remodeling was also found to be minimal fo r the optimally designed implants and, further with remodeling, the chances o f interface debonding increased only marginally. [DOI: 10.1115/1.4029061] Keywords: total hip arthroplasty, finite element analysis, bone remodeling, inteiface failure, genetic algorithms, shape optimization, multi-objective
1
Introduction
The adverse effects of implant-induced bone resorption and ex cessive interface stresses have often been reported as potential causes of failure in cementless THA [1,2]. The success of THA notwithstanding, about 10% of the hip replacements need revision surgery depending on the patient’s conditions, disorder, and the type of implants used [3,4], It is well known that the shape and
geometry of a cementless femoral implant influence periprosthetic bone remodeling and implant-bone interface stresses [5,6]. The mechanical designs of hip prosthesis, therefore, need to search for optimal implant geometry, such that the effects of stress shielding induced bone resorption and interface stresses are minimized. Shape (or geometry) optimization is a particular stage of struc tural optimization, which searches for the optimal configuration of a design domain [7]. Employing shape optimization as a design tool, stem geometries can be evaluated based on some biomechan ical cost (or objective) functions. Therefore, the formulation of suitable cost functions, representing the global effects of these failure mechanisms is necessary. These cost functions may some times be mutually conflicting in nature [2,8]. Therefore, develop ment of a multi-objective shape optimization scheme is necessary to achieve a trade-off between these two conflicting objectives, while minimizing them within the domain of feasibility. Despite a few earlier studies on single-objective optimization [9-12], the hrst major study on multi-objective shape optimization was carried out by Fernandes et al. [13] to analyse the design pros pects of a 2D cementless hip stem model. Thereafter, the hrst GA based multi-objective shape optimization strategy, to minimize proximal and distal stress in cement mantle of a cemented THA. was proposed by Ishida et al. [14], In a more recent study, Ruben et al. [15] presented a 3D shape optimization scheme based on a multicriterion objective function to design implants with an improved initial stability and reduced proximal bone resorption. Although several shape optimization schemes were employed and new stem designs were proposed, these studies were based on generic femoral bone model, consisting of a cortex layer covering a homogeneous core of cancellous bone, thereby neglecting the heterogeneous bone material property across the proximal femur [12-15], The study by Katoozian and Davy [10] attempted to model bone heterogeneity by linearly varying bone material prop erty from cancellous to cortical bone. However, such method also suffers from oversimplification of bone property. Moreover, these studies suffered from additional limitations regarding the choice of objective function for bone remodeling simulation. Most stud ies inappropriately considered strain energy density (SED) of the implanted bone only as the mechanical stimulus for bone remod eling [13,15], instead of the difference in SED between intact and implanted cases [16]. The parameterization schemes used in earlier studies for defin ing 3D implant geometry were unable to account for a wide vari ety of implant shapes [10,14,15]. The search for an optimal geometry demands extensive exploration of all possible and ad missible shapes. However, this may add to more complexities which may not be suitably handled by traditional optimization tools. Moreover, shape optimizations problems were reported to be typically nonconvex and characterized by multiple local optima [7], which requires a robust heuristic search method, such as the GA, for its complex multivariable analysis capability and gradient-independent approach. Thus, the GA seems to be a more appropriate tool for solving such multimodal and multivariable space optimization problems. Nevertheless, hardly any such study has been attempted in the field of shape optimization of cementless femoral implant. The objective of this study was to develop a novel multi objective shape optimization scheme for cementless femoral implant by integrating FE analysis and GA. The study was based on a subject-specific FE model; in which heterogeneous bone material properties were assigned element-wise using the CT-scan dataset. Geometric parameters were introduced to char acterize the femoral stem shape. The objective functions were formulated based on common failure mechanisms of cementless THA.
2 'Corresponding author. Manuscript received March 25. 2014: final manuscript received November 5, 2014; published online January 29, 2015. Assoc. Editor: David Corr.
Journal of Biomechanical Engineering
Materials and Methods
The present optimization scheme comprised of a parametric implant model generator, a script for virtual implantation, an
Copyright © 2015 by ASME
MARCH 2015, Vol. 137 / 034502-1
V
Fig. 1 Parametric scheme for evolution of key sections of the implant: ellipse (p = 2) to superellipse (p >2)
automatic mesh generator, an FE analysis script and finally an optimization controller. In order to integrate all these components, a master program was written in m atla b . The master program served as the prime mover and could invoke all components logi cally and sequentially via batch mode processing. An elaborate description of each method that collectively contributed toward development of the optimization scheme is stated in Secs. 2.1 to 2.4. 2.1 Parameterization and Constraint Specification of Implant Geometry. In a local coordinate system ary, the pair of equations that defined the geometric parameterization of stem transverse section was expressed by the parametric form as follows: .v = a cos(2/,/'' t y = b sin*2|/,,) t
0>
where the limits of t were prescribed as 0 < r < zr/2. The parame ters a, b, and p at four key sections, thus, characterized the design parameters of the femoral stem (Fig. 1). Using these four key sec tions, a virtual 3D parametric model of the implant was generated
in SolidWorks cad package (SolidWorks, Dassault Systemes S.A.) by executing an Application Programming Interface script written for the purpose. The initial design of the stem was similar to a generic design of a collarless TriLock (DePuy) prosthesis having stem-length 113 mm (Fig. 2(a)). The stem-length was kept unaltered throughout the optimization procedure. A total of 18 geometric parameters was introduced as design variables for the modeling purpose (Fig. 2(a)). These design varia bles had lower and upper bounds and geometric constraints in order to achieve clinically admissible shapes. Furthermore, it was also required that the width and thickness of cross sections would decrease from proximal to distal end in order to facilitate ease of withdrawal, if required, in case of a revision surgery. The parame ters al-a6 corresponded to the length of major axes on the medial and lateral sides of the four transverse sections, whereas b\-b4 characterized the lengths of the minor axes. The degree of the curve segments (values of p) on the medial side were designated as p 1, /;3, /;5, and p i and on the lateral side as p2, p4, p6, and p8. The complete set of geometric constraints that belonged to the admissible family is defined as follows: o2 —a 1 > 0 b2-b\ >0 a4 —a2 > 0 b3 —b2 > 0 a6 —a4 > 0 M - b3 > 0 a3 - a 2 -4 .5 > 0 a5 —a3 - 4.5 > 0 a5 —a3 —7.0 < 0, and a2 + a5 —2.0 x a3 > 0
(2)
Although the present parametric modeling and selection of con straints were, by and large, influenced by the scheme proposed by Ruben et al. [17], more shapes were explored in the present scheme by introducing an increased number of design variables in the cad routine (Fig. 3). The nature of stem transverse sections was varied from medial to lateral side by changing the degree of
Fig. 2 Design parameters used in the study: (a) eighteen design variables characterizing four key sections of the implant and ( b ) initial implant model based on a collarless TriLock (DePuy) prosthesis
034502-2 / Vol. 137, MARCH 2015
Transactions of the ASME
Fig. 3 Probable shapes of the implant key sections: a compari son between the two parameterization schemes
the cross-sectional curve (parameter p) on both sides of the stem after maintaining continuity at the intersection of both the seg ments (Fig. 2(a)). Considering all these cross-sectional features was important for the optimum design, since a combination of dif ferent features would provide the best outcome [18]. The values of p were varied between 2 and 5 and only integer values were considered. This gave rise to a wide range of shapes starting from ellipse (p = 2) to super-ellipses of different degrees (p > 2). 2.2
Development of FE Models
2.2.1 Bone Geometry and FE Meshing. The 3D FE models of right proximal femur (intact and implanted) were developed using CT scan data of a 31 yr old male patient and ansys v14.0 software (ANSYS, Inc.. PA). Images of the femur were stored in 512x512 pixels, with a pixel size of 0.799 mm, slice thickness of 0.699 mm. The CT gray value was expressed in terms of Hounsfield units (HU). The CT-scan data were processed using a medical image processing program, MIMICS (Materialise, Leuven, Belgium) to acquire the 3D surface geometry of the fe mur. Within a CT-scan slice, segmentation of an image was per formed to identify the bone from soft tissue using threshold values of CT number HU. ahductors tensorfasciae (atae iCiotihiaCtract
A ten-node tetrahedral element was used to mesh both intact and implanted models, since its convergence rate was sensibly higher than a mesh of four-node tetra-elements; and the mesh technique was reported to produce best results, especially when a solid model was available [19]. Other than the numerical accu racy, quadratic tetra mesh is particularly useful for modeling com plex femur geometry and also for geometries that offer wide ranges of surface and boundary curvatures [20,21 ], as in the pres ent parametric scheme. The element edge length for both implanted and intact FE models ranged between 0.2 and 1.0mm. The present study employed an automatic mesh generator capa ble of generating elaborate and well-conditioned volume mesh, based on fully defined solid models of the implant, while main taining topology at the interface and edges. One-to-one correspon dence between the FE models and the design variables was avoided, which is an essential requirement for solving any shape optimization problem [22]. A finer mesh was employed at the implant-bone interface for this purpose (Fig. 4). A pilot study on mesh convergence was carried out following recommendations from earlier study [23] in order to ensure that the accuracy of FE analyses remained adequate throughout the entire process of optimization. 2.2.2 Material Properties. Bone was assumed to be a linearly elastic and isotropic material. Elastic modulus of each bone ele ment was extracted from the CT-scan data using shareware soft ware bonemat v2.0, developed at the Institute Ortopedici Rizzoli, Bologna, Italy [24,25]. The apparent density (p in g cm- '1) for bone element was computed from corresponding CT number (HU) in Hounsfield Unit using the following relationship: p = p, + [(Pa —P,)/(HU2 —HU,)] x]H U -H U ,]
(3)
where (ph HU|) represented the apparent density and CT number of water, i.e., no bone condition (0.022g cm-1, 0 HU), and (p2, HU2) were the apparent density in the hardest cortical region and corre sponding CT number value (1.73 g cm"'1, 1700 HU), respectively.
flip contact
implant
femur
vastus (ate rafts
Finer mesh at the implant-bone interface
L
■vastus m afia (is
Sectional view of implanted FE model
Fig. 4 Finite element model of the implanted femur subjected to musculoskeletal loading conditions of normal walking and stair climbing. Finer mesh density was used at the implant-bone interface as indicated in the sectional view of the FE model. For magnitude of loading, see Table 1.
Journal of Biomechanical Engineering
MARCH 2015, Vol. 137 / 034502-3
The Young’s moduli (£ in MPa) of bone elements were related to the apparent density (p in g-cm~3) as follows: £ = 728 Ip 1'52
(4)
The values of constants in Eq. (4) correspond to the table men tioned in Morgan et al. [26]. Young’s modulus for the implant ma terial (Titanium alloy) was taken as I lOGPa [27], Poisson’s ratio for all the materials, including bone tissues, was considered as 0.3 [28], The implant-bone interface was considered to be perfectly bonded for the implanted FE model (Fig. 4). 2.2.3 Applied Loading Conditions. Loading conditions for the FE models included the maximum loads during stance phase of normal walking and stair climbing, which were applied as two static load cases [29,30], A combined effect of the two load cases is considered, by giving equal weightage to each case. Typical values of the force system are presented in Table 1. For the intact femur model, the joint reaction force was uniformly distributed over an approximately circular area of diameter 30 mm [31]. In case of the implanted FE model, the joint reaction force was pro jected onto the centre of the femoral head. The muscle forces were distributed over a circular patch, centred on the point of attachment of the respective muscles or group of muscles in the femur. Although there exist interindividual variations in musculo skeletal system, data on area of muscle insertions were taken from Duda et al. [32]. The muscle forces and hip contact force were cal culated assuming a moderate body weight of 70 kg, which approx imately represents the weight of average population in the same age group of the subject. Nodes located on the distal most end of the femur, about 150 mm below the lesser trochanter, were fixed in all directions (Fig. 4).
in local mechanical stimulus SED between the intact and the implanted bone. The reference stimulus (Srcf in J-g- 1) and the actual stimulus (S in J g- ) were the local (per element) elastic strain energy (U in MPa) per unit bone mass averaged over a loading history (n), for an intact and implanted femur, respectively. The remodeling stim uli, i.e„ the SED values were calculated based on average strain energy (Ua in MPa), estimated over the two (n = 2) load cases (walking and stair climbing), divided by the apparent density values (p in g em- '1), mathematically expressed as [16]
If the difference between S and Srcf, is less than the dead zone (s), the bone resorption would be initiated [34,36], The reference stim ulus Srcro f each bone element was compared with the correspond ing actual stimulus S of implanted bone element at a comparable location. It should be noted that the intact and the implanted femur models contained dissimilar meshes. In order to map the bone ele ment of the implanted model on a bone element of the intact model, a matlab script was written. Based on the centroid location of a bone element of the implanted femur model, this script deter mined the nearest element of the intact femur model using mini mum Eucledian distance criterion [37,38], The difference in mean value of centroidal distance between the intact and the implanted femur was calculated as 0.32 mm (range: 0.0024-0.57 mm), which was less than the pixel size of 0.799 mm. The validity of this tech nique depends on the bone element size, which was kept very small (maximum 1.0 mm) for both the intact and the implanted model. A function K(U,p) was introduced, the values of which were set as follows:
2.3 Formulation of the Objective Functions. Two objective functions, one based on proximal resorbed BMF, and the other based on implant-bone interface failure criterion, were formulated as described in Secs. 2.3.1 and 2.3.2. 2.3.1 Objective Function fo r Bone Resorption. The objective function, considered in the study to assess the effect of implantinduced bone resorption, was expressed in terms of percentage resorption of BMF in immediate postoperative case [8,33], Since bone resorption predominantly occurs in the proximal part of the femur [34], the calculations for BMF were confined to the bone elements located in the proximal femur (intact and implanted) upto the lesser trochanter. The calculation for the BMF was based on the adaptive bone remodeling theory [34,35]. According to this theory, the change in bone density is dependent on the difference T a b le 1
K (U ,p)
1, 0,
if S < (1 —sjSref otherwise
(6)
The dead zone (,v) was considered as ±0.75 of Srcr [35], Thus, the objective function for resorbed BMF {%) was written as
/l=i“ J M
K(U ,p)pdV
(7)
where M (gm) and V (cm3) represent the total mass and volume of the implanted proximal femur, respectively. 2.3.2 Objective Function fo r interface Stresses. Formulation of the other objective function related to interface stress was based
A p p lie d fo r c e s , e x p r e s s e d in p e r c e n ta g e o f b o d y w e ig h t c o r r e s p o n d in g to lo a d c a s e s , n o r m a l w a lk in g , a n d s t a ir c lim b in g .
F o r c e c o m p o n e n ts , F„, Fy, a n d F& c o r r e s p o n d to x , y a n d z d ir e c tio n s , r e s p e c tiv e ly , in th e g lo b a l c o o r d in a t e s y s te m a s s h o w n in F ig . 4. T h e a b b r e v ia tio n m . s ta n d s fo r m u s c le . D a ta a d o p te d fr o m B e r g m a n n e t al. [2 9 ].
Force (percentage of body weight) Load cases
Type of force
Fy
Fy
F.
Resultant F
Normal walking
hip contact m. abductor m. tensor fascia lata, proximal part m. tensor fascia lata, distal pan m. vastus lateralis
-5 4 .0 58.0 7.2 -0 .5 -0 .9
32.8 -4 .3 -1 1 .6 0.7 -18.5
-229.2 86.5 13.2 -1 9 .0 -9 2 .9
237.7 104.2 18.9 19.0 94.7
Stairs climbing
hip contact m. abductor iliotibial tract, proximal part iliotibial tract, distal part m. tensor fascia lata, proximal part m. tensor fascia lata, distal part m. vastus lateralis m. vastus medial is
-5 9 .3 70.1 10.5 -0 .5 3.1 -0 .2 -2 .2 8.8
60.6 -28.8 3.0 0.8 -4 .9 0.3 -2 2 .4 -3 9 .6
-236.3 84.9 12.8 -16.8 2.9 -6 .5 -135.1 -267.1
251.0 1 13.8 16.8 16.8 6.4 6.5 136.9 270.2
034502-4 / Vol. 137, MARCH 2015
Transactions of th e A SM E
Fig. 5
A schematic representation of the optimization procedure
on the multiaxial failure theory known as Hoffman’s criterion [39], This theory relates the interface stresses with the strength of the bone element adjacent to the interface and ascertains the area with likely chances of debonding [8,33,36,40]. The Hoffman fail ure index (FL) was calculated at each interfacial nodal location using
^
S,SCff +
s,
( 8)
where
Sanjay Gupta1 Department of Mechanical Engineering, Indian Institute of Technology Kharagpur, Kharagpur, West Bengal 721 302, India e-mail: [email protected]
Dilip Kumar Pratihar Department of Mechanical Engineering, Indian Institute of Technology Kharagpur, Kharagpur, West Bengal 721 302, India
The shape and geometry o f femoral implant influence implantinduced periprosthetic bone resorption and implant-bone inter face stresses, which are potential causes o f aseptic loosening in cementless total hip arthroplasty (THA). Development o f a shape optimization scheme is necessary to achieve a trade-off between these two conflicting objectives. The objective o f this study was to develop a novel multi-objective custom-based shape optimization scheme for cementless femoral implant by integrating finite ele ment (FE) analysis and a multi-objective genetic algorithm (GA). The FE model o f a proximal femur was based on a subjectspecific CT-scan dataset. Eighteen parameters describing the na ture o f four key sections o f the implant were identified as design variables. Two objective functions, one based on implant-bone interface failure criterion, and the other based on resorbed proximal bone mass fraction (BMF), were formulated. The results predicted by the two objective functions were found to be contra dictory; a reduction in the proximal bone resorption vras accom panied by a greater chance o f inteiface failure. The resorbed proximal BMF was found to be between 23% and 27% fo r the trade-off geometries as compared to ~39% fo r a generic implant. Moreover, the overall chances o f interface failure have been mini mized fo r the optimal designs, compared to the generic implant. The adaptive bone remodeling was also found to be minimal fo r the optimally designed implants and, further with remodeling, the chances o f interface debonding increased only marginally. [DOI: 10.1115/1.4029061] Keywords: total hip arthroplasty, finite element analysis, bone remodeling, inteiface failure, genetic algorithms, shape optimization, multi-objective
1
Introduction
The adverse effects of implant-induced bone resorption and ex cessive interface stresses have often been reported as potential causes of failure in cementless THA [1,2]. The success of THA notwithstanding, about 10% of the hip replacements need revision surgery depending on the patient’s conditions, disorder, and the type of implants used [3,4], It is well known that the shape and
geometry of a cementless femoral implant influence periprosthetic bone remodeling and implant-bone interface stresses [5,6]. The mechanical designs of hip prosthesis, therefore, need to search for optimal implant geometry, such that the effects of stress shielding induced bone resorption and interface stresses are minimized. Shape (or geometry) optimization is a particular stage of struc tural optimization, which searches for the optimal configuration of a design domain [7]. Employing shape optimization as a design tool, stem geometries can be evaluated based on some biomechan ical cost (or objective) functions. Therefore, the formulation of suitable cost functions, representing the global effects of these failure mechanisms is necessary. These cost functions may some times be mutually conflicting in nature [2,8]. Therefore, develop ment of a multi-objective shape optimization scheme is necessary to achieve a trade-off between these two conflicting objectives, while minimizing them within the domain of feasibility. Despite a few earlier studies on single-objective optimization [9-12], the hrst major study on multi-objective shape optimization was carried out by Fernandes et al. [13] to analyse the design pros pects of a 2D cementless hip stem model. Thereafter, the hrst GA based multi-objective shape optimization strategy, to minimize proximal and distal stress in cement mantle of a cemented THA. was proposed by Ishida et al. [14], In a more recent study, Ruben et al. [15] presented a 3D shape optimization scheme based on a multicriterion objective function to design implants with an improved initial stability and reduced proximal bone resorption. Although several shape optimization schemes were employed and new stem designs were proposed, these studies were based on generic femoral bone model, consisting of a cortex layer covering a homogeneous core of cancellous bone, thereby neglecting the heterogeneous bone material property across the proximal femur [12-15], The study by Katoozian and Davy [10] attempted to model bone heterogeneity by linearly varying bone material prop erty from cancellous to cortical bone. However, such method also suffers from oversimplification of bone property. Moreover, these studies suffered from additional limitations regarding the choice of objective function for bone remodeling simulation. Most stud ies inappropriately considered strain energy density (SED) of the implanted bone only as the mechanical stimulus for bone remod eling [13,15], instead of the difference in SED between intact and implanted cases [16]. The parameterization schemes used in earlier studies for defin ing 3D implant geometry were unable to account for a wide vari ety of implant shapes [10,14,15]. The search for an optimal geometry demands extensive exploration of all possible and ad missible shapes. However, this may add to more complexities which may not be suitably handled by traditional optimization tools. Moreover, shape optimizations problems were reported to be typically nonconvex and characterized by multiple local optima [7], which requires a robust heuristic search method, such as the GA, for its complex multivariable analysis capability and gradient-independent approach. Thus, the GA seems to be a more appropriate tool for solving such multimodal and multivariable space optimization problems. Nevertheless, hardly any such study has been attempted in the field of shape optimization of cementless femoral implant. The objective of this study was to develop a novel multi objective shape optimization scheme for cementless femoral implant by integrating FE analysis and GA. The study was based on a subject-specific FE model; in which heterogeneous bone material properties were assigned element-wise using the CT-scan dataset. Geometric parameters were introduced to char acterize the femoral stem shape. The objective functions were formulated based on common failure mechanisms of cementless THA.
2 'Corresponding author. Manuscript received March 25. 2014: final manuscript received November 5, 2014; published online January 29, 2015. Assoc. Editor: David Corr.
Journal of Biomechanical Engineering
Materials and Methods
The present optimization scheme comprised of a parametric implant model generator, a script for virtual implantation, an
Copyright © 2015 by ASME
MARCH 2015, Vol. 137 / 034502-1
V
Fig. 1 Parametric scheme for evolution of key sections of the implant: ellipse (p = 2) to superellipse (p >2)
automatic mesh generator, an FE analysis script and finally an optimization controller. In order to integrate all these components, a master program was written in m atla b . The master program served as the prime mover and could invoke all components logi cally and sequentially via batch mode processing. An elaborate description of each method that collectively contributed toward development of the optimization scheme is stated in Secs. 2.1 to 2.4. 2.1 Parameterization and Constraint Specification of Implant Geometry. In a local coordinate system ary, the pair of equations that defined the geometric parameterization of stem transverse section was expressed by the parametric form as follows: .v = a cos(2/,/'' t y = b sin*2|/,,) t
0>
where the limits of t were prescribed as 0 < r < zr/2. The parame ters a, b, and p at four key sections, thus, characterized the design parameters of the femoral stem (Fig. 1). Using these four key sec tions, a virtual 3D parametric model of the implant was generated
in SolidWorks cad package (SolidWorks, Dassault Systemes S.A.) by executing an Application Programming Interface script written for the purpose. The initial design of the stem was similar to a generic design of a collarless TriLock (DePuy) prosthesis having stem-length 113 mm (Fig. 2(a)). The stem-length was kept unaltered throughout the optimization procedure. A total of 18 geometric parameters was introduced as design variables for the modeling purpose (Fig. 2(a)). These design varia bles had lower and upper bounds and geometric constraints in order to achieve clinically admissible shapes. Furthermore, it was also required that the width and thickness of cross sections would decrease from proximal to distal end in order to facilitate ease of withdrawal, if required, in case of a revision surgery. The parame ters al-a6 corresponded to the length of major axes on the medial and lateral sides of the four transverse sections, whereas b\-b4 characterized the lengths of the minor axes. The degree of the curve segments (values of p) on the medial side were designated as p 1, /;3, /;5, and p i and on the lateral side as p2, p4, p6, and p8. The complete set of geometric constraints that belonged to the admissible family is defined as follows: o2 —a 1 > 0 b2-b\ >0 a4 —a2 > 0 b3 —b2 > 0 a6 —a4 > 0 M - b3 > 0 a3 - a 2 -4 .5 > 0 a5 —a3 - 4.5 > 0 a5 —a3 —7.0 < 0, and a2 + a5 —2.0 x a3 > 0
(2)
Although the present parametric modeling and selection of con straints were, by and large, influenced by the scheme proposed by Ruben et al. [17], more shapes were explored in the present scheme by introducing an increased number of design variables in the cad routine (Fig. 3). The nature of stem transverse sections was varied from medial to lateral side by changing the degree of
Fig. 2 Design parameters used in the study: (a) eighteen design variables characterizing four key sections of the implant and ( b ) initial implant model based on a collarless TriLock (DePuy) prosthesis
034502-2 / Vol. 137, MARCH 2015
Transactions of the ASME
Fig. 3 Probable shapes of the implant key sections: a compari son between the two parameterization schemes
the cross-sectional curve (parameter p) on both sides of the stem after maintaining continuity at the intersection of both the seg ments (Fig. 2(a)). Considering all these cross-sectional features was important for the optimum design, since a combination of dif ferent features would provide the best outcome [18]. The values of p were varied between 2 and 5 and only integer values were considered. This gave rise to a wide range of shapes starting from ellipse (p = 2) to super-ellipses of different degrees (p > 2). 2.2
Development of FE Models
2.2.1 Bone Geometry and FE Meshing. The 3D FE models of right proximal femur (intact and implanted) were developed using CT scan data of a 31 yr old male patient and ansys v14.0 software (ANSYS, Inc.. PA). Images of the femur were stored in 512x512 pixels, with a pixel size of 0.799 mm, slice thickness of 0.699 mm. The CT gray value was expressed in terms of Hounsfield units (HU). The CT-scan data were processed using a medical image processing program, MIMICS (Materialise, Leuven, Belgium) to acquire the 3D surface geometry of the fe mur. Within a CT-scan slice, segmentation of an image was per formed to identify the bone from soft tissue using threshold values of CT number HU. ahductors tensorfasciae (atae iCiotihiaCtract
A ten-node tetrahedral element was used to mesh both intact and implanted models, since its convergence rate was sensibly higher than a mesh of four-node tetra-elements; and the mesh technique was reported to produce best results, especially when a solid model was available [19]. Other than the numerical accu racy, quadratic tetra mesh is particularly useful for modeling com plex femur geometry and also for geometries that offer wide ranges of surface and boundary curvatures [20,21 ], as in the pres ent parametric scheme. The element edge length for both implanted and intact FE models ranged between 0.2 and 1.0mm. The present study employed an automatic mesh generator capa ble of generating elaborate and well-conditioned volume mesh, based on fully defined solid models of the implant, while main taining topology at the interface and edges. One-to-one correspon dence between the FE models and the design variables was avoided, which is an essential requirement for solving any shape optimization problem [22]. A finer mesh was employed at the implant-bone interface for this purpose (Fig. 4). A pilot study on mesh convergence was carried out following recommendations from earlier study [23] in order to ensure that the accuracy of FE analyses remained adequate throughout the entire process of optimization. 2.2.2 Material Properties. Bone was assumed to be a linearly elastic and isotropic material. Elastic modulus of each bone ele ment was extracted from the CT-scan data using shareware soft ware bonemat v2.0, developed at the Institute Ortopedici Rizzoli, Bologna, Italy [24,25]. The apparent density (p in g cm- '1) for bone element was computed from corresponding CT number (HU) in Hounsfield Unit using the following relationship: p = p, + [(Pa —P,)/(HU2 —HU,)] x]H U -H U ,]
(3)
where (ph HU|) represented the apparent density and CT number of water, i.e., no bone condition (0.022g cm-1, 0 HU), and (p2, HU2) were the apparent density in the hardest cortical region and corre sponding CT number value (1.73 g cm"'1, 1700 HU), respectively.
flip contact
implant
femur
vastus (ate rafts
Finer mesh at the implant-bone interface
L
■vastus m afia (is
Sectional view of implanted FE model
Fig. 4 Finite element model of the implanted femur subjected to musculoskeletal loading conditions of normal walking and stair climbing. Finer mesh density was used at the implant-bone interface as indicated in the sectional view of the FE model. For magnitude of loading, see Table 1.
Journal of Biomechanical Engineering
MARCH 2015, Vol. 137 / 034502-3
The Young’s moduli (£ in MPa) of bone elements were related to the apparent density (p in g-cm~3) as follows: £ = 728 Ip 1'52
(4)
The values of constants in Eq. (4) correspond to the table men tioned in Morgan et al. [26]. Young’s modulus for the implant ma terial (Titanium alloy) was taken as I lOGPa [27], Poisson’s ratio for all the materials, including bone tissues, was considered as 0.3 [28], The implant-bone interface was considered to be perfectly bonded for the implanted FE model (Fig. 4). 2.2.3 Applied Loading Conditions. Loading conditions for the FE models included the maximum loads during stance phase of normal walking and stair climbing, which were applied as two static load cases [29,30], A combined effect of the two load cases is considered, by giving equal weightage to each case. Typical values of the force system are presented in Table 1. For the intact femur model, the joint reaction force was uniformly distributed over an approximately circular area of diameter 30 mm [31]. In case of the implanted FE model, the joint reaction force was pro jected onto the centre of the femoral head. The muscle forces were distributed over a circular patch, centred on the point of attachment of the respective muscles or group of muscles in the femur. Although there exist interindividual variations in musculo skeletal system, data on area of muscle insertions were taken from Duda et al. [32]. The muscle forces and hip contact force were cal culated assuming a moderate body weight of 70 kg, which approx imately represents the weight of average population in the same age group of the subject. Nodes located on the distal most end of the femur, about 150 mm below the lesser trochanter, were fixed in all directions (Fig. 4).
in local mechanical stimulus SED between the intact and the implanted bone. The reference stimulus (Srcf in J-g- 1) and the actual stimulus (S in J g- ) were the local (per element) elastic strain energy (U in MPa) per unit bone mass averaged over a loading history (n), for an intact and implanted femur, respectively. The remodeling stim uli, i.e„ the SED values were calculated based on average strain energy (Ua in MPa), estimated over the two (n = 2) load cases (walking and stair climbing), divided by the apparent density values (p in g em- '1), mathematically expressed as [16]
If the difference between S and Srcf, is less than the dead zone (s), the bone resorption would be initiated [34,36], The reference stim ulus Srcro f each bone element was compared with the correspond ing actual stimulus S of implanted bone element at a comparable location. It should be noted that the intact and the implanted femur models contained dissimilar meshes. In order to map the bone ele ment of the implanted model on a bone element of the intact model, a matlab script was written. Based on the centroid location of a bone element of the implanted femur model, this script deter mined the nearest element of the intact femur model using mini mum Eucledian distance criterion [37,38], The difference in mean value of centroidal distance between the intact and the implanted femur was calculated as 0.32 mm (range: 0.0024-0.57 mm), which was less than the pixel size of 0.799 mm. The validity of this tech nique depends on the bone element size, which was kept very small (maximum 1.0 mm) for both the intact and the implanted model. A function K(U,p) was introduced, the values of which were set as follows:
2.3 Formulation of the Objective Functions. Two objective functions, one based on proximal resorbed BMF, and the other based on implant-bone interface failure criterion, were formulated as described in Secs. 2.3.1 and 2.3.2. 2.3.1 Objective Function fo r Bone Resorption. The objective function, considered in the study to assess the effect of implantinduced bone resorption, was expressed in terms of percentage resorption of BMF in immediate postoperative case [8,33], Since bone resorption predominantly occurs in the proximal part of the femur [34], the calculations for BMF were confined to the bone elements located in the proximal femur (intact and implanted) upto the lesser trochanter. The calculation for the BMF was based on the adaptive bone remodeling theory [34,35]. According to this theory, the change in bone density is dependent on the difference T a b le 1
K (U ,p)
1, 0,
if S < (1 —sjSref otherwise
(6)
The dead zone (,v) was considered as ±0.75 of Srcr [35], Thus, the objective function for resorbed BMF {%) was written as
/l=i“ J M
K(U ,p)pdV
(7)
where M (gm) and V (cm3) represent the total mass and volume of the implanted proximal femur, respectively. 2.3.2 Objective Function fo r interface Stresses. Formulation of the other objective function related to interface stress was based
A p p lie d fo r c e s , e x p r e s s e d in p e r c e n ta g e o f b o d y w e ig h t c o r r e s p o n d in g to lo a d c a s e s , n o r m a l w a lk in g , a n d s t a ir c lim b in g .
F o r c e c o m p o n e n ts , F„, Fy, a n d F& c o r r e s p o n d to x , y a n d z d ir e c tio n s , r e s p e c tiv e ly , in th e g lo b a l c o o r d in a t e s y s te m a s s h o w n in F ig . 4. T h e a b b r e v ia tio n m . s ta n d s fo r m u s c le . D a ta a d o p te d fr o m B e r g m a n n e t al. [2 9 ].
Force (percentage of body weight) Load cases
Type of force
Fy
Fy
F.
Resultant F
Normal walking
hip contact m. abductor m. tensor fascia lata, proximal part m. tensor fascia lata, distal pan m. vastus lateralis
-5 4 .0 58.0 7.2 -0 .5 -0 .9
32.8 -4 .3 -1 1 .6 0.7 -18.5
-229.2 86.5 13.2 -1 9 .0 -9 2 .9
237.7 104.2 18.9 19.0 94.7
Stairs climbing
hip contact m. abductor iliotibial tract, proximal part iliotibial tract, distal part m. tensor fascia lata, proximal part m. tensor fascia lata, distal part m. vastus lateralis m. vastus medial is
-5 9 .3 70.1 10.5 -0 .5 3.1 -0 .2 -2 .2 8.8
60.6 -28.8 3.0 0.8 -4 .9 0.3 -2 2 .4 -3 9 .6
-236.3 84.9 12.8 -16.8 2.9 -6 .5 -135.1 -267.1
251.0 1 13.8 16.8 16.8 6.4 6.5 136.9 270.2
034502-4 / Vol. 137, MARCH 2015
Transactions of th e A SM E
Fig. 5
A schematic representation of the optimization procedure
on the multiaxial failure theory known as Hoffman’s criterion [39], This theory relates the interface stresses with the strength of the bone element adjacent to the interface and ascertains the area with likely chances of debonding [8,33,36,40]. The Hoffman fail ure index (FL) was calculated at each interfacial nodal location using
^
S,SCff +
s,
( 8)
where
