Support Tool for Anchoring System Optimization of Titanium Craniofacial Prostheses

Titanium prostheses are artificial components used during cranioplasty, which can be made with the Incremental Sheet Forming (ISF). To repair the craniofacial defects, the implant must be fixed to the skull, so an adequate anchoring system is necessary. Because there are no specific guidelines regarding this system, the aim of this study was to present a methodology for the identification of the optimal screw shank diameter and the overlap length, considering different damage areas and accidental external loads. In this method, a numerical model for structural analysis, a statistical approach and a metamodeling technique were coupled. The results demonstrate that this methodology can be used as a decision support tool in the titanium prostheses design, ensuring the best


Introduction
Cranioplasty is a surgical procedure aimed at repairing craniofacial defects, which are changes in the anatomy and morphology of skull and facial bones due to different causes, such as traumas, congenital malformations or diseases that require neurosurgical interventions [1].
Initially, during this surgery, bone grafts were employed and the first one was implanted in the 1668, even if the use of autografts for cranioplasty became popular in the 20th century [2]. This type of intervention requires a second surgery and, so, a good amount of substitutive bone [3], which can be both considered as disadvantages.
So, in order to replace the human bone, artificial prostheses may be used to correct the defect and to fill the hole. These implants are made with substitutive materials, which must be inert, thermally non-conductive, non-corrosive, stable, resistant to infection, with good biomechanical properties malleable and easily contoured, inexpensive, readily available, and can be sterilized [4]. The most adopted material is the titanium and its alloys, in particular the bio-alloy Ti6Al4V, characterized by very good properties [5]. Furthermore, this material favors the osseointegration with the bone that indicates the direct contact at microscopic level between living bone and implant, creating a permanent anchorage that improves the long-term function of the prosthesis [6].
Furthermore, the prosthesis is anchored to the craniofacial bone using screws, which can be of different types, such as self-tapping or self-drilling (for example: Stryker Global Headquarters, Kalamazoo, Michigan, USA and Synthes, Inc., West Chester, Pennsylvania, USA).
Because the face is the most important aesthetic part of the human body [4], craniofacial prostheses have both functional and aesthetic aspects, which must not be compromised. For this reason, the anchoring system can be considered as one of the major issues for the surgical treatment success. Indeed, this system allows the fixation of prosthesis to the bone and, consequently, the correct positioning between them, and ensures the protection of intracranial structures, as the brain, from mechanical impacts and a normal and appreciable appearance [1].
Thanks to the computational advantages in medical image acquisition and to new design methodologies, recently custom-made prostheses can be made starting from the reconstruction of patient-specific 3D models obtained processing the computer tomography (CT) and the magnetic resonance imaging (MRI) scans [7].
Moreover, titanium custom-made craniofacial prostheses can be made using the Incremental Sheet Forming (ISF), a manufacturing process with significant vantages [8][9][10][11]. Generally, with this manufacture, the prosthesis has a greater area and, so, a larger perimeter in respect to the defect [11], but no information regarding the anchoring system is available, in terms of overlap length and diameter of screw shank.
Since cranioplasty is a research branch where surgical expertise has to meet engineering knowledge for a good success of craniofacial reconstruction, the aim of this study is to present a methodology for the optimal design of the anchoring system in case of patient-specific titanium prostheses made using ISF. In details, different defect areas have been analyzed in a numerical environment in order to identify the better overlap length and the diameter of screw shank. So, a "ready to use" decision support (DS) tool was developed in order to put into closer contact the players involved in the clinical process for obtaining custom and faster implants. The idea is to make easier the integration between bioengineers and surgeons during the preliminary phase by adopting the DS tool, which merges the expert knowledge and helps to converge toward the optimal solution.

Methodology
To pursue the aim previously reported, a hybrid approach has been developed, based on the coupled use of numerical model for structural analysis, statistical approach and metamodeling technique. Actually, this choice results necessary for enriching the robustness and wideness of methodology.
Hybrid approaches are usually implemented in engineering problems [12]. Focusing to the interested application field, a dataset derived by the implementation of a robust finite element model is typically analyzed by ANOVA or Response Surfaces [13] as far as used for training Neural Network or optimization tools [14]. Their use is aimed at describing and predicting the behavior of such response variable.
A similar methodology is proposed in this study, by following the sequential steps described in Figure 1. First of all, a finite element model (FEM) of the anchoring system has been implemented in order to analyze the stress and displacements distributions and to generate a base of knowledge. The outcomes of the FEM analysis represent the incomes for the Response Statistical Methodology (RSM), used both for understating of the factors significance on the investigated output and for defining a prediction function that analytically describes the input-output relationship.
After that, a multi-objective model has been set and developed for selecting the best value of different features of the cranial prosthesis for each damage area value.
All details are provided in the next sections.

Design of experiment
The first step is to define the boundaries of the analysis fixing a lower-upper value of each investigated feature. So, the range of each input was chosen as reported in Table 1. Subsequently, it was necessary to define the sampling strategy with which a set of data was collected in order to test and train the model. For the specific case, a full factorial Design of Experiment (DOE) was created, considering three levels for all the variables, except for the damage area, in which only two possible values have been tested (Table 1). The number of experiments necessary to carry out this design is N = L k , where L represents the number of levels for the investigation and k represents the number of variables, or factors (five in our case). The number of simulations for the specific case study is therefore equal to 162 (3 4 •2 1 ).

