Customized models of ocular aberrations across the visual field during accommodation

Abstract

We aimed to create individual eye models that accurately reproduce the empirical measurements of wave-front aberrations across the visual field at different accommodative states, thus providing a mechanistic explanation for the changes in the eye's aberration structure due to accommodation. Structural parameters of a generic eye model were optimized using optical design software to account for published measurements of wave-front aberrations measured for 19 individuals at 37 test locations over the central 30°-diameter visual field at eight levels of accommodative demand. Biometric data for individual eyes were used as starting values and normative data were used to constrain optimizations to anatomically reasonable values. Customizations of the accommodating eye model accurately accounted for ocular aberrations over the central 30° of visual field with an averaged root mean square fitting error typically below 0.2 μm at any given field location. Optimized structural parameters of the eye models were anatomically reasonable and changed in the expected way when accommodating. Accuracy for representing spherical aberration was significantly improved by relaxing anatomical constraints on the anterior surface of the lens to compensate for not including gradient-index media. Use of the model to compute pan-retinal image quality revealed large penalties of accommodative lag for activating photoreceptor responses to the retinal image.

Introduction

The human eye is a dynamic optical system that accommodates to changes of viewing distance by changing its refractive power, thereby maintaining a clearly focused retinal image. This change in power is accomplished by changes in the shape of the internal crystalline lens brought about by activation of the ciliary muscle, which adjusts the tension in suspensory ligaments. To produce extra refracting power needed for close objects, the anterior and posterior curvatures of the lens increase (Dubbelman, Sicam, & Van der Heijde, 2006; Rosales, Dubbelman, Marcos, & van der Heijde, 2006), the lens thickens, and the anterior chamber shortens (Dubbelman, Van der Heijde, & Weeber, 2005). In addition, the aspherical shapes of the refracting surfaces of the lens change, which changes the higher order aberrations of the eye, especially spherical aberration (Dubbelman et al., 2005; López-Gil & Fernández-Sánchez, 2010). 

Optical changes of the eye produced by accommodation occur not only along the primary line of sight (LoS) associated with foveal vision but also along the secondary LoS associated with peripheral vision. For example, changes in refractive power measured foveally occur uniformly across the central 30° diameter of the visual field in response to changes in viewing distance of foveal targets (Liu, Sreenivasan, & Thibos, 2016; G. Smith, Millodot, & McBrien, 1988). Some other aberrations, such as astigmatism, are surprisingly robust to accommodation over the central field (Liu & Thibos, 2017; Lundström, Mira-Agudelo, & Artal, 2009; Mathur, Atchison, & Charman, 2009). The goal of the current study was to develop customized optical models that can account for such changes, and lack of changes, in lower order and higher order aberrations of individual eyes during accommodation. 

One approach to developing eye models is to measure the curvatures, shapes, and positions of the eye's four refracting surfaces using optical biometry techniques such as Purkinje images or Scheimpflug slit-lamp photographic images (Atchison et al., 2008; Brown, 1973; Dubbelman et al., 2005; Koretz, Cook, & Kaufman, 1997, 2002; Rosales et al., 2006). Measurements are also required for the retinal surface, typically obtained with low-coherence tomography (Atchison & Charman, 2011). A disadvantage of this approach is that internal surfaces are imaged by more anterior surfaces, which hampers the accurate measurement of more posterior surfaces. Another disadvantage is that some properties of the lens, such as gradient refractive index (Kasthurirangan, Markwell, Atchison, & Pope, 2008) and asphericity of posterior surfaces, are hard to measure directly (Dubbelman et al., 2005). Optical measurements of retinal shape and position are also subject to uncertainties that may reduce the accuracy of the model for predicting off-axis refractive errors (Atchison, 2004; Kooijman, 1983; Pomerantzeff, Pankratov, Wang, & Dufault, 1984). Although the lack of information about the crystalline lens could be circumvented for applications such as evaluating intraocular-lens implantation (Rosales & Marcos, 2007; Tabernero, Piers, Benito, Redondo, & Artal, 2006), taken together these difficulties encountered in the ocular-biometry approach to modeling may lead to uncertainties in the eye's optical properties derived from such models. 

An alternative approach to modeling inspired by computer-assisted tomography (CAT scans of internal body structures) uses wave-front-aberration maps obtained along multiple LoSs to determine the shape and position of the major refracting elements of an eye. This computational method, called ocular wave-front tomography (OWT), uses optical design software to optimize a generic eye model in order to account for wave-front aberrations measured over a significant portion of the visual field (Wei & Thibos, 2008). Previous solutions to this inverse problem of determining a system from its behavior have been reported for foveal vision (Kong et al., 2009; Lou et al., 2015; Navarro, González, & Hernández-Matamoros, 2006) and for an expanse of visual field (Goncharov, Nowakowski, Sheehan, & Dainty, 2008; Polans, Jaeken, McNabb, Artal, & Izatt, 2015) for eyes in a relaxed state of accommodation. We have extended that work into the domain of accommodation with an improved method that uses biometric data to constrain the search space of optimization, thus yielding models that are both functionally equivalent and anatomically reasonable. Our ultimate goal was to create individual eye models that accurately reproduce the empirical measurements of wave-front aberrations of a given eye at any accommodative state, thus providing a mechanistic explanation for the changes in the eye's aberration structure that occur over the central visual field during accommodation. 

Methods

We used the Navarro eye model (relaxed accommodation) with four refracting surfaces and uniform refractive index for the crystalline lens (Navarro, Santamaría, & Bescós, 1985) as the starting point for customization by the optical design program OpticStudio (Zemax, Kirkland, WA). The first step in customization was to alter the values of parameters derived from ocular biometry of individual eyes (e.g., corneal topography, axial length), described in detail later. The second step was to use the Zemax optimization engine to find the best values of model parameters (e.g., curvatures of anterior and posterior lens surfaces, anterior-chamber depth [ACD], lens thickness [LT]), so that the model's aberrations replicate empirical measurements of wave-front aberration at 37 locations in the central 30° diameter of the visual field at eight levels of nominal accommodative demand. The methodology of the empirical experiments has been described previously (Liu et al., 2016). Anatomical similarity of the model to human eyes is maintained (Wei & Thibos, 2008) by introducing constraints on the range of parametric variation used by the Zemax optimizer (Dubbelman et al., 2005; Kong et al., 2009; Rozema, Atchison, & Tassignon, 2011; Rozema, Rodriguez, Navarro, & Tassignon, 2016). 

Customized models were constructed for the left eyes of 19 healthy subjects (nine with emmetropia—mean age: 26 years; age range: 21–34 years; and 10 with myopia—mean age: 25 years; age range: 20–31 years) who participated in our empirical survey of aberration across the central visual field (Liu et al., 2016; Liu & Thibos, 2017). In addition, two experienced subjects with presbyopia (LNT, age 68 years, and AB, age 62 years) were recruited for a control group of eyes with zero accommodative amplitude. Automatic data management was achieved by letting Zemax be controlled programmably by MATLAB (MathWorks, Natick, MA) via ZOS-API (Radiant Zemax, 2014). 

Biometry measurement

Corneal topography was measured using the Medmont corneal topographer E300 version 4.12 (Medmont International Pty, Ltd., Victoria, Australia) with the measurement axis normal to the cornea (commonly called the VK axis) while the subject fixated the center of the topographer's mires. To satisfy the data-entry requirements of Zemax, the raw elevation height data were exported and fitted by the method of least squares to an aspheric, rotationally symmetric conic surface with conic constant p (p = 1 for a sphere, with other positive values being ellipses and negative values being hyperbolas; Rabbetts, 2007) and curvature c. To this conic surface we added irregularities defined by Zernike polynomials up to the sixth order (Navarro, Rozema, & Tassignon, 2013; Schwiegerling, Greivenkamp, & Miller, 1995) that provide a least-squares fit to the corneal height data. 

Central corneal thickness, ACD, and LT were measured at the relaxed accommodative state for 11 subjects with a partial coherence interferometer (Lenstar LS 900, Haag-Streit AG, Koeniz, Switzerland) with high resolution (about 0.01 mm; Suheimat, Verkicharla, Mallen, Rozema, & Atchison, 2015). Lenstar data were not available for the other eight subjects, so for those individuals we assumed the average values from our normative database. Axial length of the whole eye for all subjects was measured with an IOL Master 500 (Carl Zeiss AG, Oberkochen, Germany), which is also a partial coherence interferometer. These biometric data were available only for foveal vision and therefore we constrained axial length of the model only along the primary LoS. 

Ocular-aberration measurement

