Adaptive model of the aging emmetropic eye and its changes with accommodation

Abstract

A general schematic model of the optical system of the emmetropic human eye is proposed, capable of adapting to changes with age and accommodation through adjustment of the optical surfaces and the internal gradient index structure of the lens. The specific models of the cornea and lens consist of minor generalizations of previous work by assuming them to be the sum of a biconic plus three higher order Zernike modes. The internal gradient index distribution adapts to the external shape so that the analytical expression is invariant with the changes with age and accommodation. The model also includes tips, tilts, and decentrations of the surfaces according to experimental findings. The specific parameters of these models are either constants or functions of age and/or accommodation. The model is polychromatic, and its optical performance was evaluated along the keratometric axis. Chromatic aberrations (longitudinal and transverse), astigmatism, coma, trefoil, and spherical aberration show good agreement with experimental averages. The change of these aberrations as a function of age or accommodation is also consistent with experimental findings (except for trefoil in eyes older than 50 years). This means that the average structure seems to predict the average performance. Nevertheless, the present model is too schematic to account for other higher order aberrations, such as tetrafoil, also present in real eyes.

Introduction

Models can help to predict function from structure and hence gain insight on the links between these two totally different aspects of the same system. In this way, realistic models are capable of gathering important information and can be useful in a variety of applications. In particular, there are a considerable number of eye models published in the literature. A full review of models is beyond the scope of this work, but specific reviews can be found in Atchison and Smith (2000) or in Navarro (2009), among others. Specific models were developed to account for the changes of the human eye as a function of age (Díaz, Pizarro, & Arasa, 2008; Goncharov & Dainty, 2007; Smith, Atchison, & Pierscionek, 1992) and accommodation (Navarro, Santamaría, & Bescós, 1985; Norrby, 2005; Popiolek-Masajada & Kasprzak, 2002). 

Custom models were also developed for a variety of goals and applications, such as modeling normal (Navarro, González, & Hernández-Matamoros, 2006a) or pseudophakic (Rosales & Marcos, 2007) eyes, or studying specific properties of the eye combining experimental measurement and ray-tracing modeling (Artal & Tabernero, 2008; Faria-Ribeiro et al., 2014). Statistical eye models (Rozema, Atchison, & Tassignon, 2011) seem an especially interesting approach as they allow implementing Monte Carlo simulations among other applications. 

The goal of the present study was to develop a generic model of the human eye, adaptive (Blaker, 1980) to the changes with age and accommodation. “Generic model” means that its structure is based on average experimental data. The specific set of average parameters can be realized as a sort of canonical proportions. Therefore, the resulting model with this sort of ideal proportions may differ from real eyes, whose specific parameters (Navarro et al., 2006a) may differ quite significantly from these averages. The second goal was to study the performance of that sort of average or “canonical” eye and then address the question as to what extent the performance of the average structure matches the average performance observed experimentally (Liang & Williams, 1997; Salmon & van de Pol, 2006; Thibos, Hong, Bradley, & Cheng, 2002)—in particular, the trends observed as a function of age (Applegate, Donnelly, Marsack, Koenig, & Pesudovs, 2007; Artal, Ferro, Miranda, & Navarro, 1993; Brunette, Bueno, Parent, Hamam, & Simonet, 2003; Calver, Cox, & Elliott, 1999; Guirao et al., 1999) and accommodation (Cheng et al., 2004; He, Burns, & Marcos, 2000; López-Gil et al., 2008). 

In the development of this type of model, one faces several limitations. A key limitation is the lack of reliable or accurate data on specific parameters that have a relevant impact on optical performance. Examples of this are the external three-dimensional topography of the crystalline lens (Ortiz, Pérez-Merino, Gambra, de Castro, & Marcos, 2012) or its internal structure (Dubbelman, Van der Heijde, Weeber, & Vrensen, 2003), for which average values are available but knowledge on their changes with age and accommodation are still scarce and/or show high uncertainties. For this reason, four parameters of the lens, chosen among those having both higher uncertainties and higher impact on optical performance, were used as variables for fine tuning of the model. In particular, their changes as a function of age were adjusted to optimize the predictions about age-related changes in optical performance. A key methodological aspect was the trade-off between minimum complexity and maximum realism of the model. In this way, the resulting model is somewhat more complex than previous schematic eye models, but the incorporation of additional parameters was supported by experimental evidence and/or by a relevant impact of these parameters on optical performance. In its current version the model is limited to the average emmetropic eye, even though its generalization to myopic eyes seems straightforward (Atchison, 2006). The present version of the model is not wide angle, as its field of view is limited mainly by the diameter of 9 mm of the cornea model. Consequently, the optical performance was evaluated for central vision only. 

Eye model

Main assumptions and design principles

The model is based on average experimental values of both structure (geometry and refractive indexes) and function (optical performance). The great majority of the parameters in the model were taken from either experimental studies or previous models. Instead of trying to average data across different studies, which could sometimes be difficult and risky due to different experimental methods and conditions, the model is based mostly on results of specific studies, especially well adapted to the goals of this work. The total number of nonzero parameters is 46, many of which change with age and/or accommodation. From these parameters, only four were used as variables for fitting the optical performance of the model: the central index of the lenticular gradient index (GRIN), n₀; the two conic constants of the isoindical surfaces of the lens, q_a and q_p, used for fine adjustment of both the power and the spherical aberration (SA); and the tip angle of the iris and lens (similar to angle κ), slightly adjusted to fit the change of coma as a function of age. The choice of these four variables obeys two main criteria: (a) their values have a critical impact on optical performance, and (b) they are especially difficult to determine experimentally, showing large uncertainties. However, special care was taken to keep the values of these four parameters within ranges consistent with experimental data. 

