The Effect of Intraocular Pressure on the Outcome of Myopic Photorefractive Keratectomy : A Numerical Approach

Photorefractive Keratectomy (PRK) is a surgical procedure widely performed to correct myopia. In this work, the effect of the intraocular pressure (IOP) on the refractive correction achieved by the PRK surgery was analyzed using a numerical model. Simulations of PRK surgery at 10, 15 and 21 mmHg of IOP were performed and the post-surgical diopters were estimated. For low and medium values of IOP (10 and 15 mmHg), the computed results were close to those used by clinicians based on experience and defined without considering the IOP, while an undercorrection was predicted for the highest value of IOP (21 mmHg). From these results, we suggest that IOP should be considered in the determination of the depth of ablation, in addition to other factors such as the level of myopia or the corneal central thickness.


INTRODUCTION
The human eye is an optical system of around 59 diopters (D) composed of four refractive media: cornea, aqueous humour, crystalline lens and vitreous.The light crossing these media reaches the retina which sends electrical impulses along the optic nerve to the brain where the image is formed.The cornea is the interface between the eye and the external world.It has both protective and refractive functions being the lens with highest optical power, about 43 D. It is a convergent lens with a diameter of around 11.5 mm and central corneal thickness (CCT) between 0.503 and 0.565 mm in adult non-pathologic eyes [1].The corneal tissue joins the sclera in a ring-shaped part known as limbus.The cornea is composed of five layers: epithelium, Bowman's membrane, stroma, Descemet's membrane and endothelium.The epithelium is the outermost layer of the cornea, comprising about 10% of its thickness [2].The most important layer to structural behaviour is the stroma, whose thickness accounts for around 90% of the total corneal thickness, being composed of several lamellae collagen fibers running parallel among them [3,4].
Laser is a technique widely used in eye refractive surgery for the correction of myopia, hyperopia and astigmatism.In particular, this work is focused on the correction of myopia, which is a refractive error due to an overgrowth of the eye resulting in an excessively long eyeball in which incoming light produces image in front of the retina.A patient with myopia sees near objects clearly but distant objects blurred.Photorefractive keratectomy (PRK) consists of the ablation of the corneal tissue by laser to reshape the cornea, modifying its optical properties and consequently those of the whole eye.The bases of this technique were introduced by the Spanish ophthalmologist J. Barraquer [5].Barraquer developed the first microkeratome to cut thin flaps in the cornea to alter its shape, calling this procedure keratomileusis.PRK was subsequently developed in 1983 at Columbia University by S. Trokel in 1983 [6], who outlined the potential benefits of using the excimer laser in refractive surgery.In PRK surgery, part of the corneal stroma tissue is removed so that the anterior stromal surface is reshaped to correct the refractive error, leaving a residual stromal bed.Although the epithelium is removed before-corneal ablation, epithelial cells located at the edges of the residual epithelium usually complete the wound healing process in 3-4 days [7].Removing the epithelium involves some clinical post-operative issues related to time of healing and ectasia risk that should be considered by the ophthalmologic surgeon to decide whether to perform this type of surgery or alternative ones like LASIK (laser in situ keratomielusis).
A complete examination of the cornea is usually accomplished before the PRK surgery in order to establish the optimal depth of corneal ablation.Sometimes this examination warns the surgeon of the risks of this refractive surgery in the patient, due to a small pachymetry of the central cornea or some other anatomical or physiological features that may cause post-surgical problems [8].
Finite-element-based biomechanical models of the eye have been presented recently as a powerful tool for a better prediction of the effects of refractive surgeries [9][10][11][12][13].One of the main problems still remaining for its systematic use corresponds to the simulation of the complex behavior of the tissues involved.These tissues present a hyperelastic, incompressible and anisotropic constitutive behavior, strongly dependent on the physiological collagen fiber distribution.Hyperelastic anisotropic models have been recently employed in several studies [14][15][16] to simulate different aspects of eye functional behaviour.
In this work, we use a biomechanical model of the human eye previously developed by ourselves [17] to investigate the effect of the intraocular pressure (IOP) on the PRK

462
The Effect of Intraocular Pressure on the Outcome of Myopic Photorefractive Keratectomy: A Numerical Approach refractive surgery outcome.IOP is a very significant parameter since it is directly related to the stress of the corneal tissue in physiological conditions and after surgery.Most previous works on refractive surgery simulation used a mean value of the IOP of about 15 mmHg [10,14,18].However, normal IOP varies between 10 and 21 mmHg [19].The hypothesis of this study is that the effect of IOP is significant and therefore should be taken into account in determining the corneal ablation parameters for the PRK surgery.To assess this hypothesis, we evaluate the effect of IOP on the outcome of PRK surgery for myopia correction.

