Modeling basics

The structural model of the cornea can be approached by creating a continuum model or shell model. Continuum models explicitly model the geometry of the cornea in three dimensions, although two-dimensional (2D) approximations are sometimes used. Shell-based models use the theory of thin shell behavior to decrease the geometrical complexity. The modeled geometry is a surface, with the thickness of the shell as a defined parameter. Tensile and bending loads on the shell create deformations in the shell in accordance with shell theory. While this approach decreases geometrical complexity and increases efficiency of the formulation, it is difficult to extend to complex constitutive models. Anderson et al.⁽²³⁾ proposed a corneal model by using shell analysis to develop a nonlinear finite element model to study the mechanical behavior of the cornea. Howland et al.⁽²⁴⁾, Pinsky and Datye⁽²⁵⁾, Li and Tighe⁽²⁶⁾, and Elsheikh and Wang⁽²⁷⁾ also based their models on this method.

Geometry

The structure of the human cornea has a nonuniform curvature with variable thickness throughout⁽²⁷⁾. It is thinner at the center and thickens towards the edges. Dubbelman et al.⁽²⁸⁾ recorded the corneal geometry for the internal and external surfaces of the human cornea using Scheimpflug photography. Roy et al.⁽¹⁴⁾ used magnetic resonance imaging to determine the geometric profile of the eye.

Some early models of the cornea, e.g.⁽²⁹⁾, use a 2D axisymmetric approximation of the cornea. This models the geometry as a surface of revolution and is efficient compared with three-dimensional models (3D). However, it cannot capture the subtle variations in the geometry of the cornea from a surface of revolution, prevalent anisotropies in the material, or any loads that occur off-center, such as an impact or incision. Figure 1 illustrates common shell elements and 3D elements used for FEA of the cornea.

Figure 1 Common shell and solid elements used for finite element analysis (FEA) of the cornea. A) Triangle and B) quadrilateral are common shell elements. Each shell element node has 5 degrees of freedom (DOF): 3 translational and 2 rotations (approximately 2 orthogonal axes on the shell surface). C) Tetrahedral and D) hexahedral are common 3D elements. Each element node has 3 translational DOF. 

Pandolfi and Holzapfel⁽²⁾ used the corneal geometry data published by Dubbelman et al.⁽²⁸⁾ to create biconic interior and exterior surfaces for a 3D corneal model. Pandolfi and Manganiello⁽¹⁰⁾ used an ellipsoidal shape to define the anterior and posterior corneal surfaces. Nguyen et al.⁽³⁰⁾ and Nguyen and Boyce⁽³¹⁾ also approximated the corneal surface to be ellipsoidal according to digital image correlation (DIC) measurements of the corneal surface profile. Buzard⁽³²⁾ and Howland et al.⁽²⁴⁾ used axisymmetry to create their corneal models.

Salimi et al.⁽³³⁾ meshed the corneal model into triangular shell elements. For 3D models, hexahedral elements are more commonly used^(2,6), although tetrahedral elements may also be used. Linear tetrahedrals and standard trilinear hexahedrals are known to perform poorly for nearly incompressible materials such as soft tissues. Trilinear hexahedral elements can be modified using techniques such as selective reduced integration⁽³⁴⁾; alternatively, higher order elements may be employed. Figure 2 shows an example of an FE mesh of the cornea.

Figure 2 An example of a finite element (FE) mesh of the cornea using hexahedral elements. 

Boundary conditions

The viscous fluid filled inside the anterior chamber of the eye exerts pressure known as IOP on the cornea. Under physiological conditions, the IOP inflates the cornea and gives it shape. A model developed for the cornea using the entire eyeball allows appropriate realistic displacement at the limbus. Uchio et al.⁽³⁵⁾ and Amini and Barocas⁽³⁶⁾ modeled the entire eye globe for the analysis. However, the whole eye model may not be efficient because it is highly time consuming and not economical for development and analysis.

Alastrue et al.⁽⁵⁾ analyzed the corneal model with the limbus constrained against displacement. Anderson et al.⁽²³⁾ compared the deformation pattern of the entire eye model to the rigid cornea-limbus boundary. They proposed an approximate boundary condition, with roller support at the edges inclined at 40° with respect to the horizontal axis, to represent the cornea-limbus behavior. The results obtained from their boundary assumption were similar to the results obtained for the whole eye. Other studies^(11,37,38) adopted similar boundary conditions for their corneal models. Roy et al.⁽¹⁴⁾ created the whole eye model and compared the displacement results with those for the corneal model with fixed sclera to show that the corneoscleral limbus plays an important part in predicting the response to refractive surgery. Pandolfi and Holzapfel⁽²⁾ proposed a corneal model allowing rotation at the limbus and restricting displacement at the edges. Rotation at the limbus plays a role in changing the curvature of the cornea. Figure 3 illustrates these boundary conditions.

Figure 3 Schematic of boundary conditions at the limbus used by researchers. A) The cornea is constrained against displacement and rotation, as proposed in⁽⁵⁾. B) The cornea is fixed against displacement, but may rotate at the edges, as in⁽²⁾. C) Springs normal to the cross-section. Researchers in⁽¹³⁾ investigated boundary conditions (a), (b), and (c). D) The cornea is allowed to expand along a given surface, as in^{(11, 23, 37, 38)}. Forty degrees from horizontal is taken as a typical angle. 

Figure 4 illustrates displacement results of FE simulation of a cornea subjected to IOP, with rotation at the limbus.

Figure 4 An example of vertical displacement mapping of a mouse cornea subjected to intraocular pressure (IOP) obtained from finite element (FE) simulations is shown in the cross-section. The wireframe on the bottom is the undeformed mesh, and vertical displacement colored by magnitude is shown on top. The corneal displacement, but not rotation, is restricted along the limbus. 

Material models

For creating a material model, the corneal material is often approximated as the stroma, which is the major contributor to mechanical strength^(2,6,31). The stromal layer may be modeled as a nonlinear, anisotropic, viscoelastic material that undergoes large deformation. Finite deformation theory is incorporated in the corneal model, which successfully captures its true nonlinear response under deformation^(2,31).

Earlier finite element models were simple linear elastic models. One of the earlier models developed by Vito et al.⁽³⁹⁾ considered the stroma as a linear, elastic, homogeneous, and isotropic material. Similarly, Bryant et al.⁽⁴⁰⁾, Hanna et al.⁽⁴¹⁾, Gefen et al.⁽¹²⁾, and Velinsky and Bryant⁽⁴²⁾ also assumed that the cornea would undergo a small deformation and used linear elasticity for their analysis. Pinsky and Datye⁽²⁵⁾ developed a linear material corneal model based on the anisotropic constitutive model using the predominant fiber directions in the cornea.

The first geometrically nonlinear models represented the stroma as a nearly incompressible, isotropic, hyperelastic material. The Neo-Hookean model is widely used for modeling tissues^(2,6). Other hyperelastic models have also been used. Niroomandi et al.⁽⁴³⁾ used a simpler Saint Venant-Kirchhoff hyperelastic model. Bryant and McDonnell⁽⁴⁴⁾ created a 2D axisymmetric model to compare various isotropic and a couple of transversely isotropic constitutive models of the cornea. Of the models they compared, an exponential nonlinear model best fit the inflation experiments.

Many researchers used a base isotropic model for the matrix of the stroma and additional terms for the effects of the fibers, as discussed below. Alastrue et al.⁽⁵⁾, Hanna et al.⁽⁴¹⁾, and Pandolfi et al.⁽¹³⁾ considered the matrix as a Mooney-Rivlin model. Nguyen et al.⁽³⁰⁾ and Pinsky et al.⁽⁶⁾ used a Neo-Hookean model as their matrix model.

Anisotropic nonlinear hyperelastic models with embedded collagen fibers have been frequently used in recent corneal models^(2,5,6,31). Alastrue et al.⁽⁵⁾ used a discrete fiber model with embedded collagen fibers in two preferred orientations, the nasal-temporal and superior-inferior directions. Most other recent models do not simulate each fiber explicitly, which can be computationally expensive; however, they include a smeared effect of oriented fibers. Cristobal et al.⁽⁴⁵⁾ and Pandolfi and Manganiello⁽¹⁰⁾ developed their models with two preferred directions at the center and circumferential at the limbus region. Pandolfi and Holzapfel⁽²⁾ used a distributed model with two preferred orientations. Pinsky et al.⁽⁶⁾ and Nguyen et al.⁽³⁰⁾ developed a continuous distributed fiber model with preferred orientation.

Nguyen et al.⁽³⁰⁾ created an anisotropic corneal model, taking into consideration the viscoelastic properties. Yoo et al.⁽⁴⁶⁾ also used a viscoelastic ocular tissue model for their study.

Parameter measurement and optimization

Corneal tissue material properties are determined by various laboratory tests of corneal specimens. The material parameters used for the cornea model are based on these studies. Wollensak et al.⁽⁴⁷⁾, Hoeltzel et al.⁽²⁹⁾, and Zeng et al.⁽⁴⁸⁾ performed uniaxial tension tests on corneal tissue to record stress-strain relation. Kohlhaas et al.⁽⁴⁹⁾ and Wollensak et al.⁽⁴⁷⁾ extended these tests to study the effects of CXL on the mechanical response of the cornea. Bryant et al.⁽⁴⁴⁾, Elsheikh et al.⁽⁵⁰⁾, and Anderson et al.⁽²³⁾ conducted inflation tests on the entire cornea to determine the apical displacement at different pressure values. Petsche et al.⁽⁵¹⁾ examined the depth-dependent material properties of the stroma. However, further experimental work to completely determine the mechanical behavior of the cornea, particularly in vivo, is required.

Inverse simulation is a mathematical approach to determine unknown model parameters to match observed or assumed physical responses. In inverse FEM, an optimization algorithm is coupled with an FEM to find a set of optimal parameters for a given model to be used in the FE simulation. Such methods are important for determining the physical properties of actual corneas.

Nguyen and Boyce⁽³¹⁾ presented an inverse FEM to estimate bovine corneal material properties using an in vitro inflation experiment and investigate the influence of variations in these properties on the bovine inflation response. They determined a set of anisotropic material properties, minimizing the error between simulation results and experimental measurements of the surface displacement field.

Summary of finite element modeling choices

Different choices in finite element models of the cornea can change the results. 2D or shell models are efficient, but they may not be able to accurately capture the complex material behavior of some applications. While modeling the entire eye may avoid approximate boundary conditions at the limbus, the computational time increases significantly. However, modeling the changes in IOP during an impact may be difficult without a full model of the eye.

Material models are also increasingly accurate and, at the same time, computationally more expensive. While linear models are very efficient, they may be inaccurate over large changes in deformation. Isotropic models may not accurately capture the behavior of the complex material in the stroma, but they may be adequate for problems that only require a rough idea of the displacements. The depth of knowledge of the properties of corneal tissue is improving. In applications with rapid load, the viscoelastic and dynamic properties of the cornea may be important. For other analyses, such as those of the long-term results of surgery, a quasistatic model is acceptable. It is important to select a model that can simulate an application with the necessary accuracy. However, a model with too much detail may computationally become very expensive.

Surgery

The curvature of the internal and external surfaces of the cornea significantly affect its refractive power. Refractive eye surgery can improve the visual acuity in patients with common refractive maladies such as astigmatism, hyperopia, and myopia. In several surgical techniques, the dioptric power of the cornea is improved by removing thin layers of biological tissues and adjusting the curvature of the cornea using excimer lasers⁽¹³⁾. Laser-assisted in-situ keratomileusis (LASIK), photorefractive keratectomy (PRK), and intrastromal photorefractive keratectomy (ISPRK) are the 3 main types of refractive surgeries. Refractive surgeries on the periphery of the cornea, such as conductive keratoplasty (CK) and intracorneal rings (ICR), have recently attracted more interest. These surgeries are primarily performed on hyperopic eyes that require steepening rather than flattening of the corneal surface⁽³⁸⁾.

FEA has potential usefulness as a simulation tool for refractive surgery planning using accurate in vivo corneal material properties and geometry data⁽¹⁴⁾. FEA has been extensively used to model deformation of the cornea for corneal surgery simulation and predict mechanical and refractive effects of surgery⁽²⁴⁾.

Pandolfi et al.⁽¹³⁾ implemented an FEM in astigmatic and myopic corneas to evaluate the mechanical and refractive outcomes of laser refractive surgery. They computed the change in the curvature and thickness of the human cornea exposed to laser ablation and evaluated the corrected dioptric power in individual patients.

Roy et al.⁽¹⁴⁾ proposed an FE model of the whole eye to give insight into refractive surgery planning. They later presented an inverse FEM to determine the undeformed state of the cornea⁽¹⁶⁾. They also developed a 3D, patient-specific, corneoscleral FE model to estimate the surgical impact on corneal shape variations and changes in corneal elastic properties.

Niroomandi et al.⁽⁴³⁾ presented a novel numerical technique for actual surgery simulation using an extended FEM (X-FEM). The X-FEM enhances standard FE by defining additional degrees of freedom, in this case, for fracture opening (those interested in further theory of the X-FEM can refer to the work of Belytschko et al.⁽⁵²⁾, among others). In their model, fracture opening induced by the incision placed during refractive surgery was reproduced using the X-FEM method. Their presented FEM was able to simulate the complex constitutive model of the cornea and reproduce the discontinuities of the scalpel incision for actual surgery simulations.

Nomograms in ophthalmology are sets of tables of corrected values or graphics used to plan the surgical procedure, particularly incisional surgeries for astigmatism and similar diseases⁽⁴⁵⁾. Nomograms can manage the incision effect by revising the cut parameters such as depth, length, and optical zone in incisional surgery. Cristóbal et al.⁽⁴⁵⁾ proposed an FE model of a limbal incision to predict patient-specific optical power that can be induced during surgery. The nomogram outcomes were compared to simulation results, and, in some cases, the nomograms were revised accordingly.

CK is a subablative thermal treatment wherein thin electrodes are inserted into the stroma to apply radiofrequency (RF) heating. The heat causes the peripheral corneal tissue to shrink. Jo et al.⁽¹⁹⁾ developed a 3D FE model, including the cornea, aqueous humor, and RF electrodes, to investigate the resultant thermal damage during RF heating in the CK procedure. Fraldi et al.⁽³⁸⁾ also presented a viscoelastic FEA to study the mechanical response of the hyperopic cornea to CK and assess the postsurgical stability of the applied refractive correction. Hameed et al.⁽⁵³⁾ created a 3D FE model to simulate the ISPRK technique for correcting myopic corneas and predicting treatment outcomes. The simulation outcomes were compared with the results from an earlier 2D FE model developed by Bryant et al.⁽⁴⁰⁾ for the same purpose. The model outputs were more strongly correlated with the clinical data when the 3D model was used than when the 2-D model was used.

Corneal diseases

Finite element analysis has been widely used to study biomechanical interactions in corneal maladies (e.g. keratoconus) and better understand the etiology and treatment of the diseases. FEM can be used to investigate the effects of variations in corneal material properties on vision power in patients with corneal diseases.

Keratoconus is a degenerative corneal disease characterized by irregular thinning and bulging of the corneal structure and progressive topographic irregularities in the cornea. This shape distortion leads to an optical aberration that can be corrected with the use of hard contact lenses or glasses in mild cases or corneal transplantation in more severe cases⁽¹²⁾. CXL is a new method for treating keratoconus. In this technique, the epithelium is usually removed from the central zone of the cornea, and riboflavin is absorbed into the stroma⁽¹⁶⁾. Under ultraviolet (UV) radiation, the riboflavin creates additional cross-linking in the collagen fibrils, stiffening the cornea.

Gefen et al.⁽¹²⁾ presented a 3D linear anisotropic FE model of normal and keratoconic corneas to investigate the mechanical behavior and optical performance. The results of analysis revealed that IOP had a considerable influence on the optical power of the keratoconic cornea, while having little effect on the refractive power of the normal cornea. Roy et al.⁽¹⁶⁾ presented a 3D patient-specific FE model of the whole eye to study keratoconus progression. The model was able to demonstrate a patient-specific procedure for investigating the locally decreased corneal elastic properties during keratoconus development.

Carvalho et al.⁽¹¹⁾ presented an FEM of a keratoconic cornea to predict the biomechanical behavior and evolution of the cornea in keratoconus. Using shell elements, they investigated how variations in IOP and material properties of the cornea (i.e., decrease in rigidity of the cornea) can lead to a localized increase in corneal curvature. Foster et al.⁽⁵⁴⁾ examined the stiffening effect of CXL in corneal tissue, fitting the material parameters to experiments on corneal strips.

Impact and trauma

Ocular trauma can cause long-term vision disorders and can be expensive to treat⁽⁵⁵⁾. FEM is a useful tool for simulating ocular trauma. It may offer possible solutions for decreasing eye impact injuries and designing protective tools against globe trauma. The impact of different projectiles on the eye has been simulated using FEA to investigate the injury potential.

Stitzel et al.⁽⁵⁵⁾ offered a nonlinear FEM of the human eye to predict eye rupture injuries caused by high-speed blunt impact. The numerical model was verified with experiments to predict the injury due to different types of impacting projectiles and loading conditions.

Weaver et al.⁽⁵⁶⁾ implemented a numerical model of eye impact that considered a variety of blunt projectiles and loading conditions and was based on many experimental impact tests. They investigated the effects of different projectile mass, size, material properties, and velocity on the response of the eye to impact. The study results gave insights for predicting eye rupture in various loading conditions and designing eye safety equipment.

Gray et al.⁽⁵⁷⁾ studied the impact of paintballs on the cornea, with the aim of designing safer paintballs. To this end, they implemented a numerical model of the human eye, orbit, and paintball into the CTH computer code (a 3D, finite-volume, and large deformation numerical hydrocode developed by Sandia National Laboratory) to study the physical mechanisms causing impact eye injuries. By varying material properties, they concluded that paintballs with less mass can be safer for the eye.

Uchio et al.⁽³⁵⁾ also presented a 3D FEA supercomputer simulation of traumatic impacts on the eyeball. They investigated the threshold of impact eye injuries caused by bodies of different sizes and velocities.

Other applications

FEM can be used as a potential noninvasive tool in clinical applications for investigating the biomechanical behavior of the cornea and obtaining characteristic measurements of the tissue. FEM may be a critical method for assessing corneal mechanical and electrical properties and for obtaining IOP measurements.

Applanation tonometry is among the most widely used contact tonometry techniques for IOP measurement. FEM has been used to examine applanation tonometry for measuring IOP⁽³³⁾. Ghaboussi et al.⁽⁵⁸⁾ developed a computationally efficient numerical method to accurately measure IOP in the cornea using a modified applanation tonometer. They used the combination of a neural network and genetic algorithm to fit clinical applanation tonometry outcomes (force and displacement history data) to the outcomes achieved from the FE model.

As the heart beats, IOP behind the cornea changes; the vibration characteristics of the eye also vary with IOP. Salimi et al.⁽³³⁾ developed a coupled fluid-structure FEM of the eye and investigated the dynamic response of the eye to changing IOPs. They validated the vibrational characteristics from the simulation model with experimental modal analysis (EMA) of a water-filled spherical shell.

Rhee et al.⁽⁵⁹⁾ developed a nonlinear, anisotropic FE model of the mouse cornea to study the formation of spiral patterns on the epithelium. They compared the obtained numerical curves of maximum shear strain with actual spiral patterns observed on X-ray images of mouse cornea. The presented work helped in providing insight into the possible development of the cornea.