Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 47, 3, pp. 573-598, Warsaw 2009 OPTIMIZATION OF DENTAL IMPLANT USING GENETIC ALGORITHM Tomasz Łodygowski Krzysztof Szajek Marcin Wierszycki Poznań University of Technology, Institute of Structural Engineering, Poznań, Poland e-mail: tomasz.lodygowski@put.poznan.pl The subject of the present work is optimization of the modern implant systemOsteoplant, which was created and is still developed by Founda- tion of University of Medical Sciences in Poznań. Clinical observations point to the occurrence of both early and late complications in the case of all two-component implant systems. In many cases, these problems are caused by mechanical fractures of the implants themselves. The ob- tained results of the previous studies focused on necessary changes of the implant mechanical behavior, which helped to achieve the required long-termstrength.However,modifications of the presentdental implant system are not obvious. In this paper, an optimization of theOsteoplant dental implant system, with the use of FEA and genetic algorithms is discussed. Key words: design optimization, genetic algorithm, dental implant 1. Introduction The use of implants is the commonly applied treatment method of dental re- storations. Themodern implant systemOsteoplant was created and has been developed by Foundation of University of Medical Sciences in Poznań since over 10 years. TheOsteoplant is a two-component implant system. It consists of the abutment and root, which are connected by a titanium screw. Clini- cal observations point to the occurrence of both early and late complications in the case of all two-component implant systems (Goodacre et al., 2003). In many cases, these problems are caused by mechanical fractures of the im- plants themselves (Hędzelek et al., 2004; Kąkol et al., 2002; Zagalak et al., 574 T. Łodygowski et al. 2005). One of the most dangerous complications is fracture and cracking of thedental implant parts.Damage of dental implant componentsmakes further treatment very difficult. The previous studies (Wierszycki, 2007; Wierszycki et al., 2006a,b,c) confirmed fatigue as a reason of implant damage. The obta- ined results of previous studies focused on necessary changes of the implant mechanical behavior, which helped to achieve the required long-term strength. However, modifications of the present dental implant system are not obvious. To find out a better design of the implant system, an optimization procedures based onFEAmust be used. In this paper, the optimization of theOsteoplant dental implant system, with the use of FEA, genetic algorithms is discussed. Fig. 1. Dental implant Osteoplant Genetic algorithms (GAs) have received wide popularity as optimization techniques during the last decades, in particular, for very complex designs and they can successfully compete with gradient-based approaches in many areas (Goldberg, 1989). GAs are stochastic search approaches which rely on the principle of survival of the fittest in natural selection. Unlike conventional optimization techniques, GAs explore simultaneously the entire design space and, therefore, are likely to reach the global minimum. An improvement of the global search process can be performed by incorporating neural networks (NN) inoptimization,which can learnandadapt changes over time. Ingeneral, GAs require a lot computations (structural analyses in our case) and hence high performance computing ideally addresses their needs, especially when combined with NN. In the created tool, the existing open source libraries have been used: Galileo (for GA) and ffnet (for NN). The whole optimization procedure was implemented with the use of Python scripting language. The FE model, nu- merical analysis and post-processing of results were performed with the use of the Abaqus Unified FEA product suite from SIMULIA. The integration of the GA and NN libraries with the FE tools was done by using the Abaqus Scripting Interface (ASI). Optimization of dental implant... 575 2. Model The learning of a genetic algorithm and neutral network needs a huge number of analyses to be carried out. For this reason, the crucial characteristic of the numerical model of implant, which is used in optimization process, is performed. In practice, the time of calculation can not exceed several dozen minutes. This limitation causes that a fully three-dimensional model of an implant and typical modeling approaches can not be used (Wierszycki et al., 2006a). For this study, a special approachwas proposed.Amodeling approach described in detail below enables us to carry out a nonlinear, 3D simulation of a dental implant in an acceptable time with a satisfactory level of accuracy of the results. Fig. 2. Numerical models of the implant: axisymmetric (a) and 3D (b) A three-dimensional model of the implant has been created by revolving an axisymmetric model about its axis of symmetry. The symmetric model generation capability of Abaqus/Standard enables one to create automatically a three-dimensional model (Abaqus..., 2007a). The nodes, elements, section definitions, material and contact definitions of the three-dimensional model are created automatically based on the axisymmetric model description. Only kinematic constraints and boundary conditions must be redefined. In order to reduce the time of calculation, the asymmetric deformation of the three- dimensional model was assumed to be symmetric with respect to the radial – symmetry axis plane at an angle equal to 0 or 180◦. The symmetric results transfer capability of the program and it was used to transfer the results from the axisymmetric simulation of the assembly to the final three-dimensional model (Abaqus..., 2007a). 576 T. Łodygowski et al. 2.1. Geometry The geometry of the numerical model was simplified to the axisymmetric description (Fig.2). The internal threads of the implant body and screwwere simplified to axisymmetric, parallel rings. Because themain goal of this study is the optimization of the screw connection, the external thread of implant bo- dywas omitted.TheparametricAbaqus/CAEmodel of the implant consists of three axisymmetric parts. The shape of each of them corresponds to the cross- section of dental implant components: abutment, body and screw. The 2D sketches of this parts have fully parametric geometry description and are fully constrained. These constraints with the dimensions and parametric equations added to the sketch enable us to automatically modify the shape of implant components (Abaqus..., 2007b). The six global independent parameters were defined: • screw head diameter, • screw head conic surface opening angle, • screw head height, • hexagonal slot diameter, • hexagonal slot height, • hexagonal slot conic surface opening angle. All parts share the same parameters, so the instances of the parts are always consistent. The parameters are further coded in GA as gens (see Sec- tion 3.2.2.). The geometry of the three-dimensional model of the implant was not defi- ned directly. The FE three-dimensional model is automatically created based on the axisymmetric model (see Section 3.3.4). 2.2. Material All components of an implant are made of medical alloys of titanium. For general stress-strain analyses, isotropic and non-linear elastic-plastic charac- teristics of material models were taken into account. The material properties were based on the Certificates of Conformity and on literature (Wang, 1996). The mechanical properties of titanium alloys are shown in Table 1. For fati- gue calculations, the model of material has to be simplified to a linear elastic description. The material models and characteristic geometry definitions are automa- tically transferred from the axisymmetric to three-dimensional model during realization of the symmetric model generation procedure (Abaqus..., 2007a). Optimization of dental implant... 577 Table 1.Mechanical properties of implant materials Implant body Abutment Screw Young’s modulus [MPa] 105200 105200 105200 Poisson ratio 0.19 0.19 0.19 Yield [MPa] 615.2 832.3 802.8 Tensile Yield [MPa] 742.4 1004.0 970.4 2.3. Assembly and contact The assembly is done at the first axisymmetric stage of creation of the implant model. The three two-dimensional instances of implant parts were positioned relative to each other in a global coordinate system. Relative posi- tion constraints were applied to align with: • conic surface of hexagonal slot and outside of abutment, • conic surface of screw head and inside of abutment. Fully relative position constraints of implant part instances make possible automatic redefinition of the implantmodel assembly (Abaqus..., 2007a). The assembling simulation involves solving a contact problem (Wierszycki et al., 2006b). For this purpose, it is necessary to define three contact areas, between: • root and abutment, • root and screw, • abutment and screw. This contact conditions produce typical assembly problems, so we decided to use a standard small-sliding contact formulation. In order to minimize the dependence onmesh density, the surface-to-surface contact discretization was used. Significant penetrations of master nodes into the slave surfaces do not occur with the used space discretization. Moreover, surface-to-surface discre- tization provides more accurate stress and pressure results, especially in the cases of contacts at corners like on the threads. In some situations, the surface- to-surface discretization approach can generate additional solution costs (Aba- qus..., 2007a). The reason for this is that more nodes per each constraint are involved. In this particular application, the surface-to-surface discretization approach significantly reduces the solution cost. For a three-dimensional mo- del with two cylindrical elements in 180◦ segment, the node-to-surface di- scretization causes a double increase in the number of increments and sever discontinue iterations. There are no significant differences in time of calcula- tion depending on themethod of contact constraint enforcement. The penalty 578 T. Łodygowski et al. method as the contact constraint enforcement method was selected for both normal and tangential directions. Tangential surface behavior was defined as the classical isotropic Coulomb friction model. The friction coefficient is the same for all contact pairs and amounts to 0.19. The definitions of contact pa- irs are automatically transferred from the axisymmetric to three-dimensional model in the symmetric model generation procedure (Abaqus..., 2007a). 2.4. Loads The loading of the implant model is a two-step process. The first step corresponds to simulation of a tightening process. In this step, bothmodel and its response are axisymmetric. This simulation can be carried outwith the use of the axisymmetric model. The second step is bending, which is caused by the worst component of service load, perpendicular to the axisymmetric axis. In this case, themodel which can describe asymmetric deformation is needed. For a two-component implant, one of the most crucial aspects of numeri- cal modeling is simulation of the mechanical assembly subject to tightening of the implant screw (Lang et al., 2003; Wierszycki et al., 2006c). For the axisymmetric or a simplified three-dimensional model, tightening simulations cannot be performed as a real physical process. A work-around of this ap- proach is necessary. For simulation of tightening, a prescribed assembly load has been used. The middle part of the screw was defined as pre-tension sec- tion. The tightening load (150N) was applied to the pre-tension section as a concentrated load. During calculation, the screw length was reduced in this pre-tension section to achieve the assumed tightening force. The implant body and abutment were tightened as results of the change of screw length. The value of axial force in the tightened screw was calculated from the empirical formulation (Bozkaya andMüfüt, 2005; Lang et al., 2003;Merz et al., 2000). It depends on the friction coefficient and torque. This assumption and procedure of tightening were verified during full simulation of screw tightening with the help of a fully three-dimensional FE model of an implant (Wierszycki et al., 2006c). The results obtained in the axisymmetric simulation were transfered into the final three-dimensional model. The second step was bending of the tightened implant. The bending force was applied to the tip of abutment by means of the concentrated load (10N) perpendicular to the axisymmetric axis of the implant. 2.5. Mesh For all parts of the model, a quad-dominated shape of elements was used. At each edge of the mesh region, the seed (approximated node locations) has Optimization of dental implant... 579 beendefinedto controlmeshdensity.Correctmeshdensitymustbeensured for different geometry configurations. The seeds have been defined by specifying average element sizes along the edges. To ensure the proper mesh density for different geometry configurations, the seed densities have been partially constrained. This approach ensured that even if the number of elements along the edge was changed its size remained such as had been defined (Abaqus..., 2007b). The mesh of the three-dimensional model was generated automatically in the symmetricmodel generation procedure. In wholemodel of implant, CCL9 and CCL12 elements were used (Abaqus..., 2007a). The cylindrical elements available in Abaqus/Standard use standard isoparametric interpolation in the radial-symmetry axis plane, combined with trigonometric interpolation func- tions with respect to the angle of revolution. Cylindrical elements can span angles between 0 and 180◦. The cost of calculation was significantly reduced by this fact itself. The number of cylindrical elements along the circumferen- tial direction is the compromise between time of calculation and accuracy of the results. Five tests were carried out to evaluate which number of elements is optimal. Themodels with three-dimensional solid element C3D8 andC3D6 were used as reference solutions. In the models, 32 and 16 elements per 180◦ segment were used. The maximum values of the equivalent Mises stress at characteristic notches and global bending stiffness of the whole implant struc- ture were used to compare the results. The comparison of computations for differentmodels are shown inFig.3. Detailed results of comparable studies are shown inTable 2. Because there are no significant differences in the results for two- and four-cylindrical elements, two elementswere used in the optimization process. This approach enabled us to describe nonlinear asymmetric deforma- tion for axisymmetric geometry and simultaneously significantly reduced size of the problem (ca. 94000 dof) in comparisonwith the fully three-dimensional model (ca. 600000 dof). Fig. 3. TheWallclock time for 3D analyses with the use of cylindrical (CCL-x) and 3D solid elements (C3D-x) 580 T. Łodygowski et al. Table 2. Selected parameters of comparable 3D simulations with the use of cylindrical (CCL-x) and 3D solid elements (C3D-x) Float.point Minimum Required Number Number Number Wall- operations memory disks- of equa- of incre- of itera- clock per itera- required diskspace tions ments tions time tion [–] [MB] [MB] [–] [–] [–] [s] CCL-1 6.51 ·109 47.31 150.73 56226 7 40 160 CCL-2 2.81 ·1010 87.25 366.96 93588 7 45 401 CCL-4 1.61 ·1011 176.93 1034.24 168312 9 62 1565 C3D-16 9.48 ·1011 447.38 2969.60 317760 13 104 8730 C3D-32 5.15 ·1012 1157.12 8427.52 616656 13 109 47685 3. Dental implant optimization The optimization of the presented dental implant is a very complex problem. The model, which is the basis (starting point) for the optimization, is well developed. It estimates the implant behavior giving consideration tomaterial nonlinearity, complex contact definition and prestress of the assembly screw. The level of complexity and large number of design parameters result inmany local minima in the design space. Thus, the optimization algorithm has to search for the optimal solution globally, based on a limited set of solutions in the design space. Moreover, a few problems with FE analysis can occur and should be taken into account. Themost frequent ones are rebuilding and no-convergence problems. The complexity of geometry and large number of design parameter combinations makes it impossible to predict all problems with the geometry. As a consequence, for some configurations of the design parameters which are in the feasible region, the proposed shape is incorrect. Additionally, the problemwith mesh generation can occur, which can be also treated as a rebuildproblem.Reasons for no-convergence areusually causedby too large plastic strain or contact problems. Both can occur for points in the feasible region of thedesign space andhave tobe expectedduringoptimization process. One of the optimization algorithmswhich can fulfill all these requirements is a genetic algorithm. In the presented optimization, a classic binary form with two genetic operators: crossover andmutation is used. Optimization of dental implant... 581 3.1. Genetic algorithm (GA) The genetic algorithm has been inspired by evolutionary biology and in- corporates techniques such as inheritance,mutation, selection and crossover to findabetter solution.Themost important advantages of the genetic algorithm over classic methods of optimization is that it works basing on the problem solution instead of analytical relations. It can optimize linguistic variables and use parallel computing by nature. Despite being a powerful optimization tool, the genetic algorithm uses simple rules, whichmakes it easy to implement. In the module presented in this work, the galileo (GPL) library was used. Fig. 4. Scheme of genetic algorithm processing A genetic algorithm operates on individuals which are abstract represen- tations of a real solution. Each individual contains a set of encoded design variables, which is called a chromosome. The binary encoding is a very classic one and is also used here. The range and number of bits was set for each de- sign parameter. The string representing a chromosome was constructed over the alphabet {0,1} and can be symbolically represented as follows A= a1a2a3 . . .ai . . .al (3.1) In equation (3.1), ai represents a particular bit, an element from alphabet, at the position i, and l denotes the total string length. In order to collect chromosomeswith similar features, we can introduce a schema. The schema is constructed over an extended alphabet of three symbols {0,1,∗}, where ∗ is do not care symbol. The schema can be symbolically represented as fallows H =h1h2h3 . . .hi . . .hl (3.2) hi denotes a particular element from the extended alphabet at the position i, and l represents the schema length. The order of schema H, denoted by o(H) 582 T. Łodygowski et al. represents the number of fixed positions, in this case, 0’s and 1’s, whereas length, denoted by δ(H) is the distance between the first and the last fixed position. The chromosome matches the schema if for i ∈ 〈1, l〉, hi = ai or hi = ∗. Each schema H can be described by the order and the length. The main idea of GA processing consists in generation of a set of initial solutions and use of evolutionary operators to improve them in successive iterations, called generations. An initial set of individuals, called the initial population, is usually randomly created. Every new population is subjected to evaluation. The evaluation process consists of assigning a fitness value to each individual. The fitness value describes solution quality and is calculated according to an objective function and results of FE analyses. Denoting the fitness value of the j-th individual matching the schema H by fHj and their number after t-th iteration by m(H,t), the average fitness value for the schema H after t-th iteration can be represented as follows f(H,t)= ∑ fHj m(H,t) (3.3) whereas the average fitness value of entire population may be given by f = ∑ fj n (3.4) The population size and j-th representation of individuals is denoted by n and fj, respectively. Fitness values are the basis for the selection process. The selectionmecha- nism watches over the improvement of the next population quality. During selection, statistically only the most fitted individuals are chosen in order to allow them to take part in creation of the next population. The probability of choosing the particular individual in a single trial using the wheel-roulette selection method can be given by the following equation pi = fj ∑ fj (3.5) Assuming that the number of a selected individual has to be equal to the po- pulation size, the expected number of individuals whichmatch the schema H in the t+1 population is as follows m(H,t+1)=m(H,t) f(H,t) f (3.6) The number of schema H representatives in the next population grows with the ratio of the average value of the schema to the average fitness of the entire Optimization of dental implant... 583 population. Thus, the schema with a higher average fitness is more strongly promoted. Moreover, one should expect that the number of individuals mat- ching the schema will grow if their average value is higher than the whole population average value. In the next stage, the selected individuals are subjected to crossover and mutation. The operators use selected chromosomes in order to create new design parameter configurations. The primary function of crossover is to mix the parent individual chromosomes and create a new couple of offsprings. The crossover proceeds with a probability defined by the crossover rate pc. The moderation of crossover prevents premature convergence getting stuck in a localminimum. In a classic form, crossover consists of sampling a chromosome partition point, dividing parent chromosomes and swapping obtained parts between them. The crossover improves the population diversity but on the other hand reduces individuals, whichmatch the schema H regardless of their fitness. The probability of survival of the scheme H, in spite of crossover, can be defined as follows ps =1−pc δ(H) l−1 (3.7) In the equation above, l denotes the chromosome length. Mutation is responsible for random distortion of chromosomes. Random distortions keep diversity of the population on a sufficient level and help to escape froma localminimum. Inbinary encoded chromosomes,mutation sam- ples a gene and changes its value to the opposite. The intensity of mutation is controled with the mutation rate pm. Similar to crossover, mutation can destroy individualsmatching the schema H. In the case ofmutation, the pro- bability of H schema representative survival can be given by the equation ps =1−o(H)pm (3.8) Summing up the effect of three operators: selection, mutation and cros- sover, the expected number of representatives of the schema H in the next population can be given as follows m(H,t+1)­m(H,t) f(H,t) f [ 1−pc δ(H) l−1 −o(H)pm ] (3.9) In the last stage, the new population replaces the previous one. Usually, the genetic algorithm ends when either the maximum fitness value for the best fitted individual is obtained or the maximum number of generations is achieved. A general chart-flow of the genetic algorithm is presented in Fig.4. 584 T. Łodygowski et al. 3.2. Definition of the optimization problem 3.2.1. The goal The present study is expected to find a new dental implant shape. The shape should lead to lower principal stresses in comparison with the curren- tly encountered values. In the presented FE model, geometry and material properties are similar to a real dental implant. The loads which were applied come from (Zagalak, 2003) and are recognized as test loads. For the presented configuration, the maximum principal stresses are greater than 625MPa and are localized in the screw corner (Fig.5). Fig. 5. Huber-Mises stress in the initial design of a dental implant In spite of looking for a newproposal of shape, an other goal is verification of the used numerical method. The study should provide an answer to the question wether this methodology can be used in further improvement of the dental implant design. The genetic algorithm tries to find the most fitted individual with the hi- ghest fitness value. In the case ofminimization problem, the objective function has to be modified in order to evaluate individuals. In the present work, the fitness function was defined according to the equation f(X)= ( 1 σmax maxPrinc )2 c (3.10) where X = [A,B,C,D,E,F] and c=109. The value of σmaxmaxPrinc represents the maximum principal stress in the upper part of the implant and the vector X denotes configuration of design Optimization of dental implant... 585 parameters. An additional multiplier cwas used in order to prevent obtaining values smaller than machine precision. 3.2.2. Design parameters Geometry of the presented two-component implantology system is quite complex and the choice of correct design parameters is crucial. The thread was excluded from consideration at this stage of the study. Six geometrical parameters were defined as real variables (Fig.6). All of them refer to the upper part of the implant. Fig. 6. Design parameters The design parameters were encoded into a binary string. There were de- termined for each design parameter range and number of bits for encoding. The ranges come from both the geometry limitations and manufacture requ- irements, and are listed in Table 3. The larger number of bits for encoding, the better accuracy. On the other hand, however, the total time of optimization substantially increases. The decision on the number of bits is always a compromise between time and required accuracy. In the present case, for angle variables four-bit stringswere constructed, whereas for all the rest only three bits were used. The choice of particular design parameters was based on the range and expected influence on final results. The obtained resolution of the design parameter d, can be calculated as follows dπ = dmax−dmin 2n−1 (3.11) 586 T. Łodygowski et al. Table 3.Range and number of bits for design parameter encoding Name Symbol Min Max Range Bits Reso- value value lution Screw head diameter A 0.95 1.6 0.65 3 0.093 Screw head conic B 0.0 80.0 80.0 4 5.33 surface opening angle Screw head height C 1.5 6.4 4.9 3 0.7 Screw diameter D 1.0 1.5 0.5 3 0.071 Hexagonal slot height E 0.05 2.65 2.6 3 0.371 Hexagonal slot conic F 11.0 90.0 79.0 4 5.27 surface opening angle where dmax, dmin denotes the maximum and minimum limit for the variable for string length equal to n. All the binary strings reference to the design variables and are joined in a chromosome. In the present work a twenty bit string is used as below chi = [g A 1 ,g A 2 ,g A 3 |g B 4 ,g B 5 ,g B 6 ,g B 7 |g C 8 ,g C 9 ,g10 C|g11 D,g12 D,g13 D| (3.12) g14 E ,g15 E ,g16 E|g17 F ,g18 F ,g19 F ,g20 F ] The superscripts denote design parameters symbols in accordance with Ta- ble 3. 3.2.3. Constraints The material used in the FE model simulates the elasticplastic behavior. To some extent, it is allowed to exceed the plastic stress limit in the dental implantbut theplastic strain cannotbe too large. InFEanalysis, it is assumed that the maximum equivalent of plastic strain can not be higher than 10%. The equivalent plastic strain can be calculated as ǫ pl = ǫpl ∣ ∣ 0 + t ∫ 0 √ 2 3 ǫ̇pl : ǫ̇pl dt (3.13) where ǫpl ∣ ∣ 0 denotes initial equivalent strain, which was equal zero, whereas ǫ̇pl denotes the equivalent plastic strain rate. The principal stress reduction, whichwas setas thegoal, directs changes in individual in theoppositedirection to this constraint. Thus, no stress or plastic strain constraints were necessary. In the case of rebuild or convergence problem, the death penalty method is applied. The individual is eliminated from further processing. Optimization of dental implant... 587 3.2.4. The genetic algorithm In the evaluation of a single individual, theFEmodel of the dental implant which is timeconsuming significantly limits the population size. On the other hand, the population has to be large enough to provide sufficient diversity and represent the design space properly. In the present work, a forty-individual population was used. The number of epochs was not limited a priori. The optimizationwas stopped, based onmonitored results. The signal to brake the procedure was sent after a few epochs without any new proposal for better solution. Thefirstpopulationwas randomly created.Moreover, because of the death penalty method constraints, individuals in the first population were genera- ted and checked as long as the whole population of correct individuals was not achieved. The procedure prevents starting from individuals, which will be eliminated in the second generation. The ranking-based method of selection was used. This type of selec- tion consists in sequential arrangement of individuals according to their fit- ness and assigning to each i-th individual in the range from rmini to r max i where rmin0 =0 r min i = r max i−1 rmaxi = r min i + pi ∑n j=0pj rmaxn =1.0 (3.14) In the above equations pi and n denote the position in the ranking and the population size, respectively. During the selection, a value between 0 and 1 is randomly generated and indicates the individual which refers to the range in which the value is included. In contrast to the classic method of selection, i.e. roulettewheel, this method sets probability of being selected according to the position in the ranking instead of thefitness.As a consequence, theprobability of choosing an individual with much lower fitness than the maximum one increases and results in higher diversity in the late stage of the optimization process. The onepoint crossover method is applied. The selected individuals are joined into couples, their chromosomes are divided into three substrings at random positions and swapped between parent individuals. A new couple of offsprings represents combinations of both parent features. In order to keep diversity on a sufficient level during the optimization process, not all couples are subjected to crossover. The intensity of crossover is determined with the crossover rate, which is constant and equal to 0.75 in the present work. 588 T. Łodygowski et al. In the last stage of recombination, random distortions are introduced to chromosomes. The probability of each gene distortion is controlled by the mutation rate. In the present work, the constant value of 0.02 during the whole optimization process was kept.When a gene is selected tomutation its value is changed to opposite one, from 0 to 1 and vice versa. The classic form of population replacement was used. The whole previous populationwas replaced with the new one completely. Thismethod allows for increasing the diversity and search through the design spacemore extensively in comparison with alternative methods, e.g. steadystate replacement. Unfor- tunately, at the same time, it moderates the convergence and increases the number of necessary FE analyses. The small population size, which was used, limited the variation of the design parameters, therefore the diversity was an important factor. 3.3. Incorporation of the genetic algorithm into Abaqus/CAE The genetic algorithm processes chromosomes which represent various combinations of design parameters. Each new chromosome has to be evalu- ated in order to assess its quality. A great advantage of the genetic algorithm is that the evaluation process can be fully elaborated in an exterior proce- dure. The procedure returns fitness values based on chromosomes. The tool employed in the evaluation procedure has to provide possibility to rebuild the FEmodel for any design parameter configuration.Moreover, some procedures like meshing has to be created automatically. Thus, the reliability and capa- city of the rebuilding tool are crucial. In the present work, the Abaqus/CAE environment in cooperation with Abaqus/Standard code is used. 3.3.1. Optimization module for Abaqus/CAE Abaqus/CAE provides an efficient environment for user scripting. The object-oriented scripting language, called python is available. The user can add newprocedures into the kernel and graphical user interface (GUI) aswell. In order to carry out the dental implant optimization, a newmodule has been created with a fully functional graphical user interface (Fig.7). The optimization driver was imported as an external library into Aba- qus/CAE environment. The open source library called galileo was used. The implemented module was divided into two main sections. The first one orga- nizes genetic optimization and employes galileo library mostly, whereas the second provides a tool for chromosome evaluation. Optimization of dental implant... 589 Fig. 7. Optimizationmodule for Abaqus/CAE 3.3.2. Individual evaluation Every chromosome created in the reproduction process is a starting point for the evaluation procedure. In the first stage, each i-th chromosome is di- vided into substrings which refer to the design parameters. Next, each binary substring is decoded into an integer value ej and i-th design parameter is calculated with the given equation dj = ej dmaxj −d min j 2n−1 (3.15) In the equation above, the symbol n denotes the substring length. Finally, a set of real design parameter values is obtained and provides the basis for model rebuilding. All design parameters refer to geometry of the dental implant. Aba- qus/CAE kernel script is created and run in order to change them in the model. The obtained geometry is controlled and if the shape is invalid, the rebuild error is raised. The weakest point of the procedure is the creation of a mesh. Themesh is generated automatically and due to the large number of analyses it is not possible to control it manually. It is assumed that any er- ror in the mesh result in higher principal stresses and eliminates the solution during selection. Based on the modified geometry, the FE analysis definition is created and sent to the Abaqus/Standard solver. In the case of the pre- sent work, an additional analysis has to be carried out to provide results for 590 T. Łodygowski et al. the axisymmertic model. The spatial dental FEmodel was built based on the results for axisymmetric one andadditional data such as lateral forcemagnitu- de. Themodule waits for analysis termination until themaximum time is not exceeded. In the case of too long analysis, the solver stops and no convergence error is raised. This limitation can also eliminate well fitted solutions but is necessary to proceed the optimization. The analysis results and all information about its proceeding is stored in an output database, which can be read with the use of the python script. All necessary information is extracted from the output database and the objective function is calculated. Finally, the fitness value returns to the optimization driver.Thewholeflowchart of the individual evaluationprocedure is presented in Fig.8. Fig. 8. Individual evaluation procedure In the case of the rebuild error, the zero fitness value is returned for an individual. The zero fitness value guarantees that invalid solution will be eli- minated during the selection process. 3.3.3. Parallelization Genetic processing generates a large number of solutions that have to be checked. Taking into account the time of a single FEanalysis, the total time of optimization (computations) is a frequent limitation.The single analysis of the presentedFEmodel took 5 to 40minutes and a large spreadwas caused by the nonlinear behavior. Because the time of analysis is significant, the algorithm was modified to support parallel analyses. During the evaluation stage, the optimization driver waits for fitness values of all individuals. The population is divided into subgroups according to available resources. In the work, an Optimization of dental implant... 591 eight-processors machine was used to determine the number of individuals, which were evaluated simultaneously. For each individual in a group, an FE analysis definition is prepared.Next, all analyses are submitted simultaneously and the longest one determines the time which is necessary for evaluation of thewhole group.The obtained results are the basis of fitness calculation. They are returned to the optimization driver. The procedure is presented in Fig.9. Fig. 9. Scheme of genetic algorithm evaluation processing 3.3.4. Procedure of building of an FE implant model during the optimization process Theprocedure of building of anFE implantmodel during the optimization process is schematically shown inFig.10.At thebeginningof theoptimization, theAbaqus/CAEfilewith the fully parametric axisymmetricmodel of the im- plant is opened. Thismodel was described in details before (Section 2).When the optimization loop is started, the axisymmetric model of the implant is modified. TheAbaqus/CAE creates an INPfile with amodified axisymmetric model of the implant and the Optimization Module submits jobs. The results of these axisymmetric simulations are assembled in the implant structure. In the next step, the optimization module generates PARAM files and simul- taneously runs three-dimensional jobs. The INP file of the three-dimensional model contain the definition of symmetric model generation procedure, and the second step of implant simulation ismanually prepared.ThePARAMfiles are created based on information from axisymmetric ODB files and options which are defined at the beginning of the optimization procedure. These files contain parameters of three-dimensional models of the implant: • segment angle through which the cross-section is to be revolved, • number of elements to be used in the segment, 592 T. Łodygowski et al. • offset for element numbering, • offset used for node numbering, • global number of node in the axisymmetric model. Fig. 10. The procedure of building of an FE implant model during the optimization process Finally, the results of three-dimensional simulations are read fromODBfi- les andtheOptimizationModulecan start theGeneticAlgorithmprocess.This proceduremust be repeated for each population. The describedmodeling ap- proachwith the use of axisymmetric geometry description and semi-analytical discretization enable us to carry out a large number of implant simulations in a realistic time period and to carry out the optimization. 3.4. Results The optimization was performed using 8 cpus in computations (Fujitsu- SiemensPRIMERGYTX300 S3). During thewhole optimization process over 1600 FEA analyses were curried out. The total time of optimization was ca. 250 hours. The evolution of the optimization process is presented in Fig.11. It shows changes of the maximum principal stresses in the upper part of the implant Optimization of dental implant... 593 screw for successively generated solutions. The average value of the principal stress is drawn by the solid curve. It can be observed that the designs are being improved - themaximumprincipal stress is decreasing. Startingwith an value of 625MPa, the principal stress is reduced down to 150MPa. The graph does not consist of the first stage when the initial population was established. More than 320 FE analyses were necessary to find forty correct individuals. Startingwith the first population of random individuals yields better solution. Moreover, the reduction of principal stress is meaningful and equals 475MPa, what ismore than 75% of the initial value. In the next generations, the design parameters from the best solution (peaks) were promoted stronger and thus they strongly influenced the population improvement. It resulted in constant decreasing of both average andminimumvalue of principal stresses. After the 20th generation, the best solution was established and no further essential reduction of stress was observed. Because of the chosen reproduction strategy, the diversity within population was kept on a sufficient level and even in the last monitored population the individuals represented a wide range of fitness value. Fig. 11. Fitness value for sucessively generated individuals Two proposals of new shapes are shown in Fig.12. The lowest principal stress was obtained for shape (b) and equaled 150MPa. The next proposal (Fig.12c) provides reduction down to 180MPa but still has a hexagonal slot instead of the most optimal solution. The haxagonal slot plays an important role in preventing screw loosening, thus solution (b) will not be considered in further analyses. There are a few clear tendencies easy to observe. In the first place, the screw head is wider than the initial one and the opening angle 594 T. Łodygowski et al. of the screw head conic surface is moderated. Together with modification of the opening angle of the hexagonal slot, the conic surface makes clamping of the abutment and the body more efficient. Moreover, moderation of the opening angle of the screw head conic surface reduces stress concentration in the interior corner of the screw. Fig. 12. Initial (a) and new (b, c) shapes of dental implant All the modifications make the joint stiffer and reduce relative rotation between the abutment and the body. The implant starts to behave as a canti- lever (Fig.13). The final implant incorporates a body to carry the lateral load. As a result, the displacement of points in the upper edge is reduced almost twice from 0.085mm down to 0.045mm. The change of dental implant behavior results also in amore uniformstress distribution (Fig.14). The contact pressure is realized on the whole contact surface in contrast to the initial implant. Consistently, the screw is less bent. All the proposals of new shapeswere verifiedwith the use of the full three- dimensional FE model. The results confirmed that the FE model built with cylindrical and solid elements provides compatible values and can be used for further optimization. 4. Conclusions The present approach successfully incorporates an FE solver into a genetic optimization procedure. The complex dental implant model was optimized based on FE analyses. The study and the final results give evidence that the Optimization of dental implant... 595 Fig. 13. Deformation of the initial (a) and final (b) dental implant (displacement scaled up 40 times) Fig. 14.Mises stress in the initial (a) and final (b) dental implant presented method is efficient and can be used for further analyses. Moreover, the FE model built with cylindrical and solid elements was validated and no meaningful differences with the full three-dimensional solid FE model were observed. A new shape was foundwhich provided a solution with substantial reduc- tion in principal stresses and displacemens. The joint between the abutment and the bodywasmuch stiffer. It made the implant behavemore as a cantile- ver than a two-element system.As a consequence, the obtained dental implant 596 T. Łodygowski et al. transmitted stress from the abutment to the bodymore smoothly and thema- ximumdisplacementwas reduced twice. In comparisonwith theoptimal dental implant it is visible that the initial one behaved more as a hinge where the screw is bent. The obtained reduction of the principal stress allowed also for assumming a significant extension of long-life strength of the implant. Fatigue damage of implant components is not a common phenomenon but it makes further treatment very difficult. One of themain final goals of dental implant optimization is the minimization of fatigue damage of the implant. The obtained level of principal stress in the key parts of the implant connec- tion enable us to reduce or even eliminate the risk of fatigue damage. It should be noted that despite the level of principal stress reduction only a part of it can be obtained in the real product. In the first place, the limitation comes frommanufacturing requirements. The dimensions have to be adopted to the given manufacture facility. The second limitation refers to the screw loose- ning. The shape of the hexagonal slot which connects the root of the implant with the abutment plays the key role in kinematic behavior of the implant. The final geometry of screw connection must protect against screw loosening phenomenon as well. The final shape of this part must also ensure easiness of in-vivo assembling of the implant. On the one hand, the hexagonal slot must safely transfer torsional load during implant fixing and abutment assembling, on the other hand the dentist should easily place the abutment in the slot. As shown above, the real implant must consider many mentioned limita- tions. The obtained results of the optimization procedure cannot be directly adopted in implant systems used at present, but can be a start point for new designs of such systems. However, the presented approach of optimization can be used in further studies and designing of modified and new dental implants. Theprimarydisadvantage of thepresentednumerical approach is thenum- ber of FE analyses, which are necessary to be carry out in the procedure. Despite the total time of the analysis can be shortened using parallel com- putations, it is still a considerable limitation. Thus, it is an important goal of the study to modify the procedure in order to reduce the number of time- consuming analyses. The first study with the use of a neural network as a surrogate model has been already elaborated, and the obtained results seem to be promising. Acknowledgment The support of Poznan University of Technology grant DS 11-034/09 as well as the possibility of using the computational environment of Poznan Supercomputing and Networking Centre are kindly acknowledged. Optimization of dental implant... 597 References 1. Abaqus Analysis User’s Manual, 2007a, SIMULIA, Pawtucket 2. Abaqus/CAEUser’s Manual, 2007b, Pawtucket 3. Bozkaya D., Müfüt S., 2005, Mechanics of the taper integrated screwed-in (TIS) abutments used in dental implants, Journal of Biomechanics, 38, 8797 4. BurczyńskiT.,DługoszA.,KuśW., 2006,Parallelevolutionaryalgorithms in shape optimization of heat radiators, Journal of Theoretical and Applied Mechanics, 44, 2, 351-366 5. Burczyński T., Kuś W., Długosz A., Orantek P., 2004, Optimization anddefect identificationusingdistributed evolutionaryalgorithms,Engineering Applications of Artificial Inteligence, 17, 4, 337-344 6. Chao H.K., Rowlands R.E., 2007, Reducing tensile stress concentration in performed hybrid laminate by genetic algorithm, Composites Science and Technology, 67, 13, 2877-2883 7. GoldbergD.E., 1989,GeneticAlgorithm in Search, Optimization andMachi- ne Learning, 1st Ed., Addison-Wesley Professional 8. GoodacreC.J., Bernal G., RungcharassaengK., Kan J.Y., 2003, Cli- nical complicationwith implants and implant prostheses, International Journal of Prosthodontics, 90, 121-129 9. HędzelekW., ZagalakR., Łodygowski T,WierszyckiM., 2004,Bada- nia biomechaniczne elementów protetycznych implantów z zastosowaniemme- tody elementów skończonych,Protetyka Stomatologiczna, 51, 1, 23-29 10. Kąkol W., Łodygowski T., Wierszycki M., 2002, Numerical analysis of dental implant fatigue,Acta of Bioengineering and Biomechanics, 4, 1, 795-796 11. Lang L.A., Kang B., Wang R.F., Lang B.R., 2003, Finite element analy- sis to determine implant preload, The Journal of Prosthetic Dentistry, 90, 6, 539-546 12. Merz B.R., Hunenbart S., Belser U.C., 2000,Mechanics of the implant- abutment connection: an 8-degree taper compared to a butt joint connection, The International Journal of Oral and Maxillofacial Implants, 15, 4, 519-526 13. Michalewicz Z., Schoenauer M., 1996, Evolutionary algorithms for con- strainedparameter optimizationproblem,EvolutionaryComputation,4, 1, 1-32 14. Natali A.N., 2003,Dental Biomechanics, Taylor & Francis 15. Szajek K., Kąkol W., Łodygowski T., Wierszycki M., 2008, Optimi- zation module for Abaqus/CAE based on Genetic Algorithm, Abaqus Users’ Conference, Newport, USA, 447-462 598 T. Łodygowski et al. 16. Wang K., 1996, The use of titanium for medical applications in the USA, Materials Science and Engineering,A213, 134-137 17. WierszyckiM., 2007,Numeryczna analiza wytrzymałościowa wszczepów uzę- bienia oraz segmentu kręgosłupa ludzkiego, PhDThesis, Poznań 18. Wierszycki M., Kąkol W., Łodygowski T., 2006a, Fatigue algorithm for dental implant, Foundations of Civil and Environmental Engineering, 7, 363-380 19. Wierszycki M., Kąkol W., Łodygowski T., 2006b, Numerical comple- xity of selected biomechanical problems, Journal of Theoretical and Applied Mechanics, 44, 4, 797-818 20. Wierszycki M., Kąkol W., Łodygowski T., 2006c, The screw loosening and fatigue analyses of three dimensional dental implant model, ABAQUS Users’ Conference 2006, BostonMA 21. Zagalak R., 2003, Evaluation of Mechanical Properties of Two Dental Im- plants Osteoplant, PhD Thesis, University of Medical Sciences in Poznań [in Polish] 22. Zagalak R., Hędzelek W., Łodygowski T., Wierszycki M., 2005, Wpływ zaniku kości i gęstości na ryzyko złamań implantów - badania metodą elementów skończonych, Implantoprotetyka, 6, 1, 3-7 Optymalizacja wszczepu stomatologicznego z wykorzystaniem algorytmu genetycznego Streszczenie Przedmiotem prezentowanej pracy jest problem optymalizacji systemu implanto- logicznegoOsteoplant, który został opracowany i wciąż jest ulepszany przez Fundację Uniwersytetu Medycznego w Poznaniu. Obserwacje kliniczne potwierdzają występo- wanie powikłań zarówno we wczesnej, jak i późnej fazie użytkowania implantu. Do- tychczas otrzymane wyniki wskazują, że wydłużenie bezawaryjnego okresu użytko- wania implantu wymaga wprowadzenia zmian w jego pracymechanicznej. Jednakże, ustalenie szczegółówmodyfikacji nie jest oczywiste.Wartykule zostałaopisanaproce- dura optymalizacji systemu implantologicznegoOsteoplant z użyciem analizymetodą elementów skończonych oraz algorytmu genetycznego. Manuscript received February 11, 2009; accepted for print April 8, 2009