1.

Introduction

The human eye, when considered as an optical system, is a separate doublet, i.e. the cornea and crystalline lens, which have an aperture stop in contact with the first surface of the latter. However, like many optical systems, the human eye is not free from aberrations that limit retinal-image quality. Furthermore, the eye does not maintain its optical geometry with aging, but rather the corneal and lens surfaces, lens thickness, and the spacing between the lens and cornea, change as the individual gets older.^1–5 Moreover, it is also known that the lens has a gradient-index (GRIN) structure,^6–11 and this also changes with age.^12–27 For all this, the ocular aberrations change over the human lifespan, as several studies have demonstrated.^18–31

A fundamental question arising from the studies regarding retinal-image quality concerns the relative contribution to the ocular aberrations by the cornea and lens. In the design of an optical system, the individual elements compensate for each other to form a system having a reasonable optical-image quality. This is an important feature for the human eye regarding ophthalmic clinical applications, such as the implant of an intraocular lens (IOL)^32–34 or refractive surgery.^35–40 Vertebrate crystalline lenses have an internal gradient of refractive index which reduces aberrations, most notably spherical aberration, and increases the mean refractive power of the lens.^41,42 This gradient reaches its maximum in the eyes of fish where the lens is spherical.^43,44

Studies addressing this aberration reduction have reported a balance between the corneal aberrations and those from the internal optics, thus resulting in a smaller amount of ocular aberrations for the whole eye,^27,45–50 although prior studies suggested that no such compensation exists.^45,51,52 Specifically, there are works demonstrating that the spherical aberration of the lens compensates for that of the cornea,^53,54 as well as the lateral (horizontal) third-order coma and horizontal/vertical (H/V) astigmatism.^54,55

Studies on ocular-aberration variations with age show this compensation is altered between the cornea and the lens in the eye.^{22,25–27,31,47,49,54,56–58} Thus, the lens cannot compensate for the positive spherical aberration of the cornea with increasing age, and the spherical aberration of the eye becomes more positive with age. This effect is attributed to a passive mechanism resulting from a genetically determined physiology.⁵⁷ A third-order coma is another ocular aberration which is partially compensated for by the internal ocular optics, and an active mechanism of compensation related with the angle kappa has been discussed.^{31,47,54,57,58} Thus, it seems a fine-tuned misalignment between the lens and the cornea is responsible for that compensation, as has been modelled in emmetropic, myopic, and hyperopic subjects,⁴⁹ as well as pseudophakic ones.⁵⁷ Regarding astigmatism, the compensation of corneal astigmatism by internal optics is well-known,⁵⁹ and some studies support the idea that this defect is compensated for individually, and possibly develops with age.^54,60

Another key question arising from these studies concerns the role the lens GRIN profile could play in that compensation. The results of modelling the effect of the compensation by using eye models with homogeneous and GRIN lens have shown different behavior.^49,54 Furthermore, there are experimental studies using magnetic resonance imaging (MRI) which show the lens GRIN profile changes with age. Thus, the lens has a plateau of constant index in its center that increases in thickness with aging, and the GRIN variation would take place in a thinner region in the periphery of the lens. ¹⁶ Thus, it would be worthwhile to study the separate contribution of the GRIN profile to the whole-lens compensation as the lens varies its geometry and GRIN profile as the human eye gets older. Experimentally, this contribution is difficult to separate from that of the lens surfaces, and thus, a theoretical approach becomes necessary.

This leads to the use of an accurate aging-eye model including a GRIN lens for the calculations. Recently, several works have been published proposing an aging GRIN lens to be incorporated in an eye model.^61–64 The one proposed by Goncharov et al.⁶¹ does not come from a systematic study with the age, and their proposed eye model is a centred optical system, which would introduce errors in the calculation of the compensation of ocular aberrations, since the misalignments between the lens and the cornea have been proven to drive it. Navarro et al.^62,63 proposed a lens GRIN profile which accurately reproduces the experimental data from MRI measurements. However, since the profile fits in vitro measurements, the model is representative for a full accommodative state of the lens, and therefore is not appropriate to model an in vivo eye. Moreover, the proposed lens is not included in a complete eye model showing misalignments between the ocular elements.

The recent eye model published by Díaz et al.⁶⁴ is, as far we know, the most up-to-date model accurately mimicking the average optical quality of the emmetropic population, and it includes a GRIN lens profile that also varies with age. This model is not only a chromatic one, but a decentered model also, which takes into account several experimental measurements regarding the relative misalignment of the different ocular surfaces, and could be useful to model the mechanism of compensation for ocular aberrations.⁶⁵

The purpose of this work is to evaluate the contribution of the lens GRIN profile to the overall compensation of the ocular aberrations between the cornea and lens throughout life. The study seeks to determine the main third-order ocular aberrations, i.e. spherical, lateral (horizontal) coma, and H/V astigmatism. This will provide better knowledge of the optical performance of the eye with aging, and improve the modelling of the proposed lens GRIN profile for a more accurate human-eye model.

2.

Methods

This section is divided in two parts. First, the eye model used is briefly described, and second, the calculations of third-order aberration coefficients is detailed, considering the misalignment of the cornea and the lens with respect to the iris, as well as the effect of the lens GRIN profile.

2.1.

Human Eye Model

The emmetropic eye model used is representative of the average population at different ages,⁶⁴ working with a 555-nm monochromatic light and a pupil diameter of 6 mm. Briefly, the cornea is an element with rotational symmetry, the iris is decentered nasally with respect to the corneal axis, and the lens is tilted. The lens refractive index is modelled by a GRIN distribution that varies both radially and axially. The thickness, the curvature of surfaces, and the lens GRIN profile vary with age between 20 and 60 years of age. The data are tabulated in Table 1, and a layout for the upper view of the right eye at the age of 25 years is shown in Fig. 1.