The model was built along the keratometric axis, which is used as the Z axis in the reference coordinate system. In a first approximation this axis can be considered approximately parallel to the line of sight (Atchison & Smith, 2000). The angle α between the visual and optical axes is greater than or equal to the angle between the keratometric and optical axes. The foveal field angle, between the nominal visual axis and the Z (keratometric) axis, is a variable of the model. The simplest zero value was assumed in most computations of optical performance, which means that the point object lies on the Z axis. 

When possible, additional simplifying assumptions were adopted. For example, the axis of corneal astigmatism (angle γ) was assumed to be “with the rule” (the cornea is flatter along the horizontal meridian and γ = 0) and constant with age. It is well known that “with the rule” is not the rule in the elderly (Gudmundsdottir et al., 2000), but the population average of γ shows only modest changes with age, with values not far from the rule (Navarro, Rozema, & Tassignon, 2013). Other simplifying assumptions are based either on the lack of experimental data or on the low reliability of those data. For example, it was assumed that the external surfaces as well as the internal isoindical surfaces of the lens (i.e., surfaces of equal refractive index) are concentric and aligned with each other. Recent data (Bahrami et al., 2014) suggest that this is not the case for individual eyes, but it seems reasonable to keep these simplifying assumptions for a canonical (average) model. In addition, the lens is assumed to be aligned with the iris so that they share a common iris–lens optical axis. Experimental evidence suggests that all these assumptions are not totally realistic (Chang, Wu, & Lin, 2007; Ortiz et al., 2012), but data are still too scarce and show large uncertainties, so choosing realistic values becomes difficult. In these cases, the model assumes the simplest configuration. 

In summary, the model design was the result of a trade-off between minimum complexity and maximum fidelity to both structure geometry and optical performance. Complexity was added to the model only when (a) experimental evidence suggested that such complexity was real and (b) it had a positive impact improving model predictions of experimental optical performance. In this way, the resulting model could yet be considered schematic in the sense that its level of complexity is lower than that of real eyes. 

The range of ages considered is from 20 to 70 years. This is a model of the emmetropic eye. Consequently, the model is assumed to be constantly emmetropic within this range of ages regardless of cross-sectional population data that suggest that, on average, younger (around 20 years old) people are more myopic than older (around 70 years old) people (Saunders, 1981). 

The maximum stimulus vergence used to study accommodation (for young eyes) was +8 diopters (D). The model is polychromatic, but both the power and the monochromatic wave aberrations were computed for a 555-nm wavelength along the axis parallel to the keratometric axis passing through the pupil center (as an approximation to the line of sight). 

Cornea

The cornea model corresponds to the result of linear regressions as a function of age of the main parameters (curvature radii, conic constants, orientation, significant Zernike deformations, and so on) of the anterior and posterior surface topographies from a recent study (Navarro et al., 2013). The model for each corneal surface is a general asphere comprising the sum of a biconic surface plus three Zernike polynomials, which were found to be significantly different from zero. In its canonical form (i.e., aligned and centered with the Z axis) the sag is given by  

 

 

Therefore, the surface shape is described by seven parameters: two curvatures (c_x = 1/R_x, c_y = 1/R_y), two conic constants (Q_x, Q_y), and three Zernike coefficients that mean higher order (aspherical) deformations modes of the surface sag. They are vertical trefoil Display Formula (Navarro, González, & Hernández, 2006b) plus fourth Display Formula and sixth Display Formula order rotationally symmetric aspheric modes. In addition to the shape, it is necessary to consider the three-dimensional position (x₀, y₀, z₀) and orientation (given by three Euler angles: α, β, γ) with respect to the reference coordinate system. As we said above, here we consider γ ≈ 0 and z₀ = 0 so that the total number of parameters describing the shape and position of each corneal surface is 11.  

Another parameter is the corneal central thickness (t_co) measured along the Z axis. Thus, the thickness associated with the posterior surface will correspond to the aqueous humor (t_aq). Most of these corneal parameters are linear functions of age (Atchison et al., 2008) and are assumed to be constant with accommodation. Although there is evidence for small temporal changes of the corneal topography with accommodation (He, Gwiazda, Thorn, Held, & Huang, 2003), their impact on the overall optical performance does not seem to be important enough to add complexity to the model. These corneal parameters and their changes as a function of age are listed in Table 1 for a corneal diameter of 9 mm. The case of vertical trefoil requires further explanation. The value −1.1 μm was the (significant) average over 113 topographies of the anterior surface measured with an Orbscan II instrument (Bausch & Lomb Surgical, Rochester, NY) (Navarro et al., 2006b), while the rest of the data listed in Table 1 correspond to a more recent study (Navarro et al., 2013) where both anterior and posterior topographies were measured with a Pentacam topographer (Oculus Optikgeräte, Wetzlar, Germany). The average results of both studies were consistent, except for Display Formula , which was found not to be significantly different from zero with the Pentacam data. The decision to keep the average value of the former study was based on experimental evidence (Read, Collins, Iskander, & Davis, 2009) that the relatively coarse angular sampling (i.e., limited number of radial slices) of the Pentacam might limit the accuracy in the reconstruction of high angular frequencies (trefoil, tetrafoil, and so on). For the anterior surface it was considered constant due to the lack of reliable data on its changes with age, and for the posterior surface it was assumed to be zero also due to the lack of data.  

Table 1

 

Model of the corneal surfaces for the central 9-mm diameter as a function of age in years. The thickness (t) is the distance to the next surface, i.e., t_co and t_aq are the axial thickness of the cornea and the aqueous humor respectively; hence, the bottom right expression corresponds to the thickness of the aqueous humor. D is accommodation stimulus in diopters.

