2.1.

Mechanical Measurement of the Elasticity Parameters

To measure directly the elastic modulus of a small sample, a deformation to a known force is usually measured. A mechanical technique was used to measure elasticity distributions in large tissue samples.^{7, 8, 9} This technique produces nearly direct measures of elasticity and is widely accepted.^{11, 12, 13, 14} The entire thickness of an ex-vivo human cornea was used $\left(857\mu \mathrm{m}\right)$ in the experimental protocol developed.^{7, 8, 9} The ex-vivo human cornea was obtained from a human donor $70\text{years}$ old. The specimen was stored in corneal storage medium (Optisol, Chiron, Irvine, California) within $3\mathrm{h}$ after patient death and was tested within $2\text{days}$ . A femtosecond surgical laser,^{15, 16} capable of cutting the tissue with high accuracy, was used to make five cuts at increasing corneal depths. Thus, a series of five circular disks or layers approximately $140\text{to}350\mu \mathrm{m}$ thick and $5\mathrm{mm}$ in diameter were obtained (see Fig. 1 ).

Fig. 1

A schematic view that shows the procedure followed to cut the corneal disks using the femtosecond laser.

We note that to produce disks for mechanical measurements, a plane was cut parallel to the applanated anterior corneal surface. If these disks assume natural curvatures in the intact globe, resulting anterior and posterior surfaces would no longer be planar when the applanating contact lens is removed after resection. In practice, an unloaded tissue block has insufficient rigidity to revert to its loaded curvature and assumes a planar shape through surface tension when placed on a flat surface.¹ The resulting disk sections have uniform thickness with parallel, planar surfaces.¹⁷ Each sample was positioned horizontally on the rigid platform of a scale ( $0.01\text{-}\mathrm{g}$ resolution). Deformations were applied with an $8\text{-}\mathrm{mm}$ -diam piston covering the entire surface of the sample. The piston was lowered onto the sample using a dc stepper motor driven positioning system (Unislide MB2506p5j-s2.5, Velmex Incorporated, East Bloomfield, New York). The piston was then slowly pulled back until a force of 0.01 Newton was reached. At fixed time intervals, the piston was lowered a single step and five consecutive force measurements were performed for each corneal button.⁷ After each trial, there was a waiting period to allow the sample to relax. For the next trial, the 0.01 Newton preload point is found as described earlier, and the measurement was repeated. Consecutive measurements were used to ensure that deformations were applied at a rate slow enough to measure the static elastic properties of the sample. Details about the reproducibility and error analysis of the mechanical measurements using this technique have been reported previously by Erkamp ⁷

The deformation measurements were then converted to an elastic modulus, where the conversion uses prior calibration of the system with plastisol samples of known Young’s modulus.⁷ A third-order polynomial was fitted to the raw data obtained experimentally. The stress-strain data for each layer and for 10% of deformation is shown in Fig. 2 . Layer 1 is the outermost layer. Since the width of the strip is sufficiently narrow ( $\mathrm{t}∕\mathrm{R}<0.1$ , assuming a radius of curvature of $7.8\mathrm{mm}$ for the intact eye), the possibility for changes in the transverse direction associated to the bending effect can be neglected.¹⁴ Moreover, because the corneal disk is thin when compared to its length, the bending stress profile should have little effect on the average stress through the corneal layer. Note that our reported values of stress and strain should, therefore, be interpreted as averages through the corneal strip. The measured Young’s elastic modulus as a function of tissue depth is shown in Fig. 3 . Elasticity variations with depth are shown for 10, 20, and 30% deformations.

Fig. 2

Stress-strain curves for each layer, and for 10% of deformation. A third-order polynomial was fitted to the raw data obtained experimentally. The solid lines are for the fitted curves, and the crosses are for the raw data.

Fig. 3

Elasticity variations with depth are shown for 10, 20, and 30% deformations. Depth is measured at the midpoint of the corneal disk.

2.2.

Computer Simulations

