Journal of Refractive Surgery

Biomechanics 

Biomechanical and Refractive Behaviors of Keratoconic Cornea Based on Three-Dimensional Anisotropic Hyperelastic Models

Zhaolong Han; Xiaohong Sui, PhD; Dai Zhou, PhD; Chuanqing Zhou, PhD; Qiushi Ren, PhD

Abstract

PURPOSE:

To investigate the biomechanical and refractive behaviors of normal and keratoconic corneas based on three-dimensional anisotropic hyperelastic corneal models with two layers.

METHODS:

Based on an anisotropic hyperelastic formula, the finite element method was employed to develop normal and keratoconic corneal models in which the fiber orientations and the biomechanical differences between corneal layers were taken into account. The displacements for normal and keratoconic corneal models were studied, as well as changes in corneal refractive power with intraocular pressure (IOP).

RESULTS:

There were different displacements for keratoconic and normal corneas. Positive correlations were found in the keratoconic cornea between IOP and apical displacement, as well as between IOP and corneal refractive power. Under normal IOP, both the corneal shape and refractive power map were affected by the stiffness distributions of the corneal layers.

CONCLUSIONS:

Finite element analysis can be used to demonstrate the biomechanical and refractive behavior of a cornea with keratoconus. From a biomechanical viewpoint, the displacement changes seen under normal IOP were due to the decreased stiffness in the keratoconic corneal tissue and local thinning disorders. Thus, the curvature and corneal refractive power map will be abnormal in keratoconus.

[J Refract Surg. 2013;29(4):282–290.]

From the School of Naval Architecture, Ocean and Civil Engineering (ZH, DZ), and the School of Biomedical Engineering (ZH, XS, CZ), Shanghai Jiao Tong University, Shanghai, China; and the Department of Biomedical Engineering, Peking University, Beijing, China (QR).

Supported by the Key Project of Fund of Science and Technology Development of Shanghai (No. 10JC1407900), and the National Natural Science Foundation of China (Projects No. 51078230, 11172174, 51278297), the Research Fund for the Doctoral Program of Higher Education (200802481029), Shanghai Jiao Tong University Innovation Fund for Postgraduates.

The authors have no financial or proprietary interest in the materials presented herein.

The authors thank Dr. T. FitzGibbon for comments on earlier drafts of the manuscript.

AUTHOR CONTRIBUTIONS

Study concept and design (ZH, QR); data collection (ZH); analysis and interpretation of data (ZH, QR, XS, CZ, DZ); drafting of the manuscript (ZH, XS); critical revision of the manuscript (CZ, DZ, XS); statistical expertise (CZ); administrative, technical, or material support (CZ, DZ, XS); supervision (CZ, DZ)

Correspondence: Dai Zhou, PhD, Department of Civil Engineering, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, No. 800 Dongchuan Road, Minhang District, Shanghai 200240, China. E-mail: zhoudai@sjtu.edu.cn

Received: December 14, 2010
Accepted: January 09, 2013

Abstract

PURPOSE:

To investigate the biomechanical and refractive behaviors of normal and keratoconic corneas based on three-dimensional anisotropic hyperelastic corneal models with two layers.

METHODS:

Based on an anisotropic hyperelastic formula, the finite element method was employed to develop normal and keratoconic corneal models in which the fiber orientations and the biomechanical differences between corneal layers were taken into account. The displacements for normal and keratoconic corneal models were studied, as well as changes in corneal refractive power with intraocular pressure (IOP).

RESULTS:

There were different displacements for keratoconic and normal corneas. Positive correlations were found in the keratoconic cornea between IOP and apical displacement, as well as between IOP and corneal refractive power. Under normal IOP, both the corneal shape and refractive power map were affected by the stiffness distributions of the corneal layers.

CONCLUSIONS:

Finite element analysis can be used to demonstrate the biomechanical and refractive behavior of a cornea with keratoconus. From a biomechanical viewpoint, the displacement changes seen under normal IOP were due to the decreased stiffness in the keratoconic corneal tissue and local thinning disorders. Thus, the curvature and corneal refractive power map will be abnormal in keratoconus.

[J Refract Surg. 2013;29(4):282–290.]

From the School of Naval Architecture, Ocean and Civil Engineering (ZH, DZ), and the School of Biomedical Engineering (ZH, XS, CZ), Shanghai Jiao Tong University, Shanghai, China; and the Department of Biomedical Engineering, Peking University, Beijing, China (QR).

Supported by the Key Project of Fund of Science and Technology Development of Shanghai (No. 10JC1407900), and the National Natural Science Foundation of China (Projects No. 51078230, 11172174, 51278297), the Research Fund for the Doctoral Program of Higher Education (200802481029), Shanghai Jiao Tong University Innovation Fund for Postgraduates.

