1. Introduction

The cornea and the lens form the refractive optical system of the human eye. One of the particularities of the lens is that it continuously grows throughout the entire lifespan [1]. This continuous growth causes age-related changes in lens shape and internal structure that play a key role in the development of the refractive state in childhood and also contributes to age-related changes in the retinal image quality later in life [2–6].

The effect of lens growth on the optics of the eye is not well understood, in part because of the lack of reliable age-dependent optical models of the lens [7]. Developing eye models that predict the role of the lens in refractive development and the optics of the ageing eye is challenging because the lens is an inhomogeneous optical element. Variations in the protein concentration within the lens produce a three-dimensional variation of the refractive index (refractive index gradient) that contributes to the refractive power and aberrations of the lens [8]. The presence of the gradient increases the complexity of optical modelling and computations.

To avoid the complexity associated with ray-tracing in inhomogeneous media, paraxial models of the crystalline lens generally assume that the lens is a homogeneous biconvex lens that has the same shape as the real lens and an “equivalent” refractive index that is calculated to give the equivalent lens a dioptric power equal to the measured lens power [9]. The use of a back-calculated homogeneous equivalent index is an artifice that helps simulate the age dependence of lens shape and lens power but fails to provide any insight regarding the relation between lens internal structure and lens optics. For instance, there is evidence that lens growth is associated with age-related changes in the refractive index distribution that could contribute to changes in lens power observed during refractive development and in the adult eye [10,11]. In addition, it is unclear if the homogeneous equivalent lens accurately predicts the position of the lens principal panes. Errors in the position of the lens principal planes could lead to clinically significant errors when the homogeneous equivalent lens is implemented in eye models.

Several inhomogeneous optical models of the lens using continuous functions or discrete shells have been proposed to assess the role of the refractive index gradient and its changes with age (see for instance [12–21]). Continuous models most generally use conic or figured conic sections to represent the lens surfaces and polynomial functions to model the spatial variation of the refractive index. The advantage of conics or figured conics is that they require only a small number of parameters that can be adjusted to generate closed contours that enable finite ray-tracing. However, these models typically require simplifying assumptions (e.g., parabolic refractive index profiles [14], thin lens [20], constant ray height [12,21] or parabolic ray path [14,17,21,22]) to produce analytical solutions for the paraxial characteristics in terms of the lens shape and refractive index gradient parameters. Alternatively, they rely on numerical ray-tracing [18,19,24] or they produce complex analytical expressions that require numerical computation [13,23]

We present a solution of the paraxial ray equation in a continuous gradient model of the human lens that produces an analytical expression for the relation between lens shape, refractive index distribution and paraxial optical properties of the lens. For this purpose, we start by solving the paraxial ray equation for inhomogeneous media using two common assumptions that define the refractive index distribution within the lens: 1) the axial iso-indicial curvature profile is linear (Fig. 1(B)) and 2) the axial refractive index profile is a power function (Fig. 1(C)). We find that this solution leads to a paraxial model of the gradient lens that accurately predicts the human lens power and accommodative response.

 

Fig. 1. Coordinate system and notation. The optical axis is taken to be the z-axis with light propagating in the positive z direction and ray height is measured along the y-axis.

Download Full Size | PDF

2. Coordinate system and paraxial lens geometry

The coordinate system and notations are defined in Fig. 1 and Table 1. The lens is assumed to be rotationally symmetric. The paraxial radii of curvature of the anterior and posterior lens surfaces are R₁ and R₂ respectively. The origin of the coordinate system is the intersection of the optical axis and the equatorial plane. We assume that the equatorial plane is the boundary between the anterior and posterior regions of the lens (i.e., the lens center). The anterior lens thickness t₁ is the distance from the anterior surface vertex to the equatorial plane. Similarly, the posterior lens thickness t₂ is the distance from the equatorial plane to the posterior surface vertex. The total axial thickness of the lens is t_L= t₁+ t₂. The lens is surrounded by an aqueous medium with refractive index n_a.

Table 1. Model input parameters. The accommodation dependent-model is calculated for age 28.5 years. The parameters for the relaxed lens (0 D of accommodation) are calculated using the age-dependent model

View Table | View all tables in this article

3. Paraxial representation of the refractive index distribution

3.1 Refractive index distribution

In paraxial approximation, the refractive index distribution can be modeled as a polynomial series in terms of surface height (y), limited to the second degree [12,22,25,26]:

3.2 Axial refractive index profile

As in previous models [13,24,29], we assume a power dependence of the refractive index profile along the optical axis in each half of the lens (Fig. 1(C)):

3.3 Paraxial iso-indicial curvature

We assume that the iso-indicial contours are scaled replicas of the outer lens contour [12,13,18,29]. This assumption produces a linear variation of the axial radius of curvature of the iso-indicial contours (Fig. 1(B)):

4. Paraxial ray height and slope

In paraxial approximation, the height, y(z), of a ray propagating in an inhomogeneous medium with a refractive index distribution given by Eq. (1) is the solution of the following second-order homogeneous linear differential equation [21,25,26]:

Combining Eqs. (2), (4b), (5) and (6) yields:

Since the coefficients of the differential equation do not have any singularity in the domain of interest, the solution of this differential equation can be expressed as a Taylor series. To simplify the expressions and derivation, it is convenient to use the lens equator (z = 0) as the reference coordinate for the Taylor series. This approach produces the following expression of the ray height y(z) and ray slope y'(z) (See derivation in Appendix A):

Equation (8) is an approximation produced by limiting the Taylor series of y(z) to terms of degree b+1 or less. In Appendix A we show that the error produced by the series truncation is negligible and that the solution is also valid when b is a fractional number greater than unity. Additionally, to evaluate the accuracy of Eq. (8), we compared the ray height and ray slope predicted using Eq. (8) with a numerical solution of the differential equation Eq. (7) obtained using MATLAB. The error depends on the values of G_R and of the coefficient b. For typical values expected from human lenses (b≥2, G_R<-2) the relative error of Eq. (8) is less than 0.04% for ray height and less than 0.62% for ray slope (Fig. 2). In summary, we find that the paraxial ray path, y(z) can be closely approximated by a polynomial function of degree b+1.

 

Fig. 2. Relative error between truncated polynomial series and numerical solutions of the paraxial ray equation. For these simulations, y(0) = 1 and y’(0) = -0.1. The error decreases as the value of the power coefficient b increases. The plot shows the error for b = 2.

Download Full Size | PDF

5. Ray transfer matrix of the gradient

In each half of the lens (anterior and posterior), the ray height and ray slope Eq. (8) can be expressed in matrix form:

We will use the condensed notation:

Equation (10) gives the ray vector, Y(z), in the anterior and posterior halves of the lens in terms of the ray-vector at the lens equation, Y(0). For ray-tracing purposes, we seek the value of Y(z) in terms of the ray vector at the anterior boundary of the gradient, Y(-t₁). For this purpose we first express Y(0) in terms of Y(-t₁):

Using Eq. (9) and the formula for the inverse of a 2×2 matrix, we find:

Equation (13) gives the ray-transfer matrix of the entire gradient in terms of the gradient parameters (n_s, n_e, b, G_R1, G_R2), and radii of curvature and thicknesses of the anterior and posterior lens. As expected from the properties of ray-transfer matrices, we have det G₁⁻¹(-t₁) ≅ n_s/n_e and det G₂ (t₂) ≅ n_e/n_s. The determinants are not exactly equal to the ratios of the refractive indices because the solution of the ray-equation (Eq. (8)) is an approximate solution truncated to the term of power b+1. A simplified expression of G_L can be derived by neglecting higher-order terms in g_n, considering that the lens has a shallow gradient (g_n<<1, see Eq. (3)). With this approximation the expression of the matrix G_L of Eq. (13) simplifies to (see Appendix B):

Equation (14b) shows that when the higher-order terms in g_n are omitted, the power of the gradient is approximately equal to the power of a thick lens of radii R₁ and R₂, thickness t_L and equivalent refractive index n_G. The equivalent index is determined entirely by the parameters of the axial refractive index profile (n_e, n_s, b). As long as 1 < b < 1/g_n, the equivalent index will be larger than the peak value of the axial refractive index distribution (n_e). Since n_e-n_s <<1 for the crystalline lens, the second degree term in Eq. (14c) makes a small contribution (less than 0.004). This term corresponds to the contribution of the lens thickness to the equivalent index. The thin lens approximation of the equivalent index is equal to the first two terms of Eq. (14c):

6. Ray-transfer matrix of the lens

The matrix, L, of the whole lens, including refraction at the anterior surface (aqueous- outer lens cortex boundary), transfer through the gradient, and refraction at the posterior surface (outer lens cortex-aqueous boundary) can be expressed as:

If we use the approximation of the gradient matrix Eq. (14), and neglect the contribution of the smaller terms, we find the following approximate expression for the matrix of the lens:

Equation (16) shows that the lens power is approximately equal to the power of a thick lens of radii R₁ and R₂, thickness t_L and equivalent refractive index n_G (see Eq. (14(c)) surrounded by aqueous. According to Eq. (16(a)), the position of the principal planes (H and H’) relative to the anterior and posterior vertices of the lens (V₁ and V₂) are

For comparison the position of the principal planes of the equivalent homogeneous lens with refractive index n_G is:

Comparing Eq. (17) and Eq. (18) shows that there are two sources of error in the predicted position of the principal planes when the lens is represented by the equivalent homogeneous lens. The contribution of the radius of curvature gradient is missing (terms in G_R in Eq. (17)) and the denominator of the surface contribution uses the equivalent index (n_G) instead of the surface index (n_s).

7. Evaluation of the model

7.1 Input parameters

To evaluate the ability of the model to predict lens paraxial optics, we used age-dependent measurements of lens radii of curvature, thickness, equivalent index and power by Atchison et al [30] and the accommodation-dependent data of Martinez-Enriquez et al [31]. For all analyses, we used the average surface and equatorial indices obtained using MRI [11]: n_s= 1.379, n_e= 1.410 calculated at a wavelength of 589 nm. We calculated the power coefficient (b) from the age-dependence of the equivalent index found by Atchison et al [30] using Eq. (14c), with n_s= 1.379 and n_e= 1.410. The calculated coefficient b was found to vary with age from 1.84 at age 15 years to 2.49 at 60 years (Fig. 3). For the ratio of anterior and posterior thickness to total lens thickness we used the values of Rosen et al [32]: t₁ = 0.4 t_L, t₂ =0.6 t_L. The data sets are summarized in Table 1.

 

Fig. 3. Power coefficient (b) calculated using Eq. (14) from the age-dependence of the equivalent index of Atchison et al [30] (black squares) using n_s= 1.410 and n_e = 1.379 [11]. For comparison, the values obtained using the equivalent index data of Dubbelman and van der Heijde [33] (red dots) and Chang et al [34] (green triangles) and from in vivo MRI imaging [11] and in vitro X-ray phase contrast tomography [35] are also shown. The three datasets produce comparable values of the power coefficient in younger eyes, but the age dependence is variable.

Download Full Size | PDF

7.2 Age-dependent lens matrix coefficients, power, and principal plane position

With the values of Table 1, the relative error between the product Eq. (15) and the approximate matrix Eq. (16) is less than 0.02% for A_L, less than 0.71% for B_L, less than 0.17% for C_L and less than 0.04% for D_L in the age range from 15 to 60 years. The error in power is less than 0.04D (Fig. 4). The approximation produces a slight shift of the principal planes towards the lens center compared to the exact model. The shift is less than 23 µm for the object principal plane and less than 11 µm for the image principal plane (Fig. 5). Overall, the error introduced by omitting the higher-degree terms in g_n is therefore negligible for practical purposes.

 

Fig. 4. Age-dependence of the exact and approximate ray-transfer matrix coefficients for the lens, using the data of Table 1. For convenience, the power is shown instead of the C-coefficient (P= - n_a * C_L).

Download Full Size | PDF

 

Fig. 5. Age-dependence of the principal plane locations for the exact gradient lens, approximate gradient lens and homogeneous equivalent lens, using the data of Table 1. The plots show the distance between the anterior vertex (V₁) and the object principal point (H) and the posterior vertex (V₂) and the image principal point (H’). Left: Absolute distance. Right: Distances relative to lens thickness. The relative position of the principal planes of the gradient lens is approximately constant with age. There is a significant difference in the position of the principal planes between the gradient model and the equivalent homogeneous model.

Download Full Size | PDF

The model predicts that the principal planes of the gradient lens shift towards the lens center with age (Fig. 5). However, when expressed relative to the lens thickness, the position of the principal planes remains approximately constant with age (Fig. 5). The distance between the principal planes of the gradient lens is less than 0.10 mm for the exact lens and less than 0.07 mm for the approximate lens. The lens can therefore be approximated as a thin lens located near the midpoint of the lens. The optimal location of the thin lens will depend on the conjugation ratio (i.e., object distance) [36].

The object and image principal planes of the homogeneous equivalent model are located closer to the posterior lens surface, when compared to their position in the gradient model (Fig. 5). This prediction error in the position of the principal planes decreases with age from 0.17 mm to 0.10 mm for the object principal plane and from 0.32 mm to 0.23 mm for the image principal plane. The separation of the principal planes is also larger in the equivalent homogeneous model (0.22 mm versus 0.10 mm).

We find that the model closely predicts the measured lens power, with an absolute error of less than 0.5D. The error progressively decreases from +0.13D at age 18 to -0.45D at age 60 years. This analysis shows that Eq. (14c) accurately predicts the relation between the gradient parameters and the equivalent index when the power coefficient b is calculated using published values of the refractive indices at the lens surface and equator. The power coefficient that is produced using the data of Table 1 is lower than the value derived using MRI [11] or X-ray tomography [35]. This difference may be due in part to limitations in the accuracy and precision of the measurement of the equivalent index and refractive index distribution (see discussion). The model also predicts that the gradient power represents a fraction of the total lens power that progressively decreases with age from 61% at age 18 to 53% at age 68 years (ratio P_G/P_L expressed in %), in good agreement with prior studies [37] (Fig. 6).

 

Fig. 6. Contribution of the lens gradient to total lens power and accommodation. Left: The model predicts that the gradient contribution to lens power progressively decreases with age from 61% at age 18 to 53% at age 68. The slope of the decrease in gradient power is larger than the slope of the increase in surface power, resulting in a decrease in total lens power with age. Right: The model predicts that for a 28.5 year subject, the gradient is responsible for 60% (= 0.779 / 1.297) of the total lens power, a fraction that remains approximately constant with accommodation.

Download Full Size | PDF

7.3 Simulation of accommodation

The change in lens power with accommodation was simulated using the average accommodation-dependent rate of change in lens thickness (0.069 mm/D), anterior radius (-0.60 mm/D) and posterior radius (+0.22mm/D) found by Martinez-Enriquez et al [31], where the change is expressed in terms of accommodative response. To simulate the lens at 0D of accommodation, we used the radii, thickness and equivalent index obtained from the age-dependent data of Table 1 at age 28.5 years (average age in the study of Martinez-Enriquez et al [31]). We assumed that the axial refractive index profile is unchanged with accommodation [11,38]. As expected, the simulation predicts an approximately linear relation between change in lens power and accommodative response, with a slope of 1.30D/D (Fig. 6). In other words, the model predicts that the accommodative response is 0.77D per diopter of lens power change. This value is in good agreement with the value found using various accommodation-dependent eye models (for instance 0.79D/D with the Navarro eye model [39]). The relative contribution of the gradient to the change in lens power is 60% (0.78D/D). This value is also in close agreement with results obtained in non-human primate lenses [37].

8. Discussion

We present a paraxial model of the gradient-index lens of the human eye that provides a direct relation between gradient parameters and paraxial characteristics (power and principal plane position). The model relies on two assumptions that are commonly used to represent the lens gradient: the iso-indicial radius of curvature varies linearly with depth and the axial refractive index distribution is a power function with single term of power b. With these assumptions, we demonstrate that the paraxial ray height can be closely approximated by the sum of the initial ray height and a polynomial function limited to three terms in z, z^b and z^b+1. The solution is also valid when the power coefficient b is a fractional number greater than unity. We find that the lens power can be closely approximated by the power of a thick homogeneous lens with an equivalent index that is determined solely by the parameters of the axial refractive index profile. However, this equivalent homogeneous lens does not accurately predict the position of the principal planes of the gradient lens. When combined with published biometric data, the gradient model accurately predicts the changes in lens power with age and accommodation.

The model predicts that an increase in the value of the power coefficient of the axial refractive index profile (i.e., a flatter profile) produces a decrease in equivalent index (Eq. (14c)). This prediction is consistent with the findings that there is progressive formation of a central plateau in the axial refractive index profile with age [3,11] (Fig. 1(C)) and at the same time a decrease in the equivalent refractive index with age [30,33,34]. The model also predicts that the formation of a plateau in the axial refractive index profile (i.e., increase in the coefficient b) with age causes a decrease in the power of the gradient (Eq. (16b)). The decrease in gradient power translates into a decrease in total lens power despite the progressive steepening in lens curvatures with age, thus providing an explanation for the so-called “lens paradox” [40].

We find that the equivalent index is dependent entirely on the parameters of the axial refractive index profile with the refractive index distribution defined in Section 3. The thin lens approximation of the equivalent index (Eq. (15b)) is consistent with formulas derived by others using a thin lens approximation of the gradient (for instance Eq. (14) of [13]). On the other hand, Sheil and Goncharov [20] derived an expression of the equivalent index (their Eq. (9)) which depends not only on the axial gradient profile, but also on the lens thickness and radius of curvature. To derive this expression, they equate the power of the thin lens approximation of the gradient lens (no thickness term) to the power of the thick homogeneous equivalent lens (thickness term included). If instead the equivalent index is calculated by equating the power of the thin lens approximations of both the gradient and homogeneous equivalent lens, Sheil and Goncharov’s [20] formula becomes the same as Eq. (14d) when the axial variation of the iso-indicial radius of curvature is linear. Overall, these comparisons provide further support for the validity of our analytical solution. Differences between models are due primarily to differences in the approach used to define or calculate model output parameters.

In the numerical application of our model (Section 7), we calculated the power coefficient (b) from Eq. (14c) using published values of the surface and equatorial indices [11] and of the equivalent index [30]. The power coefficient obtained with this approach and its age-dependence are strongly dependent on the measurements that are used (Fig. 3). For instance, the three datasets shown in Fig. 3 produce similar values of the power coefficient in young adults (b ranges from 1.9 to 2.3 at age 23), but very different changes with age. In addition, the power coefficient calculated for young adults is significantly lower than the measurements obtained using in vivo MRI and in vitro X-ray tomography (b = 4.9 at age 23 in both cases). This variability reflects the high uncertainty of the experimental measurements of the gradient parameters and equivalent index. The impact of this uncertainty on the calculated power coefficient is illustrated in Fig. 7. The analysis shows that the value of the power coefficient produced using Eq. (14c) is highly dependent on the value of the equivalent index and of the equatorial refractive index (n_e). The difference becomes more pronounced for lower values of the equivalent index (i.e., older subjects since the equivalent index decreases with age). This finding is consistent with the data of Fig. 3, which shows that the difference in the coefficient b predicted from the equivalent refractive index measured in three different studies becomes more pronounced in older eyes.

 

Fig. 7. Value of the power coefficient in terms of the equivalent index, surface index and equatorial index. Left: The coefficient b is less sensitive to variations in the surface index than variations in the equatorial index. The model predicts that the power coefficient will increase with age as the equivalent index decreases with age. The variation of the power coefficient is more pronounced for lower values of the equivalent index, consistent with the plot of Fig. 3. Right: The power coefficient b is strongly dependent on the difference n_G-n_e (Eq. (19)).

Download Full Size | PDF

Figure 7 also shows that uncertainties in the surface index have a less pronounced effect on the power coefficient calculated from the equivalent index and axial gradient parameters than uncertainties on the equatorial or equivalent index. This behavior can also be deduced from the equation relating the power coefficient to the refractive indices, derived from Eq. (14c):

Equation (19) shows that the power coefficient is dependent primarily on the differences n_G-n_e and n_e-n_s. These findings imply that inverse algorithms designed to reconstruct the refractive index gradient from optical measurements [41,42] can apply looser tolerances on the value of the surface index. Overall, even with the uncertainty regarding the absolute values of the parameters, the model provides general insight regarding the relation between the gradient parameters and the lens paraxial properties, and their age-dependence.

Clearly, the relation between gradient parameters and paraxial properties depends on the assumptions that are made regarding the refractive index distribution. In our paraxial model, the gradient is defined entirely by two independent functions: the axial variation of the refractive index and the axial variation of the curvature of the iso-indicial surfaces. The power function used to represent the axial index and the linear variation of the axial radius of curvature are commonly used to represent the lens gradient [7,13,18,24,29]. In our model, the axial gradient in iso-indicial curvatures is determined by the thickness and surface curvature in each half of the lens (G = ε R_s/t). With the values of Table 1, the axial gradient G progressively decreases (in absolute value) with age, from -8.3 to -5.3 mm/mm in the anterior lens and from -3.3 to -2.5 mm/mm in the posterior lens from age 15 to age 60. It is straight forward to apply our general approach to other models [20] of the axial dependence of the iso-indicial curvature, by replacing R(z) in the expression of n₁(z)

Previous efforts to derive a relation between gradient parameters and lens power produced complex solutions or alternatively relied on approximations to produce simpler expressions, including the thin lens approximation [7] or the assumption that the axial refractive index profiles or ray paths are parabolic or low-degree polynomials [21–23]. One of the sources of complexity is that prior models generally chose to place the origin of the coordinate system at the lens surface. The axial variation of the refractive index is then represented by the polynomial (z-t)^b, which includes b monomials of power ranging from 1 to b. This more generalized power series produces a recurrence relation for the coefficients [14,22,43] of the polynomial solution of ray height, which leads to a more complex polynomial series (see for instance the matrix solutions of Diaz [22]). We circumvented this difficulty by using the lens center (defined by the intersection of the lens axis and the equatorial plane) as the origin of the coordinate system, deriving the matrix defining the ray vector at any position in the lens in terms of the ray vector at the lens center, and then using the inverse matrix (Eq. (12) to find the relation with the ray vector at the anterior lens boundary. By choosing the lens center as the origin of the coordinate system we are able to express the axial variation of the refractive index using a single monomial term of power b, which leads to a simple closed-form analytical expression for the ray path and lens power. Bahrami and Goncharov [13] used a similar approach to derive an analytical solution of the differential equation which expresses the ray height as a sum of hypergeometric functions (their Eq. (22)). With our notation, this hypergeometric solution is a sum of monomials of power m*b and m*b+1, consistent with our findings. However, the two solutions are otherwise difficult to compare because the hypergeometric solution leads to an expression of ray height, ray slope and lens power that involves multiple distinct hypergeometric series.

The finding that the gradient contributes approximately 60% to the total lens power, a ratio that remains approximately constant during accommodation, is in good agreement with our prior experimental results obtained on non-human primate lenses [37]. The gradient contribution calculated using this approach is higher than the value found by Navarro and Lopez-Gil [7], but this difference simply reflects a difference in the definition of the gradient contribution. Navarro and Lopez-Gil use the power of a homogeneous lens with axial gradient refractive index gradient but no curvature gradient (G_R= 0) as the baseline to quantify the contribution of the gradient (their Eq. (4)). We quantify the gradient contribution directly by calculating the ratio of gradient power (Eq. (14b)) to total lens power (Eq. (16b)), corresponding to their Eq. (3) (in thin lens approximation) and to the method of Smith and Atchison [21].

Garner and Smith [44] found that the gradient is responsible for approximately 50% of the change in lens power with accommodation using a gradient model with nested elliptical iso-indicial contours and with a parabolic axial index profile (b = 2). The difference relative to our prediction of 60% is due primarily to differences in the refractive indices. If we enter in our model the index values used by Garner and Smith (n_s= 1.3859, n_e= 1.406, n_G= 1.4277), we find a power coefficient of b = 1.9 (from Eq. (14c)) and a contribution of the gradient of 46%, in good agreement with the estimate of Garner and Smith. As discussed earlier, nested shell models such as the one used by Garner and Smith are characterized by a linear variation of the axial iso-indicial radius of curvature (our Eq. (5)a). In addition, the indices used by Garner and Smith correspond to b = 1.9 in our model, close to the parabolic profile (b = 2) assumed by Garner and Smith. Finally, the lens shape data (radii and thickness) used in our numerical analysis is close to the values measured by Garner and Smith. The paraxial input parameters of the two models are therefore similar. The agreement between the two models when the input parameters are comparable confirms that our model is a paraxial equivalent of nested shell models using conic or figured conic contours.

In our simulations of the accommodative response, we assumed that the power coefficient in the axial profile of the refractive index remains constant with accommodation, following the findings of in vivo MRI studies [11]. The predicted change in lens power is in good agreement with expected values, which provides support for the assumption that the axial refractive index profile does not change significantly during accommodation. With this assumption, the change in lens power with accommodation is due entirely to changes in the curvature of the iso-indicial surfaces. In other words, the steepening and increase in lens thickness with accommodation produces a change in the slope of the iso-indicial curvatures while the axial refractive index gradient remains unchanged. This finding corresponds to Gullstrand’s intracapsular theory of accommodation, which postulates that changes in the internal structure of the lens during accommodation contribute to the changes in lens power. Our results are also in good agreement with the simulations of Navarro and Lopez-Gil [7]. In our model, the slope of the curvature gradient varies linearly from -7.3 to -4.4 mm/mm in the anterior lens (slope of +0.47 mm/mm/D) and from -3.0 to -2.2 mm/mm in the posterior lens (slope of +0.14 mm/mm/D).

Our analysis shows that the equivalent homogeneous thick lens model does not accurately simulate the position of the principal planes relative to the lens vertices and their age dependence. We evaluated the impact of this error on prediction errors of eye models. The lens model at age 20 years was implemented in an eye model with a 42D thin lens as the cornea and an anterior chamber depth of 3.2 mm (distance from cornea to anterior lens vertex). Using the exact gradient lens matrix, we calculated the vitreous depth (16.26 mm) that produces an emmetropic eye. The gradient lens was then replaced with the equivalent homogeneous lens. The model eye with the homogeneous lens had a refractive error of +0.73D. This example shows that the reliance on the homogeneous equivalent lens model in eye models can introduce clinically significant prediction errors.

The model is valid in the paraxial approximation, i.e, for small incident field angles and ray heights. Optical measurements of in vitro human and monkey lenses using laser ray-tracing [5, 45] suggest that the domain of validity of the paraxial approximation extends approximately to the central 3 mm optical zone (1.5 mm ray height) and for field angles that are within approximately +/-5° to +/-10° depending on the age. Beyond this domain, a more advanced model that takes into account lens asphericity must be implemented to account for the contribution of lens aberrations and the effect of lens tilt and decentration and elliptical pupil when the lens is implemented in an eye model. The model also assumes a single wavelength. The lens refractive index values used in our analysis were calculated for a wavelength of 589 nm [11]. To predict chromatic aberrations, a wavelength-dependent model must be developed that accounts for the dispersion of the lens. The dispersion of the equivalent index can be calculated from the dispersion of the lens refractive indices using Eq. (14c).

In summary we present a closed form analytical solution of the differential equation defining the paraxial ray path in the lens of the eye. From this solution we derive a paraxial model of the lens with gradient that closely predicts the changes in lens power with age and accommodation. We find that a decrease in equivalent index with age is associated with a flattening of the gradient profile and that the gradient contributes significantly to lens power and its changes with age and accommodation. The model predictions are consistent with the Gullstrand intracapsular theory of accommodation: Changes in iso-indicial radius of curvature gradient with accommodation contribute significantly to changes in lens power. In addition, the axial variation in iso-indicial surface curvature produces a shift in the principal plane position that is not captured by the homogeneous equivalent model. Consequently, the use of homogenous equivalent lens introduces a clinically significant source of error in eye models.

9. Appendix A: Solution of the paraxial ray equation

The differential equation of the ray Eq. (6) is solved using a Taylor series. If we use the lens equator (z = 0) as the reference point for the Taylor series, the ray height can be expressed as:

(A1)$$y(z) = \sum\limits_{k = 0}^\infty {{a_k} \times {z^k}} $$

where

(A2)$${a_k} = \frac{{{y^{(k)}}(0 )}}{{k!}}$$

where y^(k)(0) is the k-th derivative of the ray at the equator. In particular, a₀= y(0) is the ray height at the equator and a₁= y’(0) is the ray slope at the equator. For general expressions of the functions n₀(z) and n₁(z) the differential equation can be solved by calculating the expression of the higher-order derivative y^(k)(z) by successive derivation of the differential equation (Eq. (6)). This method provides a recurrence equation for y^(k)(0) in terms of the lower order derivatives.

In our case, since the functions n₀(z) and n₁(z) are polynomials, the higher order derivatives y(^k)(0) can be found directly by substituting y(z), y'(z) and y''(z) by their Taylor series in the differential equation (Eq. (6)). The substitution, and shifting the indices in the sums to produce a sum of terms in z^k-2, produces the following recurrence relation for the coefficients a_k:

(A3)$$\begin{aligned} \sum\limits_{k = 2}^\infty {k \times ({k - 1} )\times {a_k} \times {z^{k - 2}}} - \sum\limits_{k = b}^\infty {\left( {k - b + \frac{1}{{{G_R}}}} \right) \times b \times {g_n} \times \frac{1}{{{t^b}}} \times {a_{k - b}} \times {\varepsilon ^b} \times {z^{k - 2}}} + \ldots .\\ \textrm{ }\ldots - \textrm{ }\sum\limits_{k = 2b}^\infty {\left( {k - 2b + \frac{1}{{{G_R}}}} \right) \times b \times {g_n}^2 \times \frac{1}{{{t^{2b}}}} \times {a_{k - 2b}} \times {z^{k - 2}}} = 0 \end{aligned}$$

We find that all coefficients a_k from k = 2 to k = b-1 are equal to zero. The only non-zero coefficients are the coefficients of degree m×b and m×(b+1), where m is an integer (m≥1):

(A4a)$${a_{m \times b}} = \frac{1}{{m \times ({m \times b - 1} )}} \times \frac{{{g_n}}}{{{t^b}}} \times {\varepsilon ^b} \times \left[ {\left( {({m - 1} )\times b + \frac{1}{{{G_R}}}} \right) \times {a_{({m - 1} )\times b}} + \left( {({m - 2} )\times b + \frac{1}{{{G_R}}}} \right) \times \frac{{{g_n}}}{{{t^b}}} \times {\varepsilon^b} \times {a_{({m - 2} )\times b}}} \right]$$

(A4b)$${a_{m \times b + 1}} = \frac{1}{{m \times ({m \times b + 1} )}} \times \frac{{{g_n}}}{{{t^b}}} \times {\varepsilon ^b} \times \left[ {\left( {({m - 1} )\times b + 1 + \frac{1}{{{G_R}}}} \right) \times {a_{({m - 1} )\times b + 1}} + \left( {({m - 2} )\times b + 1 + \frac{1}{{{G_R}}}} \right) \times \frac{{{g_n}}}{{{t^b}}} \times {\varepsilon^b} \times {a_{({m - 2} )\times b + 1}}} \right]$$

In particular, remembering that a₀= y(0) and a₁= y’(0):

(A5)$${a_b} = \frac{1}{{b - 1}} \times \frac{{{g_n}}}{{{G_R}}} \times \frac{1}{{{t^b}}} \times {\varepsilon ^b} \times y(0)\qquad{a_{b + 1}}\textrm{ = }\frac{1}{{b + 1}} \times \frac{{{g_n}}}{{{G_R}}} \times \frac{1}{{{t^b}}} \times ({1 + {G_R}} )\times {\varepsilon ^b} \times y^{\prime}(0 )$$

The equation for the ray path (Eq. (A1)) then becomes:

(A6)$$y(z) = y(0 )+ y^{\prime}(0 )\times z + {a_b} \times {z^b} + {a_{b + 1}} \times {z^{b + 1}} + {a_{2b}} \times {z^{2b}} + {a_{2b + 1}} \times {z^{2b + 1}}\textrm{ } + \ldots + {a_{mb}} \times {z^{mb}} + {a_{mb + 1}} \times {z^{mb + 1}} + \ldots $$

This Taylor series converges very rapidly (see Tables 2, 3, 4). For practical purposes, the series can be truncated to the first four terms (up to degree b+1), which gives, from Eq. (A5):

(A7)$$y(z) = \left[ {1 + \frac{1}{{b - 1}} \times \frac{{{g_n}}}{{{G_R}}} \times {\varepsilon^b} \times \frac{{{z^b}}}{{{t^b}}}} \right] \times y(0) + \left[ {1 + \frac{1}{{b + 1}} \times \frac{{{g_n}}}{{{G_R}}} \times ({1 + {G_R}} )\times {\varepsilon^b} \times \frac{{{z^b}}}{{{t^b}}}} \right] \times y^{\prime}(0 )\times z$$

Or:

(A8)$$y(z) = \left[ {1 + \frac{1}{{b - 1}} \times \frac{{{g_n}}}{{{G_R}}} \times \frac{{{{|z |}^b}}}{{{t^b}}}} \right] \times y(0) + \left[ {1 + \frac{1}{{b + 1}} \times \frac{{{g_n}}}{{{G_R}}} \times ({1 + {G_R}} )\times \frac{{{{|z |}^b}}}{{{t^b}}}} \right] \times y^{\prime}(0 )\times z$$

Table 2. Value of the polynomial coefficients when input ray has height y(0)= 1 mm and slope y’(0) = -0.1, calculated for a 20 year old lens (see Table 1): R_s=-6mm, t = 2.2mm n_s= 1.379, n_e= 1.410

View Table | View all tables in this article

Table 3. Contribution of the first 8 polynomial terms to ray height when input ray has height y(0)= 1 mm and slope y’(0) = -0.1, calculated for a 20 year old lens (see Table 1): R_s=-6mm, t = 2.2mm n_s= 1.379, n_e= 1.410

View Table | View all tables in this article

Table 4. Contribution of the first 8 polynomial terms to ray slope when input ray has height y(0)= 1 mm and slope y’(0) = -0.1.

View Table | View all tables in this article

The derivation leading to this equation relies on the assumption that b is a positive integer, and that b≥2. However, we find that Eq. (A8) is also valid when b is a fractional number greater than unity, as shown in Fig. 8.

 

Fig. 8. Solution of the ray equation for a fractional value of the power coefficient b. The plots show the axial refractive index profile (left) and the ray slope (right) for b = 1.8. The ray slope was obtained using a numerical solution of the ray differential equation and from the derivative of Eq. (A8) for y(0) = 1 and y’(0)= -0.1. The numerical solution was obtained using exactly the same approach as for the analysis of Fig. 2.

Download Full Size | PDF

10. Appendix B: Approximation of the gradient matrix

We use the ABCD notation for the matrix G_L:

(B1)$${{\mathbf G}_{\mathbf L}} = \left[ {\begin{array}{{cc}} {{A_G}}&{{B_G}}\\ {{C_G} ={-} \frac{{{P_G}}}{{{n_s}}}}&{{D_G}} \end{array}} \right]$$

Exact values of the coefficients of the matrix G_L can be found by multiplying the matrices G₂(t₂) and G₁⁻¹(-t₁) (Eq. (13)). We can obtain a simplified expression by considering that the lens has a shallow gradient (g_n<<1, see Eq. (3)). We can therefore neglect second-degree or higher-degree terms in g_n and rely on the following approximation to simplify the expressions:

(B2)$$\frac{{{n_e}}}{{{n_s}}} = \frac{1}{{1 - {g_n}}} \cong 1 + {g_n}$$

With these approximations, and using the relation t × G_R = R_s, we obtain the following expression for the matrix coefficients and gradient power, where t_L= t₁+ t₂ is the lens thickness:

(B3a)$${A_G} = 1 - \frac{b}{{b - 1}} \times {g_n} \times \frac{{{t_L}}}{{{R_1}}} + \frac{1}{{b - 1}} \times {g_n} \times \left( {\frac{1}{{{G_{R1}}}} + \frac{1}{{{G_{R2}}}}} \right)$$

(B3b)$${B_G} = \left( {1 - \frac{b}{{b + 1}} \times {g_n}} \right) \times {t_L} + \frac{1}{{b - 1}} \times {g_n} \times \left( {\frac{{{t_1}}}{{{G_{R2}}}} + \frac{{{t_2}}}{{{G_{R1}}}}} \right) + \frac{1}{{b + 1}} \times {g_n} \times \left( {\frac{{{t_1}}}{{{G_{R1}}}} + \frac{{{t_2}}}{{{G_{R2}}}}} \right)$$

(B3c)$${P_G} = {n_s} \times \frac{b}{{b - 1}} \times {g_n} \times \left( {\frac{1}{{{R_1}}} - \frac{1}{{{R_2}}} + \frac{1}{{{R_1}}} \times \frac{1}{{{R_2}}} \times {g_n} \times {t_L}} \right)$$

(B3d)$${D_G} = 1 + \frac{b}{{b - 1}} \times {g_n} \times \frac{{{t_L}}}{{{R_2}}} + \frac{1}{{b - 1}} \times {g_n} \times \left( {\frac{1}{{{G_{R1}}}} - \frac{1}{{{G_{R2}}}}} \right)$$

The terms including the parameters G_R in the coefficient B_G are very small relative to the term in t_L. Neglecting these terms produces a relative error of less than 0.7% in the value of B_G. The expression of the coefficient B_G can therefore be further simplified:

(B3e)$${B_G} = \left( {1 - \frac{b}{{b + 1}} \times {g_n}} \right) \times {t_L}$$

The gradient power (Eq. (B3c)) can be written (using the same approximations as above):

(B4a)$${P_G} = ({{n_G} - {n_s}} )\times \left( {\frac{1}{{{R_1}}} - \frac{1}{{{R_2}}}} \right) + \frac{{{t_L}}}{{{n_G}}} \times \frac{{{{({{n_G} - {n_s}} )}^2}}}{{{R_1} \times {R_2}}} - {n_s} \times \frac{b}{{{{({b - 1} )}^2}}} \times {g_n}^2 \times \frac{{{t_L}}}{{{R_1} \times {R_2}}}$$

where:

(B4b)$${n_G} = {n_s} + \frac{b}{{b - 1}} \times {n_s} \times {g_n} = {n_s} + \frac{b}{{b - 1}} \times \frac{{{n_s}}}{{{n_e}}} \times ({{n_e} - {n_s}} )$$

or in terms of n_s, n_e, and b:

(B4c)$${n_G} = {n_e} + \frac{{{n_e} - {n_s}}}{{b - 1}} - \frac{b}{{b - 1}} \times \frac{{{{({{n_e} - {n_s}} )}^2}}}{{{n_e}}}$$

Eq. (B4a) shows that the power of the gradient is approximately equal to the power of a thick lens of radii R₁ and R₂, thickness t_L and refractive index n_G. With the values of Table 1, the contribution of the third term in Eq. (B4a) (term in g_n²) is found to be less than 0.07D from age 15 to 60. With these approximations, we obtain the following matrix for the gradient:

(B5a)$${{\mathbf G}_{{\mathbf{L - APPROX}}}} = \left[ {\begin{array}{{cc}} {1 - \frac{{{g_n}}}{{b - 1}} \times \left( {b \times \frac{{{t_L}}}{{{R_1}}} + \frac{1}{{{G_{R1}}}} - \frac{1}{{{G_{R2}}}}} \right)}&{\left( {1 - \frac{b}{{b + 1}} \times {g_n}} \right) \times {t_L}}\\ { - \frac{1}{{{n_s}}} \times {P_{G - APPROX}}}&{1 + \frac{{{g_n}}}{{b - 1}} \times \left( {b \times \frac{{{t_L}}}{{{R_2}}} + \frac{1}{{{G_{R1}}}} - \frac{1}{{{G_{R2}}}}} \right)} \end{array}} \right]$$

where

(B5b)$${P_{G - APPROX}} = ({{n_G} - {n_s}} )\times \left( {\frac{1}{{{R_1}}} - \frac{1}{{{R_2}}}} \right) + \frac{{{t_L}}}{{{n_G}}} \times \frac{{{{({{n_G} - {n_s}} )}^2}}}{{{R_1} \times {R_2}}}$$

With the values of Table 1, the relative error between the exact product and the approximate matrix is less than 0.03% for A_G, less than 0.7% for B_G, less than 0.5% for C_G and less than 0.02% for D_G in the age range from 15 to 60.

Funding

National Eye Institute (P30EY14801, R01EY014225, R01EY021834); Research to Prevent Blindness (GR004596); Beauty of Sight Foundation; Dr. Harry W. Flynn, Jr.; Dr. Raksha Urs and Mr. Aaron Furtado; Drs. KR Olsen and ME Hildebrandt.

Disclosures

There are no financial interests to disclose.

Data availability

No data were generated or analyzed in the presented research.

References

1. R. C. Augusteyn, “Growth of the lens: in vitro observations,” Clin. Exp. Optom. 91(3), 226–239 (2008). [CrossRef]  

2. P. Artal and A. Guirao, “Contributions of the cornea and the lens to the aberrations of the human eye,” Opt. Lett. 23(21), 1713–1715 (1998). [CrossRef]  

3. R. C. Augusteyn, “On the growth and internal structure of the human lens,” Exp. Eye Res. 90(6), 643–654 (2010). [CrossRef]  

4. R. Iribarren, “Crystalline lens and refractive development,” Prog. Retinal Eye Res. 47, 86–106 (2015). [CrossRef]  

5. B. Maceo Heilman, A. Mohamed, M. Ruggeri, S. Williams, A. Ho, J. M. Parel, and F. Manns, “Age-dependence of the peripheral defocus of the isolated human crystalline lens,” Invest. Ophthalmol. Visual Sci. 62(3), 15 (2021). [CrossRef]  

6. J. Rozema, S. Dankert, R. Iribarren, C. Lanca, and S-M. Saw, “Axial growth and lens power loss at myopia onset in Singaporean children,” Invest. Ophthalmol. Visual Sci. 60(8), 3091–3099 (2019). [CrossRef]  

7. R. Navarro and N. Lopez-Gil, “Impact of internal curvature gradient on the power and accommodation of the crystalline lens,” Optica 4(3), 334–340 (2017). [CrossRef]  

8. B. K. Pierscionek and J. W. Regini, “The gradient index lens of the eye: An optobiological synchrony,” Prog. Retinal Eye Res. 31(4), 332–349 (2012). [CrossRef]  

9. D. O. Mutti, K. Zadnik, and A. J. Adams, “The equivalent refractive index of the crystalline lens in childhood,” Vision Res. 35(11), 1565–1573 (1995). [CrossRef]  

10. R. C. Augusteyn, C. E. Jones, and J. M. Pope, “Age-related development of a refractive index plateau in the human lens: evidence for a distinct nucleus,” Clin. Exp. Optom. 91(3), 296–301 (2008). [CrossRef]  

11. S. Kasthurirangan, E. L. Markwell, D. A. Atchison, and J. M. Pope, “In vivo study of changes in refractive index distribution in the human crystalline lens with age and accommodation,” Invest. Ophthalmol. Visual Sci. 49(6), 2531–2540 (2008). [CrossRef]  

12. D. A. Atchison and G. Smith, “Continuous gradient index and shell models of the human lens,” Vision Res. 35(18), 2529–2538 (1995). [CrossRef]  

13. M. Bahrami and A. V. Goncharov, “Geometry-invariant gradient refractive index lens: analytical ray tracing,” J. Biomed. Opt. 17(5), 055001 (2012). [CrossRef]  

14. J. W. Blaker, “Toward an adaptive model of the human eye,” J. Opt. Soc. Am. 70(2), 220–223 (1980). [CrossRef]  

15. C. E. Campbell, “Nested shell optical model of the lens of the human eye,” J. Opt. Soc. Am. A 27(11), 2432–2441 (2010). [CrossRef]  

16. S. Giovanzana, T. Evans, and B. Pierscionek, “Lens internal curvature effects on age-related eye model and lens paradox,” Biomed. Opt. Express 8(11), 4827–4837 (2017). [CrossRef]  

17. A. V. Goncharov and C. Dainty, “Wide-field schematic eye models with gradient-index lens,” J. Opt. Soc. Am. A 24(8), 2157–2174 (2007). [CrossRef]  

18. R. Navarro, F. Palos, and L. González, “Adaptive model of the gradient index of the human lens. I. Formulation and model of aging ex vivo lenses,” J. Opt. Soc. Am. A 24(8), 2175–2185 (2007). [CrossRef]  

19. R. Navarro, F. Palos, and L. M. González, “Adaptive model of the gradient index of the human lens. II. Optics of the accommodating aging lens,” J. Opt. Soc. Am. A 24(9), 2911–2920 (2007). [CrossRef]  

20. C. J. Sheil and A. V. Goncharov, “Accommodating volume-constant age-dependent optical (AVOCADO) model of the crystalline GRIN lens,” Biomed. Opt. Express. 7(5), 1985–1999 (2016). [CrossRef]  

21. G. Smith and D. A. Atchison, “Equivalent power of the crystalline lens of the human eye: comparison of methods of calculation,” J. Opt. Soc. Am. A 14(10), 2537–2546 (1997). [CrossRef]  

22. J. A. Diaz, “ABCD matrix of the human lens gradient-index profile: applicability of the calculation methods,” Appl. Opt. 47(2), 195–205 (2008). [CrossRef]  

23. M. V. Perez, C. Bao, M. T. Flores-Arias, M. A. Rama, and C. Gomez-Reimo, “Gradient parameter and axial and field rays in the gradient-index crystalline lens model,” J. Opt. A: Pure Appl. Opt. 5(5), S293–S297 (2003). [CrossRef]  

24. G. Smith, D. A. Atchison, and B. K. Pierscionek, “Modeling the power of the aging human eye,” J. Opt. Soc. Am. A 9(12), 2111–2117 (1992). [CrossRef]  

25. R. K. Luneburg, Mathematical Theory of Optics (University of California Press, 1964); pp 243–249.

26. P. J. Sands, “Inhomogeneous lenses III. Paraxial Optics,” J. Opt. Soc. Am. 61(7), 879–895 (1971). [CrossRef]  

27. H. A. Buchdahl, Optical Aberrations Coefficients (Dover Publications, 1968); pp 305–307.

28. S. Doric, “Paraxial ray trace for rotationally symmetric homogeneous and inhomogeneous media,,” J. Opt. Soc. Am. A 1(8), 818–821 (1984). [CrossRef]  

29. A. de Castro, D. Siedlecki, D. Borja, S. Uhlhorn, J. M. Parel, F. Manns, and S. Marcos, “Age-dependent variation of the Gradient Index profile in human crystalline lenses,” J. Mod. Opt. 58(19-20), 1781–1787 (2011). [CrossRef]  

30. D. A. Atchison, E. L. Markwell, S. Kasthurirangan, J. M. Pope, G. Smith, and P. G. Swann, “Age-related changes in optical and biometric characteristics of emmetropic eyes,” J. Vis. 8(4), 29 (2008). [CrossRef]  

31. E. Martinez-Enriquez, P. Pérez-Merino, M. Velasco-Ocana, and S. Marcos, “OCT-based full crystalline lens shape change during accommodation in vivo,” Biomed. Opt. Express 8(2), 918–933 (2017). [CrossRef]  

32. A. M. Rosen, D. B. Denham, V. Fernandez, D. Borja, A. Ho, F. Manns, J. M. Parel, and R. C. Augusteyn, “In vitro dimensions and curvatures of human lenses,” Vision Res. 46(6-7), 1002–1009 (2006). [CrossRef]  

33. M. Dubbelman and G. L. Van der Heijde, “The shape of the aging human lens: curvature, equivalent refractive index and the lens paradox,” Vision Res. 41(14), 1867–1877 (2001). [CrossRef]  

34. Y. C. Chang, G. M. Mesquita, S. Williams, G. Gregori, F. Cabot, A. Ho, M. Ruggeri, S. H. Yoo, J. M. Parel, and F. Manns, “In vivo measurement of the human crystalline lens equivalent refractive index using extended-depth OCT,” Biomed. Opt. Express 10(2), 411–422 (2019). [CrossRef]  

35. B. Pierscionek, M. Bahrami, M. Hoshino, K. Uesugi, J. Regini, and N. Yagi, “The eye lens: age-related trends and individual variations in refractive index and shape parameters,” Oncotarget 6(31), 30532–30544 (2015). [CrossRef]  

36. V. M. Hernandez, F. Cabot, M. Ruggeri, C. de Freitas, A. Ho, S. Yoo, J. M. Parel, and F. Manns, “Calculation of crystalline lens power using a modification of the Bennett method,” Biomed. Opt. Express 6(11), 4501–4515 (2015). [CrossRef]  

37. B. M. Maceo, F. Manns, D. Borja, D. Nankivil, S. Uhlhorn, E. Arrieta, A. Ho, R. C. Augusteyn, and J. M. Parel, “Contribution of the crystalline lens gradient refractive index to the accommodation amplitude in non-human primates: in vitro studies,” J. Vis. 11(13), 23 (2011). [CrossRef]  

38. E. A. Hermans, M. Dubbelman, R. G. L. van der Heijde, and R. M. Heethaar, “Equivalent refractive index of the human lens upon accommodative response,,” Optom. Vis. Sci. 85(12), 1179–1184 (2008). [CrossRef]  

39. R. Navarro, J. Santamaría, and J. Bescós, “Accommodation-dependent model of the human eye with aspherics,” J. Opt. Soc. Am. A 2(8), 1273–1280 (1985). [CrossRef]  

40. J. F. Koretz, C. A. Cook, and P. L. Kaufman, “Aging of the human lens: changes in lens shape at zero-diopter accommodation,,” J. Opt. Soc. Am. A 18(2), 265–272 (2001). [CrossRef]  

41. A. de Castro, S. Ortiz, E. Gambra, D. Siedlecki, and S. Marcos, “Three-dimensional reconstruction of the crystalline lens gradient index distribution from OCT imaging,” Opt. Express 18(21), 21905–21917 (2010). [CrossRef]  

42. J. Birkenfeld, A. de Castro, and S. Marcos, “Contribution of shape and gradient refractive index to the spherical aberration of isolated human lenses,” Invest. Ophthalmol. Vis. Sci. 55(4), 2599–2607 (2014). [CrossRef]  

43. D. T. Moore, “Ray tracing in gradient-index media,” J. Opt. Soc. Am. 65(4), 451–455 (1975). [CrossRef]  

44. L. F. Garner and G. Smith, “Changes in equivalent and gradient refractive index of the crystalline lens with accommodation,” Optom. Vis. Sci. 74(2), 114–119 (1997). [CrossRef]  

45. B. Maceo Heilman, F. Manns, A. de Castro, H. Durkee, E. Arrieta, S. Marcos, and J. M. Parel, “Changes in monkey crystalline lens spherical aberration during simulated accommodation in a lens stretcher,” Invest. Ophthalmol. Visual Sci. 56(3), 1743–1750 (2015). [CrossRef]