Computer Simulation-Based Multi-objective Optimisation of Additively Manufactured Cranial Implants

1. Introduction and Brief Survey of the State of the Art

Cranioplasty is the surgical repair of damaged cranial bone resulting from congenital defects such as complex craniosynostosis, tumours, traumas caused by previous surgeries, infections, or fractures resulting from vehicular, sports, or everyday accidents [1]. Its primary goal is to restore the protective function of the neurocranium, and secondarily, to improve craniofacial aesthetics. A cranial implant or prosthesis is the biomedical device used mend the protective function of the affected portions of the skull [2].

The ideal material for the creation of cranial implants for cranioplasty should be: radiolucent, resistant to infections and biomechanical processes, malleable to adapt to defects with a complete closure, non-conductive of heat or cold, inexpensive, and easy to apply [3]. Realistically, there is currently no single material that meets all these requirements and is always applicable, it must be chosen taking into account the patient’s age, the complexity of the defect, the specialist’s judgement and the main reason for performing the cranioplasty[4].

Currently, the most commonly used materials for the construction of cranial implants are metals, polymers and ceramics. Particularly popular are the titanium alloy Ti6Al4V, and the polymers Polyether-ether-ketone (PEEK) and Polymethylmethacrylate (PMMA). Ti6Al4V is widely used to make customised cranial prosthesis [5]. However Kim et al. [6] reported that alloy-based cranial implants can provoke an increase in thermal conduction inside the cranial defect zone, because of their higher density, making them less suitable in cases with direct heat exposure. PEEK, a semi-crystalline polymer that is radiotransparent and chemically inert, offers greater strength, rigidity, and durability compared to other alloplastic materials. It has been used more recently in cranioplasty showing several advantages due to its high strength and toughness and its compatibility with the skull[7]. Current trends are directed to optimize the process of fabrication using PEEK, as is the case with the study presented by Smith et al. [8], which presents a series of configurations for improving the robustness and aesthetics, over other conventional approaches.

It has been also reported that compared to titanium alloys, PEEK implants presents lower rates of complications and failures [9]. One of the main disadvantages of using PEEK lies in its high cost of acquisition and manufacturing [10,11]. This has led to the use of PMMA as a cheaper alternative, given that the material cost is significantly less expensive and the process of manufacture much easier [12]. PMMA has less durability, strength and biocompatibility when compared to PEEK [13]. On the other hand, Moncayo-Matute et al. [14] compared the mechanical response of cranial implants fabricated with both PEEK and PMMA, and showing that there were no noticable discrepancies between personalised implants made from PMMA and PEEK, concerning the mechanical response under a load of 8 kN, while the Von Misses stress distribution performs significantly different for both materials, with PEEK implants presenting better results. Furthermore, PMMA implants exhibits higher infection and complications rate compared to PEEK ones [9]. As currently stands, PEEK is closer to being the gold standard material for the construction of personalised cranial implants.

When designing a customised cranial implant adjusted to the specifics requirements of the patients, computer aided design (CAD) techniques and methodologies give optimal results [15]. Common approaches consist in using commercial software such as Mimics, Magics and SolidWorks, with highly specialised features for the task, or employing classic design procedures such as the mirroring technique. While yielding particularly good results, these methods are usually time-consuming and require highly skilled operators, which represents a real disadvantage still to overcome [16]. Directed efforts towards automation have been made. In [17], the authors propose an automated design methodology for the construction of a customised cranial implant using image processing and open-source python libraries to obtain the implant in stereolithographic (STL) format from DICOM image data. Pimentel and co-authors introduce an automated approach for reconstructing damaged regions in the skull using a 3D statistical shape model and a 2D generative adversarial network trained unsupervised on healthy data [18]. With the rise and popularity of artificial intelligence, further advances are expected in this field [19,20,21].

