- Email: [email protected]
Stress free configuration of the human eye
Medical Engineering & Physics 35 (2013) 211–216
Contents lists available at SciVerse ScienceDirect
Medical Engineering & Physics journal homepage: www.elsevier.com/locate/medengphy
Technical Note
Stress free configuration of the human eye Ahmed Elsheikh a, Charles Whitford a, Rosti Hamarashid a, Wael Kassem b, Akram Joda a,∗, Philippe Büchler c a b c
School of Engineering, University of Liverpool, L69 3GH Liverpool, UK Division of Civil Engineering, University of El-Minia, El-Minia 61111, Egypt Institute for Surgical Technology & Biomechanics, University of Bern, Stauffacherstrasse 78, CH-3014 Bern, Switzerland
a r t i c l e
i n f o
Article history: Received 12 June 2012 Received in revised form 29 August 2012 Accepted 10 September 2012 Keywords: Ocular biomechanics Cornea Sclera Numerical simulation Unstressed-form
a b s t r a c t Numerical simulations of eye globes often rely on topographies that have been measured in vivo using devices such as the Pentacam or OCT. The topographies, which represent the form of the already stressed eye under the existing intraocular pressure, introduce approximations in the analysis. The accuracy of the simulations could be improved if either the stress state of the eye under the effect of intraocular pressure is determined, or the stress-free form of the eye estimated prior to conducting the analysis. This study reviews earlier attempts to address this problem and assesses the performance of an iterative technique proposed by Pandolfi and Holzapfel [1], which is both simple to implement and promises high accuracy in estimating the eye’s stress-free form. A parametric study has been conducted and demonstrated reliance of the error level on the level of flexibility of the eye model, especially in the cornea region. However, in all cases considered 3–4 analysis iterations were sufficient to produce a stress-free form with average errors in node location <10−6 mm and a maximal error <10−4 mm. This error level, which is similar to what has been achieved with other methods and orders of magnitude lower than the accuracy of current clinical topography systems, justifies the use of the technique as a pre-processing step in ocular numerical simulations. Crown Copyright © 2012 Published by Elsevier Ltd on behalf of IPEM. All rights reserved.
1. Introduction The human eye is a complex optical system whose threedimensional shape determines the clarity of vision. Light rays entering the eye are refracted four times before they reach the retina; twice at the two surfaces of the cornea, which accounts for about two-thirds of the refractive power, and twice at the surfaces of the crystalline lens. Small geometric changes such as imperfections of the cornea or elongation of the sclera may lead to a blurry and/or distorted vision. Nowadays, several surgical interventions induce changes in corneal [2] and scleral [3] shapes in order to improve vision. In conjunction with the advancement of surgical techniques, numerical models have been proposed as tools to better understand the effect of surgeries on the eye [4–8]. To aid these efforts, experimental investigations have been conducted to characterize the mechanical properties of the cornea [9] and sclera [10], and quantify their non-linear, anisotropic and visco-elastic behavior [11,12]. The 3D ocular topography can also be measured in vivo with several devices such as the Pentacam, Gallilei or OCT. However, since
∗ Corresponding author. Tel.: +44 1517944908. E-mail address: [email protected] (A. Joda).
the eye supports intraocular pressure (IOP), the shape measured experimentally corresponds to a deformed configuration with the tissue under primarily membrane tension. This deformed configuration is unsuitable for direct implementation in patient-specific models as the application of IOP on the models would result in a larger shape than that obtained experimentally. If on the other hand the IOP is not applied, its associated stresses will not be generated in the model, which also affects the model’s representation of in vivo conditions. For these reasons, it is necessary for accurate numerical modeling to determine the relaxed (stress-free) configuration of the tissue, or the stable eye geometry that would be adopted had the IOP been removed. Two main techniques have been proposed to solve this problem. Some studies adopted an inverse elastostatic approach based on the reverse application of IOP to explicitly estimate the stress-free configuration [13–15]. This was followed by the application of IOP on the estimated stress-free form to determine the resulting stresses in the deformed configuration through the solution of a forward problem. One of the limitations of this approach is that the thinwalled flexible ocular structure may exhibit bifurcation during the backward calculation of the zero pressure geometry [16]. To solve this problem, Gee et al. proposed a modified updated Lagrangian formulation to estimate the stress in the tissue using a forward calculation [17], and a similar technique was used by Lanchares et al. to
1350-4533/$ – see front matter. Crown Copyright © 2012 Published by Elsevier Ltd on behalf of IPEM. All rights reserved. http://dx.doi.org/10.1016/j.medengphy.2012.09.006
212
A. Elsheikh et al. / Medical Engineering & Physics 35 (2013) 211–216
pre-stress a cornea model to simulate arcuate keratotomy surgical intervention [18]. This approach, which attempts to avoid the need to determine the stress-free configuration, was later refined by Grytz and Downs to improve its accuracy [16]. The main difference between the two methods is that in inverse elastostatic approach, a complete simulation step is required to update the deformation gradient tensor while the framework proposed by Grytz and Downs continuously updates the deformation tensor during the application of the load on the eye model. Although these methods have been shown to provide accurate results, their implementations remain complicated and require a thorough understanding of continuum mechanics and finite element analysis. An alternative approach has been proposed by Pandolfi and Holzapfel [1], where general-purpose finite element software packages are combined with simple mathematical calculations to obtain stress-free configurations. The objective of this study is to assess this approach and determine if it provides a level of accuracy sufficient for ophthalmic applications. 2. Methods The method used in this study is based on the simulation of ocular deformation with finite element techniques. An iterative approach is used to gradually move the nodes of the finite element mesh to their stress-free configuration. 2.1. Calculation of the stress-free configuration An iterative approach based on work of Pandolfi and Holzapfel [1] is used to obtain the stress-free configuration. An initial numerical model is built based on a measured (stressed) shape of the tissue X0 . In a first step, the displacements u1 induced by the intraocular pressure (IOP) are calculated using non-linear finite element analysis. These displacements vectors are then subtracted from the node coordinates of the stressed shape to determine the first estimate of the stress-free form X1 , Fig. 1. u1 = X0 − x0
(1)
X1 = X0 − u1
(2)
In a second step, the model with the initial stress-free form X1 is re-inflated with IOP and the error vector u2 between the resulting (deformed) nodal positions x1 and the target nodal positions X0 (corresponding to the measured stressed configuration) is calculated at all nodes. The error vector is then subtracted from the node coordinates of the current assumed stress-free form X1 to define the initial configuration X2 for the following analysis iteration. u2 = X0 − x1
(3)
X2 = X1 − u2
(4)
This process is repeated while monitoring the values and distribution of nodal errors relative to the target (measured) configuration. uk = X0 − xk−1
(5)
Xk = Xk−1 − uk
(6)
The absolute nodal distance and the root mean square distance (L2 -norm) are used to quantify the maximum and average error and to monitor the convergence of the iterative process toward the stress-free form. The error estimate of the kth increment in then given by: ek = ||X0 − xk ||
(7)
2.2. Finite element model of the eye In order to evaluate the accuracy of the proposed approach, a parametric model of the human eye was developed. The model incorporated a number of features to closely represent in vivo conditions including the cornea’s and sclera’s non-uniform thickness, asphericity of both the cornea’s anterior and posterior surfaces, and weak inter-lamellar adhesion within corneal stroma. The finite element (FE) software ABAQUS 6.10 (Dassault Systèmes Simulia Corp., Rhode Island, USA) was used in the analyses. The model used sixnoded solid elements (C3D6H) to enable a good representation of the variable thickness. The elements were arranged in either 1, 2 or 4 layers and 120 rings, 20 of which were in the cornea, to provide smooth distributions of strains and displacements, see Fig. 2. The total elements numbers were 6962, 13,924 or 27,848 for the case with 1, 2 and 4 layers, respectively. Weak stromal inter-lamellar adhesion was assumed only in the cornea and in models with 2 or more layers, and was set to the level observed experimentally [19]. This feature could not be considered in the model with only one element layer. The models incorporated a material definition that combined both the hyperelasticity and hysteresis of the tissue as has been detailed in Elsheikh et al. [9,10,20]. With these two characteristics, the material experienced a gradually increasing stiffness under load, and should the material be unloaded, it would follow a different stress–strain path [20]. With both behavior patterns described in the material definition, the simulation would be able to select which stress–strain path each element was to follow based on its strain history. Moreover, varying the age in the model meant changing the constitutive model parameters of the cornea and sclera based on the relationships provided by Elsheikh et al. [9,10]. The constitutive model used to represent corneal and scleral mechanical behavior for all ages was the third order Ogden hyperelastic model. Further, the models were provided with support conditions designed to avoid interference with the behavior under IOP. The models were prevented from displacement in the anterior–posterior direction (z-direction) along all equator nodes and displacement in the superior–inferior and temporal–nasal directions at the corneal apex and posterior pole. Only a quarter of the eye was used in this work, for that the symmetry boundary conditions were applied at the two ends (Fig. 2).
2.3. Parametric study The accuracy of the method has been evaluated in a parametric study that considered eye numerical models with variations in corneal anterior radius, R, shape factor, p, sclera radius, Rs , intraocular pressure, IOP, and age. The parameter variations shown in Table 1 have been chosen to cover, and in some cases go beyond, the variability present in the general population. A reference benchmark model, with the parameter baseline values in gray cells in Table 1, was considered. Within each part of the study, only one parameter was allowed to vary while other parameters maintained their baseline values as in the benchmark model. All models considered a central corneal thickness, CCT, of 545 m – the average reported value in a number of earlier studies [21] – and a peripheral corneal thickness, PCT, that was larger than CCT by 150 m in line with Gullstrand’s No 1 schematic eye [22]. While the cornea was considered to be an ellipsoid with variable central radius and shape factor, the sclera assumed a spherical outer surface and a thickness that changed from PCT at the limbus, reducing gradually to 0.8 PCT at the equator before increasing to 1.2 PCT at the posterior pole, based on the results of an earlier study [10].
A. Elsheikh et al. / Medical Engineering & Physics 35 (2013) 211–216
213
Fig. 1. Schematic description of the approach used to obtain the stress free configuration of the eye.
3. Results The effects of various parameters on the accuracy of the resulting stress-free configuration of the eye have been evaluated. The parameters included the effect of boundary conditions (IOP pressure and number of model layers), geometric modeling (corneal
radius, shape factor and scleral radius) and mechanical properties, which vary with patient’s age. With successive iterations, the analysis predictions indicated fast decrease of the error between the measured target shape and the topography obtained after loading the stress-free configuration. The error decreased exponentially during the first six
Fig. 2. (a) Schematic of the model with the boundary conditions and (b)–(d) cross-sectional views of numerical models with 1, 2 and 4 element layers, respectively.
214
A. Elsheikh et al. / Medical Engineering & Physics 35 (2013) 211–216
Table 1 Morphomechanical values used in the parametric studies. Shaded cells point at values adopted in benchmark case. Each study considered variation in one parameter and kept other parameters constant at the benchmark values. Intraocular pressure (mmHg) Age (years)
10 30
15 40
Corneal anterior radius (mm) Corneal shape factor Number of model layers Sclera external diameter (mm)
7.0 0.74 1 21
7.4 0.78 2 23
20 50 7.8 0.82
25 60 8.2 0.86 4
30 70 8.6 0.90
80 0.94
90 0.98
1.02
1.06
1.1
25
iteration cycles, beyond which the rate of improvement of the shape decreased and appeared to diminish after about eight iterative cycles (Fig. 3). Following 8–12 iterations, the maximal error in individual models reduced to about 10−3 m while corresponding average model errors were approximately 10−6 m. The error distribution around the eye was not uniform, and higher errors were measured within the cornea and around the limbus than on the remaining sections of the eye (Fig. 4). In most models, the maximum errors within the cornea and limbus were approximately 7–15 times the errors within the sclera, which appeared to experience little variation from one location to another.
Variation of IOP had a significant impact on both the maximum and average errors with models under lower IOP producing highly accurate approximations of the stress-free configuration with less iterations than models with larger IOP (Fig. 3a). Further, the models with a number of layers greater than one converged faster to the stress-free configuration than the model with only one layer of elements (Fig. 3b). Having two element layers or more enabled the introduction of the weak-interlamellar adhesion between corneal layers, which reduced the flexural stiffness of the cornea and led to a relaxation of stresses in the tissue and concentrated deformations in the region of the limbus. This mechanical effect could be behind
Fig. 3. Variation of average and maximum errors with progression of analysis iterations in cases with (a) different values of IOP, (b) different numbers of model layers, (c) different ages, and (d) different central anterior radius of curvature. Heavy and light lines represent maximum and average errors, respectively.
A. Elsheikh et al. / Medical Engineering & Physics 35 (2013) 211–216
215
Fig. 5. Reduction of error in the refractive power of the reference numerical model with progression of iterative cycles. Fig. 4. Variation of nodal error along a meridional section of the eye obtained with the reference eye model and following the first iteration.
the variations observed in the convergence toward the stress-free configuration between the model with only one layer of elements and those with two layers or more. The mechanical properties of ocular tissue were allowed to vary in the study through their association with age. The results showed that stiffer mechanical behavior (with higher age) was associated with improved convergence of the model. Only three iteration cycles were needed to achieve the same accuracy with the 90-year-old properties as with four iterations in models with the other age groups (Fig. 3c). On the other hand, the geometric parameters showed limited effects on the convergence to the stress-free configuration. The results demonstrated that increased scleral radius or corneal shape factor did not cause any significant change in model convergence (results not shown). However, the convergence rate was slightly altered by the corneal radius, with larger radii producing slower convergence (Fig. 3d). The results show consistently that models with a high level of flexibility, especially in the corneal region, undergo slow convergence to the stress-free configuration. This was evident in the models that considered the weak interlamellar adhesion in corneal stroma, models with soft material behavior associated with young age, and models with larger corneal radius. The resulting slow convergence in these cases could be explained by the larger topography adjustments, caused by the larger model displacements under IOP, which needed to be considered to obtain the stress-free configuration. Furthermore, the parametric study demonstrated that four iterative cycles were sufficient to determine the stress-free configuration of the human eye with average and maximum errors below 10−6 mm and 10−4 mm, respectively. This level of inaccuracy induces changes in the patients’ visual acuity of about 10−4 diopter (Fig. 5), which is orders of magnitude lower than the changes that the human optical system is able to perceive. 4. Discussion The geometry of biological structures is usually only available in a stressed state under the physiological loading of the tissues. In order to provide accurate predictions of mechanical behavior, numerical models should account for this physiological deformation. In this study, we present a simple iterative method
to determine the stress-free configuration of biological structures and evaluate the method within a parametric study of the human eye. The parametric study, which covered a wide range of ocular shape parameters (cornea radius and shape factor, and sclera radius), boundary conditions (IOP), mechanical properties (controlled by age) and finite element mesh construction (with 1, 2 and 4 element layers), demonstrated similar rates of convergence toward the stress-free configuration in all cases. However, the results also indicate that the stiffness of the model and the IOP are associated with a slight effect on the convergence behavior. Fewer iteration cycles are needed to find the stress-free configuration if the distance between the target configuration and the initial (stressed) configuration is small. This effect explains the slow convergence in cases with higher internal pressure and softer material, and the higher errors observed in the corneal region (with lower stiffness relative to the sclera) compared to other regions of the eye. In all cases studied, the average error in determining the stress-free configuration using the iterative method is lower than 10−6 mm and the maximal error less than 10−4 mm, which is orders of magnitude lower than the accuracy of current clinical topographic systems. This level of accuracy is achievable within 3–4 cycles of analysis iterations, and there would be no added benefit in increasing the number of cycles since the target (stressed) configuration is in any case not perfectly determined experimentally. The accuracy of the results reported by Grytz and Dows [16] for an iterative forward prestressing method is higher than the accuracy achieved with the proposed approach. These authors report a maximal error decreasing to less than 10−11 m after about 10 iterative cycles, but the accuracy achieved after four iterations of the forward prestressing approach is similar to that obtained in the present work. The numerical model used in this study has a number of limitations concerning the mechanical model of the tissue. The well-known anisotropic distribution of the collagen fibers within the tissue was not taken into account in the constitutive equations. This was necessary since the age-dependent mechanical properties of the tissue (which were deemed more influential on model deformation than anisotropy) had not been evaluated using anisotropic mechanical models. The effect of omitting corneal anisotropy and local changes of the curvature on the
216
A. Elsheikh et al. / Medical Engineering & Physics 35 (2013) 211–216
efficiency of the model in acquiring the stress free form are difficult to quantify, although they are not expected to be significant. Although the model adopted the thickness variation in both the cornea and sclera as reported in the literature [9,10], the thickness variation around the optic nerve head was not considered as it was thought not to affect the assessment of the proposed method to acquire the stress-free form. Other ocular characterization that have not been considered in the model for being unlikely to produce a significant change in model performance include limbal ellipticity, difference between sagittal and tangential radii of curvature and the regional variation of mechanical properties across the sclera. Another limitation of the study concerns the idealized model of the eye, which overlooked the local change in curvature of patient-specific data. This approach has been adopted in order to test the proposed algorithm within a controlled parametric environment and enable the evaluation of the individual effects of each mechanical and geometrical parameter. 5. Conclusions In this study, a generic method has been proposed to determine the stress-free configuration of the human eye. The method does not require any modification of the finite element formulation and could be used with any type of material model. The method has been shown to converge quickly to a level of accuracy suitable for biomechanical analysis. Although the initial application was ocular biomechanics, this approach could be directly applied to other biological structures such as ligaments or blood vessels. Funding This work was funded by EPSRC grant EP/H052046/1. Ethical approval Not required. Conflict of interest statement The authors have no financial or personal interests in the work reported in this paper.
References [1] Pandolfi A, Holzapfel GA. Three-dimensional modeling and computational analysis of the human cornea considering distributed collagen fibril orientations. J Biomech Eng 2008;130(6):061006. [2] Munnerlyn C, Koons SJ, Mearshall J. Photorefractive qeratectomy: a technique for laser refractive surgery. J Cataract Refract Surg 1988;14:46–52. [3] Baino F. Scleral buckling biomaterials and implants for retinal detachment surgery. Med Eng Phys 2010;32(9):945–56. [4] Roy AS, Dupps WJ. Patient-specific modeling of corneal refractive surgery outcomes and inverse estimation of elastic property changes. J Biomech Eng 2011;133(1):011002. [5] Alastrue V, Calvo B, Pena E, Doblare M. Biomechanical modeling of refractive corneal surgery. J Biomech Eng 2006;128(1):150–60. [6] Villamarin A, Roy S, Hasballa R, Vardoulis O, Reymond P, Stergiopulos N. 3D simulation of the aqueous flow in the human eye. Med Eng Phys; 2012, http://dx.doi.org/10.1016/j.medengphy.2012.02.007 [7] Ganesan B, He S, Xu H. Modelling of pulsatile blood flow in arterial trees of retinal vasculature. Med Eng Phys 2011;33:810–23. [8] Pandolfi A, Fotia G, Manganiello F. Finite element simulations of laser refractive corneal surgery. Eng Comput 2009;25:15–24. [9] Elsheikh A, Geraghty B, Rama P, Campanelli M, Meek KM. Characterization of age-related variation in corneal biomechanical properties. J Roy Soc Interface 2010;7:1475–85. [10] Elsheikh A, Geraghty B, Alhasso D, Knappett J, Campanelli M, Rama P. Regional variation in the biomechanical properties of the human sclera. Exp Eye Res 2010;90:624–33. [11] Elsheikh A, Brown M, Alhasso D, Rama P, Campanelli M, GarwayHeath D. Experimental assessment of corneal anisotropy. J Refract Surg 2008;24(2):178–87. [12] Elsheikh A, Kassem W, Jones SW. Strain-rate sensitivity of porcine and ovine corneas. Acta Bioeng Biomech 2011;13(2):25–36. [13] Govindjee S, Mihalic PA. Computational methods for inverse deformations in quasi-incompressible finite elasticity. Int J Numer Methods Eng 1998;43(5):821–38. [14] Govindjee S, Mihalic PA. Computational methods for inverse finite elastostatics. Comput Methods Appl Mech Eng 1996;136(1–2):47–57. [15] Lu J, Zhou X, Raghavan ML. Inverse elastostatic stress analysis in pre-deformed biological structures: demonstration using abdominal aortic aneurysms. J Biomech 2007;40(3):693–6. [16] Grytz R, Downs JC. A forward incremental prestressing method with application to inverse parameter estimations and eye-specific simulations of posterior scleral shells. Comput Methods Biomech Biomed Eng 2012 [Epub ahead of print]. [17] Gee MW, Förster C, Wall WA. A computational strategy for prestressing patient-specific biomechanical problems under finite deformation. Int J Numer Methods Biomed Eng 2010;26(1):52–72. [18] Lanchares E, Calvo B, Cristóbal JA, Doblaré M. Finite element simulation of arcuates for astigmatism correction. J Biomech 2008;41(4):797–805. [19] Elsheikh A, Ross S, Alhasso D, Rama P. Numerical study of the effect of corneal layered structure on ocular biomechanics. Curr Eye Res 2009;34(1):26–35. [20] Elsheikh A, Wang D, Rama P, Campanelli M, Garway-Heath D. Experimental assessment of human corneal hysteresis. Curr Eye Res 2008;33(3):205–13. [21] Francis BA, Hsieh A, Lai MY, Chopra V, Pena F, Azen S, et al. Effects of corneal thickness, corneal curvature, and intraocular pressure level on Goldmann applanation tonometry and dynamic contour tonometry. Ophthalmology 2007;114(1):20–6. [22] Bennett G, Rabbetts RB. Clinical visual optics, vol. 2. London: Butterworths; 1989.
Contents lists available at SciVerse ScienceDirect
Medical Engineering & Physics journal homepage: www.elsevier.com/locate/medengphy
Technical Note
Stress free configuration of the human eye Ahmed Elsheikh a, Charles Whitford a, Rosti Hamarashid a, Wael Kassem b, Akram Joda a,∗, Philippe Büchler c a b c
School of Engineering, University of Liverpool, L69 3GH Liverpool, UK Division of Civil Engineering, University of El-Minia, El-Minia 61111, Egypt Institute for Surgical Technology & Biomechanics, University of Bern, Stauffacherstrasse 78, CH-3014 Bern, Switzerland
a r t i c l e
i n f o
Article history: Received 12 June 2012 Received in revised form 29 August 2012 Accepted 10 September 2012 Keywords: Ocular biomechanics Cornea Sclera Numerical simulation Unstressed-form
a b s t r a c t Numerical simulations of eye globes often rely on topographies that have been measured in vivo using devices such as the Pentacam or OCT. The topographies, which represent the form of the already stressed eye under the existing intraocular pressure, introduce approximations in the analysis. The accuracy of the simulations could be improved if either the stress state of the eye under the effect of intraocular pressure is determined, or the stress-free form of the eye estimated prior to conducting the analysis. This study reviews earlier attempts to address this problem and assesses the performance of an iterative technique proposed by Pandolfi and Holzapfel [1], which is both simple to implement and promises high accuracy in estimating the eye’s stress-free form. A parametric study has been conducted and demonstrated reliance of the error level on the level of flexibility of the eye model, especially in the cornea region. However, in all cases considered 3–4 analysis iterations were sufficient to produce a stress-free form with average errors in node location <10−6 mm and a maximal error <10−4 mm. This error level, which is similar to what has been achieved with other methods and orders of magnitude lower than the accuracy of current clinical topography systems, justifies the use of the technique as a pre-processing step in ocular numerical simulations. Crown Copyright © 2012 Published by Elsevier Ltd on behalf of IPEM. All rights reserved.
1. Introduction The human eye is a complex optical system whose threedimensional shape determines the clarity of vision. Light rays entering the eye are refracted four times before they reach the retina; twice at the two surfaces of the cornea, which accounts for about two-thirds of the refractive power, and twice at the surfaces of the crystalline lens. Small geometric changes such as imperfections of the cornea or elongation of the sclera may lead to a blurry and/or distorted vision. Nowadays, several surgical interventions induce changes in corneal [2] and scleral [3] shapes in order to improve vision. In conjunction with the advancement of surgical techniques, numerical models have been proposed as tools to better understand the effect of surgeries on the eye [4–8]. To aid these efforts, experimental investigations have been conducted to characterize the mechanical properties of the cornea [9] and sclera [10], and quantify their non-linear, anisotropic and visco-elastic behavior [11,12]. The 3D ocular topography can also be measured in vivo with several devices such as the Pentacam, Gallilei or OCT. However, since
∗ Corresponding author. Tel.: +44 1517944908. E-mail address: [email protected] (A. Joda).
the eye supports intraocular pressure (IOP), the shape measured experimentally corresponds to a deformed configuration with the tissue under primarily membrane tension. This deformed configuration is unsuitable for direct implementation in patient-specific models as the application of IOP on the models would result in a larger shape than that obtained experimentally. If on the other hand the IOP is not applied, its associated stresses will not be generated in the model, which also affects the model’s representation of in vivo conditions. For these reasons, it is necessary for accurate numerical modeling to determine the relaxed (stress-free) configuration of the tissue, or the stable eye geometry that would be adopted had the IOP been removed. Two main techniques have been proposed to solve this problem. Some studies adopted an inverse elastostatic approach based on the reverse application of IOP to explicitly estimate the stress-free configuration [13–15]. This was followed by the application of IOP on the estimated stress-free form to determine the resulting stresses in the deformed configuration through the solution of a forward problem. One of the limitations of this approach is that the thinwalled flexible ocular structure may exhibit bifurcation during the backward calculation of the zero pressure geometry [16]. To solve this problem, Gee et al. proposed a modified updated Lagrangian formulation to estimate the stress in the tissue using a forward calculation [17], and a similar technique was used by Lanchares et al. to
1350-4533/$ – see front matter. Crown Copyright © 2012 Published by Elsevier Ltd on behalf of IPEM. All rights reserved. http://dx.doi.org/10.1016/j.medengphy.2012.09.006
212
A. Elsheikh et al. / Medical Engineering & Physics 35 (2013) 211–216
pre-stress a cornea model to simulate arcuate keratotomy surgical intervention [18]. This approach, which attempts to avoid the need to determine the stress-free configuration, was later refined by Grytz and Downs to improve its accuracy [16]. The main difference between the two methods is that in inverse elastostatic approach, a complete simulation step is required to update the deformation gradient tensor while the framework proposed by Grytz and Downs continuously updates the deformation tensor during the application of the load on the eye model. Although these methods have been shown to provide accurate results, their implementations remain complicated and require a thorough understanding of continuum mechanics and finite element analysis. An alternative approach has been proposed by Pandolfi and Holzapfel [1], where general-purpose finite element software packages are combined with simple mathematical calculations to obtain stress-free configurations. The objective of this study is to assess this approach and determine if it provides a level of accuracy sufficient for ophthalmic applications. 2. Methods The method used in this study is based on the simulation of ocular deformation with finite element techniques. An iterative approach is used to gradually move the nodes of the finite element mesh to their stress-free configuration. 2.1. Calculation of the stress-free configuration An iterative approach based on work of Pandolfi and Holzapfel [1] is used to obtain the stress-free configuration. An initial numerical model is built based on a measured (stressed) shape of the tissue X0 . In a first step, the displacements u1 induced by the intraocular pressure (IOP) are calculated using non-linear finite element analysis. These displacements vectors are then subtracted from the node coordinates of the stressed shape to determine the first estimate of the stress-free form X1 , Fig. 1. u1 = X0 − x0
(1)
X1 = X0 − u1
(2)
In a second step, the model with the initial stress-free form X1 is re-inflated with IOP and the error vector u2 between the resulting (deformed) nodal positions x1 and the target nodal positions X0 (corresponding to the measured stressed configuration) is calculated at all nodes. The error vector is then subtracted from the node coordinates of the current assumed stress-free form X1 to define the initial configuration X2 for the following analysis iteration. u2 = X0 − x1
(3)
X2 = X1 − u2
(4)
This process is repeated while monitoring the values and distribution of nodal errors relative to the target (measured) configuration. uk = X0 − xk−1
(5)
Xk = Xk−1 − uk
(6)
The absolute nodal distance and the root mean square distance (L2 -norm) are used to quantify the maximum and average error and to monitor the convergence of the iterative process toward the stress-free form. The error estimate of the kth increment in then given by: ek = ||X0 − xk ||
(7)
2.2. Finite element model of the eye In order to evaluate the accuracy of the proposed approach, a parametric model of the human eye was developed. The model incorporated a number of features to closely represent in vivo conditions including the cornea’s and sclera’s non-uniform thickness, asphericity of both the cornea’s anterior and posterior surfaces, and weak inter-lamellar adhesion within corneal stroma. The finite element (FE) software ABAQUS 6.10 (Dassault Systèmes Simulia Corp., Rhode Island, USA) was used in the analyses. The model used sixnoded solid elements (C3D6H) to enable a good representation of the variable thickness. The elements were arranged in either 1, 2 or 4 layers and 120 rings, 20 of which were in the cornea, to provide smooth distributions of strains and displacements, see Fig. 2. The total elements numbers were 6962, 13,924 or 27,848 for the case with 1, 2 and 4 layers, respectively. Weak stromal inter-lamellar adhesion was assumed only in the cornea and in models with 2 or more layers, and was set to the level observed experimentally [19]. This feature could not be considered in the model with only one element layer. The models incorporated a material definition that combined both the hyperelasticity and hysteresis of the tissue as has been detailed in Elsheikh et al. [9,10,20]. With these two characteristics, the material experienced a gradually increasing stiffness under load, and should the material be unloaded, it would follow a different stress–strain path [20]. With both behavior patterns described in the material definition, the simulation would be able to select which stress–strain path each element was to follow based on its strain history. Moreover, varying the age in the model meant changing the constitutive model parameters of the cornea and sclera based on the relationships provided by Elsheikh et al. [9,10]. The constitutive model used to represent corneal and scleral mechanical behavior for all ages was the third order Ogden hyperelastic model. Further, the models were provided with support conditions designed to avoid interference with the behavior under IOP. The models were prevented from displacement in the anterior–posterior direction (z-direction) along all equator nodes and displacement in the superior–inferior and temporal–nasal directions at the corneal apex and posterior pole. Only a quarter of the eye was used in this work, for that the symmetry boundary conditions were applied at the two ends (Fig. 2).
2.3. Parametric study The accuracy of the method has been evaluated in a parametric study that considered eye numerical models with variations in corneal anterior radius, R, shape factor, p, sclera radius, Rs , intraocular pressure, IOP, and age. The parameter variations shown in Table 1 have been chosen to cover, and in some cases go beyond, the variability present in the general population. A reference benchmark model, with the parameter baseline values in gray cells in Table 1, was considered. Within each part of the study, only one parameter was allowed to vary while other parameters maintained their baseline values as in the benchmark model. All models considered a central corneal thickness, CCT, of 545 m – the average reported value in a number of earlier studies [21] – and a peripheral corneal thickness, PCT, that was larger than CCT by 150 m in line with Gullstrand’s No 1 schematic eye [22]. While the cornea was considered to be an ellipsoid with variable central radius and shape factor, the sclera assumed a spherical outer surface and a thickness that changed from PCT at the limbus, reducing gradually to 0.8 PCT at the equator before increasing to 1.2 PCT at the posterior pole, based on the results of an earlier study [10].
A. Elsheikh et al. / Medical Engineering & Physics 35 (2013) 211–216
213
Fig. 1. Schematic description of the approach used to obtain the stress free configuration of the eye.
3. Results The effects of various parameters on the accuracy of the resulting stress-free configuration of the eye have been evaluated. The parameters included the effect of boundary conditions (IOP pressure and number of model layers), geometric modeling (corneal
radius, shape factor and scleral radius) and mechanical properties, which vary with patient’s age. With successive iterations, the analysis predictions indicated fast decrease of the error between the measured target shape and the topography obtained after loading the stress-free configuration. The error decreased exponentially during the first six
Fig. 2. (a) Schematic of the model with the boundary conditions and (b)–(d) cross-sectional views of numerical models with 1, 2 and 4 element layers, respectively.
214
A. Elsheikh et al. / Medical Engineering & Physics 35 (2013) 211–216
Table 1 Morphomechanical values used in the parametric studies. Shaded cells point at values adopted in benchmark case. Each study considered variation in one parameter and kept other parameters constant at the benchmark values. Intraocular pressure (mmHg) Age (years)
10 30
15 40
Corneal anterior radius (mm) Corneal shape factor Number of model layers Sclera external diameter (mm)
7.0 0.74 1 21
7.4 0.78 2 23
20 50 7.8 0.82
25 60 8.2 0.86 4
30 70 8.6 0.90
80 0.94
90 0.98
1.02
1.06
1.1
25
iteration cycles, beyond which the rate of improvement of the shape decreased and appeared to diminish after about eight iterative cycles (Fig. 3). Following 8–12 iterations, the maximal error in individual models reduced to about 10−3 m while corresponding average model errors were approximately 10−6 m. The error distribution around the eye was not uniform, and higher errors were measured within the cornea and around the limbus than on the remaining sections of the eye (Fig. 4). In most models, the maximum errors within the cornea and limbus were approximately 7–15 times the errors within the sclera, which appeared to experience little variation from one location to another.
Variation of IOP had a significant impact on both the maximum and average errors with models under lower IOP producing highly accurate approximations of the stress-free configuration with less iterations than models with larger IOP (Fig. 3a). Further, the models with a number of layers greater than one converged faster to the stress-free configuration than the model with only one layer of elements (Fig. 3b). Having two element layers or more enabled the introduction of the weak-interlamellar adhesion between corneal layers, which reduced the flexural stiffness of the cornea and led to a relaxation of stresses in the tissue and concentrated deformations in the region of the limbus. This mechanical effect could be behind
Fig. 3. Variation of average and maximum errors with progression of analysis iterations in cases with (a) different values of IOP, (b) different numbers of model layers, (c) different ages, and (d) different central anterior radius of curvature. Heavy and light lines represent maximum and average errors, respectively.
A. Elsheikh et al. / Medical Engineering & Physics 35 (2013) 211–216
215
Fig. 5. Reduction of error in the refractive power of the reference numerical model with progression of iterative cycles. Fig. 4. Variation of nodal error along a meridional section of the eye obtained with the reference eye model and following the first iteration.
the variations observed in the convergence toward the stress-free configuration between the model with only one layer of elements and those with two layers or more. The mechanical properties of ocular tissue were allowed to vary in the study through their association with age. The results showed that stiffer mechanical behavior (with higher age) was associated with improved convergence of the model. Only three iteration cycles were needed to achieve the same accuracy with the 90-year-old properties as with four iterations in models with the other age groups (Fig. 3c). On the other hand, the geometric parameters showed limited effects on the convergence to the stress-free configuration. The results demonstrated that increased scleral radius or corneal shape factor did not cause any significant change in model convergence (results not shown). However, the convergence rate was slightly altered by the corneal radius, with larger radii producing slower convergence (Fig. 3d). The results show consistently that models with a high level of flexibility, especially in the corneal region, undergo slow convergence to the stress-free configuration. This was evident in the models that considered the weak interlamellar adhesion in corneal stroma, models with soft material behavior associated with young age, and models with larger corneal radius. The resulting slow convergence in these cases could be explained by the larger topography adjustments, caused by the larger model displacements under IOP, which needed to be considered to obtain the stress-free configuration. Furthermore, the parametric study demonstrated that four iterative cycles were sufficient to determine the stress-free configuration of the human eye with average and maximum errors below 10−6 mm and 10−4 mm, respectively. This level of inaccuracy induces changes in the patients’ visual acuity of about 10−4 diopter (Fig. 5), which is orders of magnitude lower than the changes that the human optical system is able to perceive. 4. Discussion The geometry of biological structures is usually only available in a stressed state under the physiological loading of the tissues. In order to provide accurate predictions of mechanical behavior, numerical models should account for this physiological deformation. In this study, we present a simple iterative method
to determine the stress-free configuration of biological structures and evaluate the method within a parametric study of the human eye. The parametric study, which covered a wide range of ocular shape parameters (cornea radius and shape factor, and sclera radius), boundary conditions (IOP), mechanical properties (controlled by age) and finite element mesh construction (with 1, 2 and 4 element layers), demonstrated similar rates of convergence toward the stress-free configuration in all cases. However, the results also indicate that the stiffness of the model and the IOP are associated with a slight effect on the convergence behavior. Fewer iteration cycles are needed to find the stress-free configuration if the distance between the target configuration and the initial (stressed) configuration is small. This effect explains the slow convergence in cases with higher internal pressure and softer material, and the higher errors observed in the corneal region (with lower stiffness relative to the sclera) compared to other regions of the eye. In all cases studied, the average error in determining the stress-free configuration using the iterative method is lower than 10−6 mm and the maximal error less than 10−4 mm, which is orders of magnitude lower than the accuracy of current clinical topographic systems. This level of accuracy is achievable within 3–4 cycles of analysis iterations, and there would be no added benefit in increasing the number of cycles since the target (stressed) configuration is in any case not perfectly determined experimentally. The accuracy of the results reported by Grytz and Dows [16] for an iterative forward prestressing method is higher than the accuracy achieved with the proposed approach. These authors report a maximal error decreasing to less than 10−11 m after about 10 iterative cycles, but the accuracy achieved after four iterations of the forward prestressing approach is similar to that obtained in the present work. The numerical model used in this study has a number of limitations concerning the mechanical model of the tissue. The well-known anisotropic distribution of the collagen fibers within the tissue was not taken into account in the constitutive equations. This was necessary since the age-dependent mechanical properties of the tissue (which were deemed more influential on model deformation than anisotropy) had not been evaluated using anisotropic mechanical models. The effect of omitting corneal anisotropy and local changes of the curvature on the
216
A. Elsheikh et al. / Medical Engineering & Physics 35 (2013) 211–216
efficiency of the model in acquiring the stress free form are difficult to quantify, although they are not expected to be significant. Although the model adopted the thickness variation in both the cornea and sclera as reported in the literature [9,10], the thickness variation around the optic nerve head was not considered as it was thought not to affect the assessment of the proposed method to acquire the stress-free form. Other ocular characterization that have not been considered in the model for being unlikely to produce a significant change in model performance include limbal ellipticity, difference between sagittal and tangential radii of curvature and the regional variation of mechanical properties across the sclera. Another limitation of the study concerns the idealized model of the eye, which overlooked the local change in curvature of patient-specific data. This approach has been adopted in order to test the proposed algorithm within a controlled parametric environment and enable the evaluation of the individual effects of each mechanical and geometrical parameter. 5. Conclusions In this study, a generic method has been proposed to determine the stress-free configuration of the human eye. The method does not require any modification of the finite element formulation and could be used with any type of material model. The method has been shown to converge quickly to a level of accuracy suitable for biomechanical analysis. Although the initial application was ocular biomechanics, this approach could be directly applied to other biological structures such as ligaments or blood vessels. Funding This work was funded by EPSRC grant EP/H052046/1. Ethical approval Not required. Conflict of interest statement The authors have no financial or personal interests in the work reported in this paper.
References [1] Pandolfi A, Holzapfel GA. Three-dimensional modeling and computational analysis of the human cornea considering distributed collagen fibril orientations. J Biomech Eng 2008;130(6):061006. [2] Munnerlyn C, Koons SJ, Mearshall J. Photorefractive qeratectomy: a technique for laser refractive surgery. J Cataract Refract Surg 1988;14:46–52. [3] Baino F. Scleral buckling biomaterials and implants for retinal detachment surgery. Med Eng Phys 2010;32(9):945–56. [4] Roy AS, Dupps WJ. Patient-specific modeling of corneal refractive surgery outcomes and inverse estimation of elastic property changes. J Biomech Eng 2011;133(1):011002. [5] Alastrue V, Calvo B, Pena E, Doblare M. Biomechanical modeling of refractive corneal surgery. J Biomech Eng 2006;128(1):150–60. [6] Villamarin A, Roy S, Hasballa R, Vardoulis O, Reymond P, Stergiopulos N. 3D simulation of the aqueous flow in the human eye. Med Eng Phys; 2012, http://dx.doi.org/10.1016/j.medengphy.2012.02.007 [7] Ganesan B, He S, Xu H. Modelling of pulsatile blood flow in arterial trees of retinal vasculature. Med Eng Phys 2011;33:810–23. [8] Pandolfi A, Fotia G, Manganiello F. Finite element simulations of laser refractive corneal surgery. Eng Comput 2009;25:15–24. [9] Elsheikh A, Geraghty B, Rama P, Campanelli M, Meek KM. Characterization of age-related variation in corneal biomechanical properties. J Roy Soc Interface 2010;7:1475–85. [10] Elsheikh A, Geraghty B, Alhasso D, Knappett J, Campanelli M, Rama P. Regional variation in the biomechanical properties of the human sclera. Exp Eye Res 2010;90:624–33. [11] Elsheikh A, Brown M, Alhasso D, Rama P, Campanelli M, GarwayHeath D. Experimental assessment of corneal anisotropy. J Refract Surg 2008;24(2):178–87. [12] Elsheikh A, Kassem W, Jones SW. Strain-rate sensitivity of porcine and ovine corneas. Acta Bioeng Biomech 2011;13(2):25–36. [13] Govindjee S, Mihalic PA. Computational methods for inverse deformations in quasi-incompressible finite elasticity. Int J Numer Methods Eng 1998;43(5):821–38. [14] Govindjee S, Mihalic PA. Computational methods for inverse finite elastostatics. Comput Methods Appl Mech Eng 1996;136(1–2):47–57. [15] Lu J, Zhou X, Raghavan ML. Inverse elastostatic stress analysis in pre-deformed biological structures: demonstration using abdominal aortic aneurysms. J Biomech 2007;40(3):693–6. [16] Grytz R, Downs JC. A forward incremental prestressing method with application to inverse parameter estimations and eye-specific simulations of posterior scleral shells. Comput Methods Biomech Biomed Eng 2012 [Epub ahead of print]. [17] Gee MW, Förster C, Wall WA. A computational strategy for prestressing patient-specific biomechanical problems under finite deformation. Int J Numer Methods Biomed Eng 2010;26(1):52–72. [18] Lanchares E, Calvo B, Cristóbal JA, Doblaré M. Finite element simulation of arcuates for astigmatism correction. J Biomech 2008;41(4):797–805. [19] Elsheikh A, Ross S, Alhasso D, Rama P. Numerical study of the effect of corneal layered structure on ocular biomechanics. Curr Eye Res 2009;34(1):26–35. [20] Elsheikh A, Wang D, Rama P, Campanelli M, Garway-Heath D. Experimental assessment of human corneal hysteresis. Curr Eye Res 2008;33(3):205–13. [21] Francis BA, Hsieh A, Lai MY, Chopra V, Pena F, Azen S, et al. Effects of corneal thickness, corneal curvature, and intraocular pressure level on Goldmann applanation tonometry and dynamic contour tonometry. Ophthalmology 2007;114(1):20–6. [22] Bennett G, Rabbetts RB. Clinical visual optics, vol. 2. London: Butterworths; 1989.