A custom-built instrument (Indiana Scanning Aberrometer for Wavefronts [ISAW]; Liu et al., 2016; Wei & Thibos, 2010a) was employed to measure ocular aberrations at 850 nm over a 30°-diameter field of view centered on the foveal LoS of the subject's left eye (the right eye was occluded). Zernike aberration coefficients depend on pupil diameter, which was measured by the aberrometer to monitor variations with accommodation. Measurements were obtained for a randomized sequence of 37 locations (eccentricities of 0°, 5°, 10°, and 13.5° along 12 visual meridians: 0° to 360° in 30° steps) over 20 s. This sequence was repeated for eight levels of monocular accommodative demand ranging from 1D beyond the far point to 6D in front of the far point, in 1D steps. Wave-front slope data were fitted with Zernike polynomials up to seventh order (36 terms) over the circumscribed circular domain of the fitted ellipse (Wei & Thibos, 2010b). The pupil coordinate system and nomenclature for reporting Zernike coefficients conformed to ANSI standard Z80.28 (American National Standards Institute, 2010), with the z-axis coinciding with the secondary LoS, which coincided with the measurement axis of the aberrometer. Zernike coefficients are reported for the measurement wavelength (not corrected for longitudinal chromatic aberration) and measured pupil size. For modeling and optimization purposes, aberration coefficients for second to fourth Zernike orders specified in the ophthalmic convention (ANSI standard Z80.28) were converted to the Zemax convention as described later (see Merit function) and rescaled to a fixed pupil diameter (the smallest pupil observed during ocular aberration measurements). 

Customization of a generic model with biometry data

The anterior corneal surface of the generic model was replaced with a customized surface referenced to the vertex norm located with Purkinje images. The corneal anterior surface (mean diameter = 10.32 mm, SD = 0.54 mm) was represented mathematically by an aspherical surface plus Zernike polynomials determined by least-squares fitting as already described (mean root mean square [RMS] error = 46.9 ± 11.4 μm, range = 31.7–75.7 μm), while the posterior corneal surface had a fixed radius of curvature and conic constant (r = −6.4 mm, P = 0.82; Dubbelman, Weeber, Van Der Heijde, & Völker-Dieben, 2002). To facilitate the OWT process, the z-axis of the customized cornea referenced to the VK axis was rotated and translated to become aligned with the primary LoS—that is, the reference axis of visual-field maps of ocular aberrations—as described in detail elsewhere (Liu & Thibos, 2017). The sign convention used for corneal rotation (Display Formula\(\def\upalpha{\unicode[Times]{x3B1}}\)\(\def\upbeta{\unicode[Times]{x3B2}}\)\(\def\upgamma{\unicode[Times]{x3B3}}\)\(\def\updelta{\unicode[Times]{x3B4}}\)\(\def\upvarepsilon{\unicode[Times]{x3B5}}\)\(\def\upzeta{\unicode[Times]{x3B6}}\)\(\def\upeta{\unicode[Times]{x3B7}}\)\(\def\uptheta{\unicode[Times]{x3B8}}\)\(\def\upiota{\unicode[Times]{x3B9}}\)\(\def\upkappa{\unicode[Times]{x3BA}}\)\(\def\uplambda{\unicode[Times]{x3BB}}\)\(\def\upmu{\unicode[Times]{x3BC}}\)\(\def\upnu{\unicode[Times]{x3BD}}\)\(\def\upxi{\unicode[Times]{x3BE}}\)\(\def\upomicron{\unicode[Times]{x3BF}}\)\(\def\uppi{\unicode[Times]{x3C0}}\)\(\def\uprho{\unicode[Times]{x3C1}}\)\(\def\upsigma{\unicode[Times]{x3C3}}\)\(\def\uptau{\unicode[Times]{x3C4}}\)\(\def\upupsilon{\unicode[Times]{x3C5}}\)\(\def\upphi{\unicode[Times]{x3C6}}\)\(\def\upchi{\unicode[Times]{x3C7}}\)\(\def\uppsy{\unicode[Times]{x3C8}}\)\(\def\upomega{\unicode[Times]{x3C9}}\)\(\def\bialpha{\boldsymbol{\alpha}}\)\(\def\bibeta{\boldsymbol{\beta}}\)\(\def\bigamma{\boldsymbol{\gamma}}\)\(\def\bidelta{\boldsymbol{\delta}}\)\(\def\bivarepsilon{\boldsymbol{\varepsilon}}\)\(\def\bizeta{\boldsymbol{\zeta}}\)\(\def\bieta{\boldsymbol{\eta}}\)\(\def\bitheta{\boldsymbol{\theta}}\)\(\def\biiota{\boldsymbol{\iota}}\)\(\def\bikappa{\boldsymbol{\kappa}}\)\(\def\bilambda{\boldsymbol{\lambda}}\)\(\def\bimu{\boldsymbol{\mu}}\)\(\def\binu{\boldsymbol{\nu}}\)\(\def\bixi{\boldsymbol{\xi}}\)\(\def\biomicron{\boldsymbol{\micron}}\)\(\def\bipi{\boldsymbol{\pi}}\)\(\def\birho{\boldsymbol{\rho}}\)\(\def\bisigma{\boldsymbol{\sigma}}\)\(\def\bitau{\boldsymbol{\tau}}\)\(\def\biupsilon{\boldsymbol{\upsilon}}\)\(\def\biphi{\boldsymbol{\phi}}\)\(\def\bichi{\boldsymbol{\chi}}\)\(\def\bipsy{\boldsymbol{\psy}}\)\(\def\biomega{\boldsymbol{\omega}}\)\(\def\bupalpha{\unicode[Times]{x1D6C2}}\)\(\def\bupbeta{\unicode[Times]{x1D6C3}}\)\(\def\bupgamma{\unicode[Times]{x1D6C4}}\)\(\def\bupdelta{\unicode[Times]{x1D6C5}}\)\(\def\bupepsilon{\unicode[Times]{x1D6C6}}\)\(\def\bupvarepsilon{\unicode[Times]{x1D6DC}}\)\(\def\bupzeta{\unicode[Times]{x1D6C7}}\)\(\def\bupeta{\unicode[Times]{x1D6C8}}\)\(\def\buptheta{\unicode[Times]{x1D6C9}}\)\(\def\bupiota{\unicode[Times]{x1D6CA}}\)\(\def\bupkappa{\unicode[Times]{x1D6CB}}\)\(\def\buplambda{\unicode[Times]{x1D6CC}}\)\(\def\bupmu{\unicode[Times]{x1D6CD}}\)\(\def\bupnu{\unicode[Times]{x1D6CE}}\)\(\def\bupxi{\unicode[Times]{x1D6CF}}\)\(\def\bupomicron{\unicode[Times]{x1D6D0}}\)\(\def\buppi{\unicode[Times]{x1D6D1}}\)\(\def\buprho{\unicode[Times]{x1D6D2}}\)\(\def\bupsigma{\unicode[Times]{x1D6D4}}\)\(\def\buptau{\unicode[Times]{x1D6D5}}\)\(\def\bupupsilon{\unicode[Times]{x1D6D6}}\)\(\def\bupphi{\unicode[Times]{x1D6D7}}\)\(\def\bupchi{\unicode[Times]{x1D6D8}}\)\(\def\buppsy{\unicode[Times]{x1D6D9}}\)\(\def\bupomega{\unicode[Times]{x1D6DA}}\)\(\def\bupvartheta{\unicode[Times]{x1D6DD}}\)\(\def\bGamma{\bf{\Gamma}}\)\(\def\bDelta{\bf{\Delta}}\)\(\def\bTheta{\bf{\Theta}}\)\(\def\bLambda{\bf{\Lambda}}\)\(\def\bXi{\bf{\Xi}}\)\(\def\bPi{\bf{\Pi}}\)\(\def\bSigma{\bf{\Sigma}}\)\(\def\bUpsilon{\bf{\Upsilon}}\)\(\def\bPhi{\bf{\Phi}}\)\(\def\bPsi{\bf{\Psi}}\)\(\def\bOmega{\bf{\Omega}}\)\(\def\iGamma{\unicode[Times]{x1D6E4}}\)\(\def\iDelta{\unicode[Times]{x1D6E5}}\)\(\def\iTheta{\unicode[Times]{x1D6E9}}\)\(\def\iLambda{\unicode[Times]{x1D6EC}}\)\(\def\iXi{\unicode[Times]{x1D6EF}}\)\(\def\iPi{\unicode[Times]{x1D6F1}}\)\(\def\iSigma{\unicode[Times]{x1D6F4}}\)\(\def\iUpsilon{\unicode[Times]{x1D6F6}}\)\(\def\iPhi{\unicode[Times]{x1D6F7}}\)\(\def\iPsi{\unicode[Times]{x1D6F9}}\)\(\def\iOmega{\unicode[Times]{x1D6FA}}\)\(\def\biGamma{\unicode[Times]{x1D71E}}\)\(\def\biDelta{\unicode[Times]{x1D71F}}\)\(\def\biTheta{\unicode[Times]{x1D723}}\)\(\def\biLambda{\unicode[Times]{x1D726}}\)\(\def\biXi{\unicode[Times]{x1D729}}\)\(\def\biPi{\unicode[Times]{x1D72B}}\)\(\def\biSigma{\unicode[Times]{x1D72E}}\)\(\def\biUpsilon{\unicode[Times]{x1D730}}\)\(\def\biPhi{\unicode[Times]{x1D731}}\)\(\def\biPsi{\unicode[Times]{x1D733}}\)\(\def\biOmega{\unicode[Times]{x1D734}}\)\(\alpha _x^{\prime} ,\alpha _y^{\prime} \)) and translation (Display Formula\({\rm{\Delta }}x,{\rm{\Delta }}y)\) is illustrated in Figure 1 (with the magnitude of translations and rotations exaggerated for clarity). Positive Display Formula\({\rm{\Delta }}x\) and Display Formula\(\alpha _x^{\prime} \) indicate the nasal visual field, while positive Display Formula\({\rm{\Delta }}y\) and Display Formula\(\alpha _y^{\prime} \) indicate the inferior visual field. This alignment was necessary to ensure that Zemax computed the optical thickness of the cornea in the direction that light rays travel through the eye's pupil. To allow for internal axial astigmatism (Liu & Thibos, 2017), the anterior surface of the crystalline lens was modeled by a toroidal surface with two independent conic constants associated with principal meridians. The starting point for optimization used a generic posterior lens surface with radius of curvature r_p = −6 mm and conic constant P = 1. The retinal surface was assumed to be spherical with an initial radius of curvature r_retina = −12 mm. 