Studies about optimisation and design considerations have also flourished, with focus on enhancing the mechanical features of the prosthesis and the aesthetics of the repaired skull [8,22,23]. Research by [24] explores extensively the optimal configurations for fixation mechanisms and the adequate separation between these points, suggesting four to five fixations, independently to size of the defect, with symmetrical orientation of the fixation screws where feasible. Another important line of development is to redesign implants based on natural analogies. Works such as [25] presented two implant variants for the studied case, based on dissimilar scaffold architecture that showed mechanical response adequate for implantation. Sharma et al. [26] They propose a workflow to design a cranial prosthesis that moulds to the specific characteristics of the patient, with a macrostructure of interconnected struts that mimics trabecular bone based on the Voronoi diagram, resulting in a lightweight titanium prosthesis with a good fit to the skull. Sivakumar et al. [27] suggest a biomimetic cranial prosthesis made of PMMA and optimised by incorporating a square-type porous lattice structure, considered a suitable for the development of cranial implants reducing the weight without compromising quality. Lattice structures are a type of architectural framework composed of repeating units or cells, typically arranged in a periodic or non-periodic pattern, characterised by their lightweight and high strength-to-weight ratios [28]. Furthermore, they can replicate the porous nature of natural bone, which is crucial for facilitating cell attachment, proliferation, and differentiation [29]. It could be considered that further developments can be achieved following this line of research.

Among lattice architectures, the so-called Triply Periodic Minimal Surfaces (TPMS) are intricate structures widely utilised within scaffold design to ensure optimal porosity and interconnectivity [30]. These surfaces, known for their topological complexity, have attracted considerable attention due to their ability to be entirely defined through mathematical expressions. This mathematical foundation enables precise control over design parameters such as porosity and unit cell size, essential for adapting the scaffold properties to the desired medical applications [31,32]. TPMS exhibit exceptionally smooth surfaces devoid of sharp edges or peaks, facilitating fluid flow without causing abrasion or fluid build-up. Moreover, their high degree of interconnectivity enhances their effectiveness as supportive frameworks in tissue engineering and medical implants, highlighting their versatility and importance in modern biomedical research and practice [33].

Additionally, gyroid TPMS have emerged as highly promising architectures for bone tissue engineering because of their unique structural characteristics like surface area and interconnected porosity [34]. When compared against other popular TPMS based structures such as the Diamond and Primitive, gyroid based scaffolds showed the better compromise between higher bone ingrowth and bone-to-implant contact in vivo, and mechanical strength, as shown by [32]. Numerous researches have also been conducting exploring the design parameters and its influences on mechanical properties of gyroid based scaffolds [35,36,37].

Surface area directly influences several key aspects of bone regeneration, being a key property in porous scaffold architectures. A larger surface area enhances cell attachment and proliferation by providing more sites for cell adhesion, which is vital for the initial stages of bone formation [38]. It also enhances protein adsorption, nutrient and waste exchange [39] and enables better mechanical interlocking with native bone tissue, enhancing the scaffold’s stability and integration [40]. The surface topography of a scaffold can influence stem cell differentiation into osteoblasts, with a larger surface area promoting osteogenic differentiation [41]. Highly porous structures, often associated with high surface areas, provide pathways for vascularisation and bone ingrowth [42].

Another critical property of porous scaffolds for bone applications is their permeability, which quantitatively measures a porous medium’s ability to allow fluid flow [43]. Permeability is affected by porosity, pore size, orientation, tortuosity, and interconnectivity [44,45,46]. It is essential for describing tissue growth and the diffusion of nutrients and oxygen through scaffold pores. While highly dependent on the type of bone, generally high permeability values provide more desirable conditions for in vivo bone formation, while inadequate values can lead to the formation of cartilaginous tissue instead of bone [47].

Numerical simulations, including methods such as Finite Element Method (FEM) and Computational Fluid Dynamics (CFD), are critical tools in evaluating and optimizing the properties of scaffolds for bone regeneration. FEM allows researchers to model the mechanical behaviour of scaffolds under various loading conditions, providing insights into stress distribution, deformation, and potential failure points, which are essential for ensuring that the scaffold can withstand physiological loads while promoting proper load transfer to the regenerating bone [48]. FEM is also widely used to predict the mechanical response of the cranial implant [14,49]. CFD, enables the analysis of fluid flow within the scaffold’s porous structure, assess the permeability of the medium [50].

