Iterative Methods for the Biomechanical Evaluation of Corneal Response . A Case Study in the Measurement Phase

the number of corneal surgeries steadily grew in recent years and boosted the develop‐ ment of corneal biomechanical models. These models can contribute to simulating surgery by re‐ ducing associated risks and the need for secondary interventions due to ectasias or other problems related to correcting other diseases. Biomechanical models are based on the geometry obtained with corneal topography, which is affected by intraocular pressure and material properties. Knowledge of stress distribution in the measurement phase is a key factor for improving the accuracy of in silico mechanical models. In this work, the results obtained by two different methods: prestress method and displacements method were compared to evaluate the stress and strain distribution in a general geometric model based on the Navarro eye geometry and two real corneal geometries. The results show that both methods are equivalent for the achievement of the stress distribution in the meas‐ urement phase. Stress distribution over the corneal geometry  in the measurement phase  is a key factor for accurate biomechanical simulations, and these simulations could help to develop patient‐ specific models and reduce the number of secondary interventions in clinical practice.


Introduction
In recent years, several efforts were made to study new mechanical models to analyze biological soft tissues [1]. These materials display nonlinear behavior, where the mechanical response depends on time and the previous load state. The diversity of biological soft tissues and the complexity of their structure, made up of several layers and material properties, make mechanical model development difficult.
Of all human tissues, the cornea is the only avascular tissue, and several investigations were carried out due to the number of clinical applications that include refractive surgery corrections, intracorneal rings, evaluation of keratoconus progression analyses, and cross-linking treatment optimization [17]. The in vivo corneal biomechanics study is currently based on biomechanical parameters evaluated from the results of noncontact tonometers [18].
The human cornea protects the eye from external agents, and it is a unique tissue due to its transparency and strength properties. The cornea has an aspherical shape and is usually represented by means of a biconical surface [19]. It is composed of several layers, where the stroma is the main structural layer with near 90% of total thickness [20]. The stroma is based on 200-400 lamellar layers with collagen fibers embedded in an extracellular matrix, these fibers run parallel in each lamella, but at different angles between them.
Lamellar thickness varies across corneal thickness. Lamellae are twice as thick in the posterior zone compared to anterior ones [21], and shear stiffness is in 3-fold greater in the anterior zone due to the percentage of interweaving among lamellae.
Both orientation and distribution were quantified by X-ray techniques [22]. According to the obtained results, collagen fibers are distributed mainly in the nasal-temporal and superior-inferior directions in the central cornea area and circumferentially in the limbus [23]. Daxer, A. and Fratzl, P. [24] estimated that two thirds of fibrils lie in the sectors contained at ±45° around the horizontal and vertical meridians. This distribution can be observed in the central 7 mm zone. According to secondary harmonic generation images, a clear 3D arrangement distribution appears and allows analyzing the mechanical influence of collagen fibers with thickness [25].
Collagen net damage causes loss of stiffness and a change in corneal curvature, which results in the corneal bulging and its local thinning, whose origin is natural or physical. This disease is known as keratoconus. Recent studies [26] established that the deterioration of the transversal bonds between fibers could be the reason for this weakening due to loss of fiber orientation in lamellae. When the origin is due to corrective vision surgery, it is called postsurgical ectasia, and, despite advances, the incidence of the unpredictability of surgical outcomes remains close to 2.8% [27].
The improvement of the surgery success percentage is associated with the development of patient-specific models, for which considerable efforts were made [28]. The development of patient-specific models starts with the topographic corneal measurement and mechanical modelling improvement [21], which strongly depends on the employed material model and materials' properties.
The material properties of the cornea can be evaluated by mechanical tests [17]. Although the material shows a viscous component, the elastic component dominates total behavior, and the material is usually modeled as an anisotropic hyperelastic material [21], with a behavior explained by a volumetric and an isochoric component. Corneal tissue is especially deformable, and tissue can be treated as quasi-incompressible material due to its high-water content [29].
Mechanical behavior depends on geometry, intraocular pressure (IOP), and materials' properties, and the mutual influence of these parameters makes their estimation difficult. Geometry is also affected by the stress state distribution in the measurement phase [15]. Our research team developed specific patient geometric models that allow characterizing the geometric imbalances derived from the distribution of tensions characterizing corneal morphology during the progression of the disease [30][31][32][33]. However, the calculation of stress distribution in the measurement phase is a key factor for improving the accuracy of mechanical corneal models.
The stress distribution in the measurement phase can be calculated following an iterative process based on the prestress level established with a modified updated Lagrangian formulation [34]. The calculated stress distribution in the measurement phase is in equilibrium with the intraocular pressure, and the maximum displacement obtained, when both are applied over the cornea model, is lower than a tolerance previously defined. The free stress configuration can be obtained by an inverse method if only the stress distribution in the measurement phase is considered.
The stress distribution can also be calculated using a simplified displacement method, like that described by Pandolfi, A. and Holzapfel, G.A. [35]. In this method, stress free geometry is assumed, and the geometry obtained when IOP is applied is compared with the corneal topography. If the observed differences are smaller than a previously defined tolerance, the proposed geometry is accepted as the free stress geometry, otherwise a new free stress configuration is proposed based on the correction of the previous one. The stress distribution in the measurement phase is obtained with the proposed free stress geometry and the application of the measured IOP.
Therefore, knowledge of the stress distribution in the measurement phase must be correctly evaluated to better assess the stress distribution over the cornea because the development of accurate patient mechanical models can contribute to improving the accuracy of the corneal refractive surgery procedures. In the current study, a comparison of the results obtained by the two main methods presented for the preliminary stress calculation was made on a theoretical geometry and on two real corneal geometries with different degrees of disease severity.

Materials and Methods
The material model, procedures and main hypothesis carried out for the mechanical analysis of the cornea are presented in the points that follow.

Material Definition
In short-term deformation cycles, material can be treated as an anisotropic hyperelastic model with no viscous loss. The characterization of reversible materials in an isothermal process, with no energy dissipation, is evaluated according to the continuum mechanics with a strain density function. Considering a Cartesian system represented by an orthonormal basis, the position of a material point in the initial and deformed configuration can be denoted respectively by the vectors and . The material deformation gradient is expressed with the two points tensor F = ∂x ∂X where the determinant of this tensor, expressed as J = det(F), considers the volumetric changes of the initial volume. For incompressible materials, J = 1.
The deformation gradient can be expressed as F = J 1 3 F, where J 1 3 represent the volumetric deformation and F the deviatoric deformation. The Right Cauchy Green deformation tensor C=F T F can be expressed by the tensorial product of the transpose of F, denoted as F T , and F and the modified Right Cauchy Green tensor can be expressed as C =F T F considering the deviatoric deformation tensor and its transpose.
Strain energy density was divided into a volumetric, depending on the Jacobian determinant J and isochoric component according to Lanchares, E. et al. [34], as can be observed in Equation (1).
The isochoric component was divided into an isotropic and anisotropic component. The isotropic component depends on the modified Right Cauchy Green tensor C , that considers randomly distributed fibers and the extracellular matrix. The anisotropic term, considers the structural tensors M = m 0 ⨂m 0 N= n 0 ⨂n 0 taking into account the fibers' preferential orientation, defined by the vectors m 0 and n 0 , in accordance with the results observed in the X-ray techniques [24]. According to these results, Nasal-Temporal and Superior inferior directions were considered as preferential orientations. This work did not include density variation across the thickness.
Equation (1) can be evaluated in terms of invariants as seen in Equation (2 Considering the definition of the trace of the modified Right Cauchy-Green deformation tensor tr (C), the invariants can be expressed according to the expressions contained in Equations (3)(4).
The pseudo-invariants, associated with fibers' families, can be expressed in terms of initial orientations and the modified Cauchy-Green deformation tensor according to the expressions contained in Equations (5-10) I ̅ 4 and I ̅ 6 represent the squared stretch of fibers. I ̅ 5 and I ̅ 7 are associated with transversal behavior. I ̅ 8 considers the interaction between two fibers' families. I ̅ 9 is a constant that depends on the material's configuration.
The volumetric component was evaluated according to Equation (11) where the compressibility factor d was calculated considering a Bulk factor K = 5.5 MPa and taken into account that d= 2 K as was considered in previous work [36]: The isochoric isotropic component was modeled with a 2-parametric Mooney-Rivlin, as seen in Equation (12), where a 1 and a 2 are the experimentally obtained materials.
This work did not bear in mind the interaction between fiber families, and the same shear modulus was contemplated in the direction of the fiber families and orthogonal to this direction. According to this hypothesis [1], the influence of pseudo-invariants I ̅ 5 , I ̅ 7 and I ̅ 8 was not considered. I ̅ 9 was constant during material deformation but was not included in the material deformation description. The isochoric anisotropic component was considered by the Holzapfel model [37], in which only pseudo invariants I ̅ 4 and I ̅ 6 were taken into account according to Equation (13).
In this equation, constants k 1 and k 2 are related to fiber behavior. The former considers stiffness at low extension, whereas the latter shows the influence at high strain levels.
According to the hypothesis, and considering the Clausius-Planck inequality, the second Piola-Kirchhoff tensor can be derived according to Equation (14).
From the second Piola-Kirchhoff tensor, the Cauchy Stress Tensor was evaluated by considering the following relation between the tensor components: That applied to Equation (14) can be expressed as the following equation: In these equations, p is the hydrostatic pressure, 1 is the identity tensor, is the modified Left Cauchy strain tensor defined as b=F F T and m and n are the fibers' orientation vectors in the deformed state.
Collagen fibers distribution was considered by the variation in the k 1 parameter. In this paper, corneal surface was divided into octants, and concentric zones with diameters from 0 to 10 mm in steps of 2 mm were considered. Furthermore, a transition zone between 10 and 11 mm was modeled until the limbus zone considered from 11 mm to the external diameter. Fibers were oriented in the vertical and horizontal directions with the exception of the limbus, where fibers were oriented circumferentially. The materials' properties are found in Table 1 where the parameters are obtained from the values obtained by [21,36]. The value of k 1 is varied depending on the zone to consider the fiber density variation [21]. The variation of k 1 is based on the work of Daxer, A. & Fratzl, P [24]. In their work, 66% of fibers are concentrated in horizontal and vertical octants, whereas 34% is concentrated in the oblique octants. The fiber's variation was considered taking into account a reduced value of in the central oblique zones. A transition zone was considered between zone with minimum and maximum values of k 1 .

Finite Element Model Definition
Three finite element models were developed based on the Navarro's theoretical model (N) and two real patient-specific corneal geometries that were measured by a Sirius tomographer (CSO, Italy). The first was a healthy cornea (G0), whereas the second one was affected by keratoconus in an incipient phase (G1) according to the Amsler-Krumeich criteria [32] (provided from Vissum Hospital, Miguel Hernandez University and Keratoconus IBERIA database, Alicante, Spain). This study was approved by the research ethics committee of the Technical University of Cartagena (CEI21_001).
The corneal model was extended to the limbus zone, where it was attached to sclera tissue. An almost 40° cut angle was used according to the average measured values [38]. To implement the finite element model, the ANSYS software (Swanson Analysis System Inc., Pensilvania, USA) was applied by considering the 3D elements of type SOLID 186 with four elements for thickness [21]. The mixed u-P formulation and reduced integration were considered for solving purposes. Figure 1 shows the different defined material zones. Light blue depicts the superiorinferior, nasal-temporal, and central zones. Red denotes the limbus zone and purple the zones in diagonal octants with low fibers' density, where a reduction in the stiffness parameter of fibers was considered according to previous works [24]. Blue zones denote the transition zones between zones with maximum and minimum fibers concentration. Due to the high fibers' density in the limbus zone, where a high stiffness is reached, the displacements were constrained on the exterior face, where the cornea was joined to sclera tissue. This hypothesis avoids the rotation of the limbus zone and was applied in previous works [29].
Intraocular pressure was applied as a pressure load in the cornea's posterior zone, but the aqueous humor influence was not considered. For Navarro's geometry, an IOP of 15 mmHg (1,998.4 Pa) was applied, whereas a respective IOP of 18 mm Hg (2,399.8 Pa) and 16 mm Hg (2,133.2 Pa) was used for the G0 and G1 geometries according to the measured IOP for each patient.

Calculation of the Stress Distribution in the Measurement Phase
The steps of these iterative processes are summarized below:

Displacement Method
This method was performed in an iterative process implemented in Python 3.8 (Python Software Foundation). Ansys files (Swanson Analysis System Inc., Pensilvania, USA) were read and updated following the iteration method found in Figure 2. The initial node coordinates were based on the geometrical values obtained with the topographer. From this geometry, the iterative process was based on the estimation of a free-stress geometry, which was submitted to IOP. The deformed geometry obtained is compared with the corneal topography and if the maximum difference between points is lower than a previously defined tolerance, the proposed geometry is considered the free stress geometry; elsewhere, the iterative process continued until the tolerance criteria was fulfilled or a maximum number of iterations was reached (Figure 2). The geometrical comparison considered the maximum distance between equivalent nodes placed in the intersections between sectors and circumferential divisions in the anterior and posterior corneal faces. The displacement methods allowed the stress-free geometry, stresses, and strains to be computed.

Prestress Method
This method was performed in an iterative process implemented in Python 3.8 (Python Software Foundation). Ansys files (Swanson Analysis System Inc., USA) were read and updated following the iteration method shown in Figure 3. The objective of this method is to calculate a stress distribution in the cornea that equilibrates the applied IOP. From the model based on the geometry obtained with the topographer, IOP was applied, and the deformation gradient of the previous case was applied up to the displacement level was lower than a previously established tolerance. If the displacement level evaluation was higher than tolerance, the process continued until the tolerance criteria was fulfilled or a maximum number of iterations was reached (Figure 3). The prestress method allows the stress distribution in the measurement phase to be computed. However, the stress-free geometry is not directly obtained as in the displacement method. To obtain the free stress geometry, a direct inverse calculation can be applied to the finite element model based on the corneal topography by considering only the stress distribution in the measurement phase without IOP.

Results and Discussion
The stress distributions in the measurement phase with the displacement and prestress method are observed in Figures 4-6 for the three analyzed corneas. In all the studied cases, the stress distribution and maximum values obtained were similar with the two implemented methods. Differences between corneal models, Navarro, G0 and G1, are due to the geometrical and thickness variation and the level of intraocular pressure considered (15 mmHg for Navarro model, 18 mmHg for G0 model, and 16 mm Hg for G1). The theoretical model showed a symmetric distribution that could be evaluated in a quarter finite element model however in patient-specific models' thickness and geometry are not symmetric and the whole model must be taken into account. For the Navarro theoretical model, maximum values are obtained in the posterior limbus zone where displacements are constrained, and fibers are oriented circumferentially ( Figure 4). In this zone, as was observed with X-ray scattering [21] high fibers' density was located. Excluding the limbus zone, the stress distribution in the posterior zone showed a concentric distribution with a stress increment in the center. These results are coherent with the results of other works [21] and could explain the different lamellae's dimensions in the corneal thickness where lamellae are twice as thick in the posterior zone compared to anterior one. In the anterior zone, the results were more homogeneous with lighter stress steps due to the definition of sectors. The obtained von Mises stress level is lower than in the posterior zone and this factor could explain the difference in lamellae thickness, which was thicker in the posterior zone [21].
The values obtained and spatial distribution are similar to those obtained by Pandolfi A. and Holzapfel, G.A. [35]. However, these authors only use theoretical geometries, whereas in the current research real patient-specific geometries were used besides a theoretical geometry.
For patient-specific models ( Figures 5 and 6), the stress distribution obtained with the implemented biomechanical models was coherent with the fibers' orientation and density observed by X-ray scattering [22]. An increase was noted for the fibers oriented in the superior-inferior and nasal-temporal directions in the central zone, which continued in a central zone of approximately 7 mm in diameter [24]. These results are consistent with those obtained on the theoretical geometric model by other authors [28,29,34], however, our models, in addition to the theoretical model, are applied on real cases of healthy and sick patients (G0 and G1, respectively).
The X-ray images confirmed an increment in fiber density in the internal zone of the limbus, which coincided with the zone where the von Mises stress level rose. The maximum values were obtained in the limbus zone for G0 configuration, but these values were affected by the intraocular pressure (IOP) considered in the study and the constraints imposed in these zones when displacements were constrained. These restrictions are common for this type of study [28,29,34]. An average maximum value was obtained in the central zone for the three analyzed geometries, as shown in Table 2. These values are affected by the IOP in the measurement phase (patient-specific). By means of a linear interpolation and by taking a reference intraocular pressure of 15 mmHg, the maximum values of the von Mises stress distribution in the measurement phase were obtained. As can be observed in Table 2, when the same IOP was considered, the maximum von Mises stress level rose with keratoconus evolution, and both values had higher values than those obtained with Navarro eye's geometry. These results are alienated by those presented by other authors [13,28,29], however, to our knowledge, this is the first study that compare at the same time the methods of displacements and prestress on real geometric models of patients, analyzing the influence of the stress field of both methods on healthy patients and in an early stage of the disease. Finally, and for validation purpose, the intraocular pressure (IOP) was directly applied to the free stress geometry of the theoretical Navarro model and patient-specific models, the results obtained showed that this IOP distorts the deformation field in the tissue of the corneal structure as shown in Figures 7-9. On the models, a tangent elastic module was derived in the central zone from the stress and strain values. An average value of 0.34 MPa was obtained for Navarro eye's geometry (15 mmHg), a value of 0.37 MPa for the G0 geometry (18 mmHg), and 0.37 MPa for the G1 geometry (16 mmHg).  The tangent elastic modulus obtained in the central zone of the modelled corneas is in the range of values reported by other authors in inflation tests [39]; therefore, our models are consistent with the research developed.
The values obtained are limited to the material model chosen and the parameters implemented, that depends on patient-specific tests. This factor does not affect to the results of this research: The two implemented methods, displacement method and prestress method, were equivalent, and the same stress distribution in the measurement phase was obtained. The displacement method allowed the stress-free geometry to be directly recovered and it proved to be a less complex method from the computational point of view.
Regarding the clinical application of simulation models, the main challenges are defined by a realistic simulation in real time [40,41]. In this sense, our research team developed specific real-time patient 3D geometric models [30][31][32][33]. The accuracy of these geometric models allows for a good quality meshing, which leads to the efficient generation of specific models used in preoperative analysis or for surgical planning.

Conclusions
In this research, the stress distribution in the measurement phase was calculated by two iterative methods: prestress method and displacement method. The distribution and maximum values reached were similar, and differences were attributed to the numerical methods.
The stress level distribution is coherent with the results of previous works. Maximum values are obtained in the limbus zone, where according to X-ray scattering, a high circumferential fiber density is observed. The stress level in the posterior central zone could explain the different lamellae's dimensions in the corneal thickness where lamellae are twice as thick in the posterior zone compared to that of the anterior one.
Although both methods allowed the stress distribution to be known in the measurement phase, the displacements method gave the stress-free geometry of the cornea directly, and the field of strains was obtained when real-patient IOP was applied. For the authors' knowledge, it is the first time that both methods were compared with real geometries of healthy and sick patients.
The knowledge about the stress distribution in the measurement phase is a key factor for accurate mechanical simulations. These simulations can help to develop patient-specific models and to reduce the number of second interventions in clinical practice.
For future research development, more complex material models and other algorithms to recover the free stress configuration can be analyzed to obtain faster and more accurate solutions to assist surgeons.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Conflicts of Interest:
The authors declare no conflict of interest.