The authors have no financial or proprietary interest in the materials presented herein.

The authors thank Dr. T. FitzGibbon for comments on earlier drafts of the manuscript.

AUTHOR CONTRIBUTIONS

Study concept and design (ZH, QR); data collection (ZH); analysis and interpretation of data (ZH, QR, XS, CZ, DZ); drafting of the manuscript (ZH, XS); critical revision of the manuscript (CZ, DZ, XS); statistical expertise (CZ); administrative, technical, or material support (CZ, DZ, XS); supervision (CZ, DZ)

Correspondence: Dai Zhou, PhD, Department of Civil Engineering, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, No. 800 Dongchuan Road, Minhang District, Shanghai 200240, China. E-mail: zhoudai@sjtu.edu.cn

Received: December 14, 2010
Accepted: January 09, 2013

Keratoconus is a non-inflammatory disease associated with characteristics of corneal ectasia and tissue deterioration.1,2 In keratoconic eyes, not only were localized corneal thinning disorders found,1,2 but the stiffness was also much lower than that of normal corneas.3 All of these changes lead to alterations in corneal biomechanics and are accompanied by the keratoectasia phenomenon.

The cornea is a structure composed of several layers and, of these, the stroma accounts for approximately 90% of the corneal tissue that withstands most of the pressure load, such as intraocular pressure (IOP). Previous studies have demonstrated that the corneal tissue can be characterized as nearly incompressible, viscoelastic, and anisotropic, factors that are mainly determined by the complex collagen fibril distributions in the corneal stroma.4–10 Based on these biomechanical properties, several numerical corneal models using the finite element method have simulated the biomechanical responses of the cornea.11–20 One important application of these models is to simulate the processes and results of keratoconic corneas.8,17,20,21 However, the cornea is composed of many layers and there could be differences between these layers due to distinct collagen orientation and distribution within the stroma. For example, experiments have shown that there are different biomechanical behaviors between the anterior 200 μm part and the rest of the corneal tissue.22,23 But this special feature has been generally neglected in most models. Although researchers still lack a gold standard for measuring corneal elasticity, a layered model that is associated with anisotropic hyperelastic characteristics may be more appropriate to represent the cornea.

We developed an anisotropic hyperelastic finite element model to simulate normal and keratoconic corneas. Both the hyperelastic characteristics and the mechanical differences between layers were considered in current corneal models. Localized tissue deterioration and thinning disorders were simulated in the keratoconic cornea to investigate its influence on corneal shape. In addition, the relationship between IOP and corneal refractive power maps, which are associated with clinical diagnosis and assessment for corneal diseases, was explored further.

Patients and Methods

Constitutive Law

Mechanical responses of the cornea were found to be nearly incompressible, anisotropic hyperelastic, and with finite deformation.6,7,24,25 To describe the special mechanical properties of the cornea, an anisotropic hyperelastic strain energy density function, W, was introduced in formula:

W=Wmatrix (C¯)+Wfibril (C¯, A ⊗ A, B ⊗ B)+Wvolume (J)Wmatrix (C¯)=c0 (I¯1−3)Wfibril (C¯, A ⊗ A, B ⊗ B)=∑i=13[pi (I¯4−1)2i+qi (I¯6−1)2i]Wvolume (J)=k (J−1)2

Here the corneal matrix was assumed to be isotropic and was described by the Neo-Hookean formula, Wmatrix. The Wfibril was anisotropic and represented the strain energy function caused by the collagen fibrils within the cornea, whereas the Wvolume showed the near incompressibility of the corneal tissue. The arguments in these functions were the invariants of the modified right Cauchy-Green tensor = J−2/3C, where C = FTF, and F was the deformation gradient. The invariants in formula were defined as

I¯1=trC¯,      I¯2=(trC¯)2−tr(C¯2),      J=det(F),I¯4=C¯:A ⊗ A,      I¯6=C¯:B ⊗ B
The vectors A, B were the two constitutive material directions of the fibers in the cornea, and yielded |A| = 1 and |B| = 1. The c0, pi, qi, (i = 1, 2, 3) and k were relative material constants to be determined.

In this study, A = (cosθ, sinθ, 0) and B = (cosθ, −sinθ, 0), where θ was the included angle of the fiber orientation and the x axial (Figure 1A). The detailed fiber orientations are presented in Figure 1B, which is constructed based on the x-ray scattering observation by Meek and Boote.9

(A) Schematic illustration of the angle θ of the two orthotropic fibers. (B) The fiber distributions within the corneal stroma. (C) Ellipsoid shape of the cornea. The anterior and posterior surfaces were approximated by the equation x2 + y2 + (1+Q) z2 = R2 / (1+Q).