Cornea	Anterior	Posterior
Rx (mm)	7.85 − 0.0051 × Age	6.20
Ry (mm)	7.72 − 0.0046 × Age	5.99
Qx	−0.367 − 0.0009 × Age	−0.56
Qy	−0.294 − 0.0023 × Age	−0.535 − 0.0026 × Age
(μm)	−1.1	0
(μm)	−0.026 + 0.030 × Age	6.02 − 0.046 × Age
(μm)	0.171 − 0.016 × Age	−2.37 + 0.013 × Age
x0 (μm)	106 + 4.393 × Age	317 + 4.031 × Age
y0 (μm)	−28.9 − 3.742 × Age	−295 − 3.373 × Age
α° (tilt)	−0.31 − 0.024 × Age	−2.51 − 0.024 × Age
β° (tip)	−0.71 − 0.036 × Age	−2.23 − 0.040 × Age
tco, taq (mm)	0.55	3.87 − 0.010 × Age − D × (0.048 − 0.0004 × Age)

Iris, lens, and retina

The iris was assumed to lie on a plane in front of the crystalline lens (0-mm thickness). Its image through the cornea, the entrance pupil, is decentered and tilted, and its axis forms an angle κ with the visual axis. In general, angle κ seems to have two horizontal (tip) and vertical (tilt) components. Recent studies (Berrio, Tabernero, & Artal, 2010; Chang et al., 2007) suggest an average value of around 3.7° horizontal tip and a much lower (∼0.1°) vertical tilt. Both x and y components show large experimental uncertainties (test–retest differences >1°) (Chang et al., 2007) and large intersubject variability (SD >2°). Here, we adopted two simplifying assumptions: (a) neglecting tilt and (b) considering that κ ≈ β (angle with respect to the Z axis; i.e., keratometric axis). Although an initial value of 3.7° was considered, we decided to use it as a model variable to fine tune the changes of coma as a function of age. This resulted in a quadratic increase of angle κ with age (see Table 2), but its range of values (from −3.50° to −4.65°) remained within the mean ± standard deviation of experimental data. The sign is reversed in the model because the sign of the Z axis follows light propagation toward the retina. 

Table 2

 

Geometry, position, and orientation of the iris and lens as a function of age and accommodation. D is accommodation stimulus in diopters. The axial thicknesses, t_le and t_vi correspond to lens and vitreous; t_a, t_p are the thicknesses of the anterior and posterior parts of the lens; and t_co and t_aq are the thickness of the cornea and the aqueous humor respectively. Note: β is approximately equivalent to angle κ.

	Iris–lens	Lens anterior	Lens posterior
c = 1/R(mm−1)	Infinity	(12.7 − 0.058 × Age)−1 + 0.0077 × D	−(5.9 − 0.013 × Age)−1 − 0.0043 × D
Q (external)	—	−4 − 0.5 × D	−3
q (internal)	—	−0.1480 × Age + 1.458 × 10−3 × Age2	0.422 − 0.131 × Age + 1.295 × 10−3 × Age2
x0 (μm)	0.113
β° (tip)	−4.18 + 0.046 × Age −7.5 × 10−4 × Age2
ta, tp (mm)	—	ta = 0.6 × tle	tp = 0.4 × tle
tle, tvi (mm)	0	2.93 + 0.0236 × Age + (0.058 – 0.0005 × Age) × D	23.465 − tco − taq − tle

The iris decentering was taken from the same cornea study (Navarro et al., 2013), as the topographer also provided the estimated coordinates of the pupil center. The average values were x₀ = 0.113 mm and y₀ = 0.015 mm, and they did not change significantly with age. The vertical displacement was neglected. Thus, y₀ and α were neglected so that only x₀ and β remain different from zero, and consequently only these two parameters are listed in Table 2. As explained above, lens and iris are considered to be mutually aligned so that the displacement x₀ and tip β apply to the whole iris–lens body. This approximation is not entirely realistic (Chang et al., 2007; Ortiz et al., 2012), but experimental data still seem too scarce to warrant more complex modelling. 

The crystalline lens GRIN was a rotationally symmetric (R_x = R_y = R; Q_x = Q_y = Q) adaptive model reported before (Navarro, Palos, & González, 2007) but with a slight generalization that allows a simultaneous fine tuning of power and SA. In the former model, the isoindical surfaces of the lens were assumed to be concentric with the external surface, meaning that they are scaled versions of the external surfaces. In the present version, this assumption is somewhat relaxed by allowing that the conic constants of the internal isoindical surfaces (q_a and q_p) can differ from those of the external surfaces (Q_a and Q_p). If q_a > Q_a and q_p > Q_p. This results in two thin peripheral zones (anterior and posterior) with zero axial thickness and constant refractive indexes equal to that of the lens surface. 

The external geometry of the lens model is based on a compilation (Norrby, 2005) of experimental data by Dubbelman, Van Der Heijde, and Weeber (2005) on the changes of the aging human lens with accommodation. The corresponding parameters and functions of age and accommodation are listed in Table 2. As explained above, the conic constants of the external surfaces Q are the same as in the former model, whereas the conic constants of the internal isoindical surfaces q were used as variables for fine fitting. Thus, the expressions for q_a and q_p listed in Table 2 are the results of model adjustment. The parameters t_a and t_p are the axial thicknesses of the anterior and posterior parts of the lens, respectively. 

There is a general agreement (Ooi & Grosvenor, 1995) that the axial length (i.e., the total thickness of the adult emmetropic eye) does not significantly change with age. Consequently, here we adopted a constant value of t_eye = 23.465 mm that is the average (range between 21 and 26 mm) found in a group of around 1,000 emmetropic eyes (Elabjer, Petrinović-Doresić, Durić, Busić, & Elabjer, 2007). 