Figure 1

 

Typical example of a customized eye model derived by ocular wave-front tomography. The line of sight in object space joins the foveal fixation point with the center of the entrance pupil and coincides with the z-axis. (A) Vertical (y-z) cross section. (B) Horizontal (x-z) cross section. Both cornea and crystalline lens are translated and tilted relative to the line of sight. The sign convention used for corneal and lens rotation (\(\alpha_{\rm{x}}^{\prime} ,\alpha_{\rm{y}}^{\prime}\)) and translation (\({\rm{\Delta }}x,{\rm{\Delta }}y)\) is: Positive \({\rm{\Delta }}x\) and \(\alpha _{\rm{x}}^{\prime} \) indicate the nasal visual field, while positive \({\rm{\Delta }}y\) and \(\alpha _{\rm{y}}^{\prime} \) indicate the inferior visual field.

Biometry was used to set the central thickness of the cornea for each eye, which was assumed to be invariant with accommodation. Biometry was used to set the starting values for optimizing ACD and lens parameters for each state of accommodation. The depth of the vitreous chamber is determined by axial length and position of the lens, so is not a free parameter for optimization. The model iris was juxtaposed to the anterior surface of the lens and perpendicular to the foveal LoS. The refractive indices of ocular media for a working wavelength of 850 nm were determined from the Herzberger formula (Navarro et al., 1985). 

A summary of all model parameters needed to implement a customized Zemax model is provided in Table 1. The numerical values in the table specify fixed-focus eyes of two subjects with presbyopia. Of the 30 parameters listed, one is derived from the others: Vitreous-chamber depth is calculated by subtracting central corneal thickness (a biometric measurement), ACD (an optimization variable), and lens thickness (an optimization variable) from the axial length of the whole eye (biometric measurement). Of the remaining 29 parameters in Table 1, only those 14 parameters pertaining to accommodation were subject to optimization. Statistical analysis of these accommodation-dependent model parameters for individual eyes yielded regression parameters that were averaged across individuals to reveal population trends regarding how accommodation affects the structure of the eye's optical system. 

Table 1

 

Customized eye model for eyes with presbyopia. Notes: The parameters in normal font indicate subject LNT, while the italic parameters indicate subject AB. The anterior lens is toroid with vertical and horizontal curvatures and conic constants. The thickness listed in the table (from top to bottom) refers to central corneal thickness, thickness of the anterior chamber, lens thickness, vitreous-chamber depth, and the axial length of whole eye. The iris (stop aperture) was juxtaposed to the anterior surface of the lens and perpendicular to the foveal line of sight. Superscripts indicate the origin of the parameters. Vitreous-chamber depth is a dependent variable computed from the other parameters. n is the refractive index at wavelength λ = 850 nm. \(\Delta x\) and \(\Delta y\) are translation of the cornea/lens relative to the line of sight in object space. \(\alpha _{\rm{x}}^{\rm{^{\prime} }}\) and \(\alpha _{\rm{y}}^{\rm{^{\prime} }}\) are rotation of the cornea/lens relative to the line of sight in object space. Positive \(\Delta x\) and \(\alpha _{\rm{x}}^{\prime} \) indicate the nasal visual field, while positive \(\Delta y\) and \(\alpha _{\rm{y}}^{\prime} \) indicate the inferior visual field. Positive Φ is rotation of the lens about its optical axis counterclockwise.

 

Implementation of OWT

Ocular wave-front tomography is the process of estimating the shapes and positions of the main refractive elements of an eye model based on wave-front aberrations measured along multiple LoSs across the visual field. Building a customized model typically requires a trade-off between anatomical similarity and functional equivalence. This trade-off is established by creating a suitable merit function based on relevant variables that is minimized by the optimization algorithm (Wei & Thibos, 2008). In the following, we summarize our implementation of the merit function to favor anatomical similarity as a four-step process. 

Configuring the template model with initial conditions

Since the search space for adjusting ocular parameters is large and complex in our case, a successful outcome requires suitable initial conditions that make it unlikely that the optimization engine will become trapped in a spurious local minimum. For this reason we initialized the template model with individual biometry data when available, or the mean of our population when individual measurements were unavailable. 

Merit function

The merit function embodies both goals of OWT: functional equivalence and anatomical accuracy. A merit function in Zemax consists of a list of operands, the target value of each operand, the relative weighting assigned to each operand, and a formula that uses the operands to compute a scalar metric to be optimized. In our case, the operands are the anatomical parameters and aberration coefficient of the model eye. The aberrometry measurements of a subject's eye are the target data values to be accounted for by the model. 

Our merit function had two parts. The first part, designed to ensure anatomical similarity, was a series of boundary conditions that constrained the possible values of anatomical parameters to be searched. For example, we forced the ACD to fall within the normal range of human eyes (2.91–3.65 mm; Rozema et al., 2011) by imposing a heavy penalty (operand weight = 10⁶) when the optimization engine attempted to go outside the allowed range. Lens parameters, including radii of curvature (Rosales et al., 2006), conic constants (Dubbelman et al., 2005), and decentration (Tabernero, Benito, Nourrit, & Artal, 2006), were not measured in our subjects yet are expected to change significantly with accommodation. Therefore, for lens parameters we imposed a soft constraint (operand weight = 1) proportional to the Euclidean distance between model parameters and population averages. In this way the optimizer tolerated modest departures from anatomical expectations if the reward was a significant improvement in the model's fit to the aberration measurements. 

The second part of the merit function, designed to ensure functional equivalence, assessed the agreement between measured wave-front aberrations and aberrations of the customized model. Zemax reports aberrations as Zernike coefficients over the circular physical pupil, but empirical wave-front aberrations measured off-axis exist over an elliptically foreshortened entrance pupil. Therefore, we put these two sets of Zernike coefficients on a common basis by converting the Zernike coefficients reported by ISAW to a new set of Zernike coefficients appropriate for Zemax. This conversion was implemented by a transformation matrix that depends on eccentricity and meridian (Wei & Thibos, 2008, 2010b). The vector of empirical Zernike coefficients from second to fourth order was then subtracted from the corresponding theoretical values of the model and the length of the difference vector (operand weight = 3 for defocus and astigmatism, 1 for other Zernike terms) was taken as an RMS measure of disagreement. This procedure was repeated for all 37 measured visual fields under the same accommodation state, with one exception. Measurements at 13.5° eccentricity along the horizontal meridian on the temporal side were often corrupted because of the optic disk and therefore were omitted for all subjects. 

In principle, the model parameters that optimally account for the aberrations of the relaxed eye could be taken as the starting point for finding the optimum model for the 1D accommodated eye, which in turn could be the starting point for the 2D accommodation state, and so on. However, the outcome of that sequential approach would depend strongly on the initial optimization of the relaxed eye and might lead to a propagation of errors. To protect against those undesirable possibilities, the same initial template was used to optimize model parameters for each level of accommodation tested. For similar reasons, ocular wave-front aberrations for all 36 field locations (12 Zernike terms × 36 locations = 432 operands) were used simultaneously to obtain the relatively small number of parameters needed to model the eye's aberration structure for a given level of accommodation. In addition to 14 (variables) × 2 (upper and lower boundaries) = 28 boundary constraints, a combined vector of 460 operands was included in the merit function. 