Figure 1. (A) Schematic illustration of the angle θ of the two orthotropic fibers. (B) The fiber distributions within the corneal stroma. (C) Ellipsoid shape of the cornea. The anterior and posterior surfaces were approximated by the equation x2 + y2 + (1+Q) z2 = R2 / (1+Q).

Numerical Models

The shape of the normal human cornea can be approximated as a symmetric ellipsoid.17 In three-dimensional orthogonal coordinates, where the z axis was superposed with the optical axis direction, the model was established by

x2+y2+(1+Q)z2=R21+Q
where Q was the corneal asphericity and R the apical curvature in the x and y directions. In our model, Q = −0.25, and the apical curvature of the anterior corneal surface R1 = 7.8 mm, whereas that of the posterior surface R2 = 6.8 mm (for more details, see Figure 1C).

Initially, we considered a normal cornea consisting of two main parts: an anterior 200-μm part (part a) and a posterior part (part p), as shown in Figure 2A. This two-layered model was built up based on the experimental illustration of Kohlhaas et al.22 By using the uniaxial tensile method, they found that the average Young’s modulus of the anterior 200 μm was approximately two times that of the posterior part. In our normal corneal model, the material constants for the two parts were set based on this evidence.

Material distributions for the (A) normal cornea, (B) keratoconic cornea, and (C) keratoconus region (viewed from the top). S = superior; I = inferior; N = nasal; T = temporal; a = anterior; p = posterior.

Figure 2. Material distributions for the (A) normal cornea, (B) keratoconic cornea, and (C) keratoconus region (viewed from the top). S = superior; I = inferior; N = nasal; T = temporal; a = anterior; p = posterior.

The keratoconus model was constructed as a degenerated thinner zone surrounded by healthy tissue. Previous keratoconus studies have demonstrated that the fibrous structure of the cornea was missing with a degenerative process in the keratoconic region, and that the stiffness was only 60% that of normal tissue.3,26 To describe these properties in the keratoconic cornea, four parts were built up to form the whole model, and these included most of the anterior 200 μm part (part a1), the posterior part (part p1), the localized deteriorated anterior part (part a2), and the localized deteriorated posterior part (part p2). The central thickness was 0.3 mm for the deteriorated tissue (parts b3 and b4), whereas the radius was 1.5 mm with the center point of the circle 3 mm from the corneal apex (Figures 2B and 2C).

Figure 3A illustrates the computational model of the corneas. The IOP ranged from 0 to 80 mm Hg in steps of 2 mm Hg and was perpendicular to the posterior surface (the normal IOP was set to 16 mm Hg). The peripheral connection at the cornea–sclera–limbus was represented as a hinged support.

(A) Balanced model of the cornea under intraocular pressure. (B) Finite element meshed model for the normal cornea. The model was meshed by 4 node tetrahedron element and contained 23,534 elements.

Figure 3. (A) Balanced model of the cornea under intraocular pressure. (B) Finite element meshed model for the normal cornea. The model was meshed by 4 node tetrahedron element and contained 23,534 elements.

The numerical models of normal and keratoconic corneas were created using the commercial finite element software, Ansys 10.0 (ANSYS, Inc., Canonsburg, PA). The Solid 185 element (4 node tetrahedron element), with the ability to simulate hyperelastic material, was chosen for the free mesh with a mesh dimension of 0.2 mm. The three-dimensional meshed model for normal corneas consisted of 23,534 elements (Figure 3B). For the keratoconus model, there were a total of 23,700 elements. Large displacement effects were considered in the calculations, and each IOP load step was divided into 20 sub-step loads for the purpose of computational convergence.

Corneal Refractive Power Maps

To simulate the refractive power maps of normal and keratoconic corneas, an algorithm based on Matlab software (Matlab R2007a; Mathwork Inc., Natick, MA) was developed to obtain the optimized sphere by the function

min Ri2=min1n∑j=1n [(xj−x0)2+(yj−y0)2+(zj−z0)2]
where Ri was the optimized spherical radius and x0, y0, z0 were the center point locations of this optimized sphere, whereas n was the number of nodes, and xj, yj, zj represented the locations for any chosen node j, which went through from 1 to n. This function was solved based on the least square theory of fit errors. For any node i in the anterior surface, the estimation of refractive power, Di, was shown approximately as:
Di=(n′−n)Ri
where n′ = 1.376 and n = 1 was the refractive index of the anterior corneal surface and air, respectively. The nodes within the region 1.5 mm around node i were selected for computing the radius Ri. By these means, the anterior refractive power map of the cornea could be determined.

Results

Parameter Validation

As difficult to obtain in vivo experimental data, the material parameters were acquired by comparing the numerical results with in vitro corneal inflation test results given by Elsheikh et al.24 and Bryant and McDonnell.14 Because there were obvious differences in their results, we used two sets of basic parameters (Table 1), which were also used to derive the material constants for the corneal models (Tables 2 and 3). Note that the parameters mi were related with the work of Elsheikh et al. and ni with the data of Bryant and McDonnell (i = 0 – 4).