In several earlier wide-angle eye models (Escudero-Sanz & Navarro, 1999; Goncharov & Dainty, 2007) the retina was assumed to be a 12-mm-radius spherical surface with its center lying on the Z axis. This would not be consistent with the present emmetropic eye model as it would mean that the diameter of that sphere is larger than the inner diameter of the eyeball (<23 mm). Experimental magnetic resonance imaging (MRI) measurements (Atchison et al., 2004) suggest that the retina can be approximated by a decentered, prolate, nonrevolution ellipsoid, which is in agreement with other recent studies (Faria-Ribeiro et al., 2014). For emmetropic eyes, the average radius (semiaxis) along the Z axis was found to be −11.45 mm, which seems fully consistent with the axial length adopted for the model. Therefore, this was the value assumed for the radius of the retinal surface. For the present on-axis model, the simplest assumption that the retina is a sphere aligned to the Z axis was adopted. 

Refractive indexes and chromatic dispersions

The refractive indexes and chromatic dispersions for the homogeneous ocular media were adapted from table 5 in Atchison and Smith (2005). With the exception of the lens, these values represent a minor but convenient improvement of earlier work (Navarro et al., 1985). In that table, Atchison and Smith (2005) listed the coefficients (A, B, C, and D) of the Cauchy equation for the change of refractive index as a function of wavelength:  

 

 

This reassigns the refractive indexes values for n_d = n(λ = 587.56 nm), used in classic eye models, to n₅₅₅ corresponding to the maximum spectral sensitivity. The adaptation made here was to turn back this action and consider that the values assigned to n₅₅₅ in that work should correspond to n_d instead. This adaptation affected only coefficient A, which was recomputed accordingly. The coefficients A, B, C, and D as well as the n_d and n₅₅₅ values are listed in Table 3 for the different ocular media. The indexes of the GRIN lens were also adapted from the former lens model (Navarro et al., 2007), which were based on experimental data (Jones, Atchison, Meder, & Pope, 2005). Again, the adaptation was reassigning values for n₅₅₅ to n_d. The result (see Table 3) is a slight increase of the n₅₅₅ values. 

Table 3

 

Coefficients of Equation 2, n_d, and n₅₅₅ for each ocular medium adapted from Atchison and Smith (2005). The coefficient A (n_∞) for the lens center—and hence n₅₅₅ and in general n(λ)—is a quadratic function of age in years.

	A (n∞)	B	C	D	nd	n555
Cornea	1.362994	6.009687 × 103	−6.760760 × 108	5.908450 × 1013	1.376	1.3774
Aqueous	1.323031	6.070796 × 103	−7.062305 × 108	6.147861 × 1013	1.336	1.3374
Lens center	See below	6.521218 × 103	−6.11066 × 108	5.908191 × 1013	1.4181	See below
Lens surface	1.356086	6.428455 × 103	−6.023738 × 108	5.824149 × 1013	1.3709	1.3726
Vitreous	1.323757	5.560240 × 103	−5.817391 × 108	5.036810 × 1013	1.336	1.3374
GRIN lens
Lens center	n∞ = A = 1.40965 + 3.55 × 10−4 × Age − 7.5 × 10−6 × Age2			n555 = 1.4264 + 3.55 × 10−4 × Age − 7.5 × 10−6 × Age2
Exponent	p = 2.85 + 1.1 × 10−7 × Age4

The distribution of refractive index is given by the following expressions for the anterior and posterior parts of the lens (see Navarro et al., 2007, for further details):  

 

and  

 

 

where δ_n is the difference between surface and central index n₀; ω² = x² + y²; a and b—with a² = R²/(q + 1) and b² = R²/(q + 1)²—are the semiaxes along z and ω, respectively, of the more external isoindical surface (anterior or posterior); for ellipses ε = +1 and for hyperbolas ε = −1; t is the axial lens thickness; Δ is the displacement along the Z axis required to place the origin of coordinates at the anterior apex: Δ_a = ε_aa_a and Δ_p = t – ε_pa_p; f_a = Display Formula and f_p = Display Formula are normalization parameters so that the argument (normalized radius) is always within the interval [0, 1]; and t_a and t_p are the axial thicknesses of the anterior and posterior parts of the lens, respectively: t = t_a + t_p (see Table 2). The coordinates (z_a,ω_a), (z_i,ω_i), and (z_p,ω_p) follow the equations of three conicoids corresponding to the anterior, internal interface, and posterior surfaces, respectively; hence, they represent the boundaries of validity of Equations 3a and 3b. The equation of the conicoid interface separating the anterior from the posterior parts of the lens is  

 

 

This surface is the locus formed by the intersections of the anterior and posterior parts of every isoindical surface. This adaptive lens model was copied by other authors (Chen, Jiang, Yang, & Sun, 2012). A similar model, using third-order (instead of second-order) isoindical surfaces thus permitting to keep continuity in the first derivatives, was proposed (Bahrami & Goncharov, 2012a) and is expected to provide similar optical performance. 

As explained above, the central index n₀ was one of the four variables used for fine adjustment of the model, whereas the surface index was considered constant with age. Moffat, Atchison, and Pope (2002) found a mean rate of decay of 3.4 × 10⁻⁴ per year, which predicts a total change in n₀ of −0.017 (from 1.431 to 1.414) between 20 and 70 years, which is the value assumed here. The only difference is that a finer tuning of the model was required here, considering a quadratic rather than a simple linear function. On the other hand, previous aging eye models assumed similar rates of decline of n₀ with age (Koretz & Cook, 2001; Díaz et al., 2008) or a decline of the equivalent index of the lens (Dubbelman & Van der Heijde, 2001) in order to explain the lens paradox (Brown, 1974). Nevertheless, not only do experimental data seem to show a wide variability among samples, but the results also seem to vary depending on the methods, wavelength, and so on (Pierscionek & Regini, 2012). This means that the assumed values and rates of change represent a plausible solution, but this solution might not be unique, as it also happens in custom eye models (Navarro et al., 2006a) or in other parameters of the present model, also subject to large experimental uncertainties. 