Fig. 1

Upper view layout corresponding to the right eye at the age of 25 years. The gradient-index (GRIN) profile is also shown at that age. The isoindicial contours from the edge to the core correspond to the values 1.371, 1.381, 1.391, 1.401, and 1.411 The evolution of the GRIN profile in the model eye at other ages can be seen in the Ref. 64. Numbers at each ocular surface are only for indexing purposes.

Table 1

Schematic eye model parameters as a function of age A (in years).

Surface	Medium	n (555 nm)	R (mm)	Asphericity	Thickness (mm)
1	Cornea	1.376	7.79	$-0.24+0.003·A$	0.579
2	Aqueous	1.336	6.53	$-0.006·A$	$3.291-0.01·A$
3 (iris)b	Aqueous	1.336		—	0
4c	Lens	n(z)a	1$12.7-0.058·A$	$-4.56$	$2.93+0.0236·A$
5	Vitreous	1.336	$-5.9+0.0015·A$	$-1.13$	$16.79-0.0136·A$
6 (retina)d	—	—	$R_x=-12.91$$R_y=-12.72$	$Q_x=0.27$$Q_y=0.25$	—

^aThe equation of the index is n(z)=1.371+n1(cos(n2z)−1)+n3 sin(n4z)+n5(x2+y2), and the index coefficients, n1, n2, n3, n4, and n5 depend on age A, which are tabulated in Table 2 in Ref. 64.

^bThe iris is decentered −0.3  mm along the x-axis; that is, it is decentered 0.3 mm nasally.

^cThe lens is tilted −4  deg about the y-axis containing its center; that is, the lens axis is 4 deg to the temporal space object.

^dThe retina is an ellipsoid tilted −11.5  deg about the y-axis and −3.6  deg about the x-axis.

2.2.

Calculation of Third-Order Coefficients

The third-order aberrations were calculated by means of the Seidel coefficients. On the one hand, these coefficients enable calculations of the individual contribution of the surface geometry, the GRIN media⁶⁶ as well as that of the tilts and decentrations in the ocular optical components.⁶⁷ On the other hand, the Seidel coefficients can be related to Zernike coefficients.^68,69

Therefore, their calculus has been divided in two parts. First, a paraxial raytracing of the principal and marginal rays was performed through the whole-eye model. This enabled calculation of Seidel coefficients for each aspherical ocular surface. Moreover, this enabled a calculation of the contribution by the lens GRIN distribution. This contribution can be divided into three terms: one due to the propagation of light within the lens (transfer contribution), and two terms due to the refraction of the wavefront on each surface of the lens (contribution to refraction).^66,70,71 Secondly, the effect that the relative misalignment of the cornea and lens exerts on the Seidel coefficients was calculated using expressions reported elsewhere.⁶⁷ This effect was calculated by adding a term to the Seidel coefficients which were calculated when all ocular elements were assumed to be centered.

Thus, the total ocular third-order Seidel spherical, SI, coma, SII, and astigmatism, SIII, coefficient values are given by:

Eq. (1)

$$\begin{gathered} \mathrm{S}{\mathrm{I}}_{\text{eye}}=\mathrm{S}{\mathrm{I}}_{\mathrm{GRIN}}+\mathrm{S}{\mathrm{I}}_{\text{decenter}/\mathrm{tilt}} \\ \mathrm{S}{\mathrm{II}}_{\text{eye}}=\mathrm{S}{\mathrm{II}}_{\mathrm{GRIN}}+\mathrm{S}{\mathrm{II}}_{\text{decenter}/\mathrm{tilt}} \\ \mathrm{S}{\mathrm{III}}_{\text{eye}}=\mathrm{S}{\mathrm{III}}_{\mathrm{GRIN}}+\mathrm{S}{\mathrm{III}}_{\text{decenter}/\mathrm{tilt}} \end{gathered}$$

into which the GRIN profile contribution and that of the tilted and/or decentered components are divided. These terms are calculated as detailed below.

The commercial software ZEMAX-EE⁷² was used to perform the raytracing. We used the right eye, and therefore considered the fovea was at 5 deg. in the temporal retinal side.

2.2.1.

Seidel coefficients for aspherical ocular surfaces

All ocular surfaces in the eye model have a conical shape. Therefore, each surface was considered centered in the expressions used to determine the Seidel coefficient value for each surface. These will be useful later when considering their misalignment. The expressions are:

Eq. (2)

$$\begin{gathered} \mathrm{S}{\mathrm{I}}_{k,\mathrm{cent}}=-{(n_k{i}_{k,m})}^2{y}_{k,m}(\frac{u_{k,m}^{\prime }}{n_k^{\prime }}-\frac{u_{k,m}}{n_k})+(n_k^{\prime }-n_k)Q_k{c}_{k,0}^3 \\ \mathrm{S}{\mathrm{II}}_{k,\mathrm{cent}}=-(n_k^2{i}_{k,m}{i}_{k,p})y_{k,m}(\frac{u_{k,m}^{\prime }}{n_k^{\prime }}-\frac{u_{k,m}}{n_k})+\left(\frac{y_{k,p}}{y_{k,m}}\right)(n_k^{\prime }-n_k)Q_k{c}_{k,0}^3 \\ \mathrm{S}{\mathrm{III}}_{k,\mathrm{cent}}=-{(n_k{i}_{k,p})}^2{y}_{k,m}(\frac{u_{k,m}^{\prime }}{n_k^{\prime }}-\frac{u_{k,m}}{n_k})+{\left(\frac{y_{k,p}}{y_{k,m}}\right)}^2(n_k^{\prime }-n_k)Q_k{c}_{k,0}^3, \end{gathered}$$