Anatomical variables

Corneal topography was measured with accommodation relaxed, and we assumed it applied also to other states of accommodation (J. C. He, Gwiazda, Thorn, Held, & Huang, 2003). As the profile of the posterior cornea for individual subjects was not available, and likely makes only a small contribution (Atchison, Suheimat, Mathur, Lister, & Rozema, 2016) to overall aberrations due to a small refractive-index difference between cornea and aqueous, we assumed the same surface profile for all subjects and all states of accommodation. We verified that the corneal thickness of the model increases monotonically from center to periphery for all subjects. ACD and LT were set as variables, and vitreous thickness was calculated automatically assuming that axial length is independent of accommodation. Radii of curvature and conic constants for both lenticular surfaces were set as variables to take account of the dramatic shape changes during accommodation. The toroidal lens was also allowed to rotate about its optical axis to duplicate the orientation of reported internal axial astigmatism (Liu & Thibos, 2017). It is well known that the crystalline lens is decentered and tilted about the LoS (Liu & Thibos, 2017; Rosales & Marcos, 2007; Tabernero, Benito, et al., 2006), so we treated lens, vitreous chamber, and retina as a single entity that could be decentered and tilted horizontally and vertically. The sign convention used for lens rotation and translation is illustrated in Figure 1 and is the same as for corneal decentration and tilt. We also permitted retinal curvature to vary, to allow for possible changes during accommodation (Walker & Mutti, 2002). The total of 14 anatomical variables is summarized in Table 1, with superscripts indicating their origin. 

Optimization

Like any optical design problem, the task of optimization is to improve the system performance of a reasonable initial model rather than find the prototype from scratch. Usually, more variables in the optimization process lead to better results. However, the price is increased computation time and proneness to becoming trapped in a local minimum. In addition, overparameterizing could lead to modeling the measurement errors rather than the properties of the tested eye. It is the optical designer's job to provide an initial model to guide the optimization engine and reduce the search space. For this reason, we developed two optimization stages. In the first stage, we obtained a moderate initial model within 10 min by solely incorporating defocus and astigmatism along the foveal LoS in the merit function and permitting lenticular curvatures and rotations about the z-axis as variables. The resultant model accounts for a significant portion of lower order aberrations, which dominate retinal image quality, and provided the initial conditions for Stage 2 of optimization. In the second stage, the complete merit function and full list of variables were included. The optimization engine, which used a damped least-squares algorithm (Meiron, 1965), adjusted all 14 variables simultaneously until the value of the merit function could not decrease further. This two-stage optimization procedure was separately applied for each accommodation state and repeated for each subject. The average optimization time was about 18.5 hr for one customized model eye on a 3.40-GHz Intel Core i7 processor with eight cores. Other optimization algorithms were tested, but no remarkable improvement of the fit over the damped least-squares algorithm was found. 

Pan-retinal metric of image quality

Most of the results of our study are presented in the form of maps showing how aberration coefficients of the eye vary across the visual field. To quantify how these aberrations influence the ability of cone photoreceptors to signal the retinal image, we computed cone activation for a point source at some visual-field eccentricity and meridian (Display Formula\(\varepsilon ,\theta \)) defined as the peak value of the convolution of three functions: the point source, the eye's point-spread function along the corresponding LoS, and the cone's spatial weighting function, assumed to be uniform across the cone aperture. Convolution was carried out in the spatial-frequency domain according to  

\begin{equation}\tag{1}{\rm{Cone}}\;{\rm{activation}}\left( {\varepsilon ,\theta } \right) = \int {\int {{\rm{Object}}\;{\rm{Contrast}}\left( {{f_x},{f_y}} \right)} \cdot {\rm{Optical}}\;{\rm{MTF}}\left( {{f_{\rm{x}}},{f_{\rm{y}}},\varepsilon ,\theta } \right)} \cdot {\rm{CTF}}\left( {{f_{\rm{x}}},{f_{\rm{y}}},\varepsilon ,\theta } \right)d{f_{\rm{x}}}d{f_{\rm{y}}}{\rm .}\end{equation}

 

In this equation the optical modulation transfer function (MTF) is derived from the local Zernike aberrations, and the cone transfer function (CTF) is the Fourier transform of a uniform disk the size of the cone aperture, which varies with retinal eccentricity (Curcio, Sloan, Kalina, & Hendrickson, 1990). We projected histology data (cone diameter and density) into the visual field for each customized eye model using a nonlinear projection method described previously (Drasdo & Fowler, 1974). 

Evaluating Equation 1 provides a map of how the quality of the retinal image of a point source or other localized object, as seen by cone photoreceptors, varies across the visual field. To reduce this map to a single number representing the total cone activation produced by an array of such stimuli across the visual field, we assumed the contribution of cones in any local area of retina is proportional to the number of cones in that area. Accordingly, we weighted cone activation from Equation 1 by cone density (which varies with visual-field location; Curcio et al., 1990) and then integrated across the visual field according to Equation 2 to produce an integrated measure of cone image quality (CIQ):  

\begin{equation}\tag{2}{\rm{CIQ}} = {{\int\limits_0^{2\pi } {\int\limits_{{0^o}}^{{{15}^o}} {{\rm{Cone}}\;{\rm{activation}}\left( {\varepsilon ,\theta } \right)} \cdot {\rm{Normalized}}\;{\rm{cone}}\;{\rm{density}}\left( {\varepsilon ,\theta } \right)\varepsilon d\varepsilon \;d\theta } } \over {{\rm{foveal}}\;{\rm{cone}}\;{\rm{activation}}\;{\rm{with}}\;{\rm{perfect}}\;{\rm{optics}}\left( {{\rm{point}}\;{\rm{source}}} \right)}}{\rm .}\end{equation}

 

For this purpose it was convenient to use a unitless, normalized cone-density map:  

\begin{equation}\tag{3}{\rm{Normalized}}\;{\rm{cone}}\;{\rm{density}}\left( {\varepsilon ,\theta } \right) = {{{\rm{Cone}}\;{\rm{density}}\left( {\varepsilon ,\theta } \right)} \over {\int\limits_0^{2\pi } {\int\limits_{{0^o}}^{{{15}^o}} {{\rm{Cone}}\;{\rm{density}}\left( {\varepsilon ,\theta } \right)d\varepsilon \;d\theta } } }}{\rm .}\end{equation}

 

To aid interpretation, CIQ was made unitless by normalizing (in Equation 2) the spatial integral of cone activation by the foveal cone activation expected for an optically perfect eye with fixed pupil diameter (3 mm) for a point object. 

Results

This section is organized into three parts. First, the customized models produced by OWT are presented for two subjects with presbyopia, whose ocular aberrations are not affected by accommodative fluctuations. Second, we describe in detail the modeling results for a typical subject (RH) to demonstrate how an individual eye model changes structure to account for ocular aberrations (including defocus, which reflects changes in refractive state) across the visual field as the eye accommodates. In the third section we present a population-based model of the average adult eye with emmetropia to report trends and individual variability for model parameters. Compared to the Navarro model used as an initial template, all the customized models demonstrated significant reduction of RMS error between measured ocular wave-front errors and model predictions. 

Control subjects with presbyopia

Although OWT has been successfully demonstrated using a physical model eye and for numerically simulated measurements of ocular aberrations (Wei & Thibos, 2008), it was important to verify the feasibility of our current implementation using empirical measurements on fixed-focus eyes. We conducted this test with Zernike coefficients measured for two subjects with presbyopia with zero accommodation. Visual-field maps of Zernike defocus (relative to infinity), astigmatism, coma, spherical aberration, and RMS of higher order aberration (HOA RMS, from third to seventh order) are shown for one subject (LNT) in Figure 2. The field maps are arranged in columns, with empirical data in the upper row, ray-tracing prediction in the middle row, and their difference in the bottom row (signed difference for defocus and spherical aberration). To facilitate the comparison, all Zernike coefficients were scaled to a 7-mm pupil diameter. For display purposes, the units of second-order aberrations (defocus, astigmatism) are diopters and the units of higher order aberrations (coma, spherical aberration) are microns. Astigmatism and coma are vector quantities displayed in polar format (magnitude indicated by color code, axis indicated by short line passing through each symbol) but were converted to Cartesian format (J₀, J₄₅) and (C₃⁻¹, C₃⁺¹) when computing the vector-difference maps. 