The model of the chromatic dispersion of the GRIN structure is again based on a simplifying assumption, which is that both the surface and the center indexes are functions of the wavelength n_s = n_s(λ) and n₀ = n₀(λ). Thus, δ_n = n_s(λ) − n₀(λ) = δ_n(λ). These are Cauchy equations (Equation 2), and the corresponding coefficients are listed in Table 3. This assumption means that the index distributions of Equation 3 are also functions of the wavelength n_a = n_a(z,ω,λ) and n_p = n_p(z,ω,λ). The resulting optical models for 20- and 60-year-old eyes are compared in Figure 1. Other dispersive GRIN models also used simplifying assumptions, such as considering that the constant term in a polynomial GRIN model is a function of λ (Díaz et al., 2008) or using concentric isodispersive contours (Bahrami & Goncharov, 2012b). 

Figure 1

 

Ray tracing through 20- and 60-year-old eye models.

Implementation

Ray tracing computations were carried out using ZEMAX (Radiant Zemax LLC, Redmond, WA) using two types of user-defined surfaces developed in C++ language. One surface type implements the corneal surfaces of Equation 1, and the other type for the GRIN lens was implemented before (Navarro et al., 2007). In the range of 20 to 40 years, the sampling grid parameter was set to 0.05 mm to speed up the ray tracing computations. For older lenses with higher values of exponent p, a finer sampling (0.01 mm) was used to keep accuracy, which means slowing down the ray tracing computation by a factor of five. The optical surfaces were introduced in their canonical forms, and the “tilt/decenter elements” ZEMAX tool was applied three times for (a) the anterior corneal surface, (b) the posterior corneal surface, and (c) the group of five surfaces formed by the iris (one surface) and lens (two external surfaces and two internal surfaces for the GRIN structure). Taking into account that the last surface is the retina, the model consists of eight optical surfaces: 2 (cornea) + 1 (iris) + 4 (lens) + 1 (retina). In addition, ZEMAX uses two special “coordinate break” surfaces, one before and one after the element or group of elements tilted/decentered. In this way the coordinate system is transformed just before the ray impacts the surface, and it is transformed inversely just after the ray emerges from the element or group. This method is highly useful to keep a global and invariant reference system of coordinates (especially when there is no effective optical axis), but this method (apparently) increases complexity quite significantly. In the most general case, the number of surfaces increases by a factor of three. With the “coordinate break” surfaces (3 × 2), the total number of surfaces increases from 8 to 14. ZEMAX files (including an “eyemodel” glass catalog, dll libraries [for Microsoft 32-bit Windows (Microsoft, Redmond, WA)] for the cornea and lens user-defined surface types, and an implementation example) are available under request to the author. Their use is permitted for academic purposes only. 

Finite ray tracing was used in all cases, except for the position of the principal planes. The free parameters q_a, q_p, and n₀ of the GRIN distribution were adjusted to keep (sphere) emmetropia constant with age. Nevertheless, a finer tuning was obtained under the assumption that the younger (20 years old) eye has pure hypermetropic astigmatism but that it evolves continuously with age to the opposite pure myopic astigmatism for the older (70 years old) eye. This minor age-related change (since the magnitude of astigmatism is <0.4 D) is consistent with the assumption that not only does the eye stay emmetropic but that the refractive state is adapted to the decline of accommodative response. Even though this seems arbitrary, it would be consistent with the increase of the curvature of the cornea and lens surfaces (and power) with age. This assumption means that, for the young eye, the focus with the highest power lies on the retina (emmetropic), whereas the second focus lies beyond the retina (hypermetropic). This makes sense since this allows the young eye to focus at any point within the Sturm interval (i.e., the segment along the Z axis between the two astigmatism foci, also called Sturm foci) using accommodation. In older eyes, where the accommodation response has diminished, myopic astigmatism might permit some level of pseudoaccommodation by having the complete Sturm interval fully available for effective depth of field. This is further discussed later. 

Refractive errors were computed for monochromatic light (λ = 555 nm) and a 3-mm entrance pupil with the criterion of minimum root mean square (RMS) size of the spot diagram. When the SA is positive, under this metric the best focus lies between the paraxial and minimum RMS wavefront error foci, which seems consistent with experimental findings (Martin, Vasudevan, Himebaugh, Bradley, & Thibos, 2011). The refractive state was later confirmed using the criterion of maximum visual Strehl ratio (Guirao & Williams, 2003). 

Wavefront analysis and computation of Zernike coefficients were carried out for entrance pupils of 4 mm and 6 mm for studying the changes with accommodation and aging, respectively. The entrance pupil provides a good match with Hartmann-Shack wavefront sensors, which measure the outgoing wavefront at the exit pupil. But the outgoing exit pupil becomes the entrance pupil for ingoing wavefronts propagating toward the retina. 

Optical performance

Power and astigmatism

The refractive power of the eye model as a function of age is given in Figure 2. The dashed red line is the ideal power of the eye to keep the image focused on the retina, which ranges from 62.25 D to 62.84 D. This small increase of nearly 0.6 D from 20 to 70 years is due to the continuous growth of the lens, which causes the distance from the image principal plane to the retina (i.e., the focal length) to decrease with age at a mean rate of about 4.5 μm/year. Figure 2 compares the power corresponding to the paraxial focus (black open triangles), the circle of least confusion (CLC; blue line), and the two upper and lower bounds of the Sturm interval (green circles and purple squares, respectively). The distance between these two bounds (Sturm foci) gives the geometrical astigmatism. The gap between the paraxial focus and the CLC is small (0.12 D) for young eyes, but it increases with age up to 0.4 D for the 70-year-old eye model due to the increase of higher order aberrations. As explained above, the criterion adopted for adjusting the power was based on assuming an optimal configuration for visual performance. As a result, for the younger eye the Sturm focus with maximum power matches the theoretical value (pure hypermetropic astigmatism), whereas we can observe that the theoretical value matches the opposite Sturm focal having minimum power (pure myopic astigmatism). The match with the CLC occurs for an intermediate age of about 40 years. This functional optimization of the model is further discussed later in the context of the lens paradox. 