Finite Element Model
A three-dimensional finite element model of the anterior half ocular globe previously developed by our team [17] was used to simulate the PRK surgery, considering average dimensions (see Figure 1a).Axial symmetry was assumed for the purpose of simplification, as in other published studies [14,16].
The model covered three different parts: cornea, limbus and sclera (see [17]).Since the stroma is the most important structural layer and accounts for around 90% of the total thickness of the cornea, the other layers (viz.epithelium, Bowman's membrane, Descemet's membrane and endothelium) of the cornea were not considered in the present model.The assumed dimensions of the cornea were the following: 7.5 mm in anterior radius [20], 12.0 mm in diameter, 550.0 µm in axial thickness, and 650.0 µm in edge thickness at the limbus [16].The posterior surface of the model was assumed conical, with an apical radius of 6.04 mm.The limbus was modelled as a narrow ring with an internal diameter of 12.0 mm and an external diameter of 12.5 mm, with increasing thickness from 650.0 µm to 1.0 mm.The sclera was considered spherical, with an external radius of 13.0 mm and thickness of 1.0 mm.Cubit 10 (Computational Modeling Sciences Department, Sandia National Laboratories, Albuquerque, New Mexico, US) [21] and I-DEAS v.9 (Structural Dynamics Research Corporation, Milford, Ohio, US) [22] were employed to create the geometry and to mesh the model.After a first convergence analysis, the final mesh was composed of 25,344 8-noded linear brick, hybrid elements with constant pressure (C3D8H) elements and 27,791 nodes.The central area of the cornea was composed of prismatic elements C3D6H.Five layers of elements were established in the area close to the outer surface of the central cornea, each corresponding to a 2-diopter (D) myopic correction.The apical size of the elements of each of those layers was established from the estimated incremental depth of ablation for 2D, corresponding therefore to myopic corrections of 2D to 10D, respectively (see Figure 1b).The density of the mesh was improved by increasing the number of elements until the maximum displacement changed by less than 2.5%.
The model was constrained by the following boundary conditions: the displacements of the nodes on the medial plane of the ocular globe in cylindrical coordinates were restrained both in the circumferential and longitudinal directions; only radial displacements were allowed.Three values of IOP were used for the simulations: 10 mmHg, 21mmHg and the commonly used mean value of 15 mmHg [9,23].

Constitutive Model
Human cornea in vivo is a highly porous tissue filled with biological fluid.Approximately 80% of the weight of the cornea is water [24].Among the five layers in 464 The Effect of Intraocular Pressure on the Outcome of Myopic Photorefractive Keratectomy: A Numerical Approach the cornea, the stroma forms about 90% of the thickness and is divided into 300-500 lamellae [3,4] where long collagen fibers are embedded in ground substance, mainly composed of proteoglycans and water.Collagen fibers lie parallel and run along the whole length of the lamella.
In the central region of the cornea, fibers are disposed predominantly along the superior-inferior and nasal-temporal directions while near the limbus they follow the circumferential direction to form a pseudo-annulus.There are some other fibers randomly oriented throughout the cornea.Such microstructure and the different distributions of collagen fibers result in a well-known anisotropic behaviour for the cornea [4,25].In this work, the cornea is considered as an anisotropic material with material properties associated with two preferred directions corresponding to two local families of collagen fibers, defined in terms of two local vectors at each point.Although it is known that the orientation of the collagen fibrils all over the stroma is not homogeneous, in the central part of the cornea subjected to ablation in the surgery, a strong orthogonal nasal-temporal and superior-inferior orientation is observed [4,25].Therefore, we considered the following distribution of fibre directions: one in the nasaltemporal direction and the other in the superior-inferior direction throughout the cornea except for the periphery where, following previous studies [13,15], we considered the preferential circumferential direction of collagen fibers by taking the limbus as a material similar to the cornea with only one preferential direction (see Figure 2).
Usually, the description of the constitutive behaviour of this type of material relies on the identification of an appropriate strain energy density function from which the stressstrain relation and the local elasticity tensors are derived [26,27].We postulate the existence of a unique decoupled representation of the strain-energy density function Ψ [28] to characterize isothermal processes.Because of the directional dependence on the deformation, we require that the function Ψ explicitly depends on both the right Cauchy-Green tensor C = F T F, where F is the deformation gradient tensor, and the fiber directions Two families of collagen fibers are present in the cornea, corresponding to the following preferred directions: superior-inferior (n 0 ) and nasal-temporal (m 0 ).The preferred orientation of the collagen fibers in the limbal tissue is circumferential (n 0 ).m 0 and n 0 in the reference configuration.Since the sign of m 0 and n 0 is not significant, Ψ must be an even function of m 0 and n 0 and so it may be expressed as Ψ = Ψ(C, M, N) where and are structural tensors [27].Based on the volumetric kinematic constraint, the free energy can be written in decoupled form as (1) where Ψ vol (J) describes the change in volume (volumetric) and Ψ -(C -, M, N) the change in shape (isochoric) responses of the material.Both are given scalar-valued functions of and N [29].For the purposes of this model, the elastic response of the tissue was assumed to arise from the resistance of the collagen fibers and the matrix, that is, from a unique strain energy function defined as in eqn.(1).Ψ vol (J) is expressed as the following [30]: (2) which quasi-enforces the null volumetric change depending on the value of the penalty coefficient 1/D.We used a strain energy function that combines a Mooney-Rivlin's model [15] to treat the isotropic part corresponding to the matrix behaviour and the Holzapfel's model [29] to take into account the anisotropic response of the tissue due to the collagen fibers: ( where and are the first two modified strain invariants of the symmetric modified right Cauchy-Green tensor (note that ).The pseudo-invariants and have a clear physical meaning, the square of the stretch λ along the fiber directions.To obtain the material constants for the cornea, a non-linear regression method was employed using the stress vs. strain curve reported by Hoeltzel [23], who performed uniaxial tests with strips of human corneal tissue and observed a non-linear behaviour.Since k 1 , k 2 represent one direction of fibers and k 3 , k 4 the other, they were assigned the same values in order to model all the fibers equally regardless of the orientation.The values for the sclera parameters were taken from published experimental results obtained again from tests of human scleral tissue and curve fitting by means of least squares method [31].The material constants for cornea, limbus and sclera used in this work are shown in Table 1.
The Effect of Intraocular Pressure on the Outcome of Myopic Photorefractive Keratectomy: A Numerical Approach

Surgery Simulation
Refraction is the change in direction of a wave when it passes from one medium to another.This is described by Snell's law: n 1 sin α = n 2 sin β, where n 1 , n 2 are the refractive indexes of the media, and α, β are the angles of incidence and refraction, respectively, of the wave.The orientation of the reflected wave is a function of both the incidence angle and the refractive indexes ratio of the media.Most refractive surgery techniques are based on Snell's law, modifying the curvature of the anterior surface of the cornea and therefore the optical power.For myopic PRK surgery, the depth of the ablated tissue to modify the curvature of the cornea has to be predetermined to achieve the expected correction of myopia.The depth of the material to be removed can be calculated by the expression proposed by Munnerlyn et al. [32]: where t 0 is the depth of the tissue to be removed, S is the diameter of the optical zone to be ablated by laser, as depicted in Figure 3. D is the optical power of myopia to be corrected, in dioptres, and n = 1.377 [32], is the refractive index of the cornea.Figure 3 shows a scheme of the ablation procedure that the cornea's radius of curvature is changed from R 1 to R 2 by the procedure.Abaqus v6.5 was used to perform all the simulations [33].Before simulating the surgery, the numerical model must reproduce the real physiological conditions of a human eye [27].The reference geometry is usually taken as a stress-free configuration, but the actual stress-free configuration of the cornea is not known a priori.The cornea deformed by the IOP [13,15] which causes a stress-distribution must be included in the model before actual loading.A deformation gradient F 0 n is obtained following an iterative procedure by Lanchares et al. [17].Briefly speaking, F 0 n is computed to get an initial configuration that leads to stresses in equilibrium with the physiological IOP, which has a different value for each patient.The maximum value of the initial maximal stress found was 0.03 MPa located in the circumferential periphery of the cornea.Since the limbus is stiffer than the cornea due to the circular orientation of the collagen fibers, a concentration of stress is clearly observed in such area.
In this work, five cases of laser surgery corresponding to corrections of 2, 4, 6, 8 and 10D of myopia were simulated.A value of 6 mm was chosen for the diameter of  ablation S [34,35].Values of t 0 for these five cases calculated by eqn.( 4) are shown in Table 2.They were used to mesh the model of the cornea, creating layers of elements to be successively removed in a five-step calculation.Simulation was also performed for two other values of IOP, 10 and 21 mmHg, both considered to be within the range of healthy IOP.

RESULTS
Before the simulation of the PRK surgery, the initial stress distribution was calculated to simulate the physiological conditions of the eye.Figure 4 shows the displacement distribution under an IOP of 15 mmHg in the FE model with and without introducing the initial stress.In both cases, the distribution follows a similar pattern, but the maximum displacement is three orders of magnitude higher in the case of the stress-free initial configuration than the case with initial stress, 9.663 • 10 − 2 mm versus 6.688 • 10 −5 mm.After surgery simulation, a stress-equilibrium configuration is reached and the deformed geometry of the anterior corneal surface determines the optical power.The

468
The Effect of Intraocular Pressure on the Outcome of Myopic Photorefractive Keratectomy: Parameters of the laser ablation to correct myopia [32].S is the diameter of the cornea's optical zone; t 0 is the maximum depth of laser ablation; R 1 and R 2 are initial and final radii of curvature.mean radius of curvature of the corneal anterior surface inside the optical zone was determined by curve fitting through least squares based on the deformed positions of the nodes in such zone (Figure 5).The optical power in diopters was then calculated by the following equation [32]: (5) where n = 1.377 is the corneal refractive index and R is the mean radius of curvature in meters.
Figure 5 shows the shape of the anterior corneal surface before ablation, reproducing the in vivo configuration with the initial stress distribution associated with 15 mmHg of IOP.The other five contours correspond to the stromal bed profiles at the five different depths of ablation.The markers represent the nodes used for fitting the radius of curvature R. We considered two values for the central area radius r ca , 1.5 mm and 2.5 mm, to estimate R, which varies with the number of nodes considered within the central area and therefore the computed myopic correction.These two values of r ca are used by some topographers (Orbscan II [36]) to estimate the corneal topography and pachymetry and to compute the keratometry.As shown in Figure 5, 11 nodes fall within the central area of 1.5 mm radius, and 19 nodes in an area of 2.5 mm radius.
The five-step simulation was carried out with and without the introduction of the initial stress distribution caused by the IOP.The values obtained for the dioptric change in both situations are tabulated in Table 3 and plotted in Figure 6.When the initial stress Journal of Healthcare Engineering • Vol.   which reduces accuracy and reliability.Therefore, the initial stress distribution must be considered in the model to reproduce the physiological conditions and to obtain more consistent results for the surgery simulation.Figure 7 shows the maximum stress distribution corresponding to the simulation of 2D and 10D ablations.Although the maximum principal stress occurs in the periphery of the cornea, we focus on the center of the new anterior surface (stromal bed) to evaluate the effect of laser ablation.In this optical area, the maximum principal stress increases from the 2D ablation to the 10D case where a maximum value of 2.07 × 10 −2 MPa was reached.In each step of the simulation, a central layer of elements of the central anterior surface was removed causing a modification of the anterior curvature and a decrease in the thickness of the cornea.As a result of stress concentration in this area, if the stromal bed becomes too thin or the level of stress is too high, tissue damage may occur.The same analysis was performed for two other values of IOP: 10 mmHg and 21 mmHg.In both cases, the initial stress distribution was previously introduced.For the three values of IOP considered, 10, 15 and 21 mmHg, the correction estimated using a small central area radius (1.5 mm) matches the expected result better than a larger radius of 2.5 mm.For the highest value of the IOP (21 mmHg) using a central area radius of 1.5 mm, an under-correction of the myopia occurs in the 8 and 10D cases.For the simulation using an IOP of 10 mmHg, the maximum principal stress at the center of the stromal bed were 1.23 × 10 −2 , 1.29 × 10 −2 , 1.38 × 10 −2 , 1.49 × 10 −2

472
The Effect of Intraocular Pressure on the Outcome of Myopic Photorefractive Keratectomy: A Numerical Approach

DISCUSSION
In this work, the stress distribution due to IOP was introduced in the model to reproduce the real in vivo conditions of the eye before PRK surgery.If IOP is not considered, the model does not provide the expected refractive correction and does not predict the behavior of the corneal tissue.Therefore, the stress due to IOP must be taken into account in performing any simulation.
In spite of using a hyperelastic constitutive model, the observed response of the corneal tissue in the considered range of IOP is almost linear [11], although anisotropic, with the stress level about 10 -2 MPa clearly falling in the quasi-linear region of the stress-stretch curve.The same behaviour was observed by Studer et al. for IOP up to 20 mmHg [37].
For the highest value of IOP considered (21 mmHg) in this study, a high undecorrection was obtained compared with the results predicted by the standard Munnerlyn's equation (eqn.4).This indicates that IOP does have an effect on the postsurgical refractive outcome.Therefore, in addition to the needed correction and the central corneal thickness, IOP measure should be considered to determine the depth of the ablation.IOP affects not only the final curvature of the cornea but also the initial stress of the tissue [12].That implies that the surgery applied to different eyes acts under different conditions, on tissues more or less stressed, which modifies the final deformed configuration and with this, the refractive outcome of the surgery.
Previous studies to simulate myopic PRK surgery presented different results.The results of Alastrué et al. [14] demonstrated undercorrections for all cases (from 2D to 10D), as it did not take into account either the initial stresses or the variation of IOP.Pandolfi [12] obtained a slight dependence of the myopic correction on the IOP, showing an undercorrection for IOP higher than 15 mmHg, in agreement with our results.
Simulation of ablative surgery also provides warning concerning an important clinical issue: the ablation should leave a sufficient pachymetric thickness at the central cornea to preserve a good stromal bed after surgery.Otherwise, the patient must not receive such operation.The stress concentration generated at the central cornea after laser ablation can cause ectasia, haze or other complications [12].Our results quantify the anticipated qualitative result that the magnitude of stress increases with IOP and the level of myopia to be corrected.
The model used in this work has some limitations.A more realistic geometry, based on a non-axisymmetric ellipsoid [38], together with the addition of other components which may influence the structural response of the cornea, such as the muscles or iris, would more accurately simulate the behaviour of the eye.The procedure could be improved by the estimation of the depth of ablation using the algorithms which govern the excimer laser equipments instead of the single Munnerlyn's expression [32].Moreover, the refractive outcome should be analyzed by optical methods to estimate the aberrations induced by the surgery.A validation of the model will be then required to guarantee reliable results.The validation would be established when the computational results of a number of cases reproduce post-surgery clinical outcomes in patients.
Despite those limitations, biomechanics is becoming a useful tool for surgical planning [12,18] allowing improving and/or optimizing current surgical procedures, or even developing new technologies.The methodology presented herein could be applied after complete validation to estimate the parameters of ablation for patient-specific cases using corneal thickness taken by the pachymeter and the value of IOP measured by the tonometer.In the near future, it is expected that the in vivo corneal material properties of the patient could be estimated and used to create a fully-customized model that could be implemented in the software of the excimer laser equipments.

CONCLUSION
The effect of IOP on the refractive correction achieved by the PRK surgery has been analyzed using a numerical model in the present work.The initial stress distribution due to IOP was introduced to the model to reproduce the real in vivo conditions of the eye before surgery.For a mean IOP value of 15 mmHg, an over-correction was observed for the cases of lower values of correction (2 and 4D), while a good estimation of myopia correction was achieved in the cases of higher correction (6, 8 and 10D).These results were also consistent with those estimated by the Munnerlyn's expression only if the initial stress distribution was taken into account.When the initial stress distribution was not considered, the diopters predicted by the present model were strongly dependent on the central area radius, due to the excessive deformation of the stress-free cornea.The present results suggest that the IOP should be considered to estimate the depth of ablation for PRK procedure.

Figure 1 .
Figure 1.(b) Perspective view of a quarter model of the cornea showing the five element layers of the outer part, each corresponding to 2-diopter depth of ablation.For a myopic correction of 2D, the outer layer is removed.For a 10D correction, the five layers are removed.
Figure 2.Two families of collagen fibers are present in the cornea, corresponding to the following preferred directions: superior-inferior (n 0 ) and nasal-temporal (m 0 ).The preferred orientation of the collagen fibers in the limbal tissue is circumferential (n 0 ).
Journal of Healthcare Engineering • Vol.

Figure 4 .Figure 5 .Table 3 .
Figure 4. Displacement distribution under 15 mmHg of IOP before simulation, with (left) and without (right) introducing the initial stress distribution in the eye model.All values in mm.

Figure 8
presents the obtained myopia correction.Journal of Healthcare Engineering • Vol. 1 • No. 3

Figure 6 .
Figure 6.Change in diopters of the laser refractive surgery simulation with (left) and without (right) introducing initial stresses in the model (IOP = 15 mmHg).Circular markers ( ) represent the expected correction; rhombus-shaped ( ) and square ( ) markers represent the estimated outcome for a central area radius of 1.5 and 2.5 mm, respectively.

Figure 8 .
Figure 8. Change in diopters of the laser refractive surgery simulation for two values of IOP: 10 mmHg (left) and 21 mmHg (right).Circular markers (᭺) represent the expected correction; rhombus-shaped (᭛) and square (ᮀ) markers represent the estimated outcome for a central area radius of 1.5 and 2.5 mm, respectively.

Figure 7 .
Figure 7. Transversal view of the Finite Element model of the cornea showing the maximum principal stress distribution (MPa) after a 2D-ablation (left) and 10D-ablation (right), IOP = 15 mmHg.