in which $i$ and $u$ are the incidence and the slope angles, of the paraxial raytracing, respectively, where the sub-indices $m$ and $p$ correspond to the axial marginal and principal rays; $n$ is the refractive index of the incidence medium at the pole of the surface; $Q_k$ is the asphericity of the surface; $c_{k,0}$ the curvature of the surface; and the primes are used for the corresponding values of the variables after refraction. The sub-index $k$ identifies the particular ocular surface (Fig. 1).

2.2.2.

Calculation of decentration and/or tilt contribution for the cornea and lens

The second term in Eq. (1) is related to the effect of the misalignment of the cornea and the lens relative to the iris. The cornea is decentered, and the lens is decentered and tilted with respect to the iris. Therefore, two raytracings are needed in order to establish the value for the Seidel coefficients at each surface given in Eq. (2). The first raytrace yields the slope and height of the marginal and principal rays in each surface considering the whole eye as a centred system. This was done by using ZEMAX-EE up to the anterior lens surface, and by using the ABCD matrix for the transfer of light to the posterior lens surface (next subsection). Finally, the standard paraxial equations for refraction were applied to the second lens surface. The second raytrace is a finite one that corresponds to a ray propagating through the decentered system as if it was the optical axis of the centered system.⁶⁷ This was done entirely by using ZEMAX-EE.

Regarding the cornea, it is a decentered element with respect to the iris, and for the anterior surface, the object is at infinity. Thus, the Seidel coefficients are calculated in a slightly different way, as proposed in Ref. 67. For the anterior corneal surface, a decentration does not change the angle of the principal ray, just its height. Therefore, we must only apply a shift of the stop for getting the original height of the principal ray before the decentering [Fig. 2(a)]. Thus, the Seidel coefficients for the decentered anterior surface are:

Eq. (3)

$$\begin{gathered} \mathrm{S}{\mathrm{I}}_{1,\mathrm{dec}}^{\text{cornea}}=\mathrm{S}{\mathrm{I}}_{1,\mathrm{cent}}^{\text{cornea}} \\ \mathrm{S}{\mathrm{II}}_{1,\mathrm{dec}}^{\text{cornea}}=\mathrm{S}{\mathrm{II}}_{1,cent}^{\text{cornea}}+\frac{\mathrm{\Delta }{y}_p}{y_m}\mathrm{S}{\mathrm{I}}_{1,\mathrm{cent}}^{\text{cornea}} \\ \mathrm{S}{\mathrm{II}}_{1,\mathrm{dec}}^{\text{cornea}}=\mathrm{S}{\mathrm{III}}_{1,\mathrm{cent}}^{\text{cornea}}+2\left(\frac{\mathrm{\Delta }{y}_p}{y_m}\right)\mathrm{S}{\mathrm{II}}_{1,\mathrm{cent}}^{\text{cornea}}+{\left(\frac{\mathrm{\Delta }{y}_p}{y_m}\right)}^2\mathrm{S}{\mathrm{I}}_{1,\mathrm{cent}}^{\text{cornea}}, \end{gathered}$$

in which $y_m$ is the height of the marginal ray at the first corneal surface, and $\mathrm{\Delta }{y}_p=y_p^{*}-y_p=e$ is the change in the principal ray height, due to the position shift of the entrance pupil to attain the old height [Fig. 2(a)]. In the case of the posterior corneal surface, the method reported⁶⁷ was followed, and the Seidel coefficients were calculated by:

Eq. (4)

$$\begin{gathered} \mathrm{S}{\mathrm{I}}_{2,\mathrm{dec}}^{\text{cornea}}=\mathrm{S}{\mathrm{I}}_{2,\mathrm{cent}}^{\text{cornea}} \\ \mathrm{S}{\mathrm{II}}_{2,\mathrm{dec}}^{\text{cornea}}=\left(\frac{h+e}{h_0}\right)\mathrm{S}{\mathrm{II}}_{2,\mathrm{cent}}^{\text{cornea}}+\left(\frac{e}{ul}\right)\mathrm{S}{\mathrm{I}}_{2,\mathrm{cent}}^{\text{cornea}} \\ \mathrm{S}{\mathrm{II}}_{2,\mathrm{dec}}^{\text{cornea}}={\left(\frac{h+e}{h_0}\right)}^2\mathrm{S}{\mathrm{III}}_{2,\mathrm{cent}}^{\text{cornea}}+2\left(\frac{e}{ul}\right)\left(\frac{h+e}{h_0}\right)\mathrm{S}{\mathrm{II}}_{2,\mathrm{cent}}^{\text{cornea}}+{\left(\frac{e}{ul}\right)}^2\mathrm{S}{\mathrm{I}}_{2,\mathrm{cent}}^{\text{cornea}}, \end{gathered}$$

in which $h_0$ is the field (object) position, and $h$ is that position after decentration, $e$ is the decentration, $l$ is the distance of the entrance pupil to the surface, and $u$ is the marginal axial ray angle [Fig. 2(b)].

Fig. 2

(a) The anterior corneal surface decentered. In this case the object is at infinity. (b) The posterior corneal surface decentered. In this case the object is at a finite distance, and the method in Ref. 67 can be applied.