Figure 2

 

Change of dioptric power of the eye model as a function of age. The values corresponding to the two Sturm foci (Max Sturm = maximum power; Min Sturm = minimum power), the CLC, and paraxial focus are compared with the theoretical ideal value necessary to keep the image focused on the retina.

Figure 3 compares the magnitude of geometric astigmatism (red dots) versus Zernike astigmatism (blue open circles) as a function of age. The geometrical astigmatism estimated from spot diagrams is slightly lower. The model predicts a small decline with age, possibly due to a partial compensation of corneal astigmatism by the increasing misalignment of the cornea and pupil–lens. This partial compensation seems consistent with experimental findings (Kelly, Mihashi, & Howland, 2004). The values predicted by the model (average 0.35 D) are somewhat higher than experimental averages of around 0.17 or 0.24 D (McKendrick & Brennan, 1996), but intersubject variability is higher (SD >0.6 D). 

Figure 3

 

Change of geometrical (based on spot diagram) and Zernike astigmatism of the eye model as a function of age.

Optical axis

The following analysis is restricted to horizontal tips. The analysis of vertical tilts would be equivalent, but angles are negligible in most cases. Figure 4 shows the tip angle necessary to align the anterior surface of the cornea (green diamonds) and the iris–lens (open red circles) with the Z (keratometric) axis. The tip angle estimated for the nominal optical axis of the eye model is also represented (blue line). The later was computed as the angle necessary to overlap the first and fourth Purkinje images, which is one of the usual ways to determine the angle α between the optical and visual axes (Atchison & Smith, 2000). Interestingly, the angle between the optical axis and the Z axis (“Eye” in Figure 4) seems somewhat more stable with age than the angles of the cornea and iris–lens axes that both increase monotonically as a function of age. The estimated average value of 2.2° is lower than the experimentally observed values of the horizontal component of angle α (3°–5°; dashed black lines in Figure 4). However, Mandell, Chiang, and Klein (1995) found an average shift angle between the keratometric axis and the line of sight of 2.7°. Adding this value to 2.2° we can estimate an average horizontal angle α of 4.9°, which is fully consistent with the experimental data. 

Figure 4

 

Changes of the tip angles of the optical axes of the anterior surface of the cornea (green diamonds), of the iris–lens (red open circles), and of the eye model (blue line) as a function of age.

Consequently, the optical performance of the eye model was compared for fields angles of 0° (along the Z axis; angle α ≈ 2.2°) and 2.7° (angle α ≈ 4.9°). The monochromatic performance showed little differences, which suggests that this field is within the isoplanatic patch of the model. However, the transverse chromatic aberration (TCA) showed a critical dependency with field angle (see below). 

Chromatic aberrations

The longitudinal chromatic aberration (LCA) for the 40-year-old eye and for a 3-mm pupil is shown in Figure 5. It was computed as the object vergence (inverse distance of object from the corneal vertex) providing the best focus on the retina. An intermediate age was chosen since the LCA seems to be constant with age (Howarth, Zhang, Bradley, Still, & Thibos, 1988). The model predictions (blue circles) are compared with the result of fitting experimental data to the Cauchy equation (Equation 2) (Atchison & Smith, 2005) (continuous green line) and with a chromatic reduced-eye model (Thibos, Ye, Zhang, & Bradley, 1992). The corresponding curves are slightly shifted in the vertical direction to a common origin (zero LCA) at 555 nm that is the reference wavelength in the present eye model. The agreement with average experimental data is excellent without adjustment of the chromatic dispersions of the eye model. 

Figure 5

 

LCA of the eye model compared with a Cauchy equation fit of experimental data (Atchison & Smith, 2005) and with the LCA of a reduced chromatic eye model (Thibos et al., 1992).

Contrary to the LCA, the TCA shows a large intersubject variability. Rynders, Lidkea, Chisholm, and Thibos (1995) reported values ranging from 0.05 to 2.67 arcmin, with an average of 0.8 arcmin, for the TCA between 497 nm and 605 nm. Here the TCA was estimated, for the 40-year-old eye model, as the (horizontal) offset between the retinal impacts of the chief rays for different wavelengths. As expected, the TCA changes dramatically with field angle. The model predictions agree with the experimental average, 0.8 arcmin, for a field angle of about 3° (with the Z axis of the eye model). As discussed above for the 40-year-old model, this means that the angle between the visual and optical axes would be about 5.2°. This value is not far from the value of 4.9° suggested in “Optical axis.” 

Wave aberration

The monochromatic wave aberration was computed for a 6-mm pupil (λ = 555 nm) for field angles of 0° and 2.7°. Since the differences were small, here we present the results for 0°. Figures 6, 7, and 8 compare model predictions with experimental data for fourth-order SA, coma, and total RMS higher order aberrations (HOAs), respectively. Among the various experimental studies available in the literature, the work of Applegate et al. (2007) was especially well suited for a direct comparison with the present model, and their results are fully consistent with those of similar studies. Experimental average data are represented by red open circles, and the error bars correspond to standard deviations. The red lines represent the experimental tendencies (second-order polynomial fit of the data points). Blue squares are model predictions. In the case of SA (Figure 6) the agreement is excellent. It is worth recalling that this agreement was obtained by a fine adjustment consisting of optimizing the values of the conic constants q_a and q_p of the internal isoindical surfaces of the lens. Their final values (see Table 2) do not differ much from experimental data and the differences with experimental mean are below 1 SD. 

Figure 6

 

SA as a function of age: model prediction versus average experimental data.

Figure 7

 

