# Application of particle swarm optimization in inverse finite element modeling to determine the cornea's mechanical behavior/Aplicacao da otimizacao por enxame de particulas na modelagem inversa de elementos finitos para determinar o comportamento mecanico da cornea.

IntroductionParticle Swarm Optimization (PSO), based on socio-psychological principles inspired by swarm intelligence and social behavior, has been widely applied to engineering (Marwala, 2010). PSO method was first forwarded in 1995 (Poli, Kennedy & Blackwell, 2007) and soon found applications in several areas (Liu, Liu, & Duan, 2007; Liu, Liu, & Cartes, 2008; Arumugam & Rao, 2008; Berlinet & Roland, 2008; Thakker, Patil, & Anil, 2009; Santos, Martins, & Santos, 2012).

Marwala, Boulkaibet, & Adhikari (2016) combined PSO with finite element (FE) modeling to get the best fit for numerically-predicted load-deformation behavior with experimental data obtained for a simple beam. More demanding structural mechanics applications were undertaken in later studies which demonstrated PSO reliability. Or rather, it may effectively be combined with FE methods in inverse modeling problems (Tang, Chen, & Peng, 2009; Oliveira eta al., 2011; Chen, Tang, Ge, An, & Guo, 2013; Reutlinger, Burki, Brandejsky, Ebert, & Buchler, 2014).

Other optimization methods such as Genetic Algorithms (GA) (Shabani & Mazahery, 2013) have been employed with FE (Song, Yang, Yu, & Kang, 2013) or with both (Herath, Natarajan, Prusty, & John, 2014) in optimization problems. PSO has also been compared with the Artificial Bee Colony (ABC) algorithm, demonstrating higher accuracy and better reliability in inverse modeling applications involving FE analysis (Bardsiri, Jawawi, Hashim, & Khatibi, 2013).

The determination of the mechanical properties of human organs, including the eye, using FE inverse modeling, has been quite common, specifically within Biomechanics (Wittek et al., 2013; Perez, Morris, Hart, & Liu, 2013). Research in ocular Biomechanics, focusing mainly on the cornea's and sclera's mechanical properties (Elsheikh et al., 2007; Elsheikh, Ross, Alhasso, & Rama, 2009; Elsheikh, Geraghty, Rama, Campanelli, & Meek, 2010; Abyaneh, Wildman, Ashcroft, & Ruiz, 2013) has progressed significantly, although authors used optimization techniques different from PSO.

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.

Material and methods

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.

Based on experimental data from human corneas (Elsheikh et al., 2008), the objective function has been obtained and represented by a five-order polynomial equation 1:

y = (3.33217 Ell)[x.sup.5]-(6.069929378E9)[x.sup.4] + (4.2467142 76 E7)[x.sup.3]-(145325.93 35)[x.sup.2] + (272.0315029)x + 0.00209844 2 (1)

where:

x is the intraocular pressure and [gamma] 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 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):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2)

where:

W is the strain energy;

[[mu].sub.i] and [[alpha].sub.i] are material parameters;

i = 1 to equation order N;