The lens is a decentered and tilted element. Therefore, the corresponding coefficients values for each surface ($i=3$, $4$) are calculated using the following equations,⁶⁷ from those when the lens is centered [Eq. (2)]:

Eq. (5)

$$\begin{gathered} \mathrm{S}{\mathrm{I}}_{i,d/t}^{\text{lens}}=\mathrm{S}{\mathrm{I}}_{i,\mathrm{cent}}^{\text{lens}} \\ \mathrm{S}{\mathrm{II}}_{i,d/t}^{\text{lens}}=\left[\frac{y_{p,i}-u_{p,i}D+l(\frac{y_{p,i}}{u_{p,i}}-\frac{y_{z,i}}{u_{z,i}}+D)}{-\frac{y_{p,i}}{u_{p,i}}+\frac{y_{m,i}}{u_{m,i}}}\right]\mathrm{S}{\mathrm{II}}_{i,\mathrm{cent}}^{\text{lens}}+\left[\frac{l(\frac{y_{p,i}}{u_{p,i}}-\frac{y_{z,i}}{u_{z,i}})}{u_{m,i}D}\right]\mathrm{S}{\mathrm{I}}_{i,\mathrm{cent}}^{\mathrm{ens}} \\ \mathrm{S}{\mathrm{III}}_{i,d/t}^{\text{lens}}={\left[\frac{y_{p,i}-u_{p,i}D+l(\frac{y_{p,i}}{u_{p,i}}-\frac{y_{z,i}}{u_{z,i}}+D)}{-\frac{y_{p,i}}{u_{p,i}}+\frac{y_{m,i}}{u_{m,i}}}\right]}^2\mathrm{S}{\mathrm{III}}_{i,\mathrm{cent}}^{\text{lens}}+2\left[\frac{l(\frac{y_{p,i}}{u_{p,i}}-\frac{y_{z,i}}{u_{z,i}})}{u_{m,i}D}\right] \\ \times \left[\frac{y_{p,i}-u_{p,i}D+l(\frac{y_{p,i}}{u_{p,i}}-\frac{y_{z,i}}{u_{z,i}}+D)}{-\frac{y_{p,i}}{u_{p,i}}+\frac{y_{m,i}}{u_{m,i}}}\right]S{\mathrm{II}}_{i,\mathrm{cent}}^{\text{lens}}+{\left[\frac{l(\frac{y_{p,i}}{u_{p,i}}-\frac{y_{z,i}}{u_{z,i}})}{u_{m,i}D}\right]}^{2}S{\mathrm{I}}_{i,\mathrm{cent}}^{\text{lens}}, \end{gathered}$$

in which ($u_{m,i}$, $y_{m,i}$) and ($u_{p,i}$, $y_{p,i}$) are the slope and the heights, respectively, of the marginal, principal rays at the $i$ ($=3$, $4$) surface established by paraxial raytracing; ($u_{z,i}$, $y_{z,i}$) is the slope and height, respectively, of a ray at each lens surface as if it was the optical axis determined by finite raytracing (Fig. 3); $D$ is the distance to the entrance pupil (iris) from the intersection of the tilted axis of the lens with that of the entrance pupil, and finally, $l$ is the distance from the object to the entrance pupil, (Fig. 4).

Fig. 3

Finite raytracing corresponding to a ray coinciding with the optical axis when the system is centered.⁶⁷ Only the angle and height values at each crystalline lens surfaces are needed when applying Eq. (5).

Fig. 4

Parameters defined when a lens is decentered and tilted in applying Eq. (5), following Ref. 67.

2.2.3.

Calculation of GRIN profile contribution

Since the human crystalline lens is a GRIN one, we know there is a contribution to the Seidel coefficients when the light refracts from/to a GRIN media, as well as another contribution from the propagation of light inside the GRIN element to the posterior surface.⁶⁶ The contribution to the refraction has a similar effect of working with an aspherical surface instead of a spherical one in the case of dealing with homogeneous media. Thus, we can calculate the first term of Eq. (1) as follows:

Eq. (6)

$$\begin{gathered} \mathrm{S}{\mathrm{I}}_{\mathrm{GRIN}}=\sum _{i=3}^4\mathrm{S}{\mathrm{I}}_{i,\mathrm{GRN}}^{\mathrm{refrac}}+S{\mathrm{I}}_{\mathrm{GRIN}}^{\text{transfer}} \\ \mathrm{S}{\mathrm{II}}_{\mathrm{GRIN}}=\sum _{i=3}^4\mathrm{S}{\mathrm{II}}_{i,\mathrm{GRN}}^{\mathrm{refrac}}+\mathrm{S}{\mathrm{II}}_{\mathrm{GRIN}}^{\text{transfer}} \\ \mathrm{S}{\mathrm{III}}_{\mathrm{GRIN}}=\sum _{i=3}^4\mathrm{S}{\mathrm{III}}_{i,\mathrm{GRN}}^{\mathrm{refrac}}+\mathrm{S}{\mathrm{III}}_{\mathrm{GRIN}}^{\text{transfer}}, \end{gathered}$$

in which

Eq. (7)

