Measurement of intraocular pressure (IOP) in the human eye is the primary screening method for early detection of those at risk for glaucoma. The most common method for IOP measurement is applanation (Goldmann) tonometry, in which an instrument is used to flatten the central part of the cornea. A number of studies1–4 have been carried out on human eyes in vivo. Ehlers et al2 established linear correlations between corneal thickness and the error of applanation tonometry. Johnson et al3 reported a special case of one patient with 0.9-mm corneal thickness. The patient’s actual IOP, later determined by cannulation, was 11 mmHg, whereas her applanation tonometer reading was 35 mmHg. There have been numerous studies to determine the relationship between central corneal thickness (CCT) and IOP measurement in normal, glaucomatous, and ocular hypertensive eyes. Whitacre et al4 also pointed out that measuring the corneal thickness is necessary to properly interpret the results of Goldmann applanation tonometry. Doughty and Zaman,1 in a thorough analysis of the literature published from 1968 to 1999, found a statistically significant correlation between apparent IOP and CCT.
All of these studies show systematic errors in IOP tonometer readings that correlate with increased CCT. Such errors have significant implications for glaucoma screening. However, all of the results show a high degree of scatter. The purpose of this study is to demonstrate, via a series of finite element computational simulations, that the variations in other biomechanical characteristics and properties of the cornea, that are not controlled in these various studies, are likely responsible for the observed scatter.
This study uses an improved geometrical and biomechanical model of the cornea. The corneal biomechanical model is first calibrated using experimental data.5 Corneal biomechanical and geometrical parameters and true IOPs are then randomly sampled within ranges established from published data to produce a virtual population of corneas. Numerical finite element simulation of Goldmann applanation tonometry is then carried out on each member of the virtual population to determine (simulated) Goldmann IOP. Comparisons are made with the true IOPs to determine the effect of the corneal biomechanical and geometric properties on errors in IOP.
Materials and Methods
A finite element model of the human cornea, including the effect of the aqueous humor in the anterior chamber, is used in this study. The model is briefly described here, and differences with previous simplified analytical models are noted. A more detailed description is available elsewhere.6,7
The corneal surfaces are assumed to be axisymmetric (with axis of symmetry along the optical axis), but not necessarily spherical. The geometry of the anterior and posterior surfaces is expressed by the following equation8:
where the origin is on the (anterior or posterior) surface of the cornea at the optical axis, the z axis points inward along the optical axis, and r is the radial distance from the optical axis. Q is an asphericity parameter (Q=0 for a spherical surface) and R is the radius of curvature at the apex. A value of Q<0 corresponds to a prolate ellipsoidal shape, Q>0 to an oblate ellipsoidal shape. Q (=−0.18) and R (=7.77 mm) are used for the anterior surface of the cornea and Q (=−0.60) and R (=6.40 mm) for the posterior surface.8
Once the anterior and posterior surfaces are defined, the aspheric finite element model of the cornea is developed using 4-noded axisymmetric solid elements. Depending on the parameters describing the anterior and posterior surfaces, the thickness of the cornea may vary from the optical axis to the limbus. Previous studies, particularly Liu and Roberts9 and Orssengo and Pye,10 approximated the cornea as a thin spherical shell of constant thickness.
Corneal Biomechanical Properties
The nonlinear constitutive behavior of the cornea has been documented in the literature.11–14 Anderson et al15 conducted inflation tests on porcine corneas and showed significant stiffening in a collagen-regulated phase. Bryant and McDonnell5 also noted that, in each of five representative corneas, the knee of inflation curves occurred around the normal IOP (12 to 20 mmHg).
In this study, the cornea is modeled as a nonlinear transversely isotropic material, in which the elastic properties depend on the stress level in the cornea, stiffening as stresses increase. The constitutive model is derived from Fung’s exponential strain energy function W as16
where W is the strain energy density and ɛ is the (Green) strain vector. The parameter α determines the degree of nonlinearity (higher values of α correspond to greater degrees of nonlinearity) and D0
is the initial (linear) elastic tangent stiffness matrix. The constitutive relation (secant stiffness Ds
, which relates 2nd Piola-Kirchhoff stress σ to Green strain ɛ) is obtained by differentiating the strain energy function with respect to strain.
A second differentiation yields the Jacobian matrix Dt.
The initial tangent stiffness matrix D0 contains five independent elastic constants for transverse isotropy; these have been reduced to a single independent elastic constant Einit (the initial Young’s modulus in the plane of isotropy) by using correlations proposed for the human cornea by Bryant and McDonnell.5 Then, the initial in-plane Young’s modulus Einit and the nonlinearity parameter α are determined by calibration with in vitro inflation tests.5 In Bryant and McDonnell’s in vitro inflation tests, each cornea in turn was mounted on an artificial anterior chamber and its apical displacement was measured during incremental pressure increases.5 Separate inflation tests are simulated with the finite element model of the cornea. Calibrated finite element results are shown in Figure 1 and compared with the experimental results for five of these tests.
Figure 1. Calibration of Nonlinear Material Model in Finite Element Simulations (lines) with the Experimental Results (circles) by Bryant and McDonnell.5
Note that there are substantial variations between the inflation responses of the different individual corneas. The calibrated values of Einit range from 0.12 to 0.35 MPa for the five tests; calibrated values of α ranged from 417 to 1000 (1/MPa). The results also appear to indicate that the greater the initial slope of the pressure-displacement curve, the more rapidly the slope continues to increase with increasing pressure. Therefore, it is found that the material constant α, which controls the extent of the nonlinearity, can be related to Einit by a linear interpolation function with sufficient accuracy.6 The resulting model, which is a function only of the single material constant Einit, is referred to herein as the “Extended Fung Model.”
Previous cornea models have often assumed that the cornea material is linear elastic (no stress dependence) and isotropic. In this regard, it should be noted that the constant Young’s modulus used in previous models of this type represents an average value over the applanation process, starting from the (initial) physiological stress level in the cornea corresponding to the IOP, and is not directly comparable to the parameter Einit used in the current model, which represents the stiffness of the cornea at zero physiological stress. Thus, the constant Young’s modulus values used previously in linear elastic cornea models are much higher than the parameter Einit used in the model described herein. At higher IOP values, and thus higher physiological stress, corneal stiffness is greater in response to applanating force.
Initial Physiological Stress Due to IOP
Equation () describes the geometry of the cornea when it is already in a state of initial stress produced by the IOP. To take account of this initial stress state in the simulation of applanation tonometry, a two-step numerical procedure is used herein.
In the first step of the computational procedure, the IOP is applied using the initial linear properties (D0), which are artificially increased by a suitable factor (a factor of 100 is used). The purpose of this first computational step is to determine a set of in situ stresses σ1 in the cornea, which equilibrate the IOP, without producing appreciable deformations. The deformed shape of the cornea at the end of the first computational step is therefore close to that described by Eq.(2) as desired.
Fictitious strains ɛ1 in this first step are computed using the artificial initial stiffness. The strains ɛ1iter. that would be consistent with the actual nonlinear constitutive model and the computed initial stresses σ1 are determined at each material reference point via an iterative (Newton-Raphson) numerical process.
The applanation of the cornea is then modeled in the second computational step via a series of displacement controlled load increments, now using the actual corneal constitutive model, with total strains corrected by the strains ɛ1 and ɛ1iter. determined from the first step as summarized in the following pseudo code:
If step equals 1, then
- σ1= kD0ɛ1, where k = 100
- store (ɛ1, σ1) at the end of the first step
- σ1 → ɛ1iter. by Newton-Raphson iteration.
If step equals 2, then
- ɛn+1= ɛn+ɛ1iter. −ɛ1
- σn+1= Ds(ɛn+1)ɛn+1
- dσn+1= Dt(ɛn+1)dɛn+1.
Finite Element Simulation of Goldmann Applanation Tonometry
In the full three-dimensional finite element model of the cornea6 that included the anterior and posterior chambers separated by the lens, the interaction between the aqueous humor and cornea was simulated with hydrostatic fluid elements. The aqueous humor was assumed to be incompressible with fluid mass density of 1000 kg/m3. The fluid elements provided the coupling between the deformation of the cornea and the pressure exerted by the aqueous humor on the anterior surface of the cornea. Aqueous humor secretion was simulated with the fluid link element by means of a functional relationship between the aqueous humor outflow rate (q) and the pressure difference (Δp) across the fluid link element, defined by q = (0.28μl/min/mmHg) ×Δp. Applanation simulations showed that the IOP was only slightly elevated above its initial value due to displacement of aqueous humor during the applanation. Because the IOP elevation was negligible in practical terms as Goldmann had postulated, the full three-dimensional model is simplified by using the current axisymmetric cornea where the IOP is applied through fluid elements without modeling of aqueous humor secretion.
Schwartz et al17 estimated the pressure drop across the tear meniscus between the cornea and the edge of the applanator. The following calculation is based on their equations. The pressure drop (Δd) across the tear surface is illustrated schematically in Figure 2.
where γ is the surface tension constant, 0.00005 N/mm; Rt
is the radius of the fluid surface, 1.52 mm (applanation radius); and rt
is the tear radius at the edge of the wetted area.Rc
is the radius of curvature of the anterior surface of the cornea, taken as 7.77 mm from the corneal geometry.
Figure 2. Tear-Cornea-Tonometer Interface.17 Rt = Radius of the Fluid Surface, rt = Tear Radius at the Edge of the Wetted Area, Rc = Radius of Curvature of the Anterior Surface of the Cornea
The calculated pressure drop (=4.7 mmHg) is assumed to be the same for all corneas modeled in this study because the corneal radii of curvature are not considered as variables. The applanation procedure is simulated by contact analysis with a rigid surface representing the applanator. The applied force (expressed as a pressure) is calculated when the cornea is applanated over the contact radius (1.52 mm) and the pressure drop 4.7 mmHg is subtracted to obtain the (simulated) measured IOPs. Unless otherwise stated, the values of simulated IOP in this study are corrected by the pressure drop Δd=4.7 mmHg to take into account the effect of the tear meniscus. Additional details of the modeling procedure can be found elsewhere.6,7
This contact analysis results in the central portion of the cornea being flattened as required. Previous models9,10 that have used a closed form (thin shell) analytical expression for the corneal deflection under applanation implicitly assume that the contact pressure between applanator and cornea is uniform. Velten et al18 showed that the assumption of uniform contact stress produced significantly different results from a displacement controlled contact analysis. The final deformed shape of an average cornea with CCT of 550 μm is shown in Figure 3.
Figure 3. Deformed Shape of Applanated Cornea with Undeformed Shape Superimposed (dotted Lines). Central Corneal Thickness = 550 μm, Einit = 0.23 MPa, Intraocular Pressure = 16 mmHG
Virtual Population of Cornea Finite Element Models
A virtual population of cornea finite element models is generated for simulation. It is intended that this virtual population be roughly comparable to the 181 cases reported in a clinical study by Bhan et al,19 with a reasonably similar scatter in the correlation of the apparent IOP with CCT. The virtual population of 181 different cornea models is generated by randomly sampling three variables: CCT, (true) IOP, and the material properties (Einit). It is recognized that there may be variables other than those examined herein that affect IOP measurement. For example, cornea radius of curvature is believed to have a minor effect on applanation reading as pointed out by Ehlers et al.2 It is assumed that the three variables chosen in this study are the ones with the most significant effect on tonometric IOP readings.
The statistical distribution of CCT is taken as 551.53±49 μm, corresponding to Bhan et al’s clinical measurements.19 Well-documented data are also available on the statistical variation of (measured) IOP. Pooled data from large epidemiological studies indicate that the mean population IOP is 16 mmHg, with a standard deviation (SD) of 3.0 mmHg.20 The data on IOP typically show a non-Gaussian distribution. However, a normal distribution is used here based on the assumption that it will satisfactorily represent the variation of IOP for the purposes of this study.
As described previously, Einit is the only independent material parameter in the Extended Fung Model formulated herein. Also, it is observed that Einit ranged from 0.12 to 0.35 MPa when calibrated using the experimental data of Bryant and McDonnell.5 This substantial variation is used to establish a mean and SD of Einit for the virtual population. The middle value (0.23 MPa) is selected for the mean of Einit and 0.04 MPa for the SD in a way that the range of calibrated Einit (0.12∼0.35 MPa) is covered by mean±3 SD.
The three variables are assumed to have normal distributions and to be uncorrelated; the mean CCT is 551.52±49 μm, IOP is 16±3.0 mmHg, and Einit is 0.23±0.04 MPa.
Bhan et al19 report the IOP measurements by Goldmann applanation tonometry from 181 normal corneas of 94 patients attending a general eye clinic over a 6-month period. The same size sample of 181 cornea finite element models is generated with the statistical distribution of the three variables (CCT, IOP, and Einit). The comparison of the results of finite element simulations with the clinical measurements, and the effect of each variable on the IOP measurements with Goldmann applanation tonometry, are reported. The finite element simulations of Goldmann applanation tonometry are performed using Abaqus software package 6.3.1 (Hibiit; Karlsson & Sorensen Inc, Pawtucket, Rhode Island).
IOP Statistics in the Virtual Population and Clinical Study
Figure 4 shows apparent IOP as a function of CCT for the virtual population of 181 cornea models, determined from the finite element simulation of Goldmann applanation tonometry. Figure 5 shows the digitized version of the IOP measurements with Goldmann applanation tonometry by Bhan et al19 in a clinical study of 181 normal corneas of 94 patients. The trend lines determined from regression analysis of the data are also shown in Figures 4 and 5, respectively.
Figure 4. Apparent Intraocular Pressure (IOP) via Simulation of Goldmann Applanation Tonometry for the Virtual Population of 181 Finite Element Cornea Models. The Square of the Correlation Coefficient (R2) Represents the Strength of the Linear Association Between the Two Variables (central Corneal Thickness [CCT] and IOP).
Figure 5. Measured Intraocular Pressure (IOP) via Goldmann Applanation Tonometry in 181 Normal Corneas of 94 Patients.19 The Square of the Correlation Coefficient (R2) Represents the Strength of the Linear Association Between Two Variables (central Corneal Thickness [CCT] and IOP).
Figures 4 and 5 show a reasonably similar distribution of IOP versus CCT for the virtual population and the clinical study. Considering the idealizations that have been made in the modeling of the corneal geometry and biomechanical properties and that only three variables (Einit, CCT, and true IOP) have been sampled to generate the virtual population, the similarity between Figures 4 and 5 is considered to be satisfactory. Numerical values of the statistics of the two samples are given in Table 1. The range of CCT values and the trend lines in IOP versus CCT are similar, although the clinical measurements show a somewhat higher degree of scatter. In the virtual cornea population, the scatter is due solely to unaccounted variation in true IOP and corneal material properties as represented by Einit, whereas in the clinical population, there may be additional variables21 involved.
Table 1: Statistics of Intraocular Pressure Versus Central Corneal Thickness Plots
Effect of Einit and True IOP on Measured IOP
Variation of calculated IOP as a function of the other two variables, Einit and true IOP, in the virtual population is shown in Figures 6 and 7, respectively. Both figures show a considerable scatter in the data, similar to Figures 4 and 5.
Figure 6. Apparent Intraocular Pressure (IOP) in the Virtual Population of 181 Finite Element Cornea Models as a Function of Cornea In-Plane Young’s Modulus Einit.
Figure 7. Apparent Intraocular Pressure (IOP) in the Virtual Population of 181 Finite Element Cornea Models as a Function of True IOP (IOPT).
Figure 6 demonstrates that the actual range of randomly sampled Einit values in the virtual population of 181 corneas is approximately 0.11 to 0.31 MPa. This range is roughly the mean ±2 SD, and essentially accounts for the scatter in the IOP versus CCT correlation for the virtual population (see Fig 4) because the correlation would be a single curve without any scatter when Einit is fixed as a constant. It also may be inferred from a comparison of Figure 4 with Figure 6 that the variation in corneal biomechanical properties (as represented in the virtual population by the parameter Einit) is as important as CCT in determining IOP by applanation tonometry.
To emphasize further the dependence of the scatter in the plot of IOP versus CCT on the variability in Einit and true IOP, two subsets of the virtual population are selected that include only the finite element cornea models whose Einit and true IOP fall within the mean±SD(σ) and mean±½SD(σ), respectively. The data for these two subsets are shown in Figures 8 and 9.
Figure 8. Data for the Subset of the Virtual Population with Einit and True Intraocular Pressure (IOP) Within mean±SD(σ). CCT = Central Corneal Thickness
Figure 9. Data for the Subset of the Virtual Population with Einit and True Intraocular Pressure (IOP) Within mean±½SD(σ). CCT = Central Corneal Thickness
These figures show that the scatter in measured IOP diminishes as the range of Einit and true IOP narrows. It is also interesting to note that the slope of the trend lines is reduced as the range of the variation of the other two variables narrows. Numerical values of the slope of the trend line and the SD from the trend line are given in Table 2.
Table 2: Numerical Values of the Statistics for Different Subsets of the Virtual Population
It has been demonstrated by means of realistic finite element simulations that the observed large scatter in standard correlations of clinical measurements of IOP versus CCT can be largely accounted for by reasonable individual variations in corneal biomechanical stiffness properties. Although previous studies by others have explored this claim via simple analytical models, the present study has used what is believed to be a more realistic geometrical and biomechanical model of the cornea as well as a more accurate simulation of Goldmann applanation tonometry. An important feature of the corneal biomechanical model used herein, which has been calibrated using experimental results, is that the cornea stiffens under increasing physiological stress (due to IOP). Thus, the true IOP becomes a significant variable influencing the error in measured IOP.
From the results of the finite element simulations on a virtual population of 181 corneas, it is found that the corneal biomechanical properties (as characterized in this study by the stiffness parameter Einit) are as important as the CCT in influencing measured (Goldmann) IOP. Plausible individual variations in corneal biomechanical properties may account for the observed large scatter in standard correlations of clinical measurement of IOP versus CCT. In addition, it is important to note that these findings are supported by the finite element simulations rather than experimental measurements in this study.
Unfortunately, no methods are currently available for assessing corneal biomechanical properties in vivo, and thus individual variations in these properties cannot be accounted for at present. It is evident that the determination of corneal biomechanical properties has important implications for improved design of individualized laser surgery procedures as well.
The authors have developed a new method for using the results of applanation tonometry for determining corneal biomechanical properties in vivo, as well as correcting measured IOPs. This method will be presented in a forthcoming publication.