The general aspheric shape (steeper centrally and flatter peripherally) of the cornea and the variation in the thickness is considered in all the simulations developed. Figure 4(a) shows the details of the geometric configuration used in the finite element model. The cornea was assumed to have a radius of $7.56\mathrm{mm}$ and a base diameter of $11.5\mathrm{mm}$ .¹⁸ Corneal thickness is assumed to vary linearly along the meridional direction [i.e., from $P_C$ to $P_L$ , see Fig. 4(a)]. It is taken to be $0.5\mathrm{mm}$ at the center $\left(T_C\right)$ and $0.67\mathrm{mm}$ at the limbus $\left(T_L\right)$ . The total area covered by the half of the corneal model was $3.6\mathrm{mm}$ .² The central corneal thickness $\left(T_C\right)$ is measured along the anterior-posterior axis of the eye that actually represents the axis of symmetry. The direction along this axis is used to measure any change in apical or central thickness. The radius is drawn from the geometric center of the corneal anterior center. The corneal thickness was divided in five layers to simulate the variation in stiffness through the thickness [see Fig. 4(a)]. Furthermore, the cornea is also modeled as an axisymmetric structure, with the axis of symmetry along the anterior-posterior axis of the eye (see Fig. 4). This greatly reduces the complexity of the analysis of the 3-D structure to a single, 2-D finite element analysis.

Fig. 4

Illustration is used in deriving the finite element analysis. (a) Geometric diagram of half of a corneal section. Note that $T_C$ is the apical corneal thickness $\left(0.5\mathrm{mm}\right)$ . $T_L$ is the peripheral corneal thickness $\left(0.67\mathrm{mm}\right)$ at a location $P_L$ a distance $D_C$ $\left(6\mathrm{mm}\right)$ from the corneal apex. The corneal segment is divided into five layers ( $L_1$ to $L_5$ ). $R_C$ is the apical radius $\left(7.56\mathrm{mm}\right)$ of the anterior corneal surface at the point $P_C$ . $R_L$ is the radius of the anterior corneal surface at the point $P_L$ (i.e., at the limbus). $H$ is the height of the cornea $\left(2.5\mathrm{mm}\right)$ measured from a plane at the point $P_L$ . (b) Finite element model for the intact cornea-baseline model. Note that the anterior and posterior surfaces follow concentric circular arcs defined by the standard size of the cornea. Also, the boundary conditions (arrows) are shown at the geometric axis and at the edge of the cornea (i.e., at the limbus) along with the intraocular pressure (IOP) on the posterior surface. (c) Finite element model for the removed lenticle-changed cornea model. A convex-convex lenticle simulates the removal of corneal tissue. Note that the anterior corneal segment is thicker at the edges than at the center, having a central thickness of $160\mu \mathrm{m}$ . The convex-convex lens has a diameter of $4\mathrm{mm}$ and a central thickness of $100\mu \mathrm{m}$ . The area covered by the lenticle is $0.13{\mathrm{mm}}^2$ in the 2-D cross section. The inset picture in the upper-right part shows the removed lenticle-changed cornea model for a plano-convex lenticle. Note that the anterior corneal segment has a uniform thickness of $160\mu \mathrm{m}$ . The plano-convex lenticle has a diameter of $4\mathrm{mm}$ and a central thickness of $100\mu \mathrm{m}$ . The inset picture in the left-lower part shows the details for the two-node gap elements. Note that each gap element is defined by two nodes: $i$ and $j$ . The direction of the gap in the undeformed state is defined as the line connecting node $i$ to node $j$ . The gap distance is the maximum allowable relative displacement between the two nodes along the gap direction. (d) Finite element model for the removed lenticle-changed cornea model after the removal of a convex-convex lenticle.

Since the change in curvature between the cornea and sclera suggests circumferential fibrils from the mechanical point of view, the boundary nodes at the edge of the cornea through the thickness were assumed to be fixed in displacement and free to rotate about the axis tangent to the boundary circle [see Fig. 4(b)]. Also, due to the assumption of axis symmetry, displacements at the nodes through the apex were fixed in the horizontal direction ( $x$ axis) and free in the vertical direction ( $y$ axis) [see Fig. 4(b)]. Any role played by the sclera was ignored.