RMS values for coma and total third-order Zernike aberrations as a function of age: model prediction versus average experimental data.

Figure 8

 

Total RMS HOAs as a function of age: model prediction versus average experimental data.

Figure 7 shows RMS values of coma and total third-order RMS (coma and trefoil) Zernike aberrations. Only the trend curve is represented for the experimental RMS third-order aberrations for the sake of clarity (error bars are higher than those of coma). Model prediction of coma is slightly higher than experimental values. In this case the angle κ was adjusted to keep that small difference approximately constant with age. Again, the range of variation of this angle (β in Table 2) is within experimental uncertainties, and the average across ages is similar to the experimental mean. A better fit could be possible, but this was not a primary goal in this study. The total third-order RMS (green circles) shows an excellent agreement for ages between 20 and 40 years, while for older eyes the model shows an increasing underestimation. This means that the model seems to slightly overestimate the coma increase and underestimate the trefoil increase as a function of age. This can be easily explained since the main source of trefoil in the model is the trefoil of the anterior corneal surface. It was assumed to be constant with age mainly due to the lack of reliable experimental data. For the same reason trefoil was ignored for other surfaces. 

A similar explanation can be applied to the total RMS HOA in Figure 8. The model predictions lie within the range delimited by the error bars, but they are clearly lower than average experimental data. This underestimation increases gradually with age. A possible explanation is that apart from trefoil the model ignores other possible higher order contributions to the topographies of the lens surfaces. In the case of the cornea, there is a small average nonzero residual (Navarro et al., 2013), but the individual Zernike coefficients are not significantly different from zero. This residual should contribute to fourth and higher order wave aberration. In the case of the lens, the surfaces were assumed to be revolution conicoids. There is also evidence, both experimental (Ortiz et al., 2012) and predictions from custom eye models (Navarro, 2009), that the lens topography exhibits higher order Zernike modes—in particular, trefoil and tetrafoil. The complex topography of the lens seems consistent with the presence of suture line patterns, and it was suggested that these could be the source of starbursts in retinal images of bright points (Navarro, 2009). The present model does not take into account these complex lens topographies since experimental data are still scarce. 

Figure 9 shows the gradual change of the wave aberration and its corresponding point spread function as a function of age (20-, 40-, and 70-year-old eye models). The patterns show the dominance of astigmatism, coma, and SA. 

Figure 9

 

Wave aberration and its corresponding point spread function for 20-, 40-, and 70-year-old eye models.

Accommodation

Figure 10 compares the accommodation response computed for the younger (20 years old) eye model (red line) with experimental results corresponding to two emmetropic subjects taken from figure 1 in the study of Taylor, Charman, O'Donnell, and Radhakrishnan (2009). The agreement seems quite reasonable. As in the case of the LCA, it represents the object vergence providing the best focus on the retina for a 3-mm pupil. The predicted response (red line with no symbols) deviates from the linear or ideal response (y = x; dotted black line) and can be approximated by the expression y = 0.7387x − 0.0046x². It means that the slope is lower than unity and that the response is slightly decelerated. Both features seem compatible with experimental responses. Experimental curves often show larger bends (see, e.g., figure 3 in Wang & Ciuffreda, 2006), which suggest somewhat higher decelerations. The accommodation lag of the eye model is explained by the fact that in the empirical expressions for the changes of the lens shape as a function of accommodation (Table 2) the variable D represents stimulus vergence rather than accommodation response (Dubbelman et al., 2005). 

Figure 10

 

Accommodation response for the 20-year-old eye model (red line) compared with experimental results corresponding to emmetropic subjects (A and C) for different spatial frequencies. Linear means ideal response (y = x). Figure adapted from figure 1 in Taylor, Charman, O'Donnell, and Radhakrishnan (2009).

The SA as a function of accommodation response is compared with experimental data (López-Gil et al., 2008) in Figure 11. It was computed for a 4-mm pupil diameter since accommodation studies require photopic luminance levels of stimuli so that the natural pupil can hardly reach much larger diameters. (Furthermore, accommodation induces pupil miosis.) The model reproduces the well-known fact that the SA decreases with accommodation so that it becomes overcorrected (negative) for accommodations higher than around 2 D. The overall agreement is good. The experimental line is straight since it corresponds to a linear regression of experimental data points, whereas the model predicts an accelerated decrease of SA as a function of accommodation response. As in the case of Figure 8, the model cannot accurately predict the change of total RMS HOA as a function of accommodation. Experimental data show a monotonic increase with accommodation, whereas the model predicts lower values, showing a not monotonic change with local minimum at the zero crossing of SA. 

Figure 11

 

Zernike coefficient of SA as a function of accommodation response: model prediction (green) versus experimental data (red).

Discussion

The present model provides a reasonable match with the average values of refractive errors and aberrations (astigmatism, coma, SA, and chromatic aberrations) by fine tuning only four parameters. Several simplifying assumptions helped to keep the model complexity within reasonable limits, but as a result it is unable to predict the magnitude of other HOAs present in real eyes (mainly tetrafoil, secondary astigmatism, and sixth-order SA). In this sense the present model should be understood not as a sort of finished product but rather as advancement in the optical modeling of the eye. I tried to justify most of the assumptions, choices, and solutions adopted here, but alternative choices might be possible in most cases. Furthermore, different eye models could provide equivalent predictions of optical performance (Navarro et al., 2006a). As explained before, both the mathematical models and the biometrical parameters of cornea and lens were chosen attending to certain criteria. Data from different studies were not mixed to compute new averages, and hence the value adopted for each parameter corresponded to only one study. That study was chosen among those best adapted to the goals of the present work (adaptation to the mathematical models, ranges of ages, accommodation, and so on) Nevertheless, rather than being based on data from a single population, this model is a combination of independent submodels (cornea, lens, and so on) that are each based on different cohorts. The basic optical models of the cornea (Navarro et al., 2006b, 2013) and lens (Navarro et al., 2007) were developed previously, and the most relevant aspects of those models have already been discussed. In the case of the cornea, Scheimpflug imaging permits complete and precise measurements, and hence models can attain a high fidelity and realism. The in vivo measurement of the lens is much more difficult. The changes of the aging lens geometry with accommodation were studied in vitro (Glasser & Campbell, 1998) and in vivo. In the last case, Scheimpflug imaging (Koretz, Cook, & Kaufman, 2002) provided the experimental data used here (Dubbelman et al., 2005). More recently, powerful techniques such as optical coherence tomography (Ortiz et al., 2012) or MRI (Kasthurirangan, Markwell, Atchison, & Pope, 2011) were also applied to study the lens in vivo. MRI results seem consistent with the data of Dubbelman et al. (2005) used here. Regarding optical performance, there is a good agreement among different studies on both the average and standard deviations of the Zernike aberration modes in the human eye as well as their changes with pupil size, age, or accommodation. 