$$\begin{gathered} \mathrm{S}{\mathrm{I}}_{3,\mathrm{grin}}^{\mathrm{refrac}}=y_{3,m}^4(0)\left[\right({4c_{3,0}{n}_1(0)+c_{3,0}^2\frac{\partial n_0(z)}{\partial z}|}_{z=0})+2c_{3,1}(n_0(0)-n_{\text{aqueous}})] \\ \mathrm{S}{\mathrm{II}}_{3,\mathrm{grin}}^{\mathrm{refrac}}=y_{3,m}^3(0)y_{3,p}(0)\left[\right({4c_{3,0}{n}_1(0)+c_{3,0}^2\frac{\partial n_0(z)}{\partial z}|}_{z=0})+2c_{3,1}(n_0(0)-n_{\text{aqueous}})] \\ \mathrm{S}{\mathrm{III}}_{3,\mathrm{grin}}^{\mathrm{refrac}}=y_{3,m}^2(0)y_{3,p}^2(0)\left[\right({4c_{3,0}{n}_1(0)+c_{3,0}^2\frac{\partial n_0(z)}{\partial z}|}_{z=0})+2c_{3,1}(n_0(0)-n_{\text{aqueous}})] \\ \mathrm{S}{\mathrm{I}}_{4,\mathrm{grin}}^{\mathrm{refrac}}=-y_{4,m}^4(d)\left[\right({4c_{4,0}{n}_1(d)+c_{4,0}^2\frac{\partial n_0(z)}{\partial z}|}_{z=d})+2c_{4,1}(n_{\text{vitreous}}-n_0(d))] \\ \mathrm{S}{\mathrm{II}}_{4,\mathrm{grin}}^{\mathrm{refrac}}=-y_{4,m}^3(d)y_{4,p}(d)\left[\right({4c_{4,0}{n}_1(d)+c_{4,0}^2\frac{\partial n_0(z)}{\partial z}|}_{z=d})+2c_{4,1}(n_{\text{vitreous}}-n_0(d))] \\ \mathrm{S}{\mathrm{III}}_{4,\mathrm{grin}}^{\mathrm{refrac}}=-y_{4,m}^2(d)y_{4,p}^2(d)\left[{(4c_{4,0}{n}_1(d)+c_{4,0}^2\frac{\partial n_0(z)}{\partial z}|}_{z=d}\right)+2c_{4,1}(n_{\text{vitreous}}-n_0(d))] \end{gathered}$$

are the GRIN profile contributions due to refraction at each crystalline lens surface, and

Eq. (8)

$$\begin{gathered} \mathrm{S}{\mathrm{I}}_{\mathrm{GRIN}}^{\text{transfer}}=+N_0(0)y_m(0)u_m^3(0)-N_0(d)y_m(d)u_m^3(d) \\ -2{\int }_0^{d}4N_2(z)y_m^4(z)+2{\int }_0^d{N}_1(z)y_m^2(z)u_m^2(z)\mathrm{d}z \\ -\frac{1}{2}{\int }_0^d{N}_0(z)u_m^4(z)\mathrm{d}z \\ \mathrm{S}{\mathrm{II}}_{\mathrm{GRIN}}^{\text{transfer}}=+N_0(0)y_m(0)u_m^2(0)u_p^2(0)-N_0(d)y_m(d)u_m^2(d)u_p^2(d) \\ -2{\int }_0^{d}4N_2(z)y_m^3(z)y_p(z)+{\int }_0^d{N}_1(z)u_m(z)y_p(z)(y_m(z)u_p(z)+u_m(z)y_p(z))\mathrm{d}z \\ -\frac{1}{2}{\int }_0^d{N}_0(z)u_m^3(z)u_p(z)\mathrm{d}z \\ \mathrm{S}{\mathrm{III}}_{\mathrm{GRIN}}^{\text{transfer}}=N_0(0)y_m(0)u_m(0)u_p^2(0)-N_0(d)y_m(d)u_m(d)u_p^2(d) \\ -2{\int }_0^{d}4N_2(z)y_m^2(z)y_p^2(z)\mathrm{d}z+2{\int }_0^d{N}_1(z)y_m(z)u_m^2(z)y_p(z)u_p^2(z)\mathrm{d}z \\ -\frac{1}{2}{\int }_0^d{N}_0(z)u_m^2(z)u_p^2(z)\mathrm{d}z \end{gathered}$$

are the GRIN profile contributions due to light propagating through the crystalline lens to the posterior surface.

In Eqs. (7) and (8), $d$ is the thickness of the lens; $u(z)$ and $y(z)$ are the slope and the height, respectively, of the paraxial rays at a position $z$ inside the lens, where the sub-indices $m$ and $p$ refer to the marginal and principal rays, respectively; the parameters $c_{i,1}=-1/2c_{i,0}^3(\sqrt{-Q_i}+1)\text{with}i=3,4$, are the first off-axis curvature of each hyperbolic lens surface,⁷³ $c_{i,0}$ being its curvature and $Q_i$ its asphericity; $N_0$, $N_1$, and $N_2$ are the GRIN profile coefficients as a function of position $z$, assuming a rotationally symmetric GRIN medium, i.e. $n(r,z,\lambda )$, as is the case for the lens as:

Eq. (9)

$$n(z,r,\lambda )=n_0(z,\lambda )+\sum _j{N}_j(z,\lambda )r^{2j},$$

in which, for the crystalline lens model,⁶⁴ $r^2=x^2+y^2$, $N_1$ is a constant value and $N_j=0\text{if}j\ge 2$.

It is known that, up to the date of the software version used here (December, 2010), ZEMAX-EE did not correctly calculate the paraxial raytracing when the optical system dealt with GRIN media. Hence, the slope and the height of both marginal and principal rays were wrongly traced to the posterior surface of the lens. This means that these values cannot be used in Eq. (7) and (8). However, the use of the ABCD matrix of the lens^70,71 permits them to be determined; therefore, the contribution of the GRIN profile to the refraction at the posterior lens surface, as well as the light transfer, can be calculated.

Thus, the slope and the height of a ray at any point within a rotational symmetric GRIN medium can be determined by using the parabolic approximation as follows:⁷⁰