Also note that in vivo, the cornea is continuously subjected to the intraocular pressure and deforms in response to the surgical treatment. Thus, to simulate the effects of the refractive surgery only, two benchmark models were constructed: 1. intact cornea baseline [see Fig. 4(b)], and 2. removed lenticle-changed cornea [see Fig. 4(c)]. Both models were subjected to the same intraocular pressure. Differences in the predicted deformation of these two models being the result of the refractive surgery procedure could be used to calculate a change in corneal curvature. On the other hand, three different lenticles were used for simulating the surgery and testing the models: 1. a convex-convex lenticle with a central corneal thickness of $100\mu \mathrm{m}$ and a diameter of $4\mathrm{mm}$ [see Fig. 4(c)]; 2. a convex-convex lenticle with a central corneal thickness of $120\mu \mathrm{m}$ and a diameter of $4\mathrm{mm}$ , and 3. a plano-convex lenticle with a central corneal thickness of $100\mu \mathrm{m}$ and a diameter of $4\mathrm{mm}$ [see Fig. 4(c), inset picture at right]. To simulate the intrastromal cavity arising in the place of the removed lenticle, we used gap elements.^{19, 20, 21} The gap is intended to model point-to-point contact and is applied to surfaces that may come into contact [see Fig. 4(c)]. The surfaces that will be in contact could be defined as stiff surfaces that these gap elements cannot penetrate, as the gap is incrementally closed. The collapse of the intrastromal cavity is equivalent to removing a thin lens of the same shape as the cavity.²² It is assumed that the cavity collapses completely without any gaps and with only the stretching or compression of the elements required to achieve collapse.

The anterior corneal segment, being thicker at the edges than at the center, was simulated with a central corneal thickness of $160\mu \mathrm{m}$ . Figure 4(d) shows the geometry of the model after the collapse of the cavity (i.e., after the removal of the convex-convex lenticle). The cornea was loaded by a uniform pressure distribution, which was taken to be $0.002\mathrm{N}∕{\mathrm{mm}}^2$ . This pressure is equivalent to a mercury column of $15\text{-}\mathrm{mm}$ height.²³ The intraocular pressure was applied normal to the posterior surface of the model [see Fig. 2(b)]. The loading produced by the eyelids and the extraocular muscles was ignored. Corneal tissue was assumed to behave as a nearly incompressible material with a Poisson’s ratio equal to 0.49. Furthermore, four-node isoparametric axisymmetric 2-D elements were used in the simulations, thus only two translational degrees of freedom were considered for each node. To generate the mesh, several mesh densities were tested to ensure accurate and consistent convergence. Mesh that covered half of the cornea consisted of 3600 elements and 3721 nodes; when 60 elements were used across the thickness, and about 60 elements were implemented along the meridian, it resulted in a total of 7442 equations. More elements were attributed to the anterior corneal layers for higher accuracy of the modeling.

A commercially available finite element analysis program (COSMOS/M, Structural Research and Analysis Corporation, Santa Monica, California) was used for the implementation of the numerical simulations. A Newton-Raphson procedure, an iterative solution method with high quadratic convergence rate, was used in the nonlinear analysis.^{19, 20, 21} In this procedure, the tangential stiffness matrix is formed and decomposed at each iteration within a particular solution step. A force control technique is used to control the progress of the computations along the equilibrium paths of the system.^{19, 20, 21} In this technique, applied loads are used as the prescribed variables, and the loads are incrementally applied using time curves. Convergence was reached in about 7 or 8 iterations for a total of about 100 solutions steps. The relative displacement tolerance for equilibrium iterations was $0.001\mathrm{mm}$ .

To facilitate the comparison of different constitutive laws, four different models were constructed:

2.3.

Clinical Data

In this section, we briefly describe the clinical results that were simulated using the previous models. These results were taken from a previously published report.²⁵ Specifically, we compare the outcome of the simulations with clinical results after picosecond laser keratomileusis (PLK).

3.1.

Side-by-Side Comparison between the Models

A comparison of the capabilities of the different models described before was developed. The objective was to determine the average corneal curvature change in a $3\text{-}\mathrm{mm}$ zone, and the alter in central corneal thickness based on the material properties used from both the membrane inflation experiments¹⁰ and the direct mechanical measurements.^{7, 8, 9} Table 2 shows the results obtained for the models that considered a convex-convex lenticle with a central corneal thickness of $100\mu \mathrm{m}$ and a diameter of $4\mathrm{mm}$ . The results for the models considering a convex-convex lenticle with a central corneal thickness of $120\mu \mathrm{m}$ and a diameter of $4\mathrm{mm}$ are shown in Table 3 .

Table 2

Results obtained using different constitutive laws for a convex-convex lenticle with a central thickness of 100μm and a diameter of 4mm .

Table 3

Results obtained using different constitutive laws for a convex-convex lenticle with a central thickness of 120μm and a diameter of 4mm .

As can be seen from the results shown in the tables, the smallest average curvature change in a $3\text{-}\mathrm{mm}$ zone diameter was obtained when stiffness variations with depth were considered, along with material properties nonlinearities. Furthermore, the biggest change in central corneal thickness was also obtained in this case.