This study focuses on the optimisation of design of TMPS scaffolds for cranial implants, using computer simulation techniques for modelling the relationship between decision variables (TMPS structure parameters) and optimisation objectives (surface area and permeability) and constraints (porosity and elasticity). The optimal scaffold configuration was used to construct, then, a customised cranial implant prototype.

2. Materials and Methods

2.1. Experimental Design

The design method used for creating the samples was the parametric modelling, given that, of all the existing methodologies to construct TPMS-based structures, this one allows for greater control over the scaffold properties. These parametric equations are obtained from nodal approximations to the Fourier series expansion of minimal surfaces [51]. The following expression shows the modified gyroid equation used for our experimental setup:

$$\begin{array}{c}\hfill G(x,y,z)=sin\left(\frac{2\pi x}{\delta }\right)cos\left(\frac{2\pi y}{\delta }\right)+sin\left(\frac{2\pi y}{\delta }\right)cos\left(\frac{2\pi z}{\phi \delta }\right)+...\hfill \\ \hfill ...+sin\left(\frac{2\pi z}{\phi \delta }\right)cos\left(\frac{2\pi x}{\delta }\right)-c;\end{array}$$

(1)

where c represents the isovalue, a constant for defining and manipulating specific shape and characteristics of TPMS structures. When $c=0$, the obtained structure is divided into regions of equal volumes. As the value of c increases or decreases, the thickness of the structure changes, therefore being an important parameter for the control of the porosity of the scaffolds [52].

On the other hand, $\delta$ gives the size of the unit cell and $\phi$ was a size factor, introduced for changing the ratio of the cell size on the z-axis with respect to the other two axes.

Based on this equation, its three parameters (i.e., the isovalue, c; the cell size, $\delta$; and the size factor, $\phi$) were chosen as decision variables of the considered optimisation. In order to carry out the simulations required for modelling how the scaffold properties depend on these variables, the following experimental levels were considered for each on them (see Table 1).

A Box-Behnken experimental design was used for defining the level of each variable, for the simulations. The Box-Bahnken was chosen as it is particularly efficient because it requires fewer experiments than other alternatives (like the full factorial design). Box-Bahnken also avoids extreme conditions, focusing instead on a more central, balanced experimental region. Table 2 shows the values of the decision variables at each of the 15 experimental points.

The 15 experimental samples were constructed using an in-house MATLAB script and converted to STL format with the STLWRITE function. Figure 1 shows some of the obtained models when changing the magnitude of the decision variables.

2.2. Computing the Porosity

One of the studied scaffold properties is the porosity, which plays an important role when designing cranial implants, because it facilitates bone ingrowth and the formation of new tissues within the scaffold structure. Porosity, also, ensures that the scaffold does not act as a barrier to the infiltration of cells and growth factors, which are crucial for the successful integration of the implant with the host tissue.

In this study, the porosity, P, was computed from the previously obtained geometrical models of the scaffolds, through the expression:

$$P=\left(1-\frac{V_{\mathrm{sol}}}{V_{\mathrm{tot}}}\right)\times 100\%;$$

(2)

where $V_{\mathrm{sol}}$ is the volume of the solid TPMS structure and $V_{\mathrm{tot}}$ is the total volume of the prism (i.e., solid volume plus void volume). Both volumes were obtained by using the CAD software Meshmixer v3.5, through its stability analysis tool.

2.3. Computing the Surface Area

The surface area, A, is another important property of the scaffolds used for cranial implants, because it is beneficial for enhancing cell attachment and proliferation. This property provides more sites for osseous cells to attach to, promoting osteointegration and accelerating the healing process, ensuring the stability and functionality of the implant.

The surface area of each scaffold specimen was also determined from its 3d geometric model, by using the stability analysis tool of Meshmixer v3.5.

2.4. Computing the Young Modulus

The elasticity also plays a key role in cranial implants due to the need of ensuring the scaffold mimics the mechanical properties of natural bone. This compatibility minimizes stress shielding and promotes the integration and longevity of the implant within the cranial structure.

