Establishing Standardization Guidelines For Finite-Element Optomechanical Simulations of Refractive Laser Surgeries: An Application to Photorefractive Keratectomy

Abstract

Purpose: Computational models can help clinicians plan surgeries by accounting for factors such as mechanical imbalances or testing different surgical techniques beforehand. Different levels of modeling complexity are found in the literature, and it is still not clear what aspects should be included to obtain accurate results in finite-element (FE) corneal models. This work presents a methodology to narrow down minimal requirements of modeling features to report clinical data for a refractive intervention such as PRK.

Methods: A pipeline to create FE models of a refractive surgery is presented: It tests different geometries, boundary conditions, loading, and mesh size on the optomechanical simulation output. The mechanical model for the corneal tissue accounts for the collagen fiber distribution in human corneas. Both mechanical and optical outcome are analyzed for the different models. Finally, the methodology is applied to five patient-specific models to ensure accuracy.

Results: To simulate the postsurgical corneal optomechanics, our results suggest that the most precise outcome is obtained with patient-specific models with a 100 µm mesh size, sliding boundary condition at the limbus, and intraocular pressure enforced as a distributed load.

Conclusions: A methodology for laser surgery simulation has been developed that is able to reproduce the optical target of the laser intervention while also analyzing the mechanical outcome.

Translational Relevance: The lack of standardization in modeling refractive interventions leads to different simulation strategies, making difficult to compare them against other publications. This work establishes the standardization guidelines to be followed when performing optomechanical simulations of refractive interventions.

Introduction

In the last two decades, corneal laser surgery has gained more and more popularity because it allows one to achieve spectacle independence by permanently correcting medium-low refractive defects. Among all the existent procedures, photorefractive keratectomy (PRK) was the first laser surgery to be used in clinics.¹ In PRK, a portion of corneal tissue is ablated from the anterior surface to change its curvature and, in turn, the eye's refractive power. 

With the advent of laser in situ keratomileusis (LASIK), the number of PRK procedures has decreased because of higher postsurgical pain and longer healing time. Still, PRK remains the preferred procedure in some cases, such as for patients with thin and irregular corneas, for retreatment because of residual refractive error after LASIK, and for soldiers and professional athletes who may be subjected to higher risk of postsurgical flap dislocation if LASIK is performed.² 

Because the eye is a biomechanical system in which corneal stresses are in equilibrium with the intraocular pressure (IOP),^3,4 removing a portion of corneal tissue alters this equilibrium, which can sometimes bring severe complications such as postsurgical corneal ectasia.^5,6 Corneal ectasia leads to a progressive thinning of the corneal thickness, caused by the growth of corneal tissue derived from the loss of the mechanical stiffness of the stromal tissue,⁷ because of fibers’ disarrangement.⁸ To prevent and avoid postsurgical complications, clinicians perform a presurgical evaluation to determine whether the patient is eligible to receive laser refractive surgery. Often in the past, presurgical screening only included corneal topographers (Pentacam; Oculus Optikgeräte GmbH, Wetzlar, Germany; and Sirius; CSO, Scandicci, Firenze, Italy). 

Corneal topographers measure the elevation of the corneal surfaces and its thickness, providing information on the corneal regularity and optical power and gaining insight into the overall refractive correction that the patient would require. Patients with thin or irregular corneal thickness are considered non-candidates and cannot receive laser surgery, because the risk of postsurgical complications increases.⁷ 

With the advent of new technologies, non-contact tonometers (e.g., Corvis ST; Oculus Optikgeräte GmbH) have been used to estimate corneal biomechanics indirectly. Especially in the case of Corvis ST, a measure of IOP is also obtained. In reality, these devices cannot explicitly characterize the mechanical properties of the corneal tissue, but they measure a global stiffness factor that accounts for different parameters such as corneal stiffness, IOP, and corneal thickness. In this test, the corneal deformation caused by the action of an air-puff is recorded with a high-speed camera to determine markers such as the maximum deformation amplitude or the applanation times, length, and speed. Nowadays, noncontact tonometers estimate the patient's corneal mechanical response, allowing the clinician to decide whether it is safe to proceed with the surgery.⁵ 

Proper presurgical assessment that includes patient-specific geometry and mechanics cannot be directly accessed in clinics, which would help determine the mechanical effects induced by laser surgery onto a patient's corneas. Computational models have been used to tackle this problem and to provide surgeons with an ancillary tool that helps in planning refractive surgeries. Briefly, finite-element (FE) models use a calibrated presurgical mechanical model to simulate a refractive intervention that estimates the postsurgical mechanical equilibrium. 

In 2006, Deenadayalu et al.⁹ developed the first FE model of laser surgery by creating a corneal-scleral shell model. Later that year, Alastrué et al.¹⁰ simulated a PRK surgery using an average conic cornea and an ablation profile that followed the Munnerlyn equation.¹¹ In 2009, Pandolfi et al.¹² analyzed the presurgical and postsurgical mechanical behavior of PRK using biconic models. Meanwhile, Roy et al.¹³ introduced optical concepts such as tangential curvature maps to analyze refractive simulations in axisymmetric models that disregarded the stress-free configuration of the eyeball, leading to an overestimation of the stress and stretch. In 2010, Lanchares et al.¹⁴ used rotationally symmetric corneal models that included the stress-free configuration of the eyeball to investigate the effect of IOP on the outcomes of PRK surgery. In 2011, Roy et al.¹⁵ introduced a methodology to build patient-specific (PS) corneal models to analyze the effect of different ablation profiles. In this work, optical metrics such as sagittal curvature maps and spherical aberrations¹⁶ were used to compare presurgical and postsurgical outcomes. Following up on this research line, Roy et al.¹⁷ compared the effect of two additional surgical techniques, laser in situ keratomileusis (LASIK) and small incision lenticule extraction (SMILE), by simulating the postsurgical outcomes and considering a depth-dependent anisotropic material model. From 2015 onward, research started to stress the importance of including optics as well,^18–25 when the aim is to obtain a representative computational model of the cornea. Still, research mostly focused on the mechanical effect of laser surgery on the corneal structure²⁶ rather than on the compound optomechanical effect. 