[[lambda].sub.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 [mu]i and [[alpha].sub.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:

v(t+l) = (w x v(t)) + ([c.sub.1] x [r.sub.1] x (p(t)-x(t)) + ([c.sub.2] x [r.sub.2] x (g(t)-x(t)) (3)

where:

w factor is the inertia weight;

v(t) (t) is the particle's velocity at time t;

x(t) is the particle's current position at time t;

[c.sub.1] and [c.sub.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.sub.1] and [r.sub.2] are random numbers resulting from an algorithm choice that normally vary between 0 and 1 (Marwala, 2005).

By applying Equation (3), the new velocity, v (t+1), is determined to compute the new particle's position, x (t+1) as shown in Equation (4).

x(t+l) = x(t)+v(t+l) (4)

Since in the Equation (2), [[mu].sub.i] and [[alpha].sub.i] are the material parameters to be defined, a finite element solver Abaqus/Standard 6.12 (Dassault Systemes 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.sub.1] and [c.sub.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 often 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.

Figure 5 revealed 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. Figure 6 represents the convergence of parameters [[mu].sub.1] and [[alpha].sub.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.

Doi: 10.4025/actascitechnol.v39i3.29884

Acknowledgements

Current research was funded CAPES under Grant 2623-13-7.

References

Abedrabboa, N., Worswicka, M., Mayer R., & Riemsdijkc, I.V. (2009). Optimization methods for the tube hydroforming process applied to advanced highstrength steels with experimental verification. Journal of Materials Processing Technology, 209(1), 110-123.

Abyaneh, M. H., Wildman, R. D., Ashcroft, I. A., & Ruiz, P. D. (2013). A hybrid approach to determining cornea mechanical properties in vivo using a combination of nano-indentation and inverse finite element analysis. Journal of the Mechanical Behavior of Biomedical Materials, 27(1), 239-248.

Arumugam, M. S., & Rao, M. V. C. (2008). Molecular docking with multi-objective particle swarm optimization. Applied Soft Computing, #(1), 666-675.

Barbieri, R., Barbieri, N., & Lima, K. F. (2015). Some applications of the PSO for optimization of acoustic filters. Applied Acoustics, 89(1), 62-70.

Bardsiri, V. K., Jawawi, D. N. A., Hashim, S. Z. M., & Khatibi, E. (2013). A PSO-based model to increase the accuracy of software development effort estimation. Software Quality Journal, 21(3), 501-526.

Berlinet A, & Roland C. (2008). A novel particle swarm optimization algorithm for permutation flow-shop scheduling to minimize makespan. Chaos, Solitons & Fractals, 35(5), 851-861.

Chen, J., Tang, Y., Ge, R., An, Q., & Guo, X. (2013). Reliability design optimization of composite structures based on PSO together with FEA. Chinese Journal of Aeronautics, 26(2), 343-349.

Clerc, M. (2013). Particle Swarm Optimization. New York, NY: John Wiley & Sons. Elsheikh, A. (2010). Understanding corneal biomechanics through experimental assessment and numerical simulation. In J. H. Levy (Ed.). Biomechanics: principles, trends and applications (p. 57-110). New York, NY: Nova Science Publishers Incorporated.

Elsheikh, A., Alhasso, D., & Rama, P. (2008). Biomechanical properties of human and porcine corneas. Experimental Eye Research, 86(5), 783-790.

Elsheikh, A., Geraghty, B., Rama, P., Campanelli, M., & Meek, K.M. (2010). Characterization of age-related variation in corneal biomechanical properties. Journal of the Royal Society Interface, 7(51), 1475-1485.

Elsheikh, A., Ross, S., Alhasso, D., & Rama, P. (2009). Numerical study of the effect of corneal layered structure on ocular biomechanics. Current Eye Research, 34(1), 26-35.

Elsheikh, A., Wang, D., Brown, M., Rama, P., Campanelli, M., & Pye, D. (2007). Assessment of corneal biomechanical properties and their variation with age. Current Eye Research, 32(1), 11-19.

Herath, M. T., Natarajan, S., Prusty, B. G., & John, N. S. (2014). Smoothed finite element and genetic algorithm based optimization for shape adaptive composite marine propellers. Composite Structures, 109(1), 189-197.

Lago, M. A., Ruperez, M. J., Martinez-Martinez, F., Monserrat, C., Larra, E., Guell, J. L., & PerisMartinez, C. (2015). A new methodology for the in vivo estimation of the elastic constants that characterize the patient-specific biomechanical behavior of the human cornea. Journal of Biomechanics, 48(1), 38-43.

Liu, L., Liu, W., & Cartes, D. A. (2008). Particle swarm optimization-based parameter identification applied to permanent magnet synchronous motors. Engineering Applications of Artificial Intelligence, 88(21-22), 1092-1100.

Liu, X., Liu, H., & Duan, H. (2007). Particle swarm optimization based on dynamic niche technology with applications to conceptual design. Advances in Engineering Software, 38(10), 668-676.

Lombardo, G., Serrao, S., Rosati, M., & Lombardo, M. (2014). Analysis of the Viscoelastic Properties of the Human Cornea Using Scheimpflug Imaging in Inflation Experiment of Eye Globes. Plos One, 9(11), e112169.

Marwala, T. (2005). Finite element model updating using particle swarm optimization. International Journal of Engineering Simulation, 6(2), 25-30.

Marwala, T. (2010). Finite-element-model updating using computional intelligence techniques: applications to structural dynamics. London, UK: Springer-Verlag London.

Marwala, T., Boulkaibet, I., & Adhikari, S. (2016) Model selection in finite element model updating, in probabilistic finite element model updating using bayesian statistics: applications to aeronautical and mechanical engineering. Chichester, UK: John Wiley & Sons, Ltd.

Ogden, R. W., Saccomandi, G., & Sgura, I. (2004). Fitting hyperelastic models to experimental data. Computational Mechanics, 34(6), 484-502.

Oliveira, M. E., Hasler, C. C., Zheng, G., Studer, D., Schneider, J., & Buchler, P. (2011). A multi-criteria decision support for optimal instrumentation in scoliosis spine surgery. Structural and Multidisciplinary Optimization, 45(6), 917-929.

Perez, B., Morris, H., Hart, R., & Liu, J. (2013). Finite element modeling of the viscoelastic responses of the eye during microvolumetricchanges. Journal of Biomedical Science and Engineering, 6(12A), 29-37.

Poli, R., Kennedy, J., & Blackwell, T. (2007). Particle swarm optimization: An overview. Swarm Intelligence, 1(1), 33-57.

Reutlinger, C., Burki, A., Brandejsky, V., Ebert, L., & Buchler, P. (2014). Specimen specific parameter identification of ovine lumbar intervertebral discs: On the influence of fibre-matrix and fibre-fibre shear interactions. Journal of the Mechanical Behavior of Biomedical Materials, 30(1), 279-289.

Santos, R. P. B., Martins, C. H., & Santos, F. L. (2012). Simplified particle swarm optimization algorithm. Acta Scientiarum. Technology, 34(1), 21-25.

Shabani, M. O, & Mazahery. A. (2013). Application of GA to optimize the process conditions of Al Matrix nanocomposites. Composites Part B: Engineering, 45(1), 185-191.

Song ,J. H., Yang, S. M., Yu, H. S., & Kang, H. Y. (2013). Fatigue strength under optimal conditions of spot weldment using GA and FEM. Applied Mechanics and Materials, 284-287(1), 123-126.

Tang, Y., Chen, J., & Peng, W. (2009). Probabilistic optimization of laminated composites considering both ply failure and delamination based on PSO and FEM. Tsinghua Science and Technology, 14(S2), 89-93.

Thakker, R. A., Patil, M. B., & Anil, K. G. (2009). Parameter extraction for PSP MOSFET model using hierarchical particle swarm optimization. Engineering Applications of Artificial Intelligence, 22(2), 317-328.

Trelea, I.C. (2003). The particle swarm optimization algorithm: convergence analysis and parameter selection. Information Processing Letters, 85(6), 317-325.

Tuppadung, Y., & Kurutach, W. (2011). Comparing nonlinear inertia weights and constriction factors in particle swarm optimization, International Journal of Knowledge-based and Intelligent Engineering Systems, 15(2), 65-70.

Whitford, C., Studerb, H., Booteb, C., Meekc, K. M., & Elsheikh, A. (2015). Biomechanical model of the human cornea: considering shear stiffness and regional variation of collagen anisotropy and density. Journal of the Mechanical Behavior of Biomedical Materials, 42(1), 76-87.

Wittek, A., Karatolios, K., Bihari, P., Schmitz-Rixen, T., Moosdorf, R., Vogt, S., & Blase, C. (2013). In vivo determination of elastic properties of the human aorta based on 4D ultrasound data. Journal of the Mechanical Behavior of Biomedical Materials, 27(1), 167-183.

Received on February 2, 2016.

Accepted on June 20, 2016.

Ricardo Magalhaes (1) *, Ahmed Elsheikh (2), Philippe Buchler (3), Charles Whitford (2) and Junjie Wang (2)

(1) Departamento de Engenharia, Universidade Federai de Lavras, Cx. Postal 3037, 37200-000, Lavras, Minas Gerais, Brazil. (2) Faculdade de Engenharia, Universidade de Liverpool, Liverpool, Reino Unido. (3) Instituto de Tecnologia Cirurgica e Biomecanica, Universidade de Berna, Berna, Suica. * Author for correspondence. E-mail: ricardorm@deg.ufla.br

Caption: Figure 1. Pressure-deformation of 53-year-old human corneas (Elsheikh, Alhasso, & Rama, 2008)

Caption: Figure 2. Cornea FEM model and clamping position.

Caption: Figure 3. Flow chart of PSO process for the cornea's definition of material properties.

Caption: Figure 4. PSO and HEEDS results (Boundary 1) for: a) [[mu].sub.i], PSO; b) [[mu].sub.1] HEEDS; c) [[alpha].sub.1] PSO; d) [a.sub.1] HEEDS; e) SSE PSO; f) SSE HEEDS.

Caption: Figure 5. Results of PSO vs. HEEDS SSE for boundary 2.

Caption: Figure 6. PSO parameters ([[mu].sub.1] and [[alpha].sub.1]) convergence for the six boundaries tested.

Table 1. Boundary material parameters and results. [[mu].sub.i] [[alpha].sub.i] [[mu].sub.i] Boundary results (min-max) (min-max) (PSO) 1 0.00001-1 20-400 0.054 2 0.00001-1 10-200 0.054 3 0.02-0.2 20-150 0.054 4 0.01-0.1 60-240 0.054 5 0.025-0.075 60-150 0.054 6 0.001-0.8 40-300 0.054 [[mu].sub.i] [[alpha].sub.i] [[alpha].sub.i] SSE Boundary results results results results (HEEDS) (PSO) (HEEDS) (PSO) 1 0.055 110.45 109.3 0.0016 2 0.055 110.43 109.75 0.0016 3 0.054 110.46 110.35 0.0016 4 0.054 110.45 110.40 0.0016 5 0.054 110.45 110.40 0.0016 6 0.053 110.45 111.50 0.0016 SSE Boundary results (HEEDS) 1 0.0017 2 0.0016 3 0.0016 4 0.0016 5 0.0016 6 0.0017