For computing the elasticity (Young) modulus, a Finite Element-based simulation was executed. For this purpose, the STL files were transformed into solid models, using FreeCAD v20.2. To carry out the analysis, two plates were modeled at both ends (top and bottom). With this configuration, they were exported in STEP format and imported into the static structural analysis module of ANSYS 2022R1.

The same boundary conditions were applied for all finite element models, as shown in Figure 2. The plates were scoped as rigid bodies, while the TPMS structure was defined as flexible. A force of 100 N was applied on the upper plate downward (i.e., in the negative direction of the z-axis), while the lower plate was constrained to keep it fixed. Frictional contact connections between the plates and the samples were established with a friction coefficient of 0.1, suitable for these cases [53].The scaffolds were meshed with tetrahedral elements with a minimum element size of 0.2 mm, while the two plates were meshed with hexahedral elements.

The scaffolds were assigned PEEK as material. As discussed previously, PEEK closely mimics the mechanical behaviour of natural bone, providing a stable and durable implant option. The parameters for the material selected are shown in Table 3, which are based on VICTREX PEEK 151G grade filament.

With the height deformation, $\Delta H$, obtained after the finite element simulation, the Young modulus is computed by using the equation:

$$E=\frac{F\Delta H}{V};$$

(3)

where F is applied force ($F=100\mathrm{N}$, in this study), and V is the volume of the scaffold specimen, equal to the product of its dimensions: length, L; width, B; and height, H:

$$V=LBH.$$

(4)

2.5. Computing the Permeability

Another important property of the scaffolds is the permeability, which facilitates the transport of nutrients and waste, supporting cell growth and tissue regeneration. High permeability ensures effective vascularisation, which is critical for healing and integrating the implant with surrounding biological tissues.

For the calculation of permeability, a laminar model of computational fluid dynamics (CFD) was designed and solved using the ANSYS Fluent module. The value of residual sensitivity criterion was set at ${10}^{-5}$ for an easier convergence. Water was chosen as the fluid for the analysis, with density, $\rho =1.0\mathrm{g}/{\mathrm{cm}}^3$, and dynamic viscosity, $\mu =0.001\mathrm{kg}/\mathrm{m}·\mathrm{s}$. The proposed model corresponds to an incompressible fluid with constant density and viscosity, and steady flow that obeys the Navier-Stokes equations:

$$\begin{array}{cc}& \rho \frac{\partial \mathbf{u}}{\partial t}-\mu {\nabla }^2\mathbf{u}+\rho (\mathbf{u}·\nabla )\mathbf{u}+\nabla p=\mathbf{F};\end{array}$$

(5)

$$\begin{array}{cc}& \nabla ·\mathbf{u}=0;\end{array}$$

(6)

where $\rho$ is the fluid density, $\mathbf{u}$ is the fluid velocity, t is the time, $\mu$ is the dynamic viscosity, p is the pressure and $\mathbf{F}$ represents the action of external forces (in this study, $\mathbf{F}=\mathbf{0}$). It should be remarked that this approach is appropriate to mimic the general behaviour of human body fluids [54].

For each design, the fluid domain was obtained by subtracting the 3D scaffold structure from a rectangular prism. A virtual volume was modelled at the top to extend the inlet; previous studies have shown how this helps to improve convergence and allows the fluid to fully develop [55,56]. Velocity at the inlet was set as $u_{\mathrm{in}}=1.0\mathrm{mm}/\mathrm{s}$, while the pressure at the outlet was assumed to be $p_{\mathrm{out}}=0$. No slip boundary condition was applied to the walls [57]. Figure 3 shows the final setup.

Given the height of the scaffold, H, the pressure drop, $\Delta p$, required to determine the permeability was calculated by subtracting the pressure at the top plane of the scaffold from the pressure at the outlet. In order to simplify the simulation and the necessary computational requirements, the entire domain of the model was set as the fluid domain. All models were meshed with tetrahedral 0.2 mm-sized elements.

With the values of volumetric flow, Q; and pressure drop, $\Delta p$, obtained from the CFD analysis, Darcy’s Law equation was used for the permeability calculation [58]:

$$K=-\frac{\mu QH}{A\Delta p};$$

(7)