Eq. (10)

$$\left[\begin{array}{c}y(z)\\ u(z)\end{array}\right]=\left[\begin{array}{cc}A(z)& B(z)\\ C(z)& D(z)\end{array}\right]\left[\begin{array}{c}y(0)\\ u(0)\end{array}\right]={\mathrm{A}\mathrm{B}\mathrm{C}\mathrm{D}}_{\mathrm{GRIN}}(z)\left[\begin{array}{c}y(0)\\ u(0)\end{array}\right],$$

where $u(0)$ and $y(0)$ are the slope and the height, respectively, of the ray after refraction at the first surface of the lens. Therefore, the slope and the height of a paraxial ray at the posterior human-lens surface can be calculated by using

Eq. (11)

$$\left[\begin{array}{c}{y}_4(z)\\ u_4(z)\end{array}\right]=\left[\begin{array}{cc}A(z)& B(z)\\ C(z)& D(z)\end{array}\right]\left[\begin{array}{c}{y}_3(0)\\ u_3(0)\end{array}\right]={\mathrm{A}\mathrm{B}\mathrm{C}\mathrm{D}}_{\mathrm{GRIN}}(z)\left[\begin{array}{c}{y}_3(0)\\ u_3(0)\end{array}\right]$$

with $z=d$, $d$ being the thickness of the GRIN medium, i.e. the lens thickness at a given age, and where

Eq. (12)

$${\mathrm{ABCD}}_{\mathrm{GRIN}}(z)=\left[\begin{array}{cc}1+\beta z^2& \frac{z}{2}(1+\frac{N_o(0)}{N_o(z)}+\beta z^2)\\ 2\beta z& \frac{N_o(0)}{N_o(z)}+\beta z^2\end{array}\right]$$

with

Eq. (13)

$$\beta (z)=\frac{{\int }_0^z{N}_1(\omega )\mathrm{d}\omega }{d{N}_o(z)-2{\int }_0^z{N}_1(\omega ){\omega }^2\mathrm{d}\omega }.$$

Therefore, the Seidel coefficients for the cornea are calculated adding Eqs. (3) and (4):

Eq. (14)

$$X^{\text{cornea}}=\sum _{i=1}^2{X}_{i,\mathrm{dec}}^{\text{cornea}}.$$

Seidel coefficients for the lens are calculated by adding Eqs. (5), (7), and (8), this applied to each lens surface:

Eq. (15)

$$X^{\text{lens}}=\sum _{i=3}^4(X_{i,d/t}^{\text{lens}}+X_{i,\mathrm{GRIN}}^{\mathrm{refrac}})+X_{\mathrm{GRIN}}^{\text{transfer}},$$

and finally, Eq. (1) can be written by taking into account Eqs. (14) and (15) as:

Eq. (16)

$$X_{\text{eye}}=X^{\text{cornea}}+X^{\text{lens}}.$$

In Eqs. (14)–(16), $X$ stands for SI (spherical), SII (lateral coma) and SIII (H/V astigmatism) Seidel coefficients.

3.

Results

3.1.

Optical Surfaces and GRIN Separated Contributions to Ocular Aberrations

Figure 5 shows the results corresponding to the fourth-order spherical aberration, lateral coma, and H/V astigmatism. These were calculated for both aspherical surfaces of the decentered cornea and for those of the misaligned lens, the contribution of the GRIN profile to refraction in both lens surfaces, contribution to light transfer through the GRIN, and finally, for the whole eye.

Fig. 5

Coefficient values corresponding to the fourth-order spherical aberration (top), horizontal coma (middle), and H/V astigmatism (bottom), calculated for the cornea, lens surfaces, contribution of the gradient-index (GRIN) profile to refraction in both lens surfaces, contribution to light transfer through the GRIN lens, and finally, the complete eye.

It can be seen that the coefficients values calculated for the cornea change with age and range from 0.19 to 0.28 microns for $C_{4,\text{cornea}}^0$, from 0.20 to 0.10 microns for $C_{3,\text{cornea}}^1$, and from $-0.03$ to $-0.003$ microns for $C_{2,\text{cornea}}^2$. The variation found for the spherical aberration agrees with the data previously reported.

The values calculated for the lateral coma and H/V astigmatism should be taken with caution, since the cornea assumed in the model is spherical, and the experimental coefficients determined for the cornea are measured using different reference points,⁷⁸ and are often reported for the anterior corneal surface. In fact, the calculated values from its posterior surface for all the aberrations are small (about one order of magnitude). Nevertheless, the values reported in this work can be used to estimate the compensation made by the lens.

From the same figure, several findings can be highlighted since the calculation of the coefficients values allowed the separation of the contribution from the aspherical surfaces (corneal and lenticular surfaces) as well as that from the GRIN profile. Thus, it can be seen that the main contribution to the ocular coefficients values corresponding to $C_4^0$, $C_3^1$, and $C_2^2$ is given by the GRIN profile to the refraction at the lens surfaces, i.e. $C_{4,\mathrm{GRINS}}^0$, $C_{3,\mathrm{GRINS}}^1$, $C_{2,\mathrm{GRINS}}^2$. These values for the three ocular aberrations are negative and basically independent of the age, their magnitude being about $-0.37$ microns, $-0.38$ microns, and $-0.11$ microns, for the spherical, horizontal coma, and H/V astigmatism, respectively. The contribution corresponding to the light transfer through the GRIN profile does not seem to depend on aging either. This contribution has the value of 0.042 microns for $C_{4,\mathrm{GRINT}}^0$, −0.07 microns for $C_{3,\mathrm{GRINT}}^1$, and 0.035 microns for $C_{2,\mathrm{GRINT}}^2$. Finally, the values for $C_{4,\text{lens}}^0$, $C_{3,\text{lens}}^1$, and $C_{2,\text{lens}}^2$ decrease with age. That is, the spherical aberration associated with the optical geometry of the lens surfaces decreases clearly with age from 0.20 to 0.13 microns, 0.13 to 0.08 microns for the horizontal coma, and 0.018 to 0.01 microns for the H/V astigmatism.