3.2.

Effect of the Shape of the Lenticle

The effect of the lenticle’s shape on curvature and central flattening of the cornea was also studied when stiffness variations with depth were considered, along with material properties nonlinearities. For comparison, a convex-convex lens shape [see Fig. 4(c)] and a plano-convex lenticle shape [see Fig. 4(c), inset picture in the upper-right part] were used in the modeling, both with the same central corneal thickness $\left(100\mu \mathrm{m}\right)$ , same diameter $\left(4\mathrm{mm}\right)$ , and with approximately the same area (about $0.13{\mathrm{mm}}^2$ ). The total area covered by half of the corneal model was $3.6{\mathrm{mm}}^2$ . The volume of the removed tissue was assumed the same in both models. Note that the cornea with the two different configurations for the lenticle pattern also has different configurations of its anterior segment. A plano-convex lenticle model has a uniform thickness for its anterior segment [see Fig. 4(c), inset picture in the upper-right part], and the convex-convex model has its anterior segment thicker at the edges than at the center [Fig. 4(c)]. Moreover, both cornea’s models have an anterior segment thickness that is in its central axis (i.e., at the apex) equal to $160\mu \mathrm{m}$ .

The model with a convex-convex lenticle (nonuniform anterior corneal segment thickness) leaves a stromal bed behind, which is thicker centrally than peripherally, giving a bigger correction within the central zone [see Figs. 4(c) and 5 ]. The convex-convex lenticle created more central flattening of the cornea than the plano-convex lens when stiffness variations with depth were introduced.

Fig. 5

Effect of the shape of the lenticle on corneal curvature. The horizontal line bars stand for the homogeneous and nonlinear elastic isotropy’s model, and the diagonal line bars represent the curvature changes for the inhomogeneous and nonlinear elastic isotropy’s model.

It is interesting to note how close the changes in curvature are for both models (homogeneous and inhomogeneous) when a plano-convex lenticle is used. This type of lenticle leads to the stroma being cut parallel to the anterior corneal surface, resulting in no cut ends of the fibrils to reconnect.²⁷ It seems a reasonable proposition that the cut of the layers in a more distorting manner using a convex-convex lenticle offers probably a singular advantage with respect to the distribution of the cohesive strengths between the anterior and posterior corneal layers. This might also afford better opportunities for future experiments to determine how to optimize the corneal structural changes as a function of the lack of uniformity in its underlying structure to make the refractive correction more reliable.

3.3.

Corneal Strain Analysis

In this section, meridional and depth-dependent strains predicted by the finite element model are analyzed. Only the baseline models are used for the comparison (i.e., only response to pressure loading is considered). Figure 6 shows the meridional variations of the equivalent strain for the inhomogeneous and nonlinear model, and the homogeneous and nonlinear model, respectively. Four regions of the cornea were defined: 1. central region (covering the central $3\text{-}\mathrm{mm}$ zone), 2. para-central region ( $1.5\text{to}3.5\mathrm{mm}$ from the corneal apex), 3. peripheral region ( $3.5\text{to}5.5\mathrm{mm}$ from the corneal apex), and 4. limbal region ( $5.5\text{to}6\mathrm{mm}$ from the corneal apex).

Fig. 6

The pressure-induced meridional strains at the anterior and posterior corneal surfaces.

The equivalent strain is calculated from:²⁷

The meridional, compressive, circumferential, and shear strains at the anterior and posterior corneal surface for the inhomogeneous and nonlinear model are shown in Figs. 7 and 8 . The para-central and peripheral regions deformed less in the meridional direction compared with the circumferential direction, suggesting a meridional strength of these parts.

Fig. 7

Meridional, compressive, circumferential, and shear strains at the anterior corneal surface for the inhomogeneous and nonlinear model.

Fig. 8

Meridional, compressive, circumferential, and shear strains at the posterior corneal surface for the inhomogeneous and nonlinear model.

Previous work has reported that the strain distribution on the anterior surface is highly nonuniform.⁶ Moreover, interferometric measurements of the corneal deformations at physiological pressures indicate that a dip occurs at a radial distance of $3\text{to}4\mathrm{mm}$ from the apex.²⁹ This high nonuniformity in strain was not predicted by the models tested, even when the inhomogeneities in stiffness through the thickness were considered. Instead, the models predicted a meridional strain, which decreased from a maximum at the apex. Since only the corneal mechanical properties were directly measured and no information was obtained on the mechanical properties of the sclera, this led us to hypothesize that this could be one of the reasons for the mismatch.