Basic Parameters Used for the Corneal Material Constants (Unit: MPa)

Table 1: Basic Parameters Used for the Corneal Material Constants (Unit: MPa)

Material Constants for Normal Corneas Constructed of One or Two Layers (Unit: MPa)

Table 2: Material Constants for Normal Corneas Constructed of One or Two Layers (Unit: MPa)

Material Constants for Keratoconic Corneal Models (Unit: MPa)

Table 3: Material Constants for Keratoconic Corneal Models (Unit: MPa)

Figure 4A shows the simulated IOP-displacement curves for normal corneas. Groups 1 (Nor-2L-m) and 2 (Nor-2L-n) represented two-layer models because parts a and p have the same constant value for the matrix but the stiffness of the fibers in part p were only half of that in part a (Table 2).The IOP in the simulations ranged from 0 to 80 mm Hg for group 1 (Nor-2L-m) and from 0 to 40 mm Hg for group 2 (Nor-2L-n). The two simulated curves both fell within the range of values defined by experimental inflation test data, showing that our simulations were in good agreement with the experimental data.

(A) The relationships between intraocular pressure (IOP) and apical displacement derived from experimental14,24 or simulated numerical data can be compared for normal corneas. (B) Simulated numerical comparisons of the one- and two-layer models. Groups 1 and 2 were two-layer models and groups 3 to 6 were one-layer models. (C) Simulated data for IOP versus apical displacement for normal and keratoconic corneas.

Figure 4. (A) The relationships between intraocular pressure (IOP) and apical displacement derived from experimental14,24 or simulated numerical data can be compared for normal corneas. (B) Simulated numerical comparisons of the one- and two-layer models. Groups 1 and 2 were two-layer models and groups 3 to 6 were one-layer models. (C) Simulated data for IOP versus apical displacement for normal and keratoconic corneas.

One- and Two-Layer Models

The one-layered model for the normal cornea could be achieved by setting the same material constants for parts a and p. In Table 2, four one-layered corneal models (groups 3 to 6) are shown. The material constants for groups 3 (Nor-1L-m) and 4 (Nor-1L-m/2) were respectively the same as those in parts a and p in group 1 (Nor-2L-m). Similarly, the material constants for groups 5 (Nor-1L-n) and 6 (Nor-1L-n/2) were obtained from the anterior and posterior parts of group 2 (Nor-2L-n). Note that the value of c0, representing the matrix material, was much smaller than that of the fibers.

We also chose the apical-displacement curve as an index to compare differences between the one- and two-layered models. The calculated results are shown in Figure 4B. At the same IOP, the apical displacement in group 1 (Nor-2L-m) was larger than that of group 3 (Nor-1L-m) but smaller than that of group 4 (Nor-1L-m/2). Similar outcomes were found in groups 2 (Nor-2L-n), 5 (Nor-1L-n), and 6 (Nor-1L-n/2).

Keratoconus

The keratoconus models with various material constants are shown in Table 3 (groups 7 to 10). The material constants of groups 7 (Ker-1L-m) and 8 (Ker-2L-m) were obtained from mi, and that for groups 9 (Ker-1L-n) and 10 (Ker-2L-n) from ni, The matrix part remained unchanged for the keratoconic corneal models if compared with the normal cornea, and the fibers in the region of local thinning (parts a2 and p2) were set to only one-fifth of that in the normal anterior tissue (part a). The keratoconus models could also be divided into one- and two-layer types. In the one-layer type (groups 7 [Ker-1L-m] and 9 [Ker-1L-n]), the fibers in the anterior parts (part a1) were assumed to be half of part a, but the posterior part was not altered. This is equal to construct a one-layer keratoconus model because the parameters in the normal cornea model for fiber and volume for part p are only half of those for part a (see Table 3). In the two-layer models (groups 8 [Ker-2L-m] and 10 [Ker-2L-n]), the material constants of fibers in the anterior (a1) and posterior (p1) parts were set to half of parts a and p, respectively.

Figure 4C compares the simulated relationship between IOP and apical displacements for normal and keratoconic corneas. The apical displacement for keratoconus was larger than that for the normal cornea at the same IOP value, which was in agreement with the clinically observed characteristics of corneal ectasia. Because the material constants of the fibers in part p1 in groups 8 (Ker-2L-m) and 10 (Ker-2L-n) were smaller, they had larger apical displacements.

Corneal Refractive Power