The optical performance of the model shows a reasonable agreement with experimental data for second-order aberrations (refractive errors), chromatic aberrations, coma, trefoil (for young eyes), and SA. However, the model is too schematic to account for the amounts of other HOAs also present in real eyes, such as tetrafoil, secondary astigmatism, and so on. Regarding refractive errors, the population histogram for spherical equivalent shows a sharp peak (histogram mode) at 0 D (Duke-Elder & Abrams, 1970), suggesting that more than 50% of the population is emmetropic. The present model is specific for emmetropic eyes, and consequently the model was adjusted to remain emmetropic in the range of ages (from 20 to 70 years) considered regardless of age-related changes in the prevalence of myopia and hyperopia. The generalization of the model to different refractive errors will be the subject of future work. 

The answer to the question of whether the performance of the average eye equals the average performance measured in the population is yes in this case. Since the model was built using canonic (average) proportions, we can say that the mean structure seems to predict the mean performance (with the exceptions outlined above, which can hardly be reproduced with this schematic model). However, a model based on canonic (average) proportions is not representative for individuals from the general population, who may have large variations in biometric values. In that case the generic (canonic) model could be used as customized to build personalized eye models (Navarro et al., 2006a). Another important aspect is that, like any other model of the aging eye in the literature, this model is not based on longitudinal data but rather represents a cross-section of the current population. Thus, we cannot completely exclude cohort effects, associated to gradual changes in the general population during the past 70 years, such as the increased prevalence of myopia in younger age groups. 

Another interesting question is the Brown's lens paradox: While the curvatures of the lens surfaces increase monotonically as a function of age, the eye does not seems to become more myopic. Furthermore, cross-sectional population studies suggest the opposite—that is, an age-related hyperopic shift (Saunders, 1981). Different authors offered explanations for this paradox. There seems to be a general agreement in that the equivalent (or effective) refractive index of the lens declines with age (Dubbelman & Van der Heijde, 2001). In GRIN measurements (Moffat et al., 2002) and models (Díaz et al., 2008) the paradox is explained by the decline of the central refractive index with age. The present model was kept constantly emmetropic also using the later assumption, but the so-called lens paradox seems to be a more complex multifactorial issue. On one hand, the cornea is also involved in the paradox since its apex curvature also increases as a function of age (Topuz, Ozdemir, Cinal, & Gumusalan, 2004) (see R_x and R_y in Table 1). As a result, the apical power of the anterior surface would increase more than 1 D (from 20 to 70 years) (Navarro et al., 2013). On the other hand, the total power of the eye must increase (about 0.6 D in the same range of ages) to keep the image focused on the retina. The explanation is that the distance from the image principal plane to the retina gets shorter due to the continuous increase of lens thickness. An additional factor requiring increased power was assumed here, which consists of the change from pure hypermetropic astigmatism in younger eyes to pure myopic astigmatism in older eyes. This would mean an additional increase of 0.35 D. This is only one hypothesis, but it helps to partially explain the increase in curvatures of cornea and lens with age. This assumption means that between the ages of 20 and 70 years the eye needs to increase its spherical power by 0.95 D (0.6 D + 0.35 D). Thus, from the total theoretical increase of power associated with the lens paradox (∼2.6 D), nearly 1 D could potentially be necessary (and real), whereas the remaining excess of power might be compensated for by a decline of the equivalent refractive index of the lens. 

The systematic study of the possible compensation of corneal aberrations by internal optics (Artal, Berrio, Guirao, & Piers, 2002; Artal & Tabernero, 2008; Kelly et al., 2004; Millodot & Sivak, 1979) is beyond the scope of the present work. Nevertheless, the case of SA was discussed before for both the corneal model (Navarro et al., 2013) and the lens model (Navarro et al., 2007). The cornea model predicted that corneal aberrations increase with age at a similar rate as the total ocular aberrations and that the balance between corneal and internal SA remains roughly constant with age. In the case of astigmatism, Figure 3 shows a decline of astigmatism with age, whereas both corneal and iris–lens tip angles increase. These two age-related changes suggest an increasing partial compensation of corneal astigmatism as a result of the increasing misalignment of the optical surfaces. The balance between corneal and internal optics for coma has not been studied yet. The systematic study of the sources of the main aberrations (Marcos, Burns, Prieto, Navarro, & Baraibar, 2001) in the model, as well as the balance between cornea and internal optics, will be part of future work. 

Acknowledgments

This research was supported by the Spanish Ministry of Economía y Competitividad and the European Union, grant FIS2011-22496, and by the Government of Aragón, research group E99. 

Commercial relationships: none. 

Corresponding author: Rafael Navarro. 

Email: [e-mail]. 

Address: ICMA, Consejo Superior de Investigaciones Científicas and Universidad de Zaragoza, Facultad de Ciencias, Zaragoza, Spain.