3.4.

Computer Simulation and Clinical Data

Computer simulations were compared with the clinical results $1\text{month}$ after picosecond laser keratomileusis procedures.²⁵ These clinical results were simulated using both the homogeneous and inhomogeneous nonlinear isotropy models presented before. Table 4 reproduces the achieved corrections obtained at $1\text{month}$ for two of the three treatment categories defined by Bryant and Juhasz.²⁵ The same input parameters defined by Bryant and Juhasz were used in our finite element models. The predictions provided by the inhomogeneous and nonlinear elastic isotropy model are also shown for comparison purposes (see Table 4). The changes in curvature were averaged in a $3\text{-}\mathrm{mm}$ zone diameter.

Table 4

Comparison between the predictions for myopic corrections provided by the finite element models, and the clinical data outcomes of PLK surgical procedures at 1month .

The model considering the inhomogeneities in stiffness overpredicted the curvature’s changes by a smaller margin when compared to the model that only took into account the material nonlinearities. The overprediction given by the inhomogeneous model was approximately $-1.26$ diopters for the $4\mathrm{mm}∕100\mu \mathrm{m}$ lenticle, and $-1.58$ diopters for the $4\mathrm{mm}∕120\mu \mathrm{m}$ lenticle. The average flattening overprediction was equal to $-1.42$ diopters. Moreover, the inhomogeneous model also overestimated the thickness changes by $13.4\mu \mathrm{m}$ for the $4\mathrm{mm}∕100\mu \mathrm{m}$ lenticle and $19.2\mu \mathrm{m}$ for the $4\mathrm{mm}∕120\mu \mathrm{m}$ lenticle. The average thinning overprediction was equal to $16.3\mu \mathrm{m}$ .

As can be seen, both the flattening and thinning overpredictions were reduced when the stiffness inhomogeneities were included in the finite element formulation. The predictions of the nonlinear-homogeneous and nonlinear-inhomogeneous models versus the achieved curvature changes at $1\text{month}$ , and for the two clinical groups, are shown in Fig. 9 . Data points represent the values for each group of patients. Error bars represent standard deviations of the measured curvature changes.

Fig. 9

Curvature changes at $1\text{month}$ for each patient versus changes in curvature predicted by the finite element models in a $3\text{-}\mathrm{mm}$ diameter zone.

3.5.

Sensitivity Analysis

A sensitivity study was performed for the inhomogeneous and nonlinear isotropy model. A design variable was perturbed at a time by a specified value, while the rest of the design variables were kept unchanged. The perturbed design variables were defined either by the actual values or by a perturbation ratio with respect to the initial value. To determine the sensitivity of the model, it was necessary to establish a baseline model. This model was derived from one of the treatment groups (group 1) defined in Table 4. The effect on corneal curvature of the intraocular pressure, central corneal thickness, and Poisson ratio was investigated. Table 5 shows the range of the parameters used in the sensitivity analysis.

Table 5

Parameters used in the sensitivity analysis.

Figure 10 shows the results when the intraocular pressure, central corneal thickness, and Poisson ratio were varied. Note that the central corneal curvature change is not significant in the normal physiological range $\left(15\text{to}20\mathrm{mm}\mathrm{Hg}\right)$ . The central curvature changed by as much as 10% below this range [see Fig. 10(a)]. The effect of corneal thickness variation on predicted curvature changes is shown in Fig. 10(b) and, as expected, thin corneas flatten more. The central curvature change distorted by as much as 20% in this range. In a thinner cornea, less uncut tissue is left behind once the lenticle is removed; whereas in a thicker cornea, more uncut tissue is left behind, which is less easily displaced and therefore produces less central flattening. Moreover, the effect of Poisson ratio variation on the predicted curvature changes is shown in Fig. 10(c). Note that the refractive change is not significant.

Fig. 10

Sensitivity analysis results obtained when the intraocular pressure, central corneal thickness, and Poisson ratio were varied.

Furthermore, to evaluate the sensitivity of the predicted curvature changes to the material properties, the group 1 $(4\mathrm{mm}∕100\mu \mathrm{m})$ results were recomputed considering a modulus of elasticity reduction of about 15% from the initial values used for all the layers. Thus, the finite element model was softer than the baseline model. The predicted curvature change of the group 1 procedure was $-15.32$ diopters, compared with the baseline model result of $-12.56$ diopters.