The nonlinear relationship between the IOP and refractive power for the central 3-mm corneal region are plotted in Figure 5. The refractive powers of the anterior surface of normal corneas were lower than those of the keratoconus models at the same IOP (group 1 vs 7, 8; and group 2 vs 9, 10); all refractive powers rose with an increase in IOP. The results also indicated that a cornea with lower stiffness may be more likely to have a higher refractive power. Note that in the curve of group 2 (Nor-2L-n), when IOP was greater than 24 mm Hg, it had a higher refractive power than that of groups 7 (Ker-1L-m) and 8 (Ker-2L-m), both of which were keratoconus models.

Numerical relationships between intraocular pressure and the refractive power of a 3-mm diameter anterior surface region for normal corneas and keratoconus models.

Figure 5. Numerical relationships between intraocular pressure and the refractive power of a 3-mm diameter anterior surface region for normal corneas and keratoconus models.

The refractive power maps of one- and two-layer normal and keratoconic corneas subjected to an IOP of 16 mm Hg are presented and compared in Figure 6. Groups 1 to 6 were models for the normal cornea; thus the refractive power maps were symmetrical with a center power of approximately 48 diopters, which was close to the theoretical value. But in each keratoconus model, which had a lower corneal stiffness and a local thinning region, the power map became asymmetrical and deviated to the side with the deteriorated region. The corneal stiffness in the model for group 10 (Ker-2L-n) was the lowest and also had a central refractive power as high as 55.01 diopters, which was typical of keratoconus.

Simulated refractive power maps for normal corneas (groups 1–6) and keratoconic corneas (groups 7–10) when intra-ocular pressure was 16 mm Hg. Central refractive powers for these modes were: groups 1 to 6 (48.05 diopters [D]), group 7 (Ker-1L-m) (49.20 D), group 8 (Ker-2L-m) (49.61 D), group 9 (Ker-1L-n) (52.02 D), and group 10 (Ker-2L-n) (55.01 D).

Figure 6. Simulated refractive power maps for normal corneas (groups 1–6) and keratoconic corneas (groups 7–10) when intra-ocular pressure was 16 mm Hg. Central refractive powers for these modes were: groups 1 to 6 (48.05 diopters [D]), group 7 (Ker-1L-m) (49.20 D), group 8 (Ker-2L-m) (49.61 D), group 9 (Ker-1L-n) (52.02 D), and group 10 (Ker-2L-n) (55.01 D).

Discussion

In this study an improved corneal finite element model for normal and keratoconic corneas has been developed implementing the three-dimensional, hyperelastic, anisotropic, and layer-specific properties of corneal biomechanics. The anisotropic property is one complexity in corneal biomechanics. Evidence shows that there are irregular collagen fibril distributions and orientations embedded in the matrix within the stroma. Meek and Boote9 investigated human collagen orientation by x-ray observation, and proposed that corneal fibrils went orthogonally through the superior-inferior and nasal-temporal directions and fused with circumferential collagen at the limbus. We used a strain energy function that could mathematically depict the anisotropic hyperelastic properties of the cornea. This strain energy function was in a polynomial form, which was similar to the exponential form adopted by Holzapfel et al.12 and Pandolfi et al.17–19 The matrix was soft and the fibers were strong with orientation vector A, B describing the anisotropy.

The cornea consists of several layers; thus, mechanical behavior may vary between these layers. As shown by Komai and Ushiki, corneal lamellae in the anterior one-third of the corneal stroma are more narrowly interwoven than those in the other two-thirds.27 Kohlhaas et al. illustrated that the Young’s modulus at 5% strain in the anterior 200 μm was more than two times that of the posterior part in the human stroma.22 These differences between layers will also partly contribute to the anisotropy of the cornea.

To describe this property, a two-layered model for the normal cornea was constructed. Several researchers have adopted an anisotropic one-layered model. For example, Pandolfi et al. incorporated the special corneal fibril orientations into a one-layer model,17–19 whereas Gefen et al. employed an anisotropic linear elastic description.20 Although it is hard to clarify which model proved to be a better descriptor of clinical observations, the two-layered model has some advantages over the one-layered model. First, the biomechanical differences between corneal layers can be easily incorporated into the two-layered model but cannot be realized when using the one-layered model. Second, the two-layered cornea can “become one layer” by setting the same material constants for the two layers. In other words, the two-layer model can be used as a one-layered model.

Fernández et al.23 also reported the results obtained from a multi-layered corneal model, but their model was created as a two-dimensional model. From the viewpoint of mechanics, the simplified two-dimensional model should be constructed in spherical or columniation coordinates instead of Cartesian coordinates. In our study, the three-dimensional corneal models were built using Cartesian coordinates, which are more suitable for representing the actual cornea.