Figure 2

 

Visual-field maps for the major Zernike aberrations for subject LNT, with presbyopia. The upper row shows empirical data, the middle row shows aberrations of the model, and the bottom row shows their difference (signed difference for defocus and spherical aberration). Color indicates the magnitude of each aberration interpolated from sample points (symbols). Lines drawn through symbols indicate the axis of astigmatism or coma. The mean ± 1 standard deviation of the 36 measurements are shown in the lower left corner. All maps of the visual field span 30° in diameter (i.e., maximum eccentricity = 15°).

Visual comparison of maps in the upper and middle rows confirms the model's ability to capture the most salient features of the aberration structure of LNT's eye, including the superior-inferior gradient of lower order refractive state and the bow-tie signature of ocular astigmatism due to axial astigmatism with little effect of oblique (off-axis) astigmatism in the central visual field (Liu & Thibos, 2016). The model also captures the main features of the coma maps, showing rotational symmetry and radial axis relative to a center of symmetry displaced about 7° into the temporal visual field. Maps for spherical aberration and total HOA RMS also indicate weaker aberrations in the temporal field for the tested eye and customized model. The similarity between maps of coma and HOA RMS in Figure 2 confirms previous observations (Howland & Howland, 1977; Salmon & van de Pol, 2006) that coma is the dominant high-order aberration, which is true also for the customized model. Although the salient features of subject AB's field maps (not show) were different from those in Figure 2, the features of AB's eye were also captured by his customized model. 

The qualitative impressions described are confirmed quantitatively by the small differences across the visual-field shown in the bottom row of maps. The close agreement between mean magnitude (averaged across the visual field) of empirical coma (0.5 μm) and HOA RMS (0.7 μm) occurs also for the customized model (coma RMS difference = 0.08 μm, HOA RMS difference = 0.30 μm). The largest disparity between empirical and model behavior was for spherical aberration, but even in that worst case the average difference between experiment and model prediction across the visual field was only 0.036 μm for subject LNT and 0.031 μm for subject AB. 

The structural parameters of the customized models derived by OWT (Table 1) and used to produce the middle row of maps in Figure 2 fall in the normal population range (Rozema et al., 2011), which confirms that the anatomical constraints programmed into the Zemax merit function had the desired effect of producing anatomically reasonable models. The optimized ACD values (2.91 mm for LNT, 2.30 mm for AB) and lens-thickness values (4.76 mm for LNT, 4.76 mm for AB) of the model were within 10% of the measured values (ACD: 3.29 mm for LNT, 2.45 mm for AB; LT: 4.51 mm for LNT, 4.65 mm for AB). 

Based on this initial experience with subjects with presbyopia, we concluded that OWT is an effective procedure for building customized eye models that are anatomically similar and functionally equivalent to real eyes for the purpose of describing optical aberrations and retinal image quality over the central visual field. The next step, described later, was to repeat the process for multiple states of accommodation by young adults. 

One potential application of customized eye models is to interpolate or extrapolate aberration coefficients to unmeasured visual-field locations or accommodative states. This raises the question: Can a customized optical model constructed from aberrometry measurement obtained in the central visual field be used to accurately predict aberrations in the peripheral field? We investigated this question by instructing a subject with presbyopia (AB) to fixate an eccentric point in the visual field, which allowed the ISAW instrument to measure an additional 37 new visual-field locations. About half of these new test locations were at greater than 15° eccentricity, to test extrapolation, and the other half were at eccentricities less than 15°, to test interpolation. (Equipment limitations required a fixed-target vergence for the eccentric fixation points, which precluded using this experimental design to investigate different accommodative states in younger subjects.) The eccentric fixation target was positioned sequentially at each of four possible locations (15° superior, inferior, nasal, and temporal), which expanded the array of tested locations to a total of 185 points covering a 50° field of view, as shown in Figure 3. Only the aberrometry measurements from the central 37 points were used to construct a customized optical model, which was then evaluated at the additional 148 test locations for comparison with empirical measurements. 

Figure 3

 

Extended visual-field maps for the major Zernike aberrations for subject AB with presbyopia. All maps of the visual field span 50° in diameter (i.e., maximum eccentricity = 25°). The central 30°-diameter field used for obtaining the customized model is indicated by the large black circle concentric to the field map. Data collected at the blind spot were excluded from the difference map for defocus and spherical aberration. Other graphic conventions of individual maps are the same as in Figure 2.

Experimental measurements of the eye's major aberrations (defocus, astigmatism, coma, and spherical aberration) are compared in Figure 3 with aberrations of the optical model at the same test locations. Each column of the figure shows the value of a single type of aberration over an extended field of view (±25° eccentricity in all meridians). The top row of field maps shows the empirical data, the middle row shows the aberrations computed from the model, and the bottom row shows the algebraic difference between the empirical and model values. Visual comparison of the top and middle rows confirms that the extrapolated and interpolated values obtained from the model agree closely with the extended set of measurements obtained from this subject's eye at all locations except the blind spot, which is not present in the model. Note that the scale of the color bars for the bottom row of difference maps is 1/10 that of the top and middle rows, which nevertheless is adequate for displaying the differences. Although the errors of prediction by the model are greater for extrapolation than for interpolation, the largest errors are less than about 10% of the predicted values. We conclude from this experiment that an optical model derived from biometry and aberration measurements over the central 30°-diameter field is sufficient for predicting with reasonable accuracy the aberrations over an extended 50° field of view. 

Accommodating model for a typical individual eye

Visual-field maps obtained for a typical young eye with emmetropia (RH) are shown for multiple accommodative states in Figures 4 (defocus, i.e., refractive state), 5 (astigmatism), 6 (coma), 7 (spherical aberration), and 8 (HOA RMS). Using the same conventions as in Figure 2, the empirical maps are shown in the upper row, maps for the model are shown in the middle row, and the magnitudes of corresponding Zernike vector differences (RMS error) are shown in the bottom row (except for defocus and spherical aberration, which are signed differences). Each column of maps is for a fixed accommodation demand (labeled at the upper right corners of the upper row), with the maximum accommodative demand in the leftmost column and minimum demand in the rightmost column. Although pupils tend to constrict slightly when the eye accommodates under monocular viewing conditions, we computed Zernike coefficients for a fixed (5-mm) pupil diameter for modeling purposes and to facilitate comparison across accommodative states. 

Figure 4

 

Defocus outcome for customized models of RH's eye at different accommodation states. Empirical data are in the upper row, aberrations of the customized models are in the middle row, and their signed difference is in the bottom row. Color indicates defocus interpolated from sample points shown by symbols. Accommodation demand is shown in the upper right corner of the top row. Graphic conventions of individual maps are the same as in Figure 2.

Figure 5

 

Astigmatism outcome for customized models of RH's eye at different accommodation states. Graphic conventions are the same as in Figure 3, except magnitude is displayed instead of signed value. Lines drawn through symbols indicate the axis of astigmatism.

Figure 6

 

Coma outcome for customized models of RH's eye at different accommodation states. Graphic conventions are the same as in Figure 3, except magnitude is displayed instead of signed value. Lines drawn through symbols indicate the axis of coma.

Figure 7

 

Spherical-aberration outcome for customized models of RH's eye at different accommodation states. Signed value is displayed. Graphic conventions are the same as in Figure 3.

Figure 8

 

Higher order aberration outcome for customized models of RH's eye at different accommodation states. Graphic conventions are the same as in Figure 3, except root mean square magnitude is displayed instead of signed value.

The variation of refractive state as function of accommodation is shown for subject RH in Figure 4. Empirical measurements (upper row) and the model (middle row) reveal nearly identical changes in refractive state as the eye accommodates. These changes are remarkably uniform across the visual field, which is a trait of Zernike defocus described previously for the eyes in our study population with myopia and with emmetropia (Liu et al., 2016). The averaged signed error shown in the bottom row of field maps is about 0.002D, independent of accommodative state, which confirms that the customized model accurately reproduces the eye's variation in refractive state over a large range (−1D to +6D) of accommodation demand. This conclusion is further supported by Figure 9A, which shows a very close match between the mean refractive states (averaged across the visual field) of measurement (black squares) and model (red crosses). This is a statistically powerful conclusion since the empirical data are very precise (error bars—i.e., the standard error of the mean of 37 field locations—are smaller than the symbols). 

Figure 9

 