Numerical model and structural analysis
In order to design the optimal anchoring system for a specific defect, 162 different numerical simulations were carried out using the finite element method (FEM) and five variation factors (Table 1). All structural simulations were performed by means of COMSOL 5.0 (COMSOL Inc, Stockholm, Sweden), a finite-element-based commercial software package.

Geometrical modeling
An assembly, including the skull bone, a prosthesis and a screw, was used as geometry in each structural analysis.
An easy way to obtain a volumetric model of head and of defects is to use the segmentation process of medical images (or slices), available in DICOM format, generated by CT or MRI. Using the segmentation process, a 3D virtual model of the patient skull can be obtained by means of different software [15] chosen for the 3D/volumetric reconstruction of a skull bone ( Figure 2). After that, the reverse engineering process was applied in order to obtain a 3D solid model. Since the reconstructed head has no defects, a cylindrical hole was created in the upper part (frontalparietal bone), considering two damage areas, of 2100 and of 2600 mm 2 , respectively. Moreover, the "removed bone" was used to model the prosthesis ( Figure 3A), that was obtained considering a thickness of 1.5 mm [8] and different overlap lengths (7, 10 and 13 mm). Furthermore, only a little part of the skull bone was considered ( Figure 3B). Because the defect and the prosthesis were axisymmetric, a simplified model of 20 degrees was considered ( Figure 3C). These two assumptions were necessary in order to reduce the computational cost obtaining realistic solutions anyway.
Regarding the screw, it was modeled considering a length of 4.5 mm, three different screw shanks (1.5, 2.0 and 2.5 mm) and ignoring the thread. Moreover, it was placed in the center of overlap length ( Figure 3C). The final considered model of the assembly is reported in Figure 3C.

Mechanical properties
Skull bone was approximated as a cortical one [17], with an ultimate compression stress of 145 MPa [17], whereas prosthesis and screw was set as Ti6Al4V [18]. Both materials were assumed to be homogeneous and isotropic [17], defined by linear elastic laws.

Boundary conditions
Because the exact fit between bone and implant is very important and is correlated with the osseointegration and, consequently, with the long-term function of the prosthesis [6], the screw-bone interface was modeled considering a perfect osseointegration. To numerically reproduce this condition, the identity pair was applied between bone and implant and bone and screw. Moreover, as the osseointegration was numerically reproduced, a rigid connection was assumed between prosthesis and screw.
Titanium and bone have different elastic modules (reported in [18] and [17] respectively), so, the anchoring system must be correctly designed in order to avoid the skull bone rupture in case of mechanical impacts. For this reason, a distributed compressive load was applied to the top surface of prosthesis ( Figure 4) in order to simulate an accidental load of 100, 300 and 500 N, considering also three tilt angles (0°, 45° and 90°).
Moreover, skull bone was fixed in the bottom part in order to simulate the connection with the spine.
Finally, as reported before, a simplified assembly of 20° was modeled ( Figure 3C), so symmetry boundary conditions were applied in the lateral faces of skull and prosthesis.

Simulations details
In all cases, the grid was created using fine tetrahedral mesh.
Since the anchoring system must be design to resist to accidental loads, the maximum stress on skull bone due to loaded prosthesis was evaluated by means of the Von Mises criteria [19]. An example of stress distribution is reported in Figure 5.

Response Surface Methodology
The numerical data provided by FEM analyses were used to build a model using the response surface methodology (RSM) [20]. It consists of a group of mathematical and statistical techniques adopted in the development of an adequate functional relationship between a response of interest, y, (in our case the stress of skull) and a number of associated control (or input) variables denoted by x 1 , x 2 , ... , x k .
In general, such a relationship is unknown but can be approximated by a low-degree polynomial model of the form: is a vector function that consists of powers and cross-products of powers of x 1 , x 2 , ... , x k up to a certain degree, β is a vector of unknown constant coefficients, and ε is a random experimental error assumed to have zero mean. In this case, the quantity f (x)β represents the expected value of y.
The purposes of building a model such as (1) are the following: • to establish a relationship, albeit approximate, between y and x that can be used to predict response values for given settings of the control variables. • to determine the relevance of the factors analyzed. • to determine the optimum settings of x.