The boundary conditions in our models were set as a hinged support for two reasons. One was due to experimental pictures and data demonstrating that the limbus, although a soft tissue, cannot be bent but can provide enough tensile strength to balance the IOP.6 The second reason was that our numerical results were compared with experimental data of inflation tests, in which most of the sclera was removed and the peripheral part of the cornea was strongly fixed in the test equipment.14,24 However, this simplified assumption for the flexible limbus may be a weakness of the finite element modeling and cause some unforeseeable errors. For example, it may be responsible for the inappropriate maximum refractive power position that is located around the central region instead of the keratoconic zone.

The mechanical responses of the cornea may vary between individuals. This point has been illustrated by the normal corneal inflation experiments conducted by Elsheikh et al.24 and Bryant and McDonnell.14 They performed similar experiments but the results were different, and because of this we used two sets of parameters, mi and ni (i = 0 – 4), to simulate their results. Although our results of numerical simulation in Figure 6 were in agreement with their experimental results, the experimental human corneas in the study were limited to older subjects and in vitro testing (65 to 79 years of age in Elsheikh et al.24; data were not given in the study of Bryant and McDonnell14). Thus, a wider use of the parameters needs more validation.

Although keratoconus is a common ocular disease, its etiology remains unclear.1 With the progressive ectasia and the local thinning in keratoconus,1 corneal mechanical responses become abnormal.3 Previous work9,28,29 found through x-ray observations that the decreasing collagen lamellae arrangement in keratoconus differed from the arrangement seen in normal corneas. Kohlhaas et al. discussed the biomechanical behaviors between corneal layers and gave a possible explanation such that in keratoconus only the anterior tissue, with a thickness of 200 μm, changed.22 Based on both the layered corneal structure and the hypothesis put forward by Kohlhaas et al., we set two types of keratoconus models. One corresponded to the assumption that only the anterior layer (part a1) changed (groups 7 and 9),22 whereas the second model considered the differences between the layers; both the anterior (part a1) and the posterior (part p1) layers deteriorated (groups 8 and 10).

Refractive power is important in the clinical diagnosis of ocular diseases such as keratoconus. Auffarth et al. used the Orbscan Topography System to investigate 71 examples of keratoconus and found that there was a positive relationship between the apex mean curvature and apex elevation.30 According to our numerical simulations, apical displacement increases and the central refractive power becomes larger as IOP increases (Figure 6), which was in accordance with Auffarth et al.’s clinical observations.

In addition to the local thinning disorders in keratoconus, the refractive map may be determined mainly by the corneal modulus of elasticity under normal IOP. In the actual keratoconus disease, some collagen fibrils and lamellae diminish, which causes the decrease in the elasticity modulus for the whole cornea. In this way, corneal ectasia occurs, as has been shown by our simulations. In the keratoconus groups in Figure 6, the maximum refractive power position happened around the central region instead of the weaker keratoconic zone. This may because of the combining effects of the fixed boundary conditions, fiber parameter assignments, and deformation compatibility.

More recently, Roy and Dupps31 have completed a similar and valuable work on three-dimensional finite element modeling on keratoconus. Their keratoconic models were set up based on clinical measurement and the progressions and responses after collagen cross-linking were further numerically studied. Meanwhile, they creatively adopted an exponential formula to describe the tissue decay from the center to the edge of the deteriorated weak zone, which is a good choice but not considered in the current study.

A two-layered anisotropic corneal model involving the biomechanical differences between layers was developed to investigate the biomechanical and refractive behaviors of the normal and keratoconic cornea. The anisotropy of the cornea was accomplished by defining the fiber orientations and material constants between layers. Numerical simulations showed a positive relationship between corneal refractive power and IOP. Due to the gradual decrease in stiffness of keratoconic corneas subjected to normal IOP, the changing corneal displacement (so-called corneal ectasia) was accompanied by a stress redistribution and an increase of refractive power. This model can be used to explain some changes in keratoconus. Further investigation will be done to simulate the effect of implantation of an intrastromal corneal ring in keratoconus.

However, the two-layer model is just a simplified one and is admittedly imperfect. For instance, it cannot describe the depth-dependent change in strength as experimentally shown in the work by Randleman et al.32 To give more accurate descriptions of clinical investigations such as SMILE, where the anterior stroma is left intact, and LASIK, in which some tissue is removed from the anterior stroma, the models with several element layers and model refinement are better choices.