where H is the height of the scaffold and A is its cross-section area.

2.6. Modelling

For modelling the relationship between the decision variables (i.e., cell size, isovalue and size factor) and the scaffold properties (i.e., porosity, surface area, Young modulus and permeability), linear regression models were fitted. When proper, some transformations, such as logarithms or linear combinations, were applied to the decision variables.

For all the obtained models, the goodness of fit was evaluated by computing the corresponding coefficient of determination. However, it cannot be assumed that a regression showing a high value of the coefficient of determination is implicitly a reliable model. Issues such as the distribution of the residuals, the statistical significance of the relationship between the variables and the statistical significance of the model coefficients, should be addresses to ensure the quality and usefulness of the models.

2.7. Optimisation

After the regression models were fitted, a multi-objetive optimisation was carried out for obtained the best scaffold designs. As it was previously stated, the decision variables were the TPMS parameters: cell size, $\delta$; isovalue, c; and size factor, $\phi$. All of them should be constrained to their respective experimental intervals:

$$\begin{array}{cc}& 2.0\mathrm{mm}\le \delta \le 3.0\mathrm{mm};\end{array}$$

(8)

$$\begin{array}{cc}& -0.8\le c\le 0.8;\end{array}$$

(9)

$$\begin{array}{cc}& 0.5\le \phi \le 1.5.\end{array}$$

(10)

The more relevant properties for the cranial implants were selected as optimisation targets. Maximising the surface area:

$$maxA=f_1(\delta ,c,\phi );$$

(11)

was selected as the first objective, in order to enhance the cell attachment and proliferation. The second objective selected was the maximisation of permeability:

$$maxK=f_2(\delta ,c,\phi );$$

(12)

for ensuring effective vascularisation.

In addition to the objectives, constraints were considered in the optimisation problem. Firstly, it was established that the porosity had to be in the interval from 60 to 80%:

$$60\%\le P=g_1(\delta ,c,\phi )\le 80\%;$$

(13)

to ensure bone formation and integration with the scaffold. This is because these conditions allow for better vascularisation and nutrient supply to the cells involved in bone regeneration. Scaffolds with this range of porosity facilitate bone ingrowth. This is critical for the long-term success of cranial implants, ensuring they remain integrated with the host bone and maintain their structural integrity over time.

On the other hand, the Young modulus of the scaffold had to be in the range of the replaced bone, i.e. from 0.5 to 5.0 GPa:

$$0.5\mathrm{GPa}\le E=g_2(\delta ,c,\phi )\le 5.0\mathrm{GPa}.$$

(14)

In order to homogenise the values of the constraints, equations 13 and 14 were rewritten as follows:

$$\begin{array}{cc}& {\gamma }_1=60/P-1\le 0;\end{array}$$

(15)

$$\begin{array}{cc}& {\gamma }_2=P/80-1\le 0;\end{array}$$

(16)

$$\begin{array}{cc}& {\gamma }_3=0.5/E-1\le 0;\end{array}$$

(17)

$$\begin{array}{cc}& {\gamma }_4=E/5.0-1\le 0.\end{array}$$

(18)

The total amount of unfeasibility, $\Gamma$, is then computed by the expression:

$$\Gamma =\sum _{i=1}^{4}max({\gamma }_i,0).$$

(19)

As both optimisation objective are conflictive (i.e., improving one of them worsens the other one.), an a posteriori approach was used, where the non-dominated solutions, which form the so-called Pareto front, are obtained and, then, the most convenient solution is chosen, depending on the user’s preferences.

The Non-Dominated Sorting Genetic Algorithm II (NSGA-II) [59] was chosen as optimisation technique, due to its efficiency in multi-objective optimisation, offering a blend of non-dominated sorting techniques, crowding distance techniques for diversity preservation, and elitist techniques to maintain high-quality solutions across generations. It achieves fast convergence to Pareto-optimal solutions, making it suitable for a wide range of problems.

Table 4 shows the NSGA-II parameters that were set for solving the addressed problem.

3. Results

3.1. Simulation Results

