Use of simplified geometrical models of a cornea for optimization purposes

Abstract

Particle Swarm Optimization (PSO) is an artificial intelligence technique (AI) that can be used to find approximate solutions to numerical problems of maximization and minimization. In this study, it was used a PSO algorithm to compare displacements from human cornea sample subjected to internal pressure of 45 mmHg with Results of numerical simulations were provided which identified optimized values for hyperelastic properties of the cornea (µ and α). By means of the results from numerical simulations via inverse analysis by the Finite Element Method (FEM), in conjunction with the PSO algorithm, optimized values of µ = 0.047 and α = 106.7 were found. When compared with optimized results from commercial software, errors around 0.15% were found. Results showed that, varying the values of particle inertia coefficients in the PSO algorithm, simulated displacements have improved when compared to experimental data. This demonstrates the potential use of PSO algorithm in conjunction with the FEM inverse analysis for hyperelastic materials characterization, using simplified geometrical models

Keywords:
Particle Swarm optimization; Finite element method; Cornea; Human

Resumo

Otimização por Enxame de Partículas (PSO) é uma técnica de inteligência artificial (AI), que pode ser usada para encontrar soluções aproximadas para problemas numéricos de maximização e minimização extremamente difíceis. Neste trabalho, utilizou-se um algoritmo PSO para comparar os deslocamentos sofridos por uma amostra de córnea humana submetida à uma pressão interna de 45 mmHg com resultados de simulações numéricas e identificar valores otimizados para propriedades hiperelásticas da córnea (µ e α). Por meio dos resultados das simulações via análise inversa pelo Método dos Elementos Finitos (MEF), em conjunto com o algoritmo PSO, foram encontrados valores otimizados de µ = 0,047 e α = 106,7. Quando comparado com resultados otimizados por meio de um software comercial, foram encontrados erros de aproximadamente 0,15%. Por meio dos resultados obtidos, verificou-se ainda que, variando os valores dos coeficientes de inércia da partícula no algoritmo PSO, os resultados podem sofrer ligeira melhoria, o que demonstra potencial uso do PSO em conjunto com análise inversa do MEF para caracterização de materiais hiperelásticos, utilizando modelos geométricos simplificados

Descritores:
Otimização por enxame de partículas; Método dos elementos finitos; Córnea; Humano

Introduction

The development and improvement of optimization techniques allows studying behavior in various materials, and is really important in determining their limitations, evaluating and classifying their properties. Thus, there is a great interest in the evaluation of the biomechanical properties of the cornea, being considered an important parameter to be determined, since it is related to several diagnostic procedures.

The cornea is the first and most powerful refraction surface of the optical system of the eyeball. Little is known about the behavior of the internal structure of the cornea. In most animal species, the cornea consists of four layers: epithelium, stroma, Descemet's membrane, and endothelium⁽¹⁾. A small change in the geometry of the cornea strikingly affects its optical capacity. Because of this sensitivity, biomechanical studies may reveal the performance and functionality of the cornea, which has certain mechanical characteristics that can be analyzed with a view to improving the methods for prevention of diseases and the development of more efficient equipment in the medical field.

The growing interest in the study of human biological tissues has led to the need to characterize some types of tissue with a hyperelastic mechanical behavior. Experimentally, the mechanical study of this type of tissue has been carried out by means of traditional techniques developed for the study of materials with a linear and isotropic behavior, which allows the determination of average values of some mechanical properties.

These materials have mechanical behaviors different from conventional materials, and it is necessary to use algorithms that allow the solution of these problems, one of them being Particle Swarm Optimization. Soares⁽²⁾ defines optimization as a mechanism for analyzing complex decisions involving selection of values for variables, with the simple objective of quantifying the performance and measuring the quality of decisions.

The algorithm implemented makes use of the deformation suffered by a specimen (human cornea), aiming to identify the properties of the hyperelastic material. The characterization of the material tested was performed using the Finite Element Method.

Theoretical Reference

The Particle Swarm Optimization method was developed by Kennedy et al.⁽³⁾ from the mathematical modeling that implements a metaphor of the social behavior of a group of birds in search of food or a place to build the nest, creating a technique of stochastic computation.

According to Parsopoulos et al.⁽⁴⁾, each particle is treated as a point within the search space adjusting its own "flight" according to its own experience, as well as the experience of the "flight" of other particles.

The Finite Element Method was used to simulate different mechanical properties of the materials and to provide displacement data at specific points in the geometry. Due to the fact that the subregions present finite dimensions, these subregions are called "finite elements", in contrast to the infinitesimal elements used in differential and integral calculus. From here comes the name "Finite Element Method" established by Ray Clough⁽⁵⁾.

Ogden⁽⁶⁾ proposed a model to describe the non-linear behavior of complex materials such as rubber, polymers and biological tissues. This model is described by Equation 1.

(1) $\Psi =\sum _{\mathit{i}=1}^{\mathit{N}}\frac{\mathit{mi}}{\mathit{ai}}\left({\mathit{I}}_1^{\mathit{ai}}+{\mathit{I}}_2^{\mathit{ai}}+{\mathit{I}}_3^{\mathit{ai}}-3\right)$

where Ψ is the deformation energy of the material, N is the total number of terms in the series, and µ_i and α_i are constants of the material.

The deformation energy function can be expanded by a series of real powers, and it is described as a function of the main stretches shown in Equation 2, with the following stability condition, µ_nα_n> 0:

(2) $\mathit{m}=\frac{1}{2}\sum _{\mathit{i}=1}^{\mathit{N}}{\mathit{m}}_{\mathit{n}}{\mathit{a}}_{\mathit{n}}$

where µ is the shear module and α is the deformation hardenability (hardening) capacity parameter. µ is a measure of the ability of a material to withstand transverse deformations (measured in KPa), and is also known as rigidity, which is a parameter to quantify the rigidity of the material. α is representative of the non-linear behavior due to the hardening by deformation of the material.

According to Lapper et al.⁽⁷⁾, new methodologies have now been developed for the simulation of the hyperelastic mechanical behavior of biological tissues. These methodologies generally contemplate the use of the Finite Element Method as a tool for numerical simulations to obtain the numerical behavior of biological materials.

Methods

For the experiments, we sought to combine the analysis capability of the Finite Element Method with the simplicity and robustness of the PSO. In order to obtain the input data in the PSO algorithm, experimental data from the University of Liverpool-UK were used, and it was possible to define the objective function. In these experiments, the specimens (corneas) were fixed at the border, loads were applied with internal pressure of 45 mmHg, in order to experimentally measure the displacements in the central region of the cornea (pupil) by means of digital processing images. The experimental apparatus is shown in figure 1.

Figure 1
Experimental apparatus adapted from Elsheik.⁽⁸⁾ 1- Region of fixation of the cornea; 2: Reservoir for application of pressure; 3 - Controller; 4 - Laser Emitter; 5 - CCD Cameras.

Based on the experimental data for human corneas aged approximately 50 years, the time vs. displacement graph was obtained.

From the experimentally generated curve, the objective function was defined by a polynomial of the 5th degree, according to figure 2.

Figure 2
Graph of the objective function (polynomial function of 5th degree).

The algorithm used in this work describes an "inverse analysis" problem, which will be solved via PSO. The objective of this algorithm was to generate deformation curves close to the experimental data, from maximum and minimum values of x0 (first unknown) and x₁ (second unknown), as input data of the algorithm. In general, the PSO is adequate for numerical problems with two or more unknowns, which is the case of the present work.

The PSO algorithm works on a loop of repetition. At each iteration of the loop, the particle is updated based on the search space and speed. The number of particles is defined by the user in the search for an optimal solution at each iteration. In the PSO algorithm proposed in this work, the user enters with four values, here denominated "limits", that are the maximum and minimum values of the unknowns to be obtained, in this case, minimum of µ, maximum of µ, minimum of α and maximum of α, according to the flowchart shown in figure 3.

Figure 3
Flowchart PSO/FEM⁽⁹⁾

The second step of the algorithm was the generation of a file called material.txt, which was used to substitute random values of µ and α at each iteration / simulation. Thus, it was possible to obtain several simulations and generate a new curve (time vs. displacement) at each iteration. The algorithm was able to compare the objective function (Figure 2) with the curves from the simulations, in order to obtain optimized values for µ and α.

At the end of the total number of iterations (user established value), the particle presents the optimized solution of "µ" and "α". These unknowns represent the mechanical characteristics of human corneas, taking into account hyperelastic material.

The sum of squared errors (SSE), Equation 3, was obtained from the difference between the displacement values from the simulations and the points generated by the curve of objective function (Figure 2):

(3) $\mathrm{SSE}=\sqrt{\sum _1^n\frac{{\left(\mathrm{Pp}-D\right)}^2}{\mathit{n}}}$

where P_p are the displacements from the polynomial equation of the objective function, D is the simulated displacements, and n is the number of points analyzed.

The geometry of the human eye was generated by a program from the University of Liverpool called Ocular Model Generation, but some geometries were modified according to the simulated scenarios. These modifications were made by the Abaqus^® software.

In the present work, the FEM was used to simulate different mechanical properties of the materials and to provide displacement data at specific points in the geometry. For this, the displacements in the central region of the cornea (Figure 4b) were measured as a function of a pressure applied in the inner region of the cornea, taking into account the region of fixation of the cornea in the ocular globe (Figure 4a). The displacement results were used by the PSO algorithm to compare with objective function values at each iteration.

Figure 4
Example of corneal contour data for FEM simulation.⁽¹⁰⁾

To generate the geometry, each element was modeled in the configuration of 6 nodes/element and 15 nodes/element, using a tetrahedral form.

The software Abaqus^® changes the geometry C3D6H (6 nodes) into C3D15H (15 nodes), and in this type of configuration the in the quadratic form must be considered.

The file InputFile.inp contains all the information that determines the structure of the cornea: the number of elements, the number of nodes, pressure, among others. This file works as input data for the software Abaqus^®.

The script Analysis_batch_file.bat is responsible for running the software Abaqus^® without the need to open its graphical interface, generate a file containing the output data with the extension *.odb, and execute the algorithm Out_Node..py (developed in the free language Python^®) which is to generate the simulated displacement curve. The PSO algorithm implements the whole process, working in synergy with Abaqus^®.

The PSO algorithm (developed on the platform Visual Basic^®) runs the file Analysis_batch_file.bat, which runs the Abaqus^® simulation from a file containing the input data withextension *.inp (usually called InputFile.inp) and generates a file containing output data with extension *.odb (file containing Abaqus^® results) and runs the file developed in Python^®, Out_Node.py that uses the file *.odb to generate the time (or pressure) vs. Displacement curve.

In the present work, the simulated curve (time vs. displacement) was generated by means of a file called "displacement.txt", which is a text file that is to be compared with the objective function in the program PSO.

The entire process of creating the three-dimensional geometric models was generated via commercial software Abaqus^® (version 6.14), using its graphical interface to configure and define the parameters used during the simulation. This whole process was based on the setup used by the University of Liverpool for the collection of experimental data and the Ocular Model Generation application. The first step was to import of the file InputFile.inp generated by the program Ocular Model Generation. Although the cornea is an anisotropic material, the respective one was considered isotropic in the simulations for reasons of simplification.

Figure 5 shows the finite element mesh of the eyeball obtained by the software Abaqus^®, comprised of 772 nodes and 768 elements, and with regions representing the entire structure of the human eye.

Figure 5
Finite element mesh of the eyeball.⁽¹⁰⁾

Considering that in order to obtain the corresponding region of the cornea, in Figure 6, the finite element mesh generated through the software Abaqus^® was created with the discretized model of the cornea of a human eye from geometric data from the program Ocular Model Generation

Figure 6
Discretized region of the cornea⁽¹⁰⁾

The cornea model used was considered hyperelastic, being the material defined according to the model OGDEN (Equation 1). The values of µ = 0.054 and α = 110.45 obtained by means of commercial optimization software (HEEDS) were considered reference.

An internal region of the cornea was defined (Figure 7) for the application of internal pressure. In addition, a monitoring point was defined to work as a reference in the process of collecting corneal displacement data.

Figure 7
Internal region of the cornea.⁽¹⁰⁾

After the whole process of configuration of the simulated scenario, a pressure with a magnitude of 45 mmHg was applied to the internal region of the cornea, therefore aiming to calculate the displacement suffered at the monitoring point when it undergoes the respective internal pressure.

For each geometric model generated, a different scenario was assigned based on minimum and maximum values of nodes supported by the version of Abaqus^® used, that is, from the 1st scenario to the 5th scenario (Table 1).

Table 1
Simulations carried out.

Cenários	Nº of Nodes	Nº of Elements	Node/ element
1st scenario	74	54	6 node/element
2nd scenario	291	54	15 node/element
3rd scenario	455	108	15 node/element
4th scenario	753	150	15 node/element
5th scenario	938	864	6 node/element

Each scenario had its values of nodes and elements defined by geometry. Random values of the number of nodes between the maximum and minimum were assigned to the other scenarios.

Results and Discussion

Predefined threshold values

For each scenario analyzed, predefined threshold values for µ₁ and α₁ (Table 2) were established to test the stability of the results due to the variation of the input values in the PSO algorithm.

Table 2
Simulations performed with different thresholds

PARAMETERS	µ 1	α1
	THRESHOLDS	THRESHOLDS
	(MIN-MAX)	(MIN-MAX)
1	0.00001-1	20-400
2	0.00001-1	10-200
3	0.02-0.2	20-150
4	0.01-0.1	60-240
5	0.025-0.075	60-150
6	0.001-0.8	40-300

These parameters were defined in order to test the stability of the algorithm, and the same test was performed several times, obtaining the same result, even in different computers.

Observing the curves shown in figure 8, it is possible to notice that the values begin to stabilize from the 40th iteration with optimized values close to those obtained by the commercial software HEEDS (µ₁ = 0.054).

Figure 8
Comparison of μ₁ values in the second scenario in all parameters.

This is because in the PSO algorithm the values are generated randomly in a search space where each particle, that is, each value corresponds to a possible solution that moves until reaching the best result. The values varied and converged to approximate values to those obtained by the commercial software HEEDS.

Observing the curves shown in figure 9, it is possible to notice that the values begin to stabilize after the 35th iteration,showing results close to those obtained by the software HEEDS (α₁ = 110.45). As with the values of µ₁, it was confirmed that the second scenario was the most appropriate for this situation when compared to the results of the software HEEDS.

Figure 9
Comparison of α₁ values in the second scenario in all parameters

Errors of simulated scenarios

The best result (lowest Error) was found in the 2nd scenario / 5th parameter, with 291 nodes. Figure 10 shows the results of the SSE for the 2nd scenario - 291 nodes, referring to the simulated parameters.

Figure 10
Comparison of values of Error (SSE) in the second scenario in all parameters.

For each scenario analyzed, the results were obtained after the program ran 50 iterations, thus noticing its stability after 40 iterations. By means of this analysis, it can be noticed that the best SSE value was obtained in the 2nd scenario / 5th parameter (approximately 15%), where the scenario used has 291 nodes.

Conclusions

It was verified that the stability and convergence, when varying maximum and minimum limits by the technique of inverse analysis by the Finite Element Method (FEM), along with the PSO algorithm, were satisfactory, because the results obtained demonstrated the potential use of the PSO along with the inverse analysis of the MEF for the characterization of hyperelastic materials.

Values of µ₁ = 0.054 and α₁ = 110.45 were obtained with the commercial software (HEEDS) licensed and used by the University of Liverpool.

For simulated scenarios executed from a restricted number of nodes (maximum of 1000 nodes per geometric model) depending on the version of the Abaqus software used in the simulations.

In relation to the standard PSO x Inertia, the value, even if little relevant, leads to the conclusion that, with experiments, improvements were also noticed in the standard PSO, but based on these parameters these values were not sufficient for a correct evaluation.

This work has succeeded in achieving the objectives, and the continuation of their experiments will be of great importance, implying improvements in both medical sciences and engineering. In medicine, it may contribute to the treatment of various eye diseases, more efficient diagnostic equipment, or even surgical procedures. It can also guide the improvement and creation of mathematical models, acting directly in the improvement of human life.

For future works, it is proposed to use different values of inertia, that is, c₁ ≠ c₂ and test values of µ and α for a larger number of variables, for example, Ogden model for N = 3.

Acnowledgements

To Universidade Federal de Lavras (UFLA) and CAPES for the financial support during the accomplishment of this project.

Referências

• ¹
  Samuelson DA. Ophthalmic anatomy. In: Gelatt KN, editor. Veterinary ophthalmology. Iowa: Blackwell. 2007. vol.1. p. 37-148.
• ²
  Soares RC. Otimização de seções transversais de concreto armado sujeitas à flexão: aplicação a pavimentos [dissertação]. São Carlos: Escola de Engenharia de São Carlos, Universidade de São Paulo; 1997.140p.
• ³
  Kennedy J, Eberhart RC. Particle swarm optimization. In: IEEE International Conference on Neural Networks, 1995, Perth. Proceedings. Perth: IEEE; 1995. vol. 4. p. 1942-8.
• ⁴
  Parsopoulos KE, Vrahatis MN. Particle swarm optimization method in multiobjective problems. In: ACM Symposium on Applied Computing, 2002, Madrid. Proceedings. Madrid: ACM; 2002. p. 603-7.
• ⁵
  Clough RW. The finite finite element method in plane stress analysis. In: Conference on Electronic Computation, 2nd, 1960, Pittsburg. Proceedings. Pittsburg: ASCE; 1960.
• ⁶
  Ogden R. Non-linear elastic deformations. New York: Dover; 1984. 560 p.
• ⁷
  Lapper R, Gasson P, Karri V. A hyperelastic finite-element model of human skin for interactive realtime surgical simulation. IEEE Trans Biomed Eng. 2011;58(4):1013-22.
• ⁸
  Elsheikh A. Understanding corneal biomechanics through experimental assessment and numerical simulation. New York: Nova Science; 2010. 63 p.
• ⁹
  Magalhaes RR, Elsheikh A, Buechler P, Whitford C, Wang J. Application of particle swarm optimization in inverse finite element modeling to determine the mechanical behavior of the cornea. Acta Scientiarum Technology. 2017:39(3):325-331.
• ¹⁰
  Almeida PR. Estudo do comportamento mecânico de córneas via simulações numéricas [dissertação]. Lavras:Universidade Federal de Lavras, Lavras; 2016. 80p

Publication Dates

• Publication in this collection
  Nov-Dec 2017

History

• Received
  13 Dec 2016
• Accepted
  15 Sept 2017