References

  1. Rabinowitz YS. Keratoconus. Surv Ophthalmol. 1998;42:297–319 doi:10.1016/S0039-6257(97)00119-7 [CrossRef] .
  2. Krachmer JH, Feder RS, Belin MW. Keratoconus and related noninflammatory corneal thinning disorders. Surv Ophthalmol. 1984;28:293–322 doi:10.1016/0039-6257(84)90094-8 [CrossRef] .
  3. Andreassen TT, Simonsen AH, Hans O. Biomechanical properties of keratoconus and normal corneas. Exp Eye Res. 1980;31:435–441 doi:10.1016/S0014-4835(80)80027-3 [CrossRef] .
  4. Roberts CJ. Biomechanical customization: the next generation of laser refractive surgery. J Cataract Refract Surg. 2005;31:2–5 doi:10.1016/j.jcrs.2004.11.032 [CrossRef] .
  5. Roberts CJ. Biomechanics of the cornea and wavefront-guided laser refractive surgery. J Refract Surg. 2002;18:S589–S592.
  6. Wollensak G, Spoerl E, Seiler T. Stress-strain measurements of human and porcine corneas after riboflavin-ultraviolet-A-induced cross-linking. J Cataract Refract Surg. 2003;29:1780–1785 doi:10.1016/S0886-3350(03)00407-3 [CrossRef] .
  7. Elsheikh A, Alhasso D. Mechanical anisotropy of porcine cornea and correlation with stromal microstructure. Exp Eye Res. 2009;88:1084–1091 doi:10.1016/j.exer.2009.01.010 [CrossRef] .
  8. Anderson K, Elsheikh A, Newson T. Application of structural analysis to the mechanical behaviour of the cornea. J R Soc Interface. 2004;1:3–15 doi:10.1098/rsif.2004.0002 [CrossRef] .
  9. Meek KM, Boote C. The use of X-ray scattering techniques to quantify the orientation and distribution of collagen in the corneal stroma. Prog Retin Eye Res. 2009;28:369–392 doi:10.1016/j.preteyeres.2009.06.005 [CrossRef] .
  10. Elsheikh A, Alhasso D, Kotecha A, Garway-Heath D. Assessment of the Ocular Response Analyzer as a tool for intraocular pressure measurement. J Biomech Eng. 2009;131:081010-1–081010-19 doi:10.1115/1.3148462 [CrossRef] .
  11. Elsheikh A, Wang DF. Numerical modelling of corneal biomechanical behavior. Comput Method Biomec. 2007;10:85–95 doi:10.1080/10255840600976013 [CrossRef] .
  12. Holzapfel GA, Eberlein R, Wriggers P, Weizsäcker HW. Large strain analysis of soft biological membranes: formulation and finite element analysis. Comput Method Appl M. 1996;132:45–61 doi:10.1016/0045-7825(96)00999-1 [CrossRef] .
  13. Deenadayalu C, Mobasher B, Rajan SD, Hall GW. Refractive change induced by the LASIK flap in a biomechanical finite element model. J Refract Surg. 2006;22:286–292.
  14. Bryant MR, McDonnell PJ. Constitutive laws for biomechanical modeling of refractive surgery. J Biomech Eng-T ASME. 1996;118:473–481 doi:10.1115/1.2796033 [CrossRef] .
  15. Pinsky PM, Heide D, Chernyak D. Computational modeling of mechanical anisotropy in the cornea and sclera. J Cataract Refract Surg. 2005;31:136–145 doi:10.1016/j.jcrs.2004.10.048 [CrossRef] .
  16. Djotyan GP, Kurtz R M, Fernández DC, Juhasz T. An analytically solvable model for biomechanical response of the cornea to refractive surgery. J Biomech Eng-T ASME. 2001;123:440–445 doi:10.1115/1.1388293 [CrossRef] .
  17. Pandolfi A, Manganiello F. A model for the human cornea: constitutive formulation and numerical analysis. Biomech Model Mechan. 2006;5:237–246 doi:10.1007/s10237-005-0014-x [CrossRef] .
  18. Pandolfi A, Fotia G, Manganiello F. Finite element simulations of laser refractive corneal surgery. Eng Comput-Germany. 2009;25:15–24 doi:10.1007/s00366-008-0102-5 [CrossRef] .
  19. Pandolfi A, Holzapfel GA. Three-dimensional modeling and computational analysis of the human cornea considering distributed collagen fibril orientations. J Biomech Eng-T ASME. 2008;130:061006-1–061006-12 doi:10.1115/1.2982251 [CrossRef] .
  20. Gefen A, Shalom R, Elad D, Mandel Y. Biomechanical analysis of the keratoconic cornea. J Mech Behav Biomed. 2009;2:224–236 doi:10.1016/j.jmbbm.2008.07.002 [CrossRef] .
  21. Carvalho LA, Prado M, Cunha RH, et al. Keratoconus prediction using a finite element model of the cornea with local biomechanical properties [article in Portuguese]. Arquivos Brasileiros de Oftalmologia. 2009;72:139–145 doi:10.1590/S0004-27492009000200002 [CrossRef] .
  22. Kohlhaas M, Spoerl E, Schilde T, Unger G, Wittig C, Pillunat LE. Biomechanical evidence of the distribution of cross-links in corneas treated with riboflavin and ultraviolet A light. J Cataract Refract Surg. 2006;32:279–283 doi:10.1016/j.jcrs.2005.12.092 [CrossRef] .
  23. Fernández DC, Niazy AM, Kurtz RM, Djotyan GP, Juhasz T. Finite element analysis applied to cornea reshaping. J Biomed Opt. 2005;10:064018 doi:10.1117/1.2136149 [CrossRef] .
  24. Elsheikh A, Alhasso D, Rama P. Biomechanical properties of human and porcine corneas. Exp Eye Res. 2008;86:783–790 doi:10.1016/j.exer.2008.02.006 [CrossRef] .
  25. Elsheikh A, Alhasso D, Rama P. Assessment of the epithelium’s contribution to corneal biomechanics. Exp Eye Res. 2008;86:445–451 doi:10.1016/j.exer.2007.12.002 [CrossRef] .
  26. Zimmermann DR, Fischer RW, Winterhalter KH, Witmer R, Vaughan L. Comparative studies of collagens in normal and keratoconus corneas. Exp Eye Res. 1988;46:431–422 doi:10.1016/S0014-4835(88)80031-9 [CrossRef] .
  27. Komai Y, Ushiki T. The three-dimensional organization of collagen fibrils in the human cornea and sclera. Invest Ophthalmol Vis Sci. 1991;32:2244–2258.
  28. Fratzl P, Daxer A. Structural transformation of collagen fibrils in corneal stroma during drying. Biophys J. 1993;64:1210–1214 doi:10.1016/S0006-3495(93)81487-5 [CrossRef] .
  29. Daxer A, Fratzl P. Collagen fibril orientation in the human cornea stroma and its implication in keratoconus. Invest Ophth Vis Sci. 1997;38:121–129.
  30. Auffarth GU, Wang L, Volcker HE. Keratoconus evaluation using the Orbscan topography system. J Cataract Refract Surg. 2000;26:222–228 doi:10.1016/S0886-3350(99)00355-7 [CrossRef] .
  31. Roy AS, Dupps WJ. Patient-specific computational modeling of keratoconus progression and differential responses to collagen cross-linking. Invest Ophth Vis Sci. 2011;52:9174–9187 doi:10.1167/iovs.11-7395 [CrossRef] .
  32. Randleman JB, Dawson DG, Grossniklaus HE, McCarey BE, Edelhauser HF. Depth-dependent cohesive tensile strength in human donor corneas: implications for refractive surgery. J Refract Surg. 2008;24(suppl):S85–S89.