Comparison between measured and modeled aberrations of RH's eye. Symbols show mean values averaged across the visual field, with error bars showing the standard error of the mean. (A) Effect of accommodation demand on refractive state (i.e., Zernike defocus). Black squares show empirical data, red crosses show values for a customized model derived with normal boundary constraints, and blue triangles show values for a customized model derived with an expanded range of anterior lens curvatures specified in the optimization merit function. (B) Effect of accommodation demand on spherical aberration for the same three configurations as in (A). (C) Effect of accommodation demand on residual root mean square (RMS) error of astigmatism. (D) Effect of accommodation demand on residual RMS error of coma. (E) Effect of accommodation demand on higher order aberration RMS for the same three configurations as in (A). (F) Comparison of total RMS residual error (red bars) and higher order aberration RMS residual error (black bars) at different accommodation states.

Figure 9

 

Continued

Unlike defocus (Figure 4), the visual-field map of ocular astigmatism (upper row of Figure 5) for this individual eye changes little with accommodation, as reported previously (Liu & Thibos, 2017). Nevertheless, the customized models captured not only the subtle and nonsystematic changes of astigmatism with accommodation but also the slightly displaced bow-tie pattern across the visual field produced when axial astigmatism interacts with oblique astigmatism (Liu & Thibos, 2016, 2017). The discrepancy between empirical data and the model (bottom row of Figure 5) is small and randomly distributed across the visual field. The mean RMS error (averaged across the visual field), which is the magnitude of the residual astigmatism not captured by the model, is only 0.05D and randomly changes with accommodation, as shown in Figure 9C (red crosses). 

Like oblique astigmatism, coma is a vector quantity with a magnitude that is rotationally symmetric and an axis that is radially orientated with respect to the eye's optical axis. Misalignment between the eye's optical axis and the LoS, which is the ISAW reference axis and the origin of all field maps, will displace the center of symmetry of the coma map, just as for the astigmatism map (Liu & Thibos, 2016). This displaced center of symmetry is evident in the empirical data for subject RH (top row of Figure 6) and in the coma maps for the model (middle row of Figure 6). These maps change little with accommodation, yet even these subtle changes are captured and explained by the model. Notice, for example, the wobbling of the symmetry center, which the model accounts for by tilt and decentration of the lens as the eye accommodates. Residual errors of the model were greatest in the superior temporal visual field for this eye (bottom row of Figure 6), with a mean RMS error (averaged across the visual field) of about 0.11 μm that fluctuates randomly during accommodation, as shown in Figure 9D (red crosses). 

The model's ability to capture the effect of accommodation on spherical aberration across the visual field is shown in Figure 7. Previous work has established that spherical aberration changes in the negative direction as the eye accommodates. This is true not only for foveal vision (Cheng et al., 2004; López-Gil & Fernández-Sánchez, 2010) but also for peripheral vision (Lundström et al., 2009; Mathur et al., 2009). Subject RH followed these population trends, as shown in the top row of Figure 7. Similar to refractive state (Figure 4), the visual-field map for spherical aberration changes uniformly across the visual field (uniform color in each map) and monotonically changes in the negative direction with accommodative demand (the descending sequence of colors). For this subject, the mean spherical aberration across the visual field changed by 0.16 μm, from +0.09 μm at the most relaxed states to −0.07 μm at the most accommodated state. Although the customized models also produced uniform field maps of spherical aberration that changed uniformity in the negative direction (middle row of Figure 7), they consistently overestimated the spherical aberration by about 0.1 μm, as displayed in the bottom row. As summarized in Figure 9B, the mean spherical aberration across the visual field for experimental data (black squares) and the model (red crosses) are almost parallel with each other, but separated by about 0.1 μm. Error bars representing the standard error of the mean are smaller than the symbols, indicating spatial uniformity of the visual-field maps. This systematic bias in representing spherical aberration also occurred for other subjects, a point addressed in the Discussion. 

Empirical maps of HOA RMS shown in the upper row of Figure 8 are similar to maps of coma (Figure 6), suggesting that coma is the dominant HOA at all states of accommodation. This dominance masks the systematic changes in spherical aberration (Figure 7) with accommodation, and consequently the HOA RMS maps vary little as the eye accommodates. Moreover, HOA RMS ignores the sign of spherical aberration, and consequently +0.04-μm spherical aberration at 2D accommodation demand makes the same contribution to HOA RMS as −0.04-μm spherical aberration at 5D accommodation demand. This effect is explicitly shown in Figure 9E. Since spherical aberration for this subject reverses sign, the mean HOA RMS across the visual field forms a V-shaped function (black squares). This behavior is not shared by the model, however, which has positive spherical aberration for the whole accommodation range, as shown in Figure 9B. Consequently, the predicted HOA RMS is slightly lower with accommodation (red crosses in Figure 9E). Nevertheless, the mean RMS error of HOA across the visual field was only 0.2 μm for all accommodation states, in spite of small variation within the central 30° visual field, as shown in the bottom row of Figure 8. This result demonstrates that the OWT method can produce customized models that replicate higher order wave-front aberrations of the individual eye across the central visual field with RMS error less than wavelength/4, which confirms the model's ability to represent the aberration structure of subject RH's eye. Compared to total RMS error (black bar in Figure 9F, from second- to seventh-order Zernike coefficients), the HOA (red bars) contributes approximately 87% of residual RMS error, indicating nearly perfect recovery of lower order aberrations (defocus and astigmatism) achieved by the customized model. Furthermore, the RMS error for both HOA and total Zernike vector are not influenced by accommodation, which indicates the robustness of the OWT method in achieving a level of accuracy limited only by diffraction. 

The preceding analysis confirms the first goal of OWT, which is to provide an economical representation of a large body of aberration data obtained over the visual field for various states of accommodation. The second purpose is to provide a mechanistic explanation for the origin of those optical data by showing the structural changes in the eye's optical system produced by accommodation. Those structural changes are documented in Figure 10, which shows how the parameters of the eye model for subject RH vary with accommodation. Anterior and posterior lens surfaces become more curved (Figure 10A, black and red crosses, respectively) and the lens thickens and ACD is reduced (Figure 10B, black and red crosses). This decrease in ACD is less than the increase of lens thickness, which indicates that the posterior surface of the lens drifts away from cornea during accommodation (Dubbelman et al., 2005). Notice, however, that the curvatures of both lenticular surfaces escaped from the soft anatomical constraints (shaded areas) imposed in the merit function that drives optimization. We will return to this important point in the Discussion. 

Figure 10

 

Structural parameters of customized models of RH's eye as a function of accommodation. Crosses represent the normal model; triangles represent the revised model with loosened anterior lens curvature. (A) Black symbols indicate the anterior surface of the lens, and red symbols indicate the posterior surface of the lens. The gray and red shaded areas are soft anatomical constraints for the anterior and posterior surfaces, respectively. Circles and diamonds are published anatomical measurements of human eyes. (B) Black symbols indicate anterior-chamber depth, and red symbols indicate lens thickness. The gray and red shaded areas are soft anatomical constraints for anterior-chamber depth and lens thickness, respectively. Translation (C) and tilt (D) of the lenticular optical axis relative to the line of sight are shown as a trajectory parameterized by accommodative demand. A positive x value indicates the nasal visual field, while a positive y value indicates the inferior visual field. Circles indicate 6D accommodation demand, while squares indicate −1D accommodation demand. Red symbols show the normal model, and blue symbols show the revised model.

The model also reveals a vertical decentration of lens position during accommodation by about 0.2 mm (Figure 10C, black crosses with red circle and square indicating the most accommodated and most relaxed states), perhaps due to the influence of gravity as the tension in the suspensory zonules is reduced (Atchison, Claydon, & Irwin, 1994; L. He & Applegate, 2011; Radhakrishnan & Charman, 2007). The lens also appeared to rotate slightly about its own axis in this eye (Figure 10D, black crosses), but the effect was not systematically related to accommodation, so it might be an artifact of the modeling procedure. This interpretation is consistent with our previous analysis of astigmatism based on the same data set, which indicated that the optical axis (along which oblique astigmatism of the whole eye is zero) relative to the LoS is not strongly influenced by accommodation (Liu & Thibos, 2017). 

Population trends

Customized models like those described for subject RH were constructed and inspected for the left eyes of 19 normal, healthy subjects (nine with emmetropia and 10 with myopia) to determine which of the many features observed in RH's eye represent general trends in these populations. Goodness of fit of the individual models was measured by the RMS error of Zernike vectors averaged across eight accommodation levels. The mean RMS errors of both populations are not significantly different for the full Zernike vector (0.23 ± 0.09 μm for eyes with emmetropia and 0.21 ± 0.11 μm for eyes with myopia) or the HOA Zernike vector (0.17 ± 0.06 μm for eyes with emmetropia and 0.17 ± 0.09 μm for eyes with myopia). HOA contributes about 80% of the RMS error. 