3.2.

Compensation Between the Cornea and Lens

The calculation of the different contributions to the aberration coefficients enables an investigation of the effect of the aberration compensation between the lens and the cornea. Particularly, Fig. 6 presents the corneal spherical aberration value, $C_{4,\text{cornea}}^0$, and that of the lens, $C_{4,\text{lens}}^0+C_{4,\mathrm{GRINS}}^0+C_{4,\mathrm{GRINT}}^0$. It is clear that the corneal spherical aberration increases with the age. In addition, abberration of the lens also increases, but has the opposite sign and not enough increasing ratio as with the cornea. Therefore, the spherical aberration generated by the cornea will not be fully compensated for. However, the total ocular fourth-order aberration does not augment with aging (from 0.066 to 0.089 microns, with a mean value of 0.077 microns). This agrees with experimental studies such as those of Atchison and Markwell²⁹ (mean value of 0.061 microns, at a 5-mm pupil, 0.082 microns scaled in this study), Plainins and Pallikaris (ranking form $-0.22$ to 0.26 microns at a 6-mm pupil), Cheng et al.⁷⁹ (0.132 microns at a 6-mm pupil), He et al.²¹ (0.06 microns at a 6-mm pupil), or Radhakrishnan and Charman²⁸ (mean value of 0.034 microns at a 4.5-mm pupil, 0.053 microns scaled in this study), and it is about one-half of that in other studies.^{17,18,20,22,23,25,26,54,56}

Fig. 6

Corneal aberrations compensation by the lens: spherical aberration (top), horizontal coma, and horizontal/verticle (H/V) astigmatism (bottom).

The lens does not compensate for the horizontal corneal coma, since its value, $C_{3,\text{lens}}^1+C_{3,\mathrm{GRINS}}^1+C_{3,\mathrm{GRINT}}^1$, increases in magnitude with an opposite sign and a higher ratio compared to $C_{3,\text{cornea}}^1$, which decreases (Fig. 6). Those values range from $-0.13$ to $-0.28$ microns while for the corneal horizontal coma declines with age from 0.2 to 0.09 microns. The trend of these values agrees with that reported in different studies, but the magnitude of the values is greater than some of them, even when the studies consider the third- and fifth-order root-mean-square (RMS) values.^{24,28,29,54,80}. Other studies agree with our results when considering mean values.^56,81

Finally, Fig. 6 also shows that the H/V astigmatism of the lens, $C_{2,\text{lens}}^2+C_{2,\mathrm{GRINS}}^2+C_{2,\mathrm{GRINT}}^2$, is almost constant (a value around $-0.07$ microns) with aging. Thus, the lens astigmatism does not compensate for that of the cornea, $C_{2,\text{cornea}}^2$. However, the magnitude of this aberration value compared to the others is very low.

3.3.

Analysis of the GRIN Profile Role

The effect of the GRIN profile was also analyzed by representing its relative contribution (both refraction and light transfer) compared to the joint optical ocular surfaces of the cornea and the lens (Fig. 7). Regarding the spherical aberration, the GRIN profile contribution, $C_{4,\mathrm{GRINS}}^0+C_{4,\mathrm{GRINT}}^0$, is almost constant and cannot compensate for the contribution corresponding to the optical surfaces, $C_{4,\text{cornea}}^0+C_{4,\text{lens}}^0$, which slightly increases with age due to the change in the optical geometry of the corneal and the lenticular surfaces.^2,3 However, the GRIN profile over-compensates for the effect of the horizontal coma, caused by the optical surfaces, and thus, the overall ocular coma. That is, the GRIN contribution, $C_{3,\mathrm{GRINS}}^1+C_{3,\mathrm{GRINT}}^1$, has a different sign and compensates for the horizontal coma that comes from the surfaces, $C_{3,\text{cornea}}^1+C_{3,\text{lens}}^1$, leading to an overall increase with age ($-0.13$ to $-0.28$ microns). This value is slightly higher compared to those reported in experimental studies,^28,29,49,82 although, the variation agrees with others.^57,80,83

Fig. 7

Role of the gradient-index (GRIN) profile in the compensation for the ocular aberrations induced by ocular surfaces: spherical aberration (top), horizontal coma, and horizontal/vertical (H/V) astigmatism (bottom).

In the case of the H/V astigmatism, both contributions, that of the GRIN profile, $C_{2,\mathrm{GRINS}}^2+C_{2,\mathrm{GRINT}}^2$, and that of the optical surfaces, $C_{2,\text{cornea}}^2+C_{2,\text{lens}}^2$, add up to 30 years, and then have opposite signs (Fig. 7). The value corresponding to the optical surfaces increases with age, while that corresponding to the GRIN profile is almost constant, and thus the astigmatism would seem to decrease through life. However, it can be seen that the magnitude of the value is small compared to those of the spherical and horizontal coma, supporting the claim that the astigmatism is not so important in central vision, at least for emmetropic eyes.

4.

Discussion