Optimization model
The input-output relationship defined through the RSM can be used to select the best features of the craniofacial prosthesis to optimize one or a set of criteria. In details, the prosthesis will not reach high stress if subjected to an accidental load, considering the maximum allowed stress equal to the ultimate compression stress of bone [17]. Furthermore, screws will select from an available set. So, the considered optimization problem has two objective functions. To find the efficient frontier, a multi-objective algorithm is necessary. Moreover, due to the non-linear input-output relationship, the problem is a multi-objective problem with non-linear constraints.
One of the most popular algorithms to solve this kind of problems is the NSGAII. The first version of this algorithm was proposed by [21]. It is a very effective algorithm, but has been generally criticized for its computational complexity and lack of elitism. A new version of this algorithm was implemented later with a minor complexity of O(MN 2 ), where M is the number of objective functions and N the size of the population [22]. This new algorithm incorporates also elitism ensuring the survival of the best solutions. Furthermore, the algorithm can be briefly described as follow. First of all, a random population is generated, and then the population is sorted according to non-domination in several fronts. The first front is totally non-dominated in the current population, while the second front is dominated only by the first front. Each individual in each front has a value of rank (fitness) based on the front to which it belongs. In addition to the fitness value is considered a function called crowding distance [23]. The crowding distance is a measure of the distance of an individual from its two closest neighbors. Large average crowding distance creates a population with a greater diversification. In the version of the NSGAII adopted in this study, the simulated binary crossovers were chosen for the crossover phase and the polynomial mutation for the mutation phase. One of the most important aspects of this technique is given by the choice of two parameters: the size of the population and the number of the generations. A too high value of these two parameters leads to an increase of the computational time. However, decreasing the number of generations the quality of the solutions found deteriorates while decreasing the size of the population the number of optimal solutions found is lower.

Discussion of the results
The mathematical approach described in section 5 was implemented to build a guided user interface, which can be considered as a useful decision making support tool.
The data collected through a full factorial DOE, obtained by means of a set of FEM simulations of structural analysis, were used to build a polynomial model.
The design and statistical analysis of experiments were done through the Design Expert (State Ease, Inc, trial version) software.
For the investigated factors, a reduced third order polynomial (regression) model was chosen to describe the relationship for the response and tested factors. A p value threshold of 0.05 was adopted to select the statistically relevant parameters interactions to build the model. The selected model is characterized by a high level of accuracy ( Table 2). The overall adequacy of the model in describing the numerical data may be evaluated by the coefficient of multiple determination, R 2 . It is a measure of the proportion of the overall variability in the response-bone stress in this study -, that is accounted for by the regressor variables and interactions that have been selected. However, R 2 will always increase when a variable is added to the model, whether or not the variable is significant, and thus models with large values of R 2 may yield poor predictions of future observations. A more meaningful measure is the adj R 2 (the adjusted value of R 2 ), which will only increase when a new variable is added if there is enough reduction in the residual sum of squares (SSE) to compensate for the loss of a degree of freedom. In our case, the "Predicted R-Squared" of 99.8% is in reasonable agreement with the "Adjusted R-Squared" of 99.6% ( Table 2). As a result, this model can be chosen to navigate the design space. Analyzing the measured response, the fit summary of the output indicates that the reduced cubic model is statistically significant. Furthermore, the normal plot of the residuals was employed to verify the normal distribution of data, that was confirmed ( Figure 7A). Also Gaussian probability density function was applied to the residuals to check the normality and detect outliers, which were confirmed since all data were in the interval of [-3σ,3σ] ( Figure 7B). Finally, to investigate the reliability of the prediction, the predicted versus actual diagnostic plot was analyzed ( Figure 7C). All data were properly predicted, also around the ultimate compressive stress value (145 MPa).
The response surface of bone stress is reported in Figure 8, considering a damage area of 2500 mm 2 and an accidental load of 350 N and 0°. The response decreases if the overlap length increases, regardless of the diameter of screw shank.
To select a suitable value of screw diameter and overlap length, the NSGAII algorithm was applied to the cubic model. The mathematical model is the following: (2) I=I 0 , L=L 0, D=D 0 S S where I 0 , L 0 and D 0 are the input values of tilt angle, load and damage area and S is the set of available types of screws. In the present case study, the preliminary hypothesis is that three types of screws are available, with a diameter equal to 1.5, 2 and 2.5 mm, chosen since they are the most used in neurosurgery.
The problem was to select the best screw shank diameter and the overlap length for a given damage area, ensuring that the stress due to an external load is less than the ultimate compression stress of bone.
In a first scenario, a tilt angle equal to 0° and a damage area of 2600 mm 2 were considered. The NSGAII algorithm was applied to define the efficient frontier for all the possible screw types. The results are reported in Figure 9. It is not possible to find an overlap length which allows to obtain a stress less than 145 MPa for all the analyzed screw diameters. Moreover, the three curves are really close, showing that the diameter value is not a critical variable (at least in the analyzed range of 1.5-2.5 mm).  In a second test, an inclination of 40° and a damage of 2400 mm 2 were considered. The efficient frontier of this scenario is reported in Figure 10. Also in this case the three curves are almost identical, confirming that the stress is mainly due to the level of overlap length between prosthesis and cranium. However, in this case, it is possible to respect the constraint on the maximum stress by using an overlap value greater than 10 mm.

Conclusions
During cranioplasty, patient-specific titanium prostheses, made using ISF, can be adopted to reconstruct the skull or facial defects. The presented methodology, that coupled a numerical model for structural analysis, a statistical approach and a metamodeling technique, allows to identify the best anchoring system, in terms of overlap length and diameter of screw shank, so it can be considered as a useful decision support tool for the design of titanium prostheses.