When the main goal is to perform a presurgical assessment of the effect of the laser on the corneal structure, both mechanics and optics must be considered during the modeling process. It has been hypothesized that the most accurate way to simulate and reproduce refractive interventions is by exploiting patient-specific models. However, there must be a clear pathway to performing these simulations to achieve reproducible and repeatable results. The present study aims to determine the minimal ingredients needed to model a refractive intervention such as PRK adequately. Likewise, we present a set of modeling practices and standards that can help ensure comparability for future studies. In particular, we present an optomechanical analysis of what existing geometrical models (conic, biconic, and patient-specific models) can be used when modeling the cornea. The impact of different boundary and loading conditions on both the postsurgical optics and stress are analyzed. Therefore this study focuses on determining the sufficient modeling complexity for obtaining an accurate laser refractive surgery simulation, applied to the specific case of PRK surgery. Finally, the presented methodology is applied to the eyes of five patients who previously underwent PRK surgery to ensure accuracy of the modeling procedure. 

Material and Methods

Overall, the FE pipeline to simulate a PRK surgery comprises six different steps thoroughly described in the following sections (Fig. 1). The first step of the process requires the acquisition of patient's corneal topographic data, which consist of a point cloud containing three-dimensional Cartesian coordinates (x-y-z; in mm). Some widely diffused corneal topographers (e.g., Pentacam) are not able to acquire the entire corneal surfaces, thus a surface reconstruction is needed to build the model. After data acquisition, the anterior and posterior corneal surfaces are reconstructed using the previously acquired topographic data (Step 2). Depending on the desired modeling complexity, different geometries that can range from a simple conic up to a PS model can be built by using patient-specific data. After surface reconstruction, the surgical ablation profile, which is based on conic, biconic, or wavefront-optimized profiles (Step 3), must be included. After that, the model is meshed, and the material properties are assigned to the different parts of the model. The simulation procedure includes the computation of the stress-free configuration (Step 4) and the refractive intervention (Step 5). Finally, an optomechanical analysis is conducted on the surgery outcome (Step 6) and compared to original presurgical and postsurgical patient's keratometric data. The automated tool was developed using a combination of MATLAB R2022b for surface reconstruction, ANSA pre-processor by BETA-CAE v22.0.1 for meshing, FORTRAN language for the material description of the model, and ABAQUS 6.13 for solving the FE problem. All blocks were assembled using MATLAB, constituting the entire process's launcher. 

Figure 1.

 

General pipeline of the tool for simulation of PRK surgeries.

Corneal Surface Reconstruction

The cornea's refractive power is mainly related to its shape; thus, using PS data is considered mandatory when the numerical simulation has to replicate the refractive changes induced by the surgery because of the new geometrical configuration. Although the PS model is the most accurate for reproducing the optical properties of the cornea,²¹ simpler analytical models that quickly describe significant refractive parameters (sphere and cylinder related to myopic and astigmatic defects, respectively) might be sufficient to represent overall changes after refractive surgery and are of direct interpretation by the target audience: optometrists and ophthalmologists. 

Topographic data of one PRK patient were retrieved using Pentacam (Fig. 2). 

Figure 2.

 

Pentacam topography used to build the models.

The patient signed a written informed consent, and data were anonymized by our clinical partner, Barraquer Ophtalmologic Center (Barcelona). Data acquisition followed good clinical practices and adhered to the tenets of the Declaration of Helsinki. This study was approved by the Comité de Ética de la Investigación con medicamentos (CEIm) of Centro de Oftalmología Barraquer, Spain (Ceim code: 207_Modelización_Qx_Rx). 

Three different models with growing complexity (conic, biconic, and PS models) were tested to evaluate the capacity to provide accurate mechanical and optical results in the surgery simulation. The retrieved data were used to create the patient-specific models and to approximate the conic and biconic models (Fig. 2). The parameters for the conic and biconic models retrieved from the topography are presented in Table 1. For the conic approximation, the average radius R_m of the steepest and flattest curvature radii (R_s and R_f) were used to reconstruct the surface. For the biconic approximation, the radii along the steepest and flattest curvature meridians (K₂ and K₁, respectively, with the relative radii R_s and R_f) and the asphericities (Q_s and Q_f) along both meridians were used to reconstruct the surface. 

Table 1.

 

Surface Parameters for Both Conic and Biconic Models

 

Conic Model

It is the simplest model used to represent corneal shape.²⁷ The apex is placed at the origin of the reference system, and it is characterized by rotational symmetry with respect to the optical axis (z-axis). The equation of the conicoid is the following:  

\begin{eqnarray} {x^2} + {y^2} + \left( {1 + Q} \right){z^2} - 2z{R_m} = 0\quad \end{eqnarray}

(1)

where R_m is the average corneal apical radius and Q is the asphericity parameter. This representation can be used to reproduce myopia and hyperopia. Still, it is not suitable for astigmatism, given its rotational symmetry, which does not allow to account for different curvatures on different meridians, constituting the highest limitation for this kind of model. 

Biconic Model

The biconic surface is defined using the following equation in cylindrical coordinates²⁸:  

\begin{eqnarray} z\left( {\rho ,\theta ,{R_s},{R_f},{Q_s},{Q_f},{\theta _f},{z_0}} \right) = {z_0} - \frac{{{\rho ^2}A}}{{1 + \sqrt {1 - {\rho ^2}B} }}\quad \end{eqnarray}

(2)

with  

\begin{eqnarray} A = \frac{{co{s^2}\left( {\theta - {\theta _f}} \right)}}{{{R_f}}} + \frac{{si{n^2}\left( {\theta - {\theta _f}} \right)}}{{{R_s}}}\quad \end{eqnarray}

(3)

and  

\begin{eqnarray} B = \left( {{Q_f} + 1} \right)\frac{{co{s^2}\left( {\theta - {\theta _f}} \right)}}{{R_f^2}} + \left( {{Q_s} + 1} \right)\frac{{si{n^2}\left( {\theta - {\theta _f}} \right)}}{{R_s^2}}\quad \end{eqnarray}

(4)

where R_s and R_f are the radii of curvature of the steepest and the flattest meridians, respectively, Q_s and Q_f are the corresponding asphericities, and θ_f represents the axis of astigmatism (Table 1). The biconic surface reaches the maximum value z₀ at ρ = 0. If the same values of radius and asphericity are assigned to both the steepest and flattest meridians, the biconic equation collapses into the simpler rotationally symmetric conic equation. 

Patient-Specific Model

The same methodology as described in Ariza-Gracia et al.²⁹ was followed. In particular, patient-specific data for the anterior and posterior surfaces were reconstructed using a Zernike polynomial expansion of order 7.³⁰ Where topographic data were not available, a quadric surface was used to extrapolate corneal elevation so that a minimum corneal radius of 6 mm is achieved for all corneal reconstructions (Fig. 3). Finally, a regularization algorithm was used to smooth the transition between patient-specific and extrapolated data to minimize mesh distortions. This top-down approach allows evaluation of whether the optomechanical performance is adequate for each modeling approach and to what degree each can be used to explain different surgical results. 