Basic Parameters Used for the Corneal Material Constants (Unit: MPa)

m0 m1 m2 m3 m4
0.040 0.40 4.2 × 107 2.94 × 1015 13.154
n0 n1 n2 n3 n4
0.004 0.04 1.4 × 106 3.267 × 1013 1.315

Material Constants for Normal Corneas Constructed of One or Two Layers (Unit: MPa)

Cornea Type/Group Part a
Part p
Description
Matrix (c0) Fiber [pi (qi) (i =1,2,3)] Volume (k) Matrix (c0) Fiber [pi (qi) (i =1,2,3)] Volume (k)
1 (Nor-2L-m) m0 mi m4 m0 mi /2 m4/2 2L
2 (Nor-2L-n) n0 ni n4 n0 ni /2 n4/2 2L
3 (Nor-1L-m) m0 mi m4 m0 mi m4 1L
4 (Nor-1L-m/2) m0 mi/2 m4 /2 m0 mi /2 m4/2 1L
5 (Nor-1L-n) n0 ni n4 n0 ni n4 1L
6 (Nor-1L-n/2) n0 ni /2 n4 /2 n0 ni /2 n4/2 1L

Material Constants for Keratoconic Corneal Models (Unit: MPa)

Cornea Type/Group Part a1
Part p1
Part a2(p2)
Description
Matrix (c0) Fiber [pi (qi) (i =1,2,3)] Volume (k) Matrix (c0) Fiber [pi (qi) (i =1,2,3)] Volume (k) Matrix (c0) Fiber [pi (qi) (i =1,2,3)] Volume (k)
Keratoconus
  7 (Ker-1L-m) m0 mi /2 m4 /2 m0 mi /2 m4 /2 m0 mi /5 m4 /5 1L
  8 (Ker-2L-m) m0 mi /2 m4 /2 m0 mi /4 m4 /4 m0 mi /5 m4 /5 2L
  9 (Ker-1L-n) n0 ni /2 n4 /2 n0 ni /2 n4 /2 n0 ni /5 n4 /5 1L
  10 (Ker-2L-n) n0 ni /2 n4 /2 n0 ni /4 n4 /4 n0 ni /5 n4 /5 2L

10.3928/1081597X-20130318-08

Sign up to receive

Journal E-contents