From these results we constructed a population-based model of a hypothetical “average eye” whose parameters change with accommodation in the same way as the averaged value of the population of customized models. Since the vitreous-chamber thickness has large variability within the population with myopia (17.58 ± 0.46 mm), we chose to base our average model on the population with emmetropia using parameters specified in Table 2. 

Table 2

 

Accommodation-dependent schematic eye parameters of hypothetical mean adult with emmetropia. Notes: Acc = accommodation demand (in diopters). All the surfaces are rotationally symmetric. The conventions are the same as Table 1.

 

Variation of selected anatomical parameters of the average-eye model with accommodation demand is displayed in Figure 11. Population means of model parameters were well fit by linear functions, excluding −1D demand, which typically produced inappropriate accommodation to targets beyond the eye's far point (Hennessy, 1975). Consistent with previous reports (Dubbelman et al., 2005; Rosales et al., 2006), the curvatures of both lens surfaces increase linearly as the average eye accommodates (adjusted R² = 0.998 for the anterior surface and 0.938 for the posterior surface), but more so for the anterior surface (as shown in Figure 11A). As shown in Figure 11B, the ACD declines linearly with accommodation (adjusted R² = 0.942), while LT increases (adjusted R² = 0.764) with a nearly identical rate, indicating that the lenticular posterior pole of the average-eye model based on our study population maintains its axial position regardless of accommodation state. The error estimates of these linear regression parameters are listed in Table 3. The lens is slightly decentered in a downward temporal direction (relative to the foveal LoS) with moderate variation during accommodation (Display Formula\({\rm{\Delta }}x\) = −0.098 ± 0.02 mm, Display Formula\({\rm{\Delta }}y\) = 0.12 ± 0.03 mm), as displayed in Figure 11C. The variation of torsional rotation of the lens is negligible on average, and the averaged angle between the optical axis of the lens and the LoS is 4.3° temporally and 1.1° inferiorly, as shown in Figure 11D. (See Figure 1 for our sign conventions for translations and rotations of the lens.) 

Figure 11

 

Structural parameters of a hypothetical, average emmetropic-eye model as a function of accommodation demand. Graphic conventions are the same as in Figure 9. Dashed lines are linear regression lines labeled with fitted equations.

Table 3

 

Error estimates of the linear regression parameters (y = p1 × Acc + p2). Notes: Acc is accommodation demand (in diopters).

 

Visual-field maps of optical aberrations of the average model eye are shown in Figure 12, with each column a different accommodation state and each row a different Zernike aberration. As we reported previously (Liu et al., 2016; Liu & Thibos, 2017) and have demonstrated in this article, defocus and spherical aberration are relatively uniform across the central visual field but vary systematically with accommodation. To the contrary, astigmatism and coma are largely invariant with accommodation, which permits the use of a higher resolution color code in Figure 12 to reveal patterns of nonuniformity caused by misalignment of the model's lens with respect to LoS. Null points for astigmatism and coma are displaced from the field-map center (LoS) but in opposite directions. Additional analysis showed that 4.3° of temporal lenticular tilt accounts for mean decentration of the coma null point, and 0.212 mm of nasal corneal decentration accounts for mean decentration of the astigmatism null point. 

Figure 12

 

Predicted field maps of major aberrations for the average-eye model at different accommodation states: (A) defocus, (B) astigmatism, (C) coma, (D) spherical aberration. Graphic conventions are the same as in Figure 3.

Discussion

This study demonstrates the feasibility of using OWT to construct a customized, accommodating eye model to account for ocular aberrations over the central 30° of visual field. These optical models successfully captured the salient changes of lower and higher order aberrations that occur during accommodation. For example, refractive state (i.e., Zernike defocus) and spherical aberration both change uniformly over the visual field as the eye accommodates, while astigmatism and coma are almost invariant with accommodation. Our simplified models employed a lens with homogeneous refractive index, yet the averaged RMS fitting errors of the models (typically about 0.2 μm) were almost the same as other customized models using a gradient-index (GRIN) lens (Goncharov et al., 2008; Polans et al., 2015). 

The structural parameters of the resulting eye models are realistic, and they changed in the expected way when accommodating. For example, during accommodation the curvature of both surfaces of the lens increases, the ACD shortens, and the lens thickens, all of which agree with physiological observations (Dubbelman et al., 2005; Rosales et al., 2006). By capturing these structural changes in a functional model we achieve a highly efficient representation of large amounts of optical data. In the present experiments, for example, nearly 300 wave-front aberration maps (37 locations × 8 levels of accommodative demand) measured for each subject were summarized by just 30 model parameters plus a few regression parameters describing how the eye's structure changes with accommodation. When those parameters are used to configure a Zemax model of the accommodating eye, we not only accurately recover the measured aberrations, but the model also serves as an interpolation tool to describe other visual-field locations and accommodative states not measured (Figure 4). 

Optical models derived by the OWT method are not unique (Wei & Thibos, 2008). For example, an accommodating model lens with uniform or gradient refractive index can yield valid solutions with similar RMS fitting error. The implication for practical applications is that including the complexity of GRIN optics in the model is not necessary for reproducing adequately the effect of GRIN optics on wave-front aberrations for various locations in the visual field or for different states of accommodation. What our study demonstrates is that the physical dimensions and overall layout of eye models produced by OWT are consistent with known anatomy of eyes, and the optical aberrations of the model are equivalent to those of the individual eye used to create the model. 

Although the customized eye models captured the major features of accommodation, they fell short of accounting for some optical features of eyes described later. Further analysis revealed that these shortcomings could be rectified by relaxing anatomical constraints to favor functional equivalence over anatomical similarity. 

Balancing the dual goals of OWT

The greatest inaccuracy of our customized models was the consistent overestimation of spherical aberration in the positive direction (about 0.1 μm for subject RH, Figure 9B) for all young subjects (0.076 μm on average for our population). Accuracy was not improved by including lens refractive index as an optimization parameter in a sample of five eyes with emmetropia and three with myopia. Since the amount of overestimation was independent of accommodation, correcting for this mean offset would be an expedient way to restore functional equivalence that might be acceptable for some applications of the model. However, we would prefer to correct the underlying cause of this systematic error, which we suspected to be a consequence of our emphasis on anatomical similarity in the OWT merit function used to optimize model parameters. 

As mentioned in the Results, the curvatures of both lenticular surfaces shown in Figure 10A escaped from the soft anatomical constraints (shaded areas) imposed by the merit function. This was an anticipated consequence of using a lens with homogeneous refractive index because greater curvature is needed to match the power of a GRIN lens (Navarro et al., 2006; Wei & Thibos, 2008), and that increased curvature also affects spherical aberration. Additional, exploratory modeling revealed that anterior lens curvature influences ocular spherical aberration more than posterior lens curvature. After loosening the constraints on anterior lens curvature in the merit function, we repeated the optimization for subject RH; the results are shown by blue triangle symbols in Figure 9. The refractive state of this revised model (Figure 9A) again matched the measured values closely. More importantly, the new model solved the problem of overestimation of spherical aberration with a negligible mean underestimation of 0.013 μm. Moreover, the averaged RMS fitting errors of astigmatism and coma diminished by 15% and 53%, respectively, for the revised model (Figure 9C and 9D). As a result, the variation of HOA RMS with accommodation for the revised model (Figure 9E) more closely matched the experimental data, with a constant offset of 0.12 μm. 

This example demonstrates that the cost of improving functional equivalence is reduced anatomical similarity for curvature of the anterior lens surface. In the revised model (black triangles, Figure 10A), curvature increased by 30% (which increased paraxial refractive power by 3.3D) and is greater than expected for human eyes (gray shaded area). However, to maintain the correct overall ocular power, the curvature of the posterior lens decreased in the revised model (red triangles), thereby becoming more anatomically reasonable (i.e., closer to the red shaded area). Lens thickness and ACD of the revised model remained anatomically reasonable (Figure 10B). The lens of the revised model was slightly more decentered (Figure 10C) and tilted (Figure 10D) slightly more toward the inferior side. 

In summary, by allowing the anterior lens surface to become more curved than human lenses, we improved the functional equivalence of the model for subject RH, especially for spherical aberration, and all the remaining structural parameters were either unchanged or moved closer to anatomical norms. We did not construct revised eye models for other subjects in our population, but that effort would appear to be justified for future applications that prioritize functional equivalence over anatomical similarity in eye models. 

It is worth noting that spherical aberration across the visual field was accurately represented by the OWT models derived for our two subjects with presbyopia (Figure 2, underestimated spherical aberration = 0.036 μm for LNT and 0.031 μm for AB). We suspect this good agreement might arise from the change of the GRIN distribution with age (Dubbelman & Van der Heijde, 2001; Hemenger, Garner, & Ooi, 1995). 

Accounting for accommodation error