Figure 3.

 

Surface reconstruction starting from Pentacam point cloud.

Surface Ablation Profile

Derived from the geometry complexity, laser ablation profiles can be personalized and tailored to each previously selected geometry: from simple conic (or biconic) ablations that can correct myopia (and astigmatism) to wavefront-optimized ablation profiles for patient-specific surgeries. When building the ablation profile, the target refractive correction was set as the one indicated in the surgery treatment plan, provided by our clinical partner: the refractive goal aimed at correcting −4D of myopia and −1D of astigmatism at an angle of 170° for the selected patient. The treatment plan was exported directly from the software associated to the employed laser for the surgery, the excimer laser Wavelight EX500 by Alcon (Forth Worth, TX, USA). Because the software of the excimer laser Wavelight EX500 is proprietary of the company, direct access to the point-wise ablation profile was not possible, but the color map of the ablation profile in the Wavelight EX500 treatment report was reproduced by means of a MATLAB code. 

To assess the central ablation depth, we compared the theoretical values proposed in the literature^11,16,31 with the values that the software of the excimer laser Wavelight EX500 set by changing the dioptric target (Table 2). The theoretical approach brings to an underestimation of the central ablation depth with respect to the laser ablation, which is currently used in clinics; thus the ablation profiles in the three models (conic, biconic, and PS) were built by forcing the central ablation depth to be the same as indicated in the treatment plan (74.09 µm for the patient of interest) with the aim of comparing the outcomes with postsurgical keratometric data. 

Table 2.

 

Theoretical Versus Excimer Laser Ablation Depths

 

As we can see from Table 2, no theoretical ablation depth matches the one computed by the laser software, and the difference between the theoretical and the laser central ablation depths increases at higher dioptric targets. All the algorithms for the point cloud reconstruction and surface manipulation were developed in MATLAB R2022b. 

Conic Profile for Myopia

For the conic model we applied a conic ablation profile, given by Equation 1, where the apical radius of the ablation surface was obtained from Lensmaker's equation for thin lenses:  

\begin{eqnarray} {R_{ablation}} = \left( {\frac{D}{{n - {n_0}}} \cdot \frac{1}{{{R_m}}}} \right)\quad \end{eqnarray}

(5)

where D is the desired dioptric correction, R_m is the apical radius of the corneal anterior surface, as defined above. n₀ = 1 is the refractive index of air, and n = 1.3771 is the corneal refractive index. Note that this refractive index value must be used to compute the ablation profile and not the keratometric index of refraction n_k = 1.3375, commonly used to represent corneal keratometric values.³² The desired dioptric correction was taken directly from the treatment plan, that indicated a correction of −4D of myopia and −1D of astigmatism. Given the impossibility of correcting astigmatism with a conic ablation, we computed the ablation apical radius for the myopic component of the target (−4D), and we assigned the ablation depth as indicated in the treatment plan. 

The asphericity value for the ablation surface was obtained with the following equation³³:  

\begin{eqnarray} {Q_{ablation}} = \frac{{R_{ablation}^3}}{{R_m^3}}\left( {1 + Q} \right) - 1\quad \end{eqnarray}

(6)

where R_m is again the mean apical radius (mm) and Q is the asphericity of the corneal anterior surface. 

Biconic Profile for Myopia and Astigmatism

In this case, the biconic parameters (R_{s ablation}, R_{f ablation}, Q_{s ablation}, Q_{f ablation}) were computed with Equations 5 and 6 for the steepest and flattest meridians of the biconic model. The values of radii and asphericities were then substituted in Equation 2 to determine the biconic ablation profile. 

Standard Treatment: Wavefront-Optimized Ablation Profile

As already mentioned, the algorithm for the computation of the ablation profile is proprietary of the manufacturing company and cannot be accessed. The user manual of the excimer laser Wavelight EX500 by Alcon describes the different ablation options, available to be used by the surgeon, by presenting multiple examples. The most used treatment is the wavefront-optimized or “standard” treatment for the calculation of the ablation map, where the user is required to insert the target refraction and the laser's software establishes the ablation profile. This kind of treatment is not patient specific: the software will propose the same ablation map for patients with the equal refraction, disregarding individual characteristics. As it can be seen in the treatment plan of a patient who received this treatment, the ablation map shows a pattern that follows an aspheric shape. The ablation profile shown in the treatment plan of the patient was reproduced by using a biconic (Equation 2) and by enforcing the same central ablation depth as the one of the real treatment. The standard or wavefront-optimized treatment must not be confused with the wavefront-guided treatment or customized treatment, which uses corneal wavefront aberrations to determine the customized ablation profile needed to correct patient's visual acuity in a personalized fashion. 

Material Model

Corneal tissue is characterized by five different layers: epithelium, Bowman's layer, stroma, Descemet's membrane, and endothelium. The stroma constitutes 90% of the corneal thickness³⁴ and is mainly responsible for the mechanical response of the corneal tissue. For this reason, in many works about corneal biomechanics,^12,18,19,35 the assumption of representing corneal material with a monolayer model, made up only by the stroma, has been used. In this work, we made the same assumption with a slight modification: we reduced the CCT of 50 µm to account for the presence of the epithelium, which does not contribute to the mechanics of the cornea, to avoid an overestimation of the mechanical contribution of the stroma. This assumption was made necessary by the current lack of information on the mechanical properties of the other corneal layers and on the processes that drive the regrowth of the epithelial layer. 

In human corneas, the stroma is characterized by the presence of collagen fibers, which show a preferred orientation, being orthogonally arranged along the nasal-temporal (N-T) and inferior-superior (I-S) directions within the central region and circumferentially oriented in the limbal region (Fig. 4a).^36–40 Because the fibers are not perfectly aligned, in-plane and out-of-plane dispersion terms must be considered. The degree of out-of-plane dispersion varies in depth, meaning that the fibers are more aligned with N-T and I-S directions within the two posterior thirds and are more isotropically oriented within the anterior third. Consequently, collagen fibers have a strong contribution to the corneal mechanical behavior, which varies depending on the location. 

Figure 4.

 

Material model characteristics. (a) Schematic of theoretical collagen fibers distribution.³⁹ (b) Apical displacement - IOP curves comparison between experimental data⁴⁰ and the material model implemented with a UMAT subroutine of Abaqus. (c) Variation of the in-plane dispersion k_ip in the corneal geometry. (d) Variation of the out-of-plane dispersion k_op as a function of the local coordinate s.