Porosity and surface area were numerically computed from the geometric model, as explained in Section 2.2 to Section 2.3. Their respective obtained values are shown in the second and third columns of Table 5. As can be seen, their values range from 22.73 to 79.45%, for the porosity, and from 1,003.1 to 2,570.0 ${\mathrm{mm}}^2$, for the surface area.

On the other hand, the Young modulus was computed from the outcomes of the finite element-based simulations, as explained in Section 2.4. Figure 4 shows and example of the obtained deformation and stress profiles, for one of the scaffolds specimens. The values of Young modulus, computed from the simulated mechanical behaviour of the scaffold specimens are shown in the fourth column of Table 5 and range from 0.038 to 2.357 GPa.

Finally, the permeability was computed, as explained in Section 2.5, from the pressure profile, obtained by CFD-based simulation (see Figure 5), by using the Darcy’s equation. The obtained values are shown in the fifth column of Table 5 and range from $0.637\times {10}^{-9}$ to $85.859\times {10}^{-9}{\mathrm{m}}^2$.

The permeability values obtained are in the order of ${10}^{-8}$ to ${10}^{-10}$ which align with the results reported in the consulted scientific literature [55,57,60]. The scaffolds with the lower values of the isovalue constant were the ones that presented higher permeability values. It can be observed how cell size factor of 0.5 can impact negatively on the permeability results, while the opposite results in an increase of this property.

3.2. Regression Models

After applied the regression technique, as explained in Section 2.6, linear models were obtained for each studied dependent variable.

For the porosity, the model was

$$P=50.18-34.14C.$$

(20)

This model has a determination coefficient equal to $99.9\%$ and a standard error of the estimate of 0.6313. It should be remarked that only the C variable remains in the model as all the other independent variables were removed, through a backward stepwise regression, due to, after the t-Student test, they were no statistically significant. Figure 6-a shows the graphical representation of the obtained model.

For the surface area, the model was obtained after applying a logarithmic transformation was applied to all variables, both independent and dependent. In the case of the isovalue, C, as part of the transformation, was increased by one unit, to make all its values positive. The model is given by the expression:

$$A=\frac{3,681{(C+1)}^{0.1138}}{{\delta }^{1.006}{\phi }^{0.3650}};$$

(21)

and has a determination coefficient of 94.6% and a standard error of the estimate equal to 0.06750. In this model, all the independent variables are significant. Consequently, they were kept into the model. Figure 6-b shows the graphical representation of the computed model.

On the other hand, the model for the Young modulus was obtained by considering a second order interaction between the variables and applying a backward stepwise regression, for removing the non-significant variables. The obtained model has, then, the expression:

$$E=0.8143\phi +0.9436C\phi ;$$

(22)

with a determination coefficient of 90.9% and standard error of the estimate 0.3546. Figure 6-c shows the model.

Finally, the model of the permeability was fitted after a logarithmic transformation, obtaining the equation:

$$K=\frac{1.401{\delta }^{1.887}{\phi }^{1.367}}{{(C+1)}^{1.281}};$$

(23)

with a determination coefficient of 84.6% and a standard error of the estimate of 0.6254. All the variables are significant and, consequently, they were kept into the model. Model is depicted in Figure 6-d.

In all these model there is a 99% confidence level-significant relationship between the variables, given from the ANOVA table.

3.3. Optimisation Results

After the optimisation process, carried out as explained in Section 2.7, the corresponding Pareto front and set were obtained. The Pareto front (i.e., the values of the objective for the non-dominated solutions – see Figure 7-a) is composed by two curves: A first segment going from point A to point B, and a second one, going from point C to point D. Nevertheless, for practical purposes, points B and C are close enough to be considered the same point.

It should be noted that the Pareto portfolio values include only a part of the target ranges. In this sense, the optimal surface area ranges from 969 to 1,817 ${\mathrm{mm}}^2$, with their respective values as predicted by the model ranging from 875 to 2,500 ${\mathrm{mm}}^2$. With respect to the permeability, its optimal values range from $7.2\times {10}^{-9}$ to $48.4\times {10}^{-9}{\mathrm{m}}^2$, while the corresponding model gives values from $0.4\times {10}^{-9}$ to $152.4\times {10}^{-9}{\mathrm{m}}^2$.