All of our aberration measurements and model parameters have been described in terms of accommodation demand, which is the dioptric difference between the eye's far point and the accommodation stimulus. For some applications it may be more useful to characterize the model in terms of accommodative response, which is the change in refractive state induced by the stimulus. Response is typically less than demand (López-Gil et al., 2013), and this lag of accommodation is known for all of the subjects in our study populations (Liu et al., 2016). By compensating for foveal lag we reexamined the effect of accommodation on surface curvature of the lens in the population average model (Figure 11A). The linear relationships described earlier persisted, but the regression slope for anterior curvature increased from 0.009 to 0.013 mm⁻¹ per diopter of accommodation response, while the slope for posterior curvature increased from 0.006 to 0.009 mm⁻¹/D. Regression slopes for ACD (Figure 11B) changed from −0.047 to −0.061 mm/D, while regression slopes for lens thickness changed from 0.047 to 0.054 mm/D. Thus in every case, the rate of change for model parameters is greater when expressed in terms of accommodation response rather than accommodation demand. 

Uniformity of accommodation predicted by customized models

One of the remarkable features of the accommodative response maps displayed in Figures 4 (subject RH) and 12 (population average) is their uniformity of color, which indicates that variation of refractive state across the central visual field is small compared to the range of accommodation investigated. Uniformity of refractive state, in turn, indicates that the retinal conjugate surface for the macular region is spherical for all states of accommodation (Liu et al., 2016). Therefore, to visualize the subtle variations present in these maps requires subtracting some reference values to produce differential maps. For example, local accommodation maps reveal the change in refractive state at each location in the visual field, whereas accommodation error maps reference refractive state at each visual-field location relative to a spherical target surface of constant vergence. Both types of relative maps are shown for subject RH in Figure 13, with experimental data in the top row and model predictions in the bottom row. The local accommodation map (left column) shows that an irregular pattern of about 0.5D variation in the measurements is represented by a smooth variation in the model. Similarly, the irregular pattern of accommodative error (right column) evident in the measurements is represented by a smooth variation in the model. This failure of the model to represent fine details in these relative accommodative maps is not unexpected, since the optical components of the model are smooth, idealized elements that contain none of the local irregularities that might be expected in biological optical elements. 

Figure 13

 

Local accommodation and accommodation error field maps for subject RH at 6D accommodation demand. The experimental data are in the top row and model predictions are in the bottom row. Graphic conventions are the same as in Figure 3. The local accommodation map (left column) reveals the change in refractive state at each location in the visual field, while the accommodation error map references refractive state to a spherical target surface of constant vergence.

Pan-retinal measurement of image quality

One practical application of customized eye models is to compute the optical limits to vision. Psychophysical evaluations are slow and tedious by comparison to the rapid data acquisition by a wave-front aberrometer. Thus an accurate optical model of ocular aberrations across the visual field derived by OWT may be especially useful for predicting the visual effects of those aberrations, especially in the peripheral field, where behavioral psychophysics is challenging even for trained subjects. Practical attempts to determine the optical limits to vision typically employ a scalar metric of image quality that is associated with some measure of visual performance (Thibos, Hong, Bradley, & Applegate, 2004). Although these metrics can be calculated from empirical wave-front measurements, an optical model of an individual's eye makes it possible to interpolate and extrapolate to visual-field locations not measured. Image-quality metrics are also useful for designing and evaluating ophthalmic devices for clinical treatment of the eye's optical flaws. The models we report under Results describe a person's aberration so accurately that the model itself is essentially a “prescription for perfection” in the sense that a hypothetical, ideal treatment that perfectly corrects the model's aberrations would, when perfectly implemented, also render the human eye well corrected even for higher order aberrations across the central visual field. 

Historically, foveal vision has been given highest priority for assessment of optical limits to function, but current theories of emmetropization and myopia progression have redirected attention to peripheral refractive errors (E. L. Smith, 2011). Animal experiments have demonstrated that hyperopic and myopic blur have opposite effects on the development of refractive error, which implies that the sign of defocus must be sensed by visual mechanisms responsible for regulating eye growth. Researchers (Wallman & Winawer, 2004) have hypothesized that retinal activity inhibits eye growth, which implicates image contrast as the key feature of retinal images relevant to myopia development (Thibos, Bradley, Liu, & López-Gil, 2013). Other experiments (Norton, Essinger, & Mcbrien, 1994; Raviola & Wiesel, 1990; Troilo & Wallman, 1991) point to the outer retina as the critical site for regulatory retinal activity, which in turn suggests the intriguing possibility that cone photoreceptors may signal the presence of image contrast directly to the mechanisms regulating eye growth (Crewther, 2000). 

To illustrate how a wide-field optical model may be used to approach the problem of specifying image quality as seen by cone photoreceptors over a significant portion of the visual field, we computed field maps of cone activation (see Equation 1) for the average emmetropic-eye model viewing an idealized point source. The resulting maps of cone activation displayed in Figure 14 for a sequence of accommodative demands reveal how image quality, as seen by cone photoreceptors, varies across the visual field. For any given accommodative state, cone activation is higher in the fovea than the surrounding region because most aberrations increase with eccentricity and peripheral cones have larger diameter. As the eye accommodates, the variation across the visual field becomes less noticeable and cone activation maps become more uniform, indicating that image quality is dominated by a uniform accommodative error (defined as the difference between absolute refractive state and foveal target vergence—e.g., map D defined previously; Liu et al., 2016), with negligible contribution from other aberrations or cone-diameter variations. 

Figure 14

 

Visual-field maps of cone activation at different accommodation states computed using the average emmetropic-eye model. Graphic conventions are the same as in Figure 3, except that the −1D accommodation demand field map was replaced with results calculated with a diffraction-limited optical system.

Integrating cone activation over the whole visual space provides a cumulative, pan-retinal measure of cone image quality (Equation 2) that declines more than 1 log unit over the 6D range of accommodative demand imposed by our experiments. As shown in Figure 15, this decline is even greater for foveal CIQ, calculated by restricting the area of integration in Equations 2 and 3 to a small region in the foveal center. These results indicate that normal levels of accommodative lag under monocular viewing conditions, coupled with other aberrations and the normal sign reversal of spherical aberration with accommodation, reduce CIQ by increasing amounts as the eye accommodates. Under binocular viewing conditions, CIQ would be expected to decline less than is shown in Figure 5, because accommodative lag is typically less for binocular compared to monocular viewing (López-Gil et al., 2013). 

Figure 15

 

Comparison of cone image quality of an average emmetropic eye at different accommodation states while viewing an ideal point source (black symbols) or a filtered point source with a 1/f spectrum typically found in natural scenes (red symbols). Filled symbols show integrated cone image quality across the central 30° of the visual field. Empty symbols show foveal cone image quality computed by the same method, but integrated over a small, central field of view.

For comparison, we repeated the calculation of CIQ for a filtered point source with a 1/f spectrum (i.e., contrast varying inversely with spatial frequency), as frequently occurs in natural scenes (Ruderman & Bialek, 1994). CIQ is an order of magnitude lower for the 1/f spectrum (red symbols) compared to the ideal point source (black symbols) at all accommodative states because object contrast at high spatial frequencies is attenuated in a natural scene. Nevertheless, for both stimuli, CIQ over the central 30° field of view (filled symbols) is less than (or equal to) foveal CIQ (open symbols) because of peripheral aberrations (mainly coma and oblique astigmatism). As accommodative demand increases, defocus due to accommodative lag dominates the retinal image everywhere in the visual field, causing CIQ to fall for both types of visual stimulus. This result suggests that mitigating accommodation error to improve image quality across the visual field may be a better strategy for increasing retinal activity and retarding myopia progression than the common strategy of deliberately introducing peripheral relative myopia (Wildsoet et al., 2019), which will leave image quality and cone image quality suboptimal. 

Conclusions

Customized optical models of the eye that accurately reproduce the wave-front aberrations of individual eyes across the visual field at different accommodative states provide a mechanistic explanation for the changes in the eye's aberration structure due to accommodation. Applications of such models include computing retinal image quality and the optical limits to vision, designing devices to modify or correct the eye's optical flaws, and understanding the role of vision in the eye's control of its own growth. 

Acknowledgments

Commercial relationships: none. 

Corresponding author: Tao Liu. 

Email: . 

Address: School of Optometry, Indiana University, Bloomington, IN, USA. 

Appendix

The authors supply Supplementary File S1 of the Zemax configuration files required to implement the hypothetical, average-eye model for seven levels of accommodative demand (indicated by the last number in each file name). An additional MATLAB file gives corresponding aberration coefficients of the model reported by Zemax and converted to the ANSI Z80.28 standard convention for a fixed pupil diameter of 6 mm.