The purpose of this manuscript is to introduce a new age-dependent model of the human lens with two GRIN power distributions (axial and radial) that allow decoupling of its refractive power and axial optical path length. The aspect ratio of the lens core can be held constant under accommodation, as well as the lens volume by varying the asphericity of the lens external surfaces. The spherical aberration calculated by exact raytracing is shown to be in line with experimental data. The proposed model is compared to previous GRIN models from the literature, and it is concluded that the features of the new model will be useful for GRIN reconstruction in future experimental studies; in particular, studies of the accommodation-dependent properties of the ageing human eye. A proposed logarithmic model of the lens core enables decoupling of three fundamental optical characteristics of the lens, namely axial optical path length, optical power and third-order spherical aberration, without changing the external shape of the lens. Conversely, the near-surface GRIN structure conforms to the external shape of the lens, which is necessary for accommodation modelling.
© 2016 Optical Society of America
Efforts to model the crystalline lens of the human eye are centred around faithfully representing the lens anatomically, biomechanically and optically. Anatomical and biomechanical accuracy is facilitated by the use of higher-order polynomials to describe the lens surfaces, which feature a small number of terms to prevent unnecessary complexity when fitting experimental data. The lens volume, for example, can be specified and used as an important physical constraint for modelling accommodation of the lens , while its growth with age is another important factor to consider. Regarding the optical properties of the lens, the gradient index (GRIN) nature of the lens medium has been represented analytically, with efforts to match its experimentally measured properties both qualitatively—e.g. replicating the shape of the iso-indicial contours— and quantitatively—e.g. reconstructing refractive power. The aim of this paper is to depart from geometry-invariant models and provide a representative model of the anatomy and optics of the ageing human lens.6], ζ has the same functional form (Eq. (1)) extending from the centre of the nucleus to all points on the lens external surface, thus resulting in invariance of geometry.
The GRIN distribution of refractive index in younger eyes has a relatively smooth increase from periphery to nucleus, and is suitably approximated by lower-order polynomials in z and r . However, the older eye has a central plateau of refractive index, with a steep gradient at the periphery. To model this more abrupt distribution, higher-order polynomials have been proposed [5,6,8]. In addition to steepening of the axial refractive index profile with age, there is an age-related relative change in parameter P between radial and axial GRIN profiles. In the older eye, the lens iso-indicial contours approach the GIGL model. In contrast, only the sub-surface region of the younger lens has the same geometry and aspect ratio as the exterior. The internal iso-indicial contours of the young lens depart from geometrical invariance as they approach the nucleus of the lens and exhibit more rapid increase in curvature [4, 9–11]. If we consider Eq. (1), the young lens has different values of P in the axial and radial directions, respectively. This effect is perhaps most clearly demonstrated using magnetic resonance imaging (MRI)  and Talbot interferometry .
A model of the ageing human GRIN lens with different axial and radial refractive index profiles has been proposed by Bahrami et al. . In this “adjustable internal structure” (AIS) model, the external surfaces of the lens are conicoids of revolution. The radii of the lens surfaces control the refractive power of the lens, whereas the conic constants are used to ensure a smooth joining of anterior and posterior segments at the equator. This ties up the conics for that purpose alone, and prevents their use in aberration matching and accommodation modelling. Following on from the conclusion of Bahrami et al. , we see that it is desirable to produce a higher-order description of the iso-indicial contours of the lens. An obvious application of this is the creation of a model which, in addition to predicting optical refractive power and represent qualitatively the shape of the iso-indicial contours, can also match optical aberrations—in particular spherical aberration—of an anatomically realistic, age-dependent lens.
Similarly, the 4-variable GRIN model of Manns, de Castro and co-workers [8, 14] contains different values of P in the axial and radial directions. Also using Pierscionek’s refractive index profile of Eq. (1), Manns and de Castro’s model has a constant refractive index at the surfaces, which are represented by conicoids. This model is included in Table 1, in addition to the AIS model, and other recent models of the GRIN lens. The disadvantage of this model is that the surface asphericity order of 2 leads to a sharp join of the anterior and posterior lens segments at the equator, unless one constrains the surfaces to ellipses, as done in the elliptical model of Smith et al. . This lack of smooth join compromises the bio-mechanical structure of the lens, and the physical meaning of the lens volume is lost. Hence, the volume of the lens cannot be used as a physical parameter for modelling accommodation. Furthermore, in this model, the value of P (Eq. (1)) is a smooth function of the angle θ to the optical axis in the tangential plane, and this function is found through optimisation against experimental data. As a result, the internal structure of the lens is not explicitly defined, and therefore the optical characteristics of the GRIN bulk cannot be defined a priori. Consequently, developing analytical raytracing through the GRIN structure is not possible; analytical raytracing is useful as a starting point for optimisation procedures and especially for solving inverse problems. On the other hand, the functional dependence of P on θ offers a high degree of flexibility of the internal GRIN structure, where the internal iso-indicial contours might take an arbitrary shape. This flexibility is gained at the expense of stability of the optimisation process.
Both the AIS model and that of Manns and de Castro, in spite of their limitations, allow one to independently model the ageing of the GRIN profile in the axial and radial directions. Furthermore, this freedom of having two adjustable (axial and radial) GRIN profiles for the internal lens structure allows decoupling of the optical path length (OPL) and refractive power of the lens, as shown later in this manuscript. In addition, one could also decouple the refractive power and spherical aberration (SA), since the latter is primarily affected by the radial GRIN profile. The effect of radial GRIN profile on SA is seen perhaps most clearly in the 2007 paper by Goncharov and Dainty , where the fourth-order radial term (n2r4) can be used to alter SA independently of power. A sixth-order polynomial representation of the GRIN structure was adopted by Smith et al. , also allowing decoupling of power or OPL and SA. In general, if polynomials are to be used, one needs a radial profile of order at least 4 to adjust SA of the bulk independently of power. For this reason, the models of Liou and Brennan  and Díaz et al.  do not offer the decoupling. The axial profile of Díaz et al. contains trigonometric functions which can be expanded in a power series in z, resulting in an equivalent polynomial of order 3.
In principle, one can decouple the power, OPL and SA by departing from the constraint that the GRIN iso-indicial contours must be concentric with the external surface of the lens. This approach was adopted in the 2014 model of Navarro , where the conic constant of the GRIN bulk is different to the conic constant of the external lens surface, resulting in a peripheral zone of zero axial thickness and constant refractive index. This approach limits the accommodative ability of the model, since the changes in the sub-surface region during accommodation have no concrete physical basis. In terms of accurately modelling the ageing of the GRIN medium, one needs to be able to adjust the structure such that the ratio between the axial and radial profile exponents (Pz and Py) can assume not only an integer value, but any rational value. The polynomial representation of refractive index profile is not ideal because of the discrete nature of the power terms involved; since, for small ray heights, the lowest-order polynomial term is responsible for the lens optical properties. This flexibility is only available in models with decoupled axial and radial profiles, such as those of Manns, de Castro et al. [8, 14], the AIS model of Bahrami et al.  and the proposed AVOCADO model. For instance, in a recent paper by Pierscionek et al. , the representative examples of GRIN profile fits for the 16 year old lens have values of Pz = 2.88 and Py = 1.78, whereas the 91 year old eye has values of Pz = 2.40 and Py = 2.53.
In terms of matching the GRIN structure to experimental data, the stability of the optimisation process depends on the degrees of freedom in the mathematical representation of the GRIN structure. Choosing the lens external shape as a basis to confine the GRIN distribution during optimisation is probably the most sensible way to ensure optimisation convergence, because it limits the number of free parameters to the practical minimum—for example, the anterior and posterior lens radii and conic constants, and central lens thickness. Such use of the external shape to define the GRIN structure is also convenient for performing raytracing from back-to-front. This invertibility for raytracing is not available in the models of Liou and Brennan  or Diaz et al. , since the GRIN is not explicitly given in terms of the external lens shape parameters. Furthermore, confining the GRIN structure to the external shape is essential for modelling accommodation. In this manuscript, we show that the proposed AVOCADO model offers simplicity and flexibility to supersede existing models.
2. Mathematical description of the lens internal iso-indicial contours
In the new lens model, the radii are scaled non-linearly according to an appropriate power, m, of ζ. Another crucial component is that the conic constant, K, of an iso-indicial contour varies with ζ. To be more precise, the shape factor Q = 1 + K of the surface is made to scale with ζ to the power of 2m. This gives the following representation for an internal contour height:Fig. 1(a). Using positive values of m enables us to describe the GRIN distribution in younger eyes, with m varying from roughly 1 to 0 as the eye ages. Eq. (2) simplifies to a geometry-invariant model for m = 0, shown in Fig. 1(b). A negative value of m produces a lens with a larger value of P in the radial profile, as observed in recent studies [4,12], and shown in Fig. 1(c). Note that in Fig. 1, the three lenses have the same external geometry and size; only m is different (m = 1, 0 and −0.5).
The new parameter m allows us to decouple the axial and radial refractive index profiles of the lens. Eq. (1) describes how the refractive index changes as a function of the normalised distance ζ in any direction from the nucleus within the lens. We can also look at how the index varies with distance z along the axis of the lens (ρ = 0); and distance ρ along the radial line z = 0. For the axial index profile—noting that Δn = nc − ns and ζ = z/Tp for the posterior portion of the lens— Eq. (1) becomes:Eq. (2) and re-write it for the posterior surface and set z = 0, giving: Eq. (2) becomes Eq. (3) can be simplified to express any height along this ordinate as a function of ρ0 according to the relation . This can be rearranged for ζ to give the radial profile: Eq. (1) as before, we see that the refractive index now varies radially as: Fig. 2. It is worth noting how the vertical profile of the GRIN structure affects the lens optical power; in the paraxial region, the ρ2P/(m+1) exponent is related to optical power by the scaling of radii: rp(ζ) = ζ2m+1Rp.
3. Average and equivalent refractive indices
The horizontal and vertical profiles of the GRIN can be used to account for the differences in the average refractive index (nav) and the equivalent refractive index (neq) of the lens. It is well known that nav and neq are two phenomena originating from simplification of the GRIN; nav accounts for axial path length measurements made using OCT, and neq for the lens refractive power. A recent paper discusses the idea of separate values for nav and neq, and their difference with age . The average refractive index is related to the time of flight of light rays passing along the optical axis of the lens. In OCT, the information provided to the user corresponds to optical path length rather than physical distance. Thus, we can see that the measurement depends on both refractive index and physical distance. In the case of the GRIN lens, the total axial optical path is given by:
Similarly, neq is used to account for another property of the lens; namely, its optical refractive power. Since refractive power involves only the paraxial properties of the lens, it is convenient to replace the GRIN of the lens with a constant equivalent refractive index, neq, thus greatly simplifying raytracing through the eye model.
The thin-lens power of the AVOCADO model can be approximated when the lens is immersed in a medium of refractive index n as:
In terms of neq, we can also express the power of the homogeneous lens using the thick-lens equation to give:Fig. 1), where the two lenses (Fig. 1(a) & 1(b)) have the same external geometry and values of P, nc and ns, and hence nav; but different values of m. From Eqs. (8) & (9), it follows that by altering only parameter m from 0 to 1, it is possible to induce large changes of ≃ 0.04 in neq, and ≃ 10 D in the refractive power of the lens. In the GIGL model, this is not possible because m = 0. Figure 3 below is a plot of lens power (F) vs m for different values of P, showing that changing m from 0 to 1 induces large changes in F. Shown in this figure are the powers calculated using the thin-lens formula of Eq. (8), and exact raytracing introduced in the following section. This figure shows that while exact raytracing is preferred, Eq. (8) provides a convenient means of predicting lens power, with a typical over-estimate of approximately 0.5 D or less. A more accurate expression for lens power can be derived and used in place of this thin-lens formula .
It can be argued that, in the GIGL model, the extra parameter required to distinguish between nav and neq could be provided by varying either nc or ns. While it is true that an age-related change in these parameters will alter the values of both nav and neq, it changes both simultaneously—they cannot be decoupled in the same way as with the new parameter m. To see this, we note from Eq. (7) that nav is affected if we alter either ns, nc, Δn or P. If we keep nc constant, and assuming we have a value of nav that we wish to match to experimental data, is there a way to alter Δn and P to keep nav constant while changing neq sufficiently? It can be shown that typical values of Δn = 0.04–0.06 require a range of P = 3.9–6.1, resulting in a maximum change in neq of only 0.0005. Thus, this approach of decoupling neq and nav does not give a satisfactory result.
Given the external shape of the lens, the GRIN medium can be optimised to account for the experimentally observed measurements of refractive power and hence neq, as well as matching the axial OCT measurements of optical path length using nav, as discussed above.
Decoupling the axial and radial refractive index profiles of the lens is an important aspect of lens modelling. For example, the 2008 study by Kasthurirangan et al.  and the 2015 study by Pierscionek et al.  show that both young (accommodated and unaccommodated) and old lenses can have larger values of P in the radial section. This result can be reproduced with the extra flexibility afforded by decoupling the axial and radial exponents in the new AVOCADO lens model. The flexibility of the new model allows the generation of a GRIN structure with any particular axial and radial refractive index profiles. The above results show that the new model is capable of producing a lens that has a larger equivalent value of P in either the axial profile (m > 0) or radial profile (m < 0).
Regarding the equivalent refractive index of the lens, neq, we emphasise the role that it plays: for a given ocular geometry (curvatures of the ocular surfaces and intraocular distances) and refractive indices of the cornea, aqueous and vitreous, neq makes up the remainder of the optical power of the eye. neq becomes particularly useful with personalised ocular modelling, where the curvatures and conic constants of the ocular surfaces can be measured, for example, and used to provide the framework on which optimization of the GRIN lens model can be based. In this type of personalised modelling, we note that it will be beneficial to use both the radius and conic constant of a surface to define the properties of the optical zone for image formation. The flexibility of the new AVOCADO lens model is important here since it allows us to use the radius and conic constant for that purpose, while the B-coefficient ensures a smooth join at the boundary.
4. Exact raytracing for power and SA analysis
In addition to the simplified paraxial analysis above, we can use exact raytracing to calculate the refractive power and aberrations of the lens, by solving the differential ray equation [24,25]. Already mentioned is the possibility of decoupling the axial optical path length and refractive power of the lens, F. We can show this by example below, where the value of P is fixed so that OPL is unchanged, but m is altered. The result is that the refractive power changes considerably, shown in Fig. 3; hence, nav and neq are decoupled.
From this figure, it becomes apparent that the AVOCADO model offers an additional degree of freedom when reconstructing the internal structure of the lens. OCT measurements provide information on time-of-flight of the light through the lens, and for a given lens thickness Δz, one can determine nav (Eq. (6)); also, for given values of Δn and nc, P can be determined from Eq. (7). If one can also measure F and the external shape of the lens, the value of m can be found using exact raytracing—see Fig. 3; alternatively, the approximate value of m can be found from the thin-lens approximation (Eq. (8)). Since refractive power and time-of-flight are decoupled, as outlined above, the method is feasible as a practical way of reconstructing the GRIN profile.
Now, we can analyse spherical aberration (SA) of the lens as a function of m. From population studies of ocular aberrations, SA is the only classical (rotationally symmetric) aberration term to have an on-axis average that is significantly different from zero ; hence we can use the value of SA as a useful metric for verification of the model. The AVOCADO model reduces to the GIGL model for m = 0, which is viable in terms of SA prediction . Note that SA is also strongly affected by the choice of the lens conic constants. Figure 4 shows a series of plots of the longitudinal spherical aberration (LSA) of the lens for different values of m, with all other lens parameters fixed, with P = 3. The data from numerical raytracing are shown in red, while 8th-order (even coefficients only) polynomial fits are shown in blue. The quadratic term (ρ2) in this polynomial fit of LSA is converted to transverse SA and then integrated to give the wave-front aberration polynomial coefficient (W4,0)  for a pupil diameter of 4 mm. Assuming the system has zero defocus, this term can be converted to Zernike SA using the formula:Figure 5 shows the change in as a function of m for several values of P. Comparing Figs. 3 & 5, one can see that with increasing m, the lens becomes more powerful and its SA is more pronounced. In particular, for m = 0.4 and P = 3, , which is comparable to the GIGL case where m = 0 and —see Fig. 4. For values of P smaller than 3, the value of becomes rapidly more negative with increasing m—see Fig. 5. It can be seen from Figs. 3 & 5 that for a given value of m, one obtains unique values of F and . To decouple optical power (F) from third-order SA ( ) for a given value of P, we will introduce the logarithmic description of the lens core in the following section.
5. Volume and aspect ratio of an internal iso-indicial contour
Accommodative changes in lens radii and thickness can be accompanied with changes in conic constants to maintain constancy of volume. It can be shown that the volume of the entire lens is given by the formula:1]. This formula demonstrates the role of conic constants in volume specification. We can obtain an expression for the volume of an internal contour of the lens by taking the equation above and replacing the constants by their scaled counterparts, such that the volume as a function of the normalised distance ζ is given by:
The aspect ratio of the external shape of the lens is given by A0 = 2ρ0/T, where ρ0 is found from Eq. (4). From Eq. (5), ρ(ζ) = ζm+1ρ0, and since the axial thickness t(ζ) = ζT, we see that the aspect ratio of an internal contour is given by:
The aspect ratio of an internal core can be plotted as a function of ζ, shown in Fig. 6. This figure highlights the geometrical significance of w, as it is seen to alter both the rate of change of aspect ratio and the limiting value in the lens core, as ζ → 0. Equation (11) can be solved to find the value of w that gives a desired aspect ratio at a given iso-indicial contour depth ζ.
Equation (14) illustrates the optical significance of parameters m and w, as follows. The radius Ra is scaled by ζϕ(ζ), and is directly responsible for the lens power. This effect can be seen in Fig. 7, where m produces a large change in power. Furthermore, in Eq. (14), the shape factor (1 + Ka) is scaled by the factor ϕ(ζ); hence, in addition to the geometrical significance of w outlined above, w has optical significance in altering SA according to this scaling of the shape factor.
Now having an additional parameter w, we can easily decouple SA and optical power. Figure 8 is a plot of SA vs m for P = 3, showing that altering w affects the change of SA with m. The LSA of a selection of lenses from Fig. 8 is shown in Fig. 9. These lenses all have the same optical power F = 28.6 D; this can be achieved by adjusting m for a given w. As a result, each lens has a different amount of SA, while we have constant values for optical power, F; and P, which corresponds to constant axial optical path length. It is evident that SA is decoupled from both P and refractive power, and hence SA is decoupled from nav and neq. Having w as an additional free parameter not only allows us to change the sign of SA, but also helps to change the higher-order SA, so that the LSA curves develop the opposite trend (bending towards the lens) at intermediate pupil heights—see the case of w = −0.3 in Fig. 9. In this case, the higher-order SA counteracts lower-order SA, and the choice of m and w will affect this SA balancing.
The logarithmic model of the lens core with P, m and w enables decoupling of three fundamental optical characteristics of the lens, to wit optical power, axial optical path length and third-order spherical aberration, without changing the external shape of the lens. On the other hand, the near-surface GRIN structure conforms to the external shape of the lens, which is necessary for accommodation modelling.
6. Discussion and conclusion
The use of a cubic profile is perhaps the simplest way to describe the lens external shape in an anatomically realistic way, while the conic constants can be used for aberration matching. The inclusion of parameter m is a convenient way to attain a particular value for the ratio nav : neq for a given external lens geometry. Furthermore, the age-dependence of m and P provides a useful way to account for the experimentally observed change in iso-indicial contour shape. Using the three parameters P, m and w, this new model can be used as a fundamental tool for experimental research, so that optical path, refractive power and SA measurements will allow reconstruction of the internal structure of the lens in a way that was not possible previously.
When compared to second-order conic descriptions of the lens, the cubic profile of this model imparts an anatomical accuracy to the model. With the resulting smooth equatorial join, the lens volume becomes physically meaningful and is a powerful constraint for accommodation . This is essential for study by finite element methods and analysis of the accommodation-dependent properties of the ageing human eye. Therefore, we are one step closer to the ultimate goal of having a biomechanically accurate age-dependent accommodative model of the human lens, capable of predicting optical path length, refractive power and spherical aberration.
One interesting application of this core flexibility could be in the study of accommodation, where the external aspect ratio A0 of the lens is changing. If the older eye has a stiff core in the very centre of the lens, then the aspect ratio of this core will not change with accommodation. Hence, we can apply Eq. (13) to the two separate accommodation states of the lens, with the condition that the aspect ratio A of the core does not change. If we wish to keep m0 constant with accommodation, we can use this condition to solve for an accommodation-dependent value of w. Alternatively, we need not use the core in the limiting case ζ → 0; we could use any particular value of ζ—e.g. a value of might represent the lens nucleus. Together with an analysis of volume using Eq. (10), this could form the basis of a mathematical description of presbyopia.
This research was funded in part by The Irish Research Council, application RS/2012/351.
References and links
1. C. J. Sheil, M. Bahrami, and A. V. Goncharov, “An analytical method for predicting the geometrical and optical properties of the human lens under accommodation,” Biomed. Opt. Express 5, 1649–1663 (2014). [CrossRef] [PubMed]
2. B. K. Pierscionek, “Presbyopia—effect of refractive index,” Clin. Exp. Optom. 73, 23–30 (1990). [CrossRef]
4. 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. Vis. Sci. 49, 2531–2540 (2008). [CrossRef] [PubMed]
5. 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, 2175–2185 (2007). [CrossRef]
7. A. V. Goncharov and C. Dainty, “Wide-field schematic eye models with gradient-index lens,” J. Opt. Soc. Am. A 24, 2157–2174 (2007). [CrossRef]
8. 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, 21905–21917 (2010). [CrossRef] [PubMed]
10. C. Jones, D. Atchison, R. Meder, and J. Pope, “Refractive index distribution and optical properties of the isolated human lens measured using magnetic resonance imaging (MRI),” Vision Res. 45, 2352–2366 (2005). [CrossRef] [PubMed]
12. 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, 1–13 (2015).
13. M. Bahrami, A. V. Goncharov, and B. K. Pierscionek, “Adjustable internal structure for reconstructing gradient index profile of crystalline lens,” Opt. Lett. 39, 1310–1313 (2014). [CrossRef] [PubMed]
14. F. Manns, A. Ho, D. Borja, and J.-M. Parel, “Comparison of uniform and gradient paraxial models of the crystalline lens,” Invest. Ophthalmol. Vis. Sci. 51, 789 (2010).
16. H.-L. Liou and N. A. Brennan, “Anatomically accurate, finite model eye for optical modeling,” J. Opt. Soc. Am. A 14, 1684–1695 (1997). [CrossRef]
17. J. A. Díaz, C. Pizarro, and J. Arasa, “Single dispersive gradient-index profile for the aging human lens,” J. Opt. Soc. Am. A 25, 250–261 (2008). [CrossRef]
20. S. Giovanzana, R. A. Schachar, S. Talu, R. D. Kirby, E. Yan, and B. K. Pierscionek, “Evaluation of equations for describing the human crystalline lens,” J. Mod. Opt. 60, 406–413 (2013). [CrossRef]
23. Alexander V. Goncharov and Conor J. Sheil, Applied Optics Group, National University of Ireland, Galway, are preparing a manuscript on raytracing through GRIN media defined as an infinite number of thin shells with a known shape.
25. B. D. Stone and G. W. Forbes, “Optimal interpolants for runge–kutta ray tracing in inhomogeneous media,” J. Opt. Soc. Am. A 7, 248–254 (1990). [CrossRef]
26. L. N. Thibos, X. Hong, A. Bradley, and X. Cheng, “Statistical variation of aberration structure and image quality in a normal population of healthy eyes,” J. Opt. Soc. Am. A 19, 2329–2348 (2002). [CrossRef]
27. J. C. Wyant and K. Creath, Applied Optics and Optical Engineering (Academic Press, 1992), vol. 11, chap. 1: Basic Wavefront Aberration Theory for Optical Metrology.