Application of particle swarm optimization in inverse finite element modeling to determine the cornea ́ s mechanical behavior

Particle Swarm Optimization (PSO) was foregrounded by finite element (FE) modeling to predict the material properties of the human cornea through inverse analysis. Experimental displacements have been obtained for corneas of a donor approximately 50 years old, and loaded by intraocular pressure (IOP). FE inverse analysis based on PSO determined the material parameters of the corneas with reference to first-order, Ogden hyperelastic model. FE analysis was repeated while using the commonly-used commercial optimization software HEEDS, and the rates of the same material parameters were used to validate PSO outcome. In addition, the number of optimization iterations required for PSO and HEEDS were compared to assess the speed of conversion onto a global-optimum solution. Since PSO-based analyses produced similar results with little iteration to HEEDS inverse analyses, PSO capacity in controlling the inverse analysis process to determine the cornea material properties via finite element modeling was demonstrated.

Current study assessed a simple inversemodeling implementation of PSO in the application of Biomechanics and compared to the commercial software HEEDS.Commercial software HEEDS has a strong record for the optimization of solutions in mechanics problems by employing a hybrid and adaptive algorithm.HEEDS employs multiple search strategies and adapts to the problem.In fact, it requires significantly fewer model evaluations to identify optimized solution for the first time.
When compared to commercial package HEEDS, PSO is simple to implement and may be coded by research groups with limited resources and fine-tuned to suit their specific needs.The comparisons presented in current paper are related to an issue in which the mechanical properties of the corneas of the 50-year-old donor are estimated in an inverse modeling exercise.

Objective function definition
In previous studies, cornea samples were tested under inflation conditions (Elsheikh et al, 2007).Figure 01 shows the hyperelastic behavior of human corneas under increasing posterior pressure.(Elsheikh, Alhasso, & Rama, 2008) Based on experimental data from human corneas (Elsheikh et al., 2008) where: x is the intraocular pressure and y is the displacement.Equation 1 was used for inverse modeling procedures for PSO and HEEDS methods to compare the rates of material parameters obtained and the speed of the loading process in general.

Finite element model
An FE model of a human cornea involving 12,168 quadratic, triangular, prismatic, hybrid elements, arranged in 3 radial layers, was generated and used in current study (Figure 2).The model was maintained along the edge nodes to simulate the conditions experienced by the donor´s corneas in the experimental program.The model included fluid elements representing the aqueous humor of the anterior chamber, taking an internal pressure of up to 45 mm Hg.The model had an anterior geometry that was sinusoidal, with a central radius of curvature of the anterior surface, Rant = 7.8 mm, a shape factor, p = 0.82, central corneal thickness, CCT = 0.545 mm and a peripheral corneal thickness, PCT = 0.695 mm.The above rates represented the average results arrived at experimentally for the tested corneas and made Acta Scientiarum.Technology Maringá, v. 39, n. 3, p. 325-331, July-Sept., 2017 possible the average values in the numerical modeling conducted in current study.

Inverse modeling
For simplification, corneal tissue was presumed to have an isotropic and hyperelastic behavior, represented by Ogden's material model (Ogden, Saccomandi, & Sgura, 2004): where: W is the strain energy; μ i and  i are material parameters; i = 1 to equation order N; j  are the stretches in the Cartesian system, j = x, y, z.
The material parameters to be defined from an inverse modeling procedure involve the finite element analysis, in which the input includes the cornea's geometry and the pressure-deformation behavior, while the expected output is the material behavior represented by μ i and  i .The solution comprises an optimization process that identifies the global minimum of the difference between the experimentally-obtained and the numericallypredicted, load-deformation behavior at various points on the corneal surface.This process commonly utilizes commercial tools such as HEEDS or, as in this case, optimization methods such as PSO.In current study, the order of the strain energy material model was assumed as 1 since, in previous studies, this order was adequate to describe experimental behavior of biological tissue (Elsheikh, 2010).Limiting the order to 1 meant the optimizing of rates of only two parameters to reduce cost of analysis.

Particle Swarm Optimization
PSO is an optimization technique based on populations with m particles (m individuals) that evolve within the hyperspace defined by the design´s variable bounds following some random criteria towards the particle with the best performance (usually the particle that is closest to the optimum point) (Barbieri, Barbieri, & Lima, 2015).The PSO algorithm is an alternative which describes each particle (material parameters) by its vector position in the search space as: where: w factor is the inertia weight; v (t) is the particle's velocity at time t; x (t) is the particle's current position at time t; c 1 and c 2 are weight constants; p (t) is the particle's best position; g (t) is the best known position found by any particle in the swarm; r 1 and r 2 are random numbers resulting from an algorithm choice that normally vary between 0 and 1 (Marwala, 2005).

1) +
Since in the Equation (2), μ i and  i are the material parameters to be defined, a finite element solver Abaqus/Standard 6.12 (Dassault Systèmes Simulia Corp., Rhode Island, USA) was used with PSO.For the numerical analysis, Abaqus was used with a PSO code written in Visual Basic (VB) to provide the best fitness between FE results and experimental data.A Python code was written to extract the displacement curve in z direction of the apical node of the cornea, using FE.The PSO code used this target curve to calculate the cost function as the sum of squared error (SSE) between FE displacements from the cornea apex and experimental data (Equation 1).The algorithm takes the discrete data points of a target curve, fit the curve to the data and compare the equation of the curve with FE results.A summary of the PSO process, specific for this kind of problem, is shown in a flow chart (Figure 3).Boundary rates were defined to constrain the search space between maximum and minimum rates for each PSO iteration.This limitation is required since rates that exceed maximum and minimum boundaries lead to unrealistic material parameters and convergence issue for the FE solver.
Ten particles, parameters c 1 and c 2 equal to 1.494, w equal to 0.729 and population size equal to 30 were set based on previous studies.Rates were recommended by Clerc (2013) and tested by Tuppadung & Kurutach (2011). Trelea (2003) used them in a large number of simulation experiments for convergence and sensitivity analysis, with good results.In case of non-converged simulation in the proposed algorithm, a specific line was added in the PSO code to discard the iteration in which errors occurred.

Results and discussion
Material parameters obtained from PSO have been compared with those by HEEDS.Analysis was performed with six random boundaries situations (Table 1), justified by the variation of the limits (maximum and minimum) of the parameters to certify the stability of the proposed methodology.
Figure 4 shows graphs from PSO and HEEDS results for the first boundary condition.Similarly, the other boundary conditions were also evaluated (Table 1).Figure 4 shows fewer iterations for PSO when compared with those for HEEDS.It should be underscored that each iteration displayed on the graphs represents the position of ten particles.This means that the number of FE simulations tested was the same for both methods.Taking the above into consideration coupled to HEEDS rates in ascending order, PSO achieved the optimized rate first.Figure 5 compares SSE rates for the two methods with boundary 2 (Table 1) as example.In other words, each method took advantage of the best attributes and solutions found from other parallel searches (Abedrabboa, Worswicka, Mayer, & Riemsdijkc, 2009).
Regarding to PSO convergence, it may be observed that running the same boundary condition more than once, all rates after each iteration remained the same even in different PC.Consequently, six different boundaries were implemented to test the convergence.Figure 6 represents the convergence of parameters 1  and 1  for the six boundaries tested (from Table 1).
Even though the benefits of PSO extend to other cases from a specific case (healthy 53-yearold human corneas modeled with Ogden constitutive law and loaded to 45 mm Hg), it may be noted that PSO is more stable than HEEDS.On the other hand, HEEDS software provided a graphical user-friendly interface easy to use and did not need any programming knowledge, or rather, assets benefitted HEEDS.
Ogden model is highly useful in the cornea for material characterization (Elsheikh et al., 2007;Lago et al., 2015), even though it has certain disadvantages, such as isotropy.One option is the use of PSO optimization scheme for other models (e.g., Holzapfel model, Fung model) not used here.It is common knowledge that the hyperelastic material model is not the only choice to represent the corneal mechanical response, whilst viscoelasticity (Lombardo, Serrao, Rosati, & Lombardo, 2014) and anisotropy (Whitford, Studerb, Booteb, Meekc, & Elsheikh, 2015) should be considered in future research works with PSO.

Conclusion
Current research focused on a novel application for Particle swarm optimization coupled to finite element modeling to predict the material properties of the human cornea through inverse analysis.One advantage of the PSO code is that it does not need the initial guess since it chooses random rates between the boundary values established by the user.However, HEEDS runs simulations in parallel with several different optimization algorithms.At the end of all simulations, it chooses the best solution from the parallel searches.PSO depends on the previous particle's position, generating a limitation in the way that the proposed application has been implemented in terms of time consuming.
It is known that results´ stability with inverse analysis used is also an advantage.Stability is mostly a function of the change in results from different initial guesses.In current case, results demonstrated no concern with PSO with regard to stability since six different boundary conditions have been used with the same results after optimization.
PSO has also the advantage that it converged first, considering the same number of iterations as HEEDS.Further, PSO is an open source code whereas HEEDS is a commercial software.Taking all the above aspects into consideration, the capacity of PSO in controlling the inverse analysis process to determine the cornea´s material properties via finite element modeling should be underscored.

Figure 3 .
Figure 3. Flow chart of PSO process for the cornea´s definition of material properties.

Figure 5
Figure5revealed that, by SSE values, PSO converged after 400 simulations, or rather, prior to HEEDS which converged after 500 simulations.However, HEEDS ran simulations in parallel using a variety of optimization algorithms such as Multiobjective SHERPA, Genetic algorithm, Quadratic programming, Simulated annealing, Response surface, Multi-start local search and Nelder-Mead simplex.In other words, each method took advantage of the best attributes and solutions found from other parallel searches(Abedrabboa, Worswicka, Mayer, & Riemsdijkc, 2009).Regarding to PSO convergence, it may be observed that running the same boundary condition more than once, all rates after each iteration remained the same even in different PC.Consequently, six different boundaries were implemented to test the convergence.Figure6represents the convergence of parameters 1  and 1  for the six boundaries tested (from Table1).Even though the benefits of PSO extend to other cases from a specific case (healthy 53-yearold human corneas modeled with Ogden constitutive law and loaded to 45 mm Hg), it may be noted that PSO is more stable than HEEDS.On the other hand, HEEDS software provided a graphical user-friendly interface easy to use and did not need any programming knowledge, or rather, assets benefitted HEEDS.Ogden model is highly useful in the cornea for material characterization(Elsheikh et al., 2007;Lago et al., 2015), even though it has certain disadvantages, such as isotropy.One option is the use of PSO optimization scheme for other models (e.g., Holzapfel model, Fung model) not used here.It is common knowledge that the hyperelastic material model is not the only choice to represent

Table 1 .
Boundary material parameters and results.