On the other hand, the Pareto set (i.e, the values of the decision variables corresponding to the non-dominated solutions) is also composed by two segment, the first one, from point A to point B, included solutions with a constant cell size, $\delta =2.0\mathrm{mm}$, and isovalues decreasing from $C=-0.29$ to $C=-0.50$, and values of shape factor increasing from $\phi =0.42$ to $\phi =1.49$.

The second segment, from point C to point D, on the other hand, included solutions with constant isovalue, $C=-0.51$, and shape factor, $\phi =1.50$, while the cell size increases from $\delta =2.0$ to $\delta =3.0$. In can be noted that, while the optimal cell sizes ranges over its whole interval, $\delta =2.0\mathrm{mm}$ to $\delta =3.0\mathrm{mm}$, the shape factor only ranges from $\phi =0.92$ to $\phi =1.50$, and the isovales barely covers a small subinterval, from $C=-0.29$ to $C=-0.51$.

4. Discussion

The first issue to be remarked when discussing the obtained results is that all the points shown in Figure 7 corresponds to solutions which are optimal in the wide sense that no other solution, in the feasible search space, can improve one of the objective without worsening the other one. Consequently, any of these points can be used for designing the optimal scaffold structure for a cranial implant. The finally chosen solution should depend on which goal (or goal combination) will be preferred by the decision-maker.

Under this point of view, there are three specially interesting points (see Figure Table 6): point A represents the best permeability values and (as they are contradictory) the worst surface area. It should be chosen if enhancing cell attachment and proliferation is the prioritised criterion. On the contrary, if facilitating the transport of nutrients and wastes, and ensuring effective vascularisation is the most important goal, point D should be chosen, as it gives the higher value of permeability. Finally, any other point (including B and C) can be considered as trade-off solutions, that can be chosen if both goals are to be considered together.

As an example, it was designed a cranial implant for the skull defect shown in Figure 8-a. The optimized custom cranial implant was designed using the optimal parameters corresponding to point B: unit cell size, $\delta =2.0\mathrm{mm}$; isovalue, $C=-0.5$; and shape factor, $\phi =1.49$. Point B was selected as it is the most balanced option for the cranial implant design due to its combination of high permeability, adequate surface area, and favourable porosity, which promotes better biological integration and mechanical stability.

As an example, it was designed a cranial implant for the skull damage shown in Figure 8-a. The optimised custom cranial implant was designed in Ntopology v4.20.2 software using a custom infill volume box, unit cell adjustments and thickness were selected according to the optimal scaffold settings. Figure 8-b shows the final results for the designed prosthesis and Figure 8-c, the location of the implant showing the fitting to the skull.

An implant prototype was finally additively manufactured (see Figure 9) for validating the technological feasibility of the processes.

5. Conclusions

The main outcome of the work is the methodology to optimise the design of cranial implants based on TPMS scaffolds. These approach considers two optimisation targets: the surface area and the permeability. Both are key issues in cranial implants as the first one promotes the cellular adhesion and the mechanical stability of the implant, while the second one enhances the cellular adhesion and bone growth. Nevertheless, they are conflictive and cannot be simultaneously optimised.

The proposed a posteriori approach, where a set of non-dominated solutions (also known as Pareto front) are obtained through the NSGA-II, showed to be effective, as it offers a set of solutions that can be considered optimal as they are those for which no other solution can be found, in the feasible search space, that simultaneously improve both objectives. It can be remarked that the proposed optimisation approach guaranties the fulfilment of the considered constraints: the elasticity, which minimises stress shielding, and porosity, which promotes the formation of new bone tissue within the implant structure.

The presented case study allowed not only to demonstrate the presented approach but also its validation through the design of an implant for a given case.

As a future development of this work, firstly, the practical validation of the results obtained can be considered, based on experimental measurements of stiffness, porosity, permeability and surface area, in order to verify the accuracy of the simulations carried out. The inclusion of manufacturing parameters as decision variables can also be considered, in order to study their influence on the variables studied and, consequently, on the results of the optimisation process.