Several experimental studies have demonstrated the robustness of the design of the human eye to compensate for ocular aberrations between the different elements. Thus, corneal aberrations, particularly spherical aberration and horizontal coma, are partially compensated for by the decentered and tilted lens, as well as decentration of the pupil. Moreover, the change of that compensation between cornea and internal optics has been postulated to be partially dependent on the variation in the lens geometry, in the alignment with respect to the cornea, as well as in the GRIN profile throughout life.^{31,47,54,57,58}

A number of studies have contended that the fourth-order spherical aberration does not vary with age^22,29 while others hold that it does.^28,31,80 This study also demonstrates that the lens cannot fully compensate for the spherical aberration of the cornea. However, it makes the total value retain its value with age. Notably, the change of the GRIN profile with aging is responsible for this effect, since the optical surfaces of the lens by themselves could not do so. The modelling shows that the corneal and the lens surfaces add to the spherical aberration with the age. However, since the lens-surface contribution decreases, this cannot account for the variation of the ocular spherical aberration.⁵⁸ In addition, the fourth-order spherical aberration is not affected by the decentration and tilt between the cornea and lens with respect to the iris, as shown in the Eqs. (3)–(5). This could support the hypothesis that the role of the GRIN lens profile is that of a passive physiological mechanism, having a genetic origin, as has been postulated by other authors.⁵⁴ Moreover, the most relevant term in the compensation is that of the GRIN profile contribution to the refraction, in contrast to the transfer term, which has a small one. This could indicate the minor effect on the spherical aberration of the refractive-index axial variation within the lens. However, when Eqs. (7) and (8) are taken into account, this axial variation is found to be responsible not only for the aberration contribution to the wavefront through the lens, but also for the refraction. Thus, a coupling effect of the geometry of the lens surfaces and the GRIN profile, given by Eq. (7), appears to be determinant in the compensation of the spherical ocular aberration.^62,63 Therefore, since the asphericities of the surfaces drive the contribution, reliable measurements are needed, as has been pointed out previously.⁵⁸

Regarding the corneal horizontal coma, there are studies reporting no significant variation with age.^29,31 The modelling reported in the present work indicates that the corneal horizontal coma does indeed diminish slightly throughout life, with a mean value of 0.17 microns. Nevertheless, it should be taken into account how these values are calculated compared to those from experimental studies. In the latter, the keratometric axis sets the anterior corneal pole in a different position from that of the corneal model. The rotation for correcting this reference point adds changes in the measured values for the corneal horizontal coma, which are not considered in the model eye. Internal optics does help to compensate slightly for the corneal coma, and thus, an increase arises with age. Again, the GRIN profile contribution to the horizontal coma is almost constant and it has an opposite sign compared to that from the optical surfaces. Horizontal coma is affected by the kappa angle, and this angle affects the coupling with the spherical aberration, as reflected by Eq. (5). Thus, it is conceivable that a fine tuning in the compensation between the cornea and lens may occur in individual subjects by varying the angle kappa, as several studies have shown.^31,49 Therefore, an active mechanism could take place in order to balance the spherical aberration with the ocular horizontal coma.

Finally, regarding the H/V astigmatism, although the trend is consistent with that reported in the literature, it should be discussed from a more conservative viewpoint. Moreover, this ocular aberration is of less relevance, since spherical and coma are the most important ocular aberrations limiting image quality in humans.⁴⁷ In this regard, it should even be noted that our model has a rotationally symmetric cornea, which does not take into account the corneal astigmatism. This may overestimate the value of the compensation of astigmatism by the lens. However, if only the contribution of the latter is considered, again the relevant term is the contribution of the GRIN profile to the refraction. In addition, the coupling of spherical aberration and coma through the offset and rotation must be taken into account, as shown by Eq. (5).

The notable weight that the GRIN profile has on the compensation of the aberrations between the cornea and the lens reveals that its role cannot be neglected in that compensation. Thus, optical surfaces themselves cannot explain the failure in the balance for the ocular aberrations, particularly for the horizontal coma. A recent modelling has been made along this line, considering both the cornea and lens as spherical thin lenses,^31,49,58 and it has been concluded that better knowledge of the changes in the lens-surface asphericities and GRIN profile could support a better understanding of that decoupling. Although those studies could have considered the aspherical thin lenses,⁸⁴ for which the horizontal coma is not linear with the shape factor because of coupling with the spherical term [Eq. (5)], this study agrees that the role played by the GRIN profile appears to be fundamental.

In summary, current eye models are based on experimental data regarding the geometrical surface parameters as well as the GRIN profile. Then, optimisation methods are performed based on the aberration data gained experimentally as well. Several degrees of freedom can be set in the eye models to fit the experimental data. However, more reliable measurements regarding the misalignments between the lens and the cornea, the lens-surface geometry and, particularly, the GRIN profile, are needed to get an accurate eye model. The aging GRIN profile equation proposed recently for the human lens⁶⁴ seems robust enough to formulate a better aging-eye model,⁶⁵ including misalignments between ocular elements as well as non-rationally symmetric surfaces. Moreover, the modelling of this study has demonstrated the coupling of the asphericities of the lens surfaces with the GRIN profile, providing the main factor in balancing the ocular aberrations.

5.

Conclusion

This paper presents a modelling study, based on an aging-eye model, of the time course of the third-order spherical aberration, horizontal coma, and H/V astigmatism, by including misalignments between ocular elements and a GRIN lens. The role of a GRIN lens profile is revealed from its main contribution to all the ocular aberrations studied. Moreover, the GRIN profile should be taken into account because the contribution of the lens optical surfaces alone fails to explain the decoupling of the aberrations between the cornea and the lens.