In this section, large deformation kinematics theory is used to describe the hyperelastic anisotropic behavior of corneal tissue.⁴¹ Given the reference configuration of a body B_ref and a point x_p belonging to the body, the reference body can undergo a motion x = χ(x_p,t) toward its deformed configuration B_def through the deformation gradient F = ∇χ and its determinant J = det(F), which constitutes the jacobian. 

The invariants are:  

\begin{eqnarray} {\bar I_1} = tr{{\boldsymbol{C}}_{dis}}\quad \end{eqnarray}

(7)

 

\begin{eqnarray} {\bar I_i} = {{\boldsymbol{C}}_{dis}}:{\boldsymbol{a}}_0^i\, \otimes \,{\boldsymbol{a}}_0^i\,for\,i = 4,\,6\quad \end{eqnarray}

(8)

 

\begin{eqnarray} {\bar I_n} = {{\boldsymbol{C}}_{dis}}:{{\boldsymbol{a}}_n} \otimes {{\boldsymbol{a}}_n}\quad \end{eqnarray}

(9)

C_dis represents the distorsional right Cauchy-Green tensor:  

\begin{eqnarray} {{\boldsymbol{C}}_{dis}} = {\boldsymbol{F}}_{dis}^T{{\boldsymbol{F}}_{dis}} = {J^{ - 2/3}}{\boldsymbol{C}}\quad \end{eqnarray}

(10)

where C = F^TF is the right Cauchy-Green tensor. 

\({\boldsymbol{a}}_0^4\) and \({\boldsymbol{a}}_0^6\) are the unit vectors of the mean referential directions of the two families of fibers that characterize the corneal stroma, and \({\boldsymbol{a}}_n\) is the unit vector normal to the plane, that identifies the out-of-plane direction. 

The following strain energy density function for a nearly-incompressible material has been used to mimic the corneal tissue behavior³⁸:  

\begin{eqnarray} \psi = {\psi ^{matrix}}\left( {{{\bar I}_1}} \right) + \mathop \sum\limits_{i = 4,6} {\psi ^{fibers}}\left( {{{\boldsymbol{C}}_{dis}},{{\boldsymbol{H}}_i}} \right) + \psi \left( J \right)\quad \end{eqnarray}

(11)

where \({\psi ^{matrix}}( {{{\bar I}_1}} )\) where represents the isotropic contribution, ψ^fibers(C_dis, H_i) accounts for the fibers of the model and ψ(J) is the volumetric term. 

A Neo-Hookean model was chosen to describe the behavior of the extracellular matrix component of the tissue:  

\begin{eqnarray} {\psi ^{matrix}}\left( {{{\bar I}_1}} \right) = {C_{10}}\left( {{{\bar I}_1} - 3} \right)\quad \end{eqnarray}

(12)

where C₁₀ is a material constant that controls matrix's behavior. 

To model the anisotropic behavior of the collagen fibers (Fig. 4b), the hyperelastic Holzapfel-Gasser-Odgen model with dispersion parameters was used:^36,38 

\begin{eqnarray} {\psi ^{fibers}}\left( {{{\boldsymbol{C}}_{dis}},{{\boldsymbol{H}}_i}} \right) = \frac{{{k_1}}}{{2{k_2}}}\left( {{e^{{k_2}{{\left( {\bar I_i^* - 1} \right)}^2}}} - 1} \right)\quad \end{eqnarray}

(13)

where k₁ and k₂ denote the fibers stiffness and fibers non-linearity, respectively. 

The distorsional generalized invariant \(\bar I_i^*\) has the following form:  

\begin{eqnarray} \bar I_i^* &=& tr\left( {{{\boldsymbol{H}}_i}{{\boldsymbol{C}}_{dis}}} \right) = 2{k_{ip}}{k_{op}}{\bar I_1} + 2{k_{op}}\left( {1 - 2{k_{ip}}} \right){\bar I_1}\nonumber\\ && + \left( {1 - 6{k_{ip}}{k_{op}} - 2{k_{op}}\left( {1 - 2{k_{ip}}} \right)} \right){\bar I_n}\quad \end{eqnarray}

(14)

where H_i is the structure tensor, that quantifies the dispersion of the fibers:  

\begin{eqnarray} {{\boldsymbol{H}}_i} = A1 + B{\boldsymbol{a}}_0^i \otimes {\boldsymbol{a}}_0^i + \left( {1 - 3A - B} \right){{\boldsymbol{a}}_n} \otimes {{\boldsymbol{a}}_n}\quad \end{eqnarray}

(15)

with constants with constants A = 2k_ipk_op and B = 2k_op(1 − 2k_ip), where in-plane and out-of-plane dispersion are now introduced. 

The following equation for in-plane dispersion (k_ip) was used (Fig. 4c):  

\begin{eqnarray} {k_{ip}}\left( {\theta ,r} \right) = k_{ip}^{min} + \frac{1}{2}\left( {{k_{ip}}\left( \theta \right) - k_{ip}^{min}} \right)\left( {1 - cos\frac{{2\pi r}}{{{R_{TZ}}}}} \right)\quad \end{eqnarray}

(16)

where k_ip ∈ [0.1, 0.5], \(k_{ip}^{min} = 0.1\) and R_TZ = 5.5 mm is the radius of the transition zone from the cornea to the limbus. A homogeneous in-plane dispersion k_ip = 0.5 was assigned to the limbus.³⁸ 

To determine the out-of-plane dispersion term (k_op), it is necessary to define a local coordinate s ∈ [0, 1], parallel to the normal unit vectors. This local coordinate is equal to 0 at the anterior surface and to 1 at the posterior surface, changing throughout the corneal thickness (Fig. 4d). The equation for the out-of-plane dispersion is the following:  

\begin{eqnarray} {k_{op}}\left( s \right) = k_{op}^{min} + \left( {k_{op}^{max} - k_{op}^{min}} \right)\left( {1 - {e^{ - {\gamma _d}s}}} \right)\quad \end{eqnarray}

(17)

Where \(k_{op}^{min} = 1/3\) and \(k_{op}^{max} = 1/2\) and the constant γ_d controls the nonlinearity of Equation 17. Given that in the range of interest of our simulations (IOP = 0–20 mmHg), the non-linear behavior of Equation 17 is similar for any γ_d (Fig. 4b), a value of γ_d = 1 was used. 

Finally, the volumetric term is:  

\begin{eqnarray} \psi \left( J \right) = \frac{1}{D}{\left( {lnJ} \right)^2}\quad \end{eqnarray}

(18)

where D is the volumetric constant. 

The material model has been implemented by means of a User Material (UMAT) subroutine in the FORTRAN language. 

As it will be further explained later, one of the tested models was also built with the sclera to analyze the effect of different boundary conditions on the simulation outcome. Given the isotropic behavior of the sclera,⁴² the Neo-Hookean model for nearly incompressible material was assigned to this structure. 

In Table 3, we present the material constants used in our model. These material constants are suitable for an average healthy population, but they must be adjusted to describe the behavior of patient-specific models,⁴³ especially when dealing with pathological cases. As will be explained in the Clinical Application section, the selected material properties may influence the final optical result of the simulation because of the variability of mechanical properties in patients’ corneal tissue. 

Table 3.

 

Material Constants Used in Our Model

 

FE Model

Meshes for the FE model are automatically generated using previously acquired surface point clouds (see the section on Corneal Surface Reconstruction) with the software ANSA version 22.0.1 by BETA CAE Systems. Raw topographic data are provided, and no surface interpolation is performed to preserve the optical properties of the initial point clouds, given that interpolation algorithms of the meshing software are proprietary and might include settings that could interfere with the initial optical quality of the corneal surface. Thus the points constituting the cloud were directly used as nodes for creating the elements of the mesh. The mesh element type was set to quadratic tetrahedrons for two reasons: (i) using hexahedrons would lead to distorted elements because of the sharp angle between the anterior and ablation surfaces and (ii) quadratic elements avoid shear locking and provide improved accuracy. 

A mesh sensitivity analysis was conducted to choose the best compromise between the accuracy of the optomechanical performance of the model and the computational cost by analyzing the convergence of the displacement field and the corneal curvature (i.e., mean curvature; see Appendix). After choosing second-order tetrahedral elements, six different element sizes were tested (Table 4). 

Table 4.

 

Dimensions of the Five Meshes Tested in the Mesh Sensitivity Analysis

 

FE Simulation Pipeline, Boundary Conditions, and Loads

Three different boundary conditions were tested to determine whether they affected the optical outcome of the simulation (Fig. 5, left): 

• a. Fixed boundary condition at the limbal surface.
• b. Sliding boundary condition at the limbal surface. Only radial displacements are allowed, while polar and azimuthal movements are blocked.
• c. Symmetric boundary conditions at the scleral base, allowing only inplane radial displacement. Out-of-plane displacements were restricted.

Figure 5.

 

FE simulation characteristics. Left: Boundary conditions tested in this work. Right: Simulations steps.

Physiologically, the eyeball maintains its shape thanks to hydrostatic IOP exerted by the humor on the tissue. Although IOP can be modeled in different ways, two were selected in this study: 

(i) as a constant distributed surface load or (ii) as an equation of state simulating a fluid-filled cavity that allows for pressure changes within the eyeball due to external forces.^44,45 In both cases, IOP was set to patient's physiological value, obtained from CORVIS non-contact tonometry test; in particular, the value of the biomechanical corrected IOP (bIOP = 19.44 mm Hg) was used. 

Given that the patient's eye is subjected to a physiological IOP at the moment of the topographic acquisition of the data (i.e., a deformed geometry is measured without knowledge of the current deformation field), the first step before running a biomechanical simulation of any human organ is to recover the initial stress-free configuration. To do so, an in-house iterative algorithm was used to recover the unpressurized (reference) geometry.²⁹ 

Once the reference configuration is retrieved, a refractive intervention can be simulated. In the present study, the PRK simulation spans two steps (Fig. 5, right). 

Pressurization Step

IOP is applied to recover the topographic geometry (including both the deformation field and measured patient geometry). 

Surgical Step

PRK laser surgery is simulated by removing the ablation portion (i.e., an element set) from the anterior section of the cornea model. All of the simulations were run with the proprietary software ABAQUS version 6.13-1. 

Optical Analysis

In clinics, axial/sagittal and tangential curvature maps are the most used to assess the preoperative and postoperative status of any patient,⁴⁶ given their simple calculation. Curvature is usually measured in diopters (D) because its value accounts for the change in refractive index between air and the corneal surface ((n − n₀), where n = 1.3375 is the corneal keratometric refractive index, and n₀ = 1 is the refractive index of air). 

Both curvatures are used to evaluate patient's corneal surface state. More specifically, sagittal curvature is able to highlight refractive defects, especially when astigmatism is present; on the other hand, tangential curvature is more suitable for highlighting the presence of surface irregularities, like ectasia. 

A more sensitive method for surface analysis and detection of eventual irregularities is by means of differential geometry,^46,47 where Gaussian and mean curvatures can be computed by performing non-trivial calculations. These two curvatures represent the first and the second invariant of a quadric surface and can better identify surface irregularities with respect to axial and tangential maps because they are independent of the optical axis and are not influenced by the astigmatic defect, when present, highlighting only surface shape (see Appendix). 

Given that healthy corneas were considered in this study, that is no other pathology was present in addition to the refractive defect, only sagittal and mean curvature maps were used to analyze curvature changes. Moreover, to evaluate the methodology's optical accuracy, the keratometric values of the objective refraction given by the presurgical and postsurgical Pentacam topographies of the patient were compared with the results given by our models in terms of sphere (S) and cylinder (C) powers, which are commonly used by ophthalmologists to quantify the refractive change. To evaluate objective refraction in our models, the minimum and maximum curvatures (K₁ and K₂) were calculated by means of an ellipsoidal fitting.⁴⁸ Subsequently, the sphere was computed as  

\begin{eqnarray} S = \left( {{n_k} - {n_0}} \right) \cdot \frac{{\left( {{K_1} + {K_2}} \right)}}{2}\quad \end{eqnarray}

(19)

and the cylinder as⁴⁷ 

\begin{eqnarray} C = \left( {{n_k} - {n_0}} \right) \cdot \left( {{K_1} - {K_2}} \right)\quad \end{eqnarray}

(20)

 

Results

This section presents the results that are of interest to standardize the simulation of refractive interventions. First, the mesh sensitivity analysis is shown to determine the minimum mesh size that is accurate enough to ensure both mechanical and optical convergence. Second, different boundary conditions and loads are analyzed to gauge the possible impact on optomechanical outputs. Then, the comparison of different geometrical approaches (from conic to patient-specific reconstructions) is analyzed to determine the minimum amount of information required to capture general predictive trends in refractive interventions. Finally, the proposed methodology is applied to other eyes of four patients who previously underwent PRK to ensure accuracy of the modeling procedure. 

Mesh Sensitivity Analysis

For the sensitivity analysis (Table 4), the full FE pipeline was run using the conic geometrical approximation of the cornea. The displacement field and the mean curvature of the anterior corneal surface were used to control the tradeoff between accuracy and computational effort. 

Curvature convergence results are more restrictive than displacement fields as they represent the second derivative of the corneal geometry. Hence, for the sake of simplicity, only mean curvature results are presented below (Fig. 6). Overall, element sizes between 0.1 and 0.2 mm represent a good trade-off between accuracy, both optical and mechanical, and computation time. Increasing beyond 0.2 mm would increase the optical error up to 0.3D, which is in the sensitivity range of current clinical topographers and, therefore, not suitable for optomechanical simulations. Changes and convergence of the displacement fields agreed with those of the optical fields. 

Figure 6.

 

Percentage error in dioptric power εH of the optical zone of the anterior surface (dashed line) of the tested meshes with respect to the finest mesh.

Boundary Conditions and Loads

In this section, the main optical differences related to changes in boundary conditions and loads are presented. For the sake of simplicity, only patient-specific results are presented in this section. 

As each boundary condition presents advantages and disadvantages, the pairwise difference in terms of mean curvature between each simulated postsurgical intervention was calculated to determine the impact of different simulation choices (Fig. 7). The greatest postsurgical optical difference (up to −0.35 D within a 6 mm optical region) occurs when choosing a fixed boundary condition at the limbal region over a sliding boundary condition at the scleral plane of symmetry. When considering a sliding boundary versus a fixed boundary at the limbal region, the maximum optical error reached −0.3 D at the periphery of the 6 mm central optical region. Finally, the difference of considering a sliding boundary condition at the limbus over a sliding boundary condition at the scleral plane of symmetry was minimal, with an average difference close to 0 D over a 4 mm central optical region with a maximum of 0.15 D, which is below the resolution of the topographer. 

Figure 7.

 

Comparison among the three tested boundary conditions (BC). Left: mean curvature difference ΔD, measured in diopters, between the models with fixed BC and symmetrical BC at the base of the sclera; center: mean curvature difference between the models with fixed BC and sliding BC. Right: Mean curvature difference between the models with symmetrical BC at the base of the sclera and sliding BC.

From a loading perspective, considering a constant distributed pressure load or a fluid cavity interaction did not show a relevant difference. In particular, when considering the IOP as a fluid-filled cavity controlled by an equation of state, there should exist a drop in pressure because of the higher postsurgical corneal compliance. In the present case, the postsurgical drop in pressure was negligible (i.e., <0.05 mmHg). Therefore it should not be necessary to include the complexity of including a fluid-filled cavity when simulating refractive interventions. 

Optomechanical Analysis for Geometrical Complexity

To analyze the difference in the optomechanical behavior for increasing geometrical complexity, three different FE models were created, including conic, biconic, and patient-specific models. They all presented the same boundary conditions and loads (sliding boundary condition; bIOP = 19.4 mm Hg as a distributed load pressure) and mechanical properties (anisotropic material model for the stroma, see Table 3). 

From a purely mechanical perspective, we considered the increase in terms of stress and strain values due to the ablation, that is the difference between the ablation and the pressurization steps. Just observing maximum stress and strain values, the difference between considering the most complex modeling (patient-specific) or the simplest was low: both models presented similar maximum principal stresses (i.e., 3.49 kPa vs. 3.70 kPa; see Table 5), with postsurgical stress increments around 17%. However, strain patterns presented differences (see Fig. 8) that are related to the highest fidelity of geometrical features for biconic and patient-specific models and the use of different ablation profiles. For all cases, strain distributions followed a regular pattern in the optical zone of interest (Ø_OZ = 6.5 mm), given the regularity of the geometries and of the ablation profiles. The strain pattern in the conic model presented a rotational symmetry, while for the biconic and patient-specific models the pattern was rotated and aligned with the astigmatic axis. Regardless, the highest mechanical effect always occurred in the same central region where the ablation was performed, and the corneal thickness decreased. In all three models, the highest strains were observed at the ablated anterior surface, where the surgery mostly acts; however, strains also concentrated at the boundary of the optical zone, due to the sudden change in the curvature caused by the ablation cut. The stresses mostly concentrated in the central portion of the posterior surface. As it can be observed in Figure 8, biconic and patient-specific models showed a quite similar mechanical behavior in terms of stress and strain distribution. 

Table 5.

 

Maximum Values of Maximum Principal Stress and Strain Changes Affecting the Three Geometries and the Percentage Increase Because of the Laser Ablation

 

Figure 8.

 

Maximum principal stress (top) and strain (bottom) changes induced by the ablation in the three models considered with respect to the pressurization step.

From an optical perspective, however, different geometrical complexity led to different outcomes. By mathematical definition, conic geometries cannot capture astigmatism, whereas biconic and patient-specific geometries can. By comparing our models with presurgical and postsurgical Pentacam data (see Table 6), we can confirm that all three models achieved the target spherical power correction; moreover, the biconic and PS models achieved a close cylinder correction to the real correction that the patient received. 

Table 6.

 

Comparison of the Objective Refraction of the Three Model With Respect to Presurgical and Postsurgical Pentacam Data, Calculated in the Optical Zone

 

Sagittal curvature maps (Fig. 9) extend and support the concepts introduced by objective refraction. Presurgical conic geometries present rotational symmetry, and thus astigmatism is disregarded. Biconic can introduce astigmatism as shown by its bow-tie pattern (Fig. 9, center). Finally, patient-specific patterns present more complex distributions that can be used to assess whether the surgery regularized the surface or not. Our results show that, after surgery, the conic model still presents a fairly rotational symmetric curvature related to its inability to represent astigmatism. Likewise, although the biconic model constitutes an analytical representation of corneal geometry and has a symmetry axis that does not allow to replicate patient's detailed surface features, the sagittal curvature pattern turned out to be very close to the PS one. Overall, as expected, the PS model was the most accurate in terms of reproducing patient's surface features: the sagittal curvature maps of the PS model turned out to accurately reproduce the pattern of the presurgical and postsurgical maps given by the Pentacam topographies. 

Figure 9.

 

Presurgical and postsurgical sagittal curvature maps for the conic, biconic, and patient-specific cases (central optical zone; R_OZ = 3 mm).

Clinical Application

In order to check the accuracy of the proposed methodology, the surgery was simulated on other four eyes of patients who previously underwent PRK. When considering different patients, unique material properties (Table 3) may not be representative enough of the variability of mechanical properties among patients, as already underlined in the Material Model section, and may affect the refractive result of the simulation. Because no device available in the market is currently able to directly extract PS material properties, we considered two different cases for the four eyes: two values of constant C10 were considered (30 kPa and 15 kPa) to account for two different stiffness values with a 50% variation with respect to the average value found in literature.³⁶ We chose to change only the constant C10, because it controls the extracellular matrix behavior, which mostly affects corneal tissue's mechanical response. The results of the objective refraction for the four patient-specific models are shown in Table 7. 

Table 7.

 

Objective Refraction of the Treatment Plan, Pentacam Topographies and Patient-Specific Models

 

Discussion

Nowadays, clinicians evaluate patient eligibility for laser surgery using corneal topography (e.g., Pentacam by Oculus or Sirius by CSO) and noncontact tonometry (e.g., CORVIS ST by OCULUS). The former allows to determine corneal geometry and quantify the refractive errors due to corneal shape abnormalities, when present (for instance, the astigmatic defect is caused by two different principal curvatures, that cause corneal shape to be similar to a rugby ball). On the other hand, CORVIS ST is used to provide only a qualitative assessment of corneal mechanical response.^42,49 When the cornea is deemed healthy in terms of geometry (central corneal thickness [CCT] > 490 µm), regularity, and biomechanics, clinicians proceed with the surgery. Otherwise, clinicians consider that the postsurgical residual stromal bed could not be sufficient to withstand the biomechanical stresses. Therefore the patient is considered unsuitable for refractive surgery because of the high risk of developing postsurgical corneal ectasia.⁷ 

This interplay between geometry and mechanics cannot be assessed with any clinical tool currently available in clinics. Hence, developing a predictive tool that uncouples both the cornea's mechanics and optics is helpful in supporting clinical diagnosis and minimizing postsurgical complications after a refractive surgery procedure. In particular, FE models can be useful for investigating aspects that cannot be assessed in daily clinical practice, such as mechanical changes that will lead to unforeseen postsurgical optical changes, or to help plan surgeries in which the mechanical dependence is so strong that they still remain an open clinical problem, such as astigmatism correction.⁵⁰ Unfortunately, reproducibility remains an issue unless all the modeling details are disclosed. 

The present work describes a methodology to determine the minimal requirements to simulate a refractive intervention such as PRK. To create reproducible models, we advocate that standards must be set to avoid the divergence in results and lack of reproducible methods. A modular FE methodology was developed to determine and disclose the minimal features required to simulate a refractive surgery. First, the geometry is provided as input to the FE pipeline, which will calculate the laser ablation profile of the appropriate dioptric correction, according to patient's treatment plan. Second, the FE mesh of the model, the boundary conditions, loads, and material definition are created automatically. Third, the surgery simulation is performed, including the pre-stretch of the eyeball due to the IOP and the laser ablation. Fourth, the cornea's optics in terms of objective refraction (Table 6) and curvature maps is calculated using an in-house optical algorithm. 

First and foremost, our analysis suggests that for presenting a good optomechanical performance, a FE model for a refractive surgery should have a mesh size of around 0.1 mm in the optical region. Even if mechanical variables present a good convergence for a larger mesh size (between 0.2 and 0.3 mm), optical variables perform better for a smaller mesh size. Commercial topographers present similar grid resolutions for clinical devices, which aligns with our analysis. A word of caution is necessary at this point: commercial meshing software might microscopically change the original raw data acquired by the topographer. Even if this change seemed negligible, a micrometric change in the elevation values of the point cloud results in an altered curvature and, therefore, presurgical corneal optics would not be preserved. 

Boundary conditions are more complex to select. In general, it does seem clear, and well accepted in the literature,^15,51 that restricting the displacements at the limbal region is too restrictive: it would be equivalent to assuming that the scleral tissue is infinitely rigid and that the cornea cannot rotate or slide in that region. Choosing between considering a sliding boundary condition at the limbal region or including a portion of the sclera with a sliding boundary poses interesting modeling decisions with different underlying hypotheses. Assuming sliding at the limbus implies that the sclera is infinitely rigid against rotations while fully compliant to radial displacement along the sliding boundary condition. Consideration of a scleral portion allows for a more realistic transition at the limbal region where radial displacements and rotations are accounted for. However, including the scleral tissue poses the uncertainty of introducing additional mechanical variables related to the tissue material properties. Unless a good mechanical characterization and modeling of the scleral material properties are performed, including the sclera introduces an additional hypothesis to be dealt with. As both present small relative differences in the optical zone (Fig. 7), we would advise avoiding including the scleral portion unless a proper mechanical characterization has been performed. This approach reduces computational time and restricts uncertainties, resulting in an acceptable relative difference of up to 0.15 D (i.e., below the topographer's resolution of 0.25 D). 

Regarding the load conditions to simulate the IOP, our results suggest that a simple approach of considering the aqueous humor as a distributed surface load applied on the inner eyeball surfaces should suffice for the aim of the present simulation. Increasing the complexity of the modeling for the IOP does not provide any additional information. As outlined in the literature,^44,45 reported clinical differences between presurgical and postsurgical IOP are probably related to noncontact tonometry measurements and their correlation to corneal thickness. This seems to be why CORVIS ST provides a corrected IOP measurement that considers diverse patient biomarkers such as CCT or age.^52,53 

To be noticed that the assumption of considering the IOP as a distributed surface load can be made only in cases where changes in the volume of the anterior chamber of the eye are negligible, as for PRK surgery. When there is a change in volume within the eyeball cavity, this assumption would not be valid anymore and the fluid cavity setting could be used. Moreover, circadian oscillation and ions transport inside the eye cavity, which may affect IOP value, are here disregarded. 

To assess the degree of geometrical complexity required, patient-specific and theoretical (conic and biconic) geometries were tested. From a mechanical standpoint, differences between the presurgical and postsurgical mechanical stress were below 18% and focused around the central region in which the surgery was performed (Fig. 8). In absolute values, mechanical changes remained low and similar between corneal models (approx. 3.5 kPa), regardless of the geometrical complexity or the ablation pattern. However, strains patterns included different features with increased geometrical complexity. In particular, not including complex features results in losing spatial information about the postsurgical strain distribution, which could be valuable to assess vulnerable corneal regions with altered mechanical properties which are thought to be linked to postsurgical ectasia. Regarding the stress distribution, it is interesting to highlight that stresses induced by the surgery mostly concentrated at the posterior surface. Further investigation would be required to assess whether these concentrations could cause a risk in surgery's safety on myopic patients, where the surgeon is unsure if they are eligible for receiving the surgery. 

From an optical standpoint, geometry does play a role, and patient-specific models gave the most accurate result. In particular, conic models should not be used when corneal astigmatism is present as they completely disregard it, while biconic models can include it. The spherical power of both conic and biconic models (which approximated the patient-specific geometry) matched quite well the postsurgical refractive target of the patient, as well as the patient-specific model. Taking into account the lower precision of the analytical models (conic and biconic) in reproducing surface features with respect to the PS one, we believe that they still represent a good compromise between low computational cost and optomechanical outcome for a fast presurgical evaluation. In particular, the conic model should be employed only for reproducing regular myopic corneas without astigmatism, while the biconic model can be used to reproduce corneas with both myopia and regular astigmatism. When patients have a more irregular surface, for instance in the presence of irregular astigmatism, that means that the angle between K₁ and K₂ is not 90°, the biconic model is not anymore able to reproduce the specific shape of patient's corneas due to the fact that this model is symmetric along the astigmatic axis. In that case, a patient-specific model is mandatory. Still, we emphasize that patient-specific model remains the most accurate in terms of optomechanical performance. 

Finally, the proposed methodology was applied to other four eyes of patients who underwent PRK, to check the accuracy of the modeling procedure for the patient-specific case. As it can be observed in Table 7, a mean error of 9.4% in the spheric correction with respect to Pentacam keratometry is obtained, depending on the selected material properties. For each patient, a very low error (0%–2%) was obtained by changing the mechanical properties of the model. Choosing appropriate material properties is crucial for obtaining an accurate refractive result. 

This study is not exempt from limitations. First and foremost, a single refractive surgery (i.e., PRK) was considered. All laser refractive surgeries work by the same principle: removing a small portion of the corneal stroma to change the corneal refractive power. Even if mechanical deformation (resp. stresses) distributions are expected to be different between different refractive procedures (i.e., LASIK, SMILE or PRK), optical requirements should be in the same range from a modeling procedure point of view. Moreover, the proposed methodology is flexible and could be adapted for additional surgical procedures. 

Second, average PS material properties are used. As already anticipated, this is a limitation for all the current studies on corneal biomechanics, given that PS material properties are still unknown, but research is advancing fast in this field.^{4,42,43,51,54} Third, even if we accurately reproduced the ablation profile from patient's treatment plan in a central optical region of 6.5 mm, we did not have direct access to the specific point-wise pattern, given that laser ablation algorithm is proprietary of the manufacturing company. 

Moreover, Equations 1 and 2, used for determining the ablation profiles, present a limitation, that is that they do not consider the biomechanical response of the corneal tissue when it receives the surgery, treating the cornea as infinitely rigid. We suppose that the standard treatment computed by the laser machine takes into account the biomechanical response of the tissue by means of a correction factor, that would explain the discrepancy between the theoretical and the real ablation depths (Table 2). For this reason, we chose to apply the ablation depth as indicated by the laser machine and not the theoretical one. 

Finally, five patient-specific geometries were included in this study. Accounting for a larger cohort of patients does not change the basic principles of this study but would provide a population-based clinical validation of the proposed procedure, which will constitute the future development of the current work. Once a standardized protocol is obtained, it may be used to isolate the 5% of patient who miss the postsurgical refractive target.⁵⁵ Also, because healthy patients were considered (i.e., when the cornea is characterized by a regular thickness and average corneal tissue properties),^36,38 no mechanical abnormalities were introduced (e.g., mechanical weakening that could lead to ectasia) and cannot be observed in the subsequent mechanical analysis. 

Conclusions

To summarize, when interested in developing a numerical model with good optomechanical performance, patient-specific models are required. Care must be taken when using biconic models to represent common optical parameters (both spherical and cylindrical powers), mostly used by optometrists and ophthalmologists, to describe the state of a patient's refractive state. Additionally, (i) mesh size should be around 0.1 mm to capture optical details properly; (ii) boundary conditions is a complex topic but, in the absence of a detailed mechanical characterization of the sclera, we would advise using a sliding boundary condition at the limbal region; and (iii) IOP can be simplified as a distributed surface load. In the future, and to ensure reproducibility, we encourage all authors to report all the details and settings of the algorithms. 

Acknowledgments

Supported by the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No 956720 and the Department of Industry and Innovation (Government of Aragon) through the research group Grant T24–20R (cofinanced with Feder 2014-2020: Construyendo Europa desde Aragón). Part of the work was performed by the ICTS “NANBIOSIS” specifically by the High Performance Computing Unit (U27), of the CIBER in Bioengineering, Biomaterials and Nanomedicine (CIBER-BBN at the University of Zaragoza). 

Disclosure: B. Fantaci, None; B. Calvo, None; R. Barraquer, None; A. Picó, None; M.Á. Ariza-Gracia, None 

Appendix

Different curvature maps can be considered when analyzing a refractive surface. Axial and tangential curvature maps are the most used by clinicians, given their easy computation and interpretation.⁴⁶ Axial/sagittal curvature is given by the inverse of the distance of a point on the anterior surface to the optical axis along its normal; tangential curvature is defined as the inverse of the instantaneous radius of curvature, obtained by fitting a sphere at each location of the anterior surface. Gaussian and mean curvatures⁴⁷ are two other useful maps, which represent the first and the second invariant of a quadric surface. At each surface point, it is possible to identify the flattest and the steepest meridians, defined as principal curvatures (K1 and K2, respectively). Mean curvature (H) and Gaussian curvature (K) are the arithmetic and geometric means of the principal curvatures, respectively. 

They can be calculated as:  

\begin{eqnarray} H = \frac{{\left( {1 + h_x^2} \right){h_{yy}} - 2{h_x}{h_y}{h_{xy}} + \left( {1 + h_y^2} \right){h_{xx}}}}{{2{{\left( {1 + h_x^2 + h_y^2} \right)}^{3/2}}}}\quad \end{eqnarray}

(21)

 

\begin{eqnarray} K = \frac{{{h_{xx}}{h_{yy}} - h_{xy}^2}}{{{{\left( {1 + h_x^2 + h_y^2} \right)}^2}}}\quad \end{eqnarray}

(22)

where h(x) represents the elevation of the corneal anterior surface (i.e. the height of the cornea over its surface, h_x = ∂h/∂x, h_y = ∂h/∂y, h_xx = ∂²h/∂x², h_yy = ∂²h/∂y² are the first and the second partial derivatives, respectively). For a complete clinical evaluation, the best option is to analyze the four curvature maps, to have a comprehensive view of the clinical status of the patient's corneas.