Single function crystalline lens capable of mimicking ciliary body accommodation

: The lens is a complex optical component of the human eye because of its physiological structure: the surface is aspherical and the structural entities create a gradient refractive index (GRIN). Most existent models of the lens deal with its external shape independently of the refractive index and, subsequently, through optimization processes, adjust the imaging properties. In this paper, we propose a physiologically realistic GRIN model of the lens based on a single function for the whole lens that accurately describes diﬀerent accommodative states simultaneously providing the corresponding refractive index distribution and the external shape of the lens by changing a single parameter that we associate with the function of the ciliary body. This simple, but highly accurate model, is incorporated into a schematic eye constructed with reported experimental biometric data and accommodation is simulated over a range of 0 to 6 diopters to select the optimum levels of image quality. Changes with accommodation in equatorial and total axial lens thicknesses, as well as aberrations, are found to lie within reported biometric data ranges.


Introduction
In spite of the apparent simplicity of the human eye, which incorporates two lenses, it is an optical system with extraordinary complexity. In particular, the refractive index of the inner lens, known as the crystalline lens is not homogeneous but distributed as a gradient index (GRIN) [1][2][3][4][5]. When an observed object is placed at different distances from the eye, the lens changes its shape to maintain the focus of the image on the retina, a process known as accommodation [6]. As the lens accommodates, its shape changes and its GRIN is modified [7]; the exact nature of the latter is not known.
Over the last four centuries, different schematic eye models have been studied and their accuracy and complexity have continued to evolve. Earlier schematic models incorporated spherical refractive surfaces whose geometrical centres were placed along the optical axis, considered elements with homogeneous refractive indices, small pupils and were designed for objects close to the optical axis. These are known as paraxial schematic eye models [8,9].
The simplest of these is the Emsley reduced model in which the eye is represented by a single refractive surface; the anterior surface of the cornea [8][9][10]. The Gullstrand-Emsley model increases the number of refracting surfaces to three; one for the cornea and two for the lens.
The refractive indices of the ocular media, the aqueous and vitreous, were up to 1.416 which is substantially higher than the real values of 1.336 [9,10]. The Le Grand full theoretical model is represented by four surfaces; two for the cornea and two for the lens [11].
Robust paraxial schematic eyes have included an inhomogeneous refractive index for the crystalline lens, as this feature has been recognized for over a hundred years [12]. Gullstrand was the first to propose a discrete model of this type using four iso-indicial surfaces for the lens. This model is known as Gullstrand's No. 1 "exact" eye [13,14]. In this model, the core of the lens has a high refractive index, limited by its interior surfaces, while the outer section has a lower refractive index limited by the external surfaces. This was the first attempt to represent the inhomogeneity of the lens using discrete increments of refractive index and including separate core and cortical sections. The next evolutionary step was to consider a lens with a continuous gradient refractive index. An example of this is the Maxwell fish-eye spherical lens [15] or the Luneburg lens [16]. However, these lenses have a spherical shape that does not correspond to the shape of the human lens.
Schematic eye models can be put into two different categories, the "encyclopaedic" type, which consists of a thoroughly compact description of knowledge about geometric and mechanical ocular properties, and "the toy train" type, which represents a pragmatic tool that models the functioning of real eyes without attempting to accurately reproduce the anatomical nor the mechanical properties of the eye [14]. Nevertheless, the best choice will depend on the particular problem that is being addressed. The simplest one covering the needs of a particular application is a good alternative [6]. An example of this is the recent work by Liu and Thibos, in which no GRIN structure is used to describe the lens [17].
Current imaging technologies allow precise in vivo and in vitro measurements of the physiological parameters of the lens, such as the anterior and posterior radii of curvature of the lens surfaces as well as its thickness [3][4][5][18][19][20][21]. Optical coherence tomography (OCT) is a leading modality in imaging ocular parameters [22][23][24][25][26] and this technique has been used to estimate the refractive index of the lens [27,28].
Due to the asymmetric shape of the lens, the GRIN has commonly been represented by two different sections separated by a plane or curved surface that intersects the lens capsule at its equator [6,[29][30][31][32][33][34][35][36][37]. To be physiologically correct, the corresponding GRIN functions of each segment must change smoothly from one to the other side of the separating surface, i.e., they must obey continuity conditions. The lens capsule is described in a similar way, early GRIN models incorporate two surfaces, an anterior and a posterior, that usually do not match with the inner iso-indicial surfaces of the GRIN structure. Nevertheless, there are few exceptions such as those of Navarro and Goncharov that accurately match the GRIN distribution with both external surfaces [31][32][33][34][35].
The first accommodative models were defined only for discrete accommodative values, so they could not be used to simulate a dynamical process [11,13,38,39]. Later improved models were defined for continuous accommodative values taking into account changes in the radius of curvature of optical surfaces, equatorial and axial lens thickness and index of refraction [34][35][36][40][41][42].
As mentioned above, most of the continuous accommodative models have incorporated GRIN distributions described by two functions. However, there are a limited number of models that describe the geometry and/or the GRIN distribution of the lens with a single function [42][43][44]. The advantage of single function models over the others resides in the fact that the smoothness of the refractive index profile of the lens is guaranteed within the whole volume without having to impose any boundary conditions.
Kasprzak introduced an elegant accommodative lens model with a single function based on a combination of hyperbolic trigonometric functions. His model introduced a set of parameters taken from biometrical data at different accommodative states to represent the capsule and refractive index of the lens [43]. The refractive index of the lens was obtained by feeding the parameters of the hyperbolic function with the corresponding biometrical data, radii of curvature and the lens thickness. Some of the parameters were obtained by optimization processes imposing the restriction that the volume of the lens remained constant during accommodation. Later, this model was applied by Popiolek-Masajada and Kasprzak to obtain an alternative GRIN model [45]. However, the single function introduced by Kasprzak in Ref. [43] is only used to describe the capsule of the lens, and the refractive index distribution is described by a different function whose iso-indicial surfaces do not match with the external shape of the lens. Popiolek-Masajada and Kasprzak [46], used the same hyperbolic single function to describe the lens but considered a homogeneous refractive index. The accommodative process was modelled by changing the curvature and thickness taken from biometrical data. These geometrical changes were not enough to obtain the expected increase in refractive power so that the homogeneous index of refraction of the lens had to be increased to generate an effective refractive power [46].
Subsequently, Huang and Moore also used a single continuous function to describe the total GRIN structure of the lens [42]. They generalized the model of Liou-Brennan by incorporating terms into the GRIN single function that take into account the asymmetry of the lens. The anterior and posterior aspheric surfaces for the different degrees of accommodation were obtained using a numerical optimization method that minimizes transverse aberrations. Díaz et al. (2008) presented an age-dependent human lens model using a single function for the GRIN profile, modelling the posterior and anterior aspheric faces using conicoid surfaces but did not consider accommodation [44].
In this paper, we present an accurate accommodative crystalline lens model with a single GRIN function that delivers the GRIN distribution and external shape of the capsule, both of which are in good agreement with recently published biometrical data [47]. With this lens model we constructed a schematic eye for which, with the exception of one, the parameters of the GRIN function were obtained through an optimization process to model a 26-year-old eye. Our investigations show that with the free parameter we can mimic the behavior of the ciliary body that provides the force on the crystalline required for accommodation. By simply changing the value of that parameter, the accommodative process can be modelled. It simultaneously provides the continuous changes in the geometry of the anterior and posterior aspheric faces of the lens capsule and of the internal GRIN distribution in order to enable focusing over a wide range of object distances while maintaining constancy of the image distance, as occurs in the human eye.
The new schematic eye model can simulate accommodation from 0 to 6 diopters, providing high image quality with aberrations that are within the range of biometric data.

Poisson-Gauss function and its characteristics
In statistics and probability, the Poisson probability distribution for the occurrence of discrete events, usually of low likelihood, occurring in a given spatial interval is defined as where m is the number of events and z is the interval for z>0 [48,49]. Fig. 1(a) shows the skewed bell shape of the Poisson distribution as a function of z in dimensionless units for several values of m ∈ N, where N is the set of natural numbers. The GRIN in the human lens is asymmetric along the optical axis and in the process of accommodation it changes its width. The Poisson function, as can be observed in Fig. 1(a), also presents a similar asymmetric behavior that can be associated with the axial GRIN of the human lens. For this purpose we propose a function where all of the variables and parameters are real numbers and it will be constructed by merging the Poisson function with a two dimensional Gaussian function, namely: where r = x 2 + y 2 is the cylindrical radial coordinate, z>0 is the longitudinal axial coordinate and m, b, a z , a r are numerical parameters to be determined from biometrical data of a real crystalline lens. Contrary to its former definition in Eq. (1), m is not an integer anymore but will take real values. This function will be referred to as Poisson-Gauss (PG) function. To simplify the visualization and notation of Eq. (2) the 3D PG function will be projected to the plane y − z, showing only the corresponding 2D spatial coordinates, namely PG(y, z), and considering as implicit the rest of the parameters of the PG function. Fig. 1 shows the steps in the construction of the PG function: Fig. 1(a) corresponds to the 1D Poisson distribution along the z-axis, Fig. 1(b) depicts a 2D Poisson distribution in the y − z plane as a function only on z, Fig. 1(c) shows a Gaussian distribution and finally the PG function density distribution resulting from the product of the two previous distribution patterns is shown as a 2D density plot in Fig. 1(d). The PG function is unbound extending over the entire half plane z>0 whereas the lens must have a finite size, then it is required a finite domain. For that purpose we delimit the domain of the PG function by cutting its surface with an intersecting flat surface parallel to the y − z plane at a given height (to be defined later) as shown in Fig. 2(a). The intersection creates a contour curve that projects onto the y − z plane and delimits the domain where the lens will be defined. In Fig. 2(b) shows the resulting 2D density plot of the PG function in the new domain. It should be noted that the latter resembles a human crystalline lens.
Next, we show that it is possible to model the human lens with the PG function by relating its geometry with that of the GRIN properties of the lens. We can associate the refractive index at the centre of the human lens with the location of the maximum of the PG function, and show that in the plane y − z this is given by the coordinates The use of the subindex e will be clear in the subsection below. From its definition in Eq. (2) and Eq. (3) we observe that the parameters a z , a r , b and m allow manipulation of the geometric properties of the PG function. The maximum value of PG(y, z) is denoted by This value is illustrated in Fig. 3 which also shows the PG function 3D contour plot starting from the aforementioned intersecting flat surface at height h. The inset, shows the projected 2D contour plot. The GRIN of the lens can also be associated with the statistical properties of the PG function noting that the square root of twice the variance σ is a good measure of the extent of a bell shaped function [50]. Along the z−direction, the PG function extends over a value of σ z from an average value located at z = µ z , y = µ y = 0 as shown in Appendix A. Both µ z and σ z are functions of the parameters of the PG function, that is, It should be noted that neither µ z nor σ z depend on the transverse behavior (x, y) of the PG function. From Fig. 3, the height of the intersecting plane is given by evaluating the PG function at y = µ y = 0, z = µ z + σ z , namely The corresponding intersection results in a level curve given by PG(y, z) = h, that can be rewritten as z 2 The family of the level curves, shown in Fig. 3, is therefore given by where c parametrizes the family and lies in the interval where h and l are the minimum and maximum values of the intersected PG function given by Eqs. (4) and (5), respectively, as shown in Fig. 3. Therefore, the finite domain D 0 of the lens in the z − y plane is limited by Eq. (6) and is defined by and shown in Fig. 2(b). Considering the lens as a solid of revolution of the above region along the z-axis, the expression of the finite 3D domain D of the refractive index of the lens is The external contour of the lens is determined by Eq. (6) with c = h, which when solved for y can be written as This equation can be related to the physiological parameters of the lens: equatorial radius R e , and the anterior and posterior thicknesses relative to the equatorial axis, z a and z p respectively. Fig. 4(a) shows the surface of revolution generated by the above equation that encloses the domain D and Fig. 4(b) shows how the external contour curve changes with increasing the value of m in Eq. (10). The increase in curvature of the anterior part of the curve, results in an increase of its refractive power and this variation in shape of the external contour is incorporated in the model as part of the accommodative process. The values used correspond to experimental data from 0 to 6 diopters. Fig. 4(c) illustrates the main physiological parameters of a radially symmetric crystalline lens that will be used in the next section.

Poisson-Gauss GRIN crystalline lens
Applying all of the above, the geometric and physical properties (GRIN) of the lens can be completed. The equatorial radius is given by the maximum of the equation that defines the capsule, Eq. (10), and is located at z = z e , namely The width parameters of the lens z a and z p can be obtained setting y + (z; h) = 0. In accordance with Eq. (5), it can be seen that z = µ z + σ z is one of the solutions of this equation: The second solution to y(z, h) = 0 gives z a ; this is done numerically since it is a transcendental equation. As mentioned above, the values of µ z and σ z are calculated as described in Apendix A.  Table 3) Taking the above equations into account, the 3D GRIN profile of the crystalline lens can be described as: for z>0, r ≥ 0, where l and h are given by Eqs. (4) and (5) respectively and n c and n s are the values of the central and peripheral refractive indices of the lens, respectively. When the PG function reaches its maximum l, the numerator in Eq. (13) equals l − h, so that the maximum value of n(r, z) is n c . Analogously, when the PG function reaches it minimum h, the numerator equals zero and the minimum of n(r, z) is n s , that corresponds to the surface of the capsule. Between these two values, h, l, the proposed PG function will determine the GRIN distribution of the lens within the domain D defined in Eq. (9). Eq. (13) will be referred to as PG GRIN lens. This last function in Eq. (13) is one of the main contributions of this work; it provides the refractive index and the shape of the capsule of the human lens. In the next section the results presented above will be related to the physiological parameters used in the specialized literature.

Accommodative properties of the PG GRIN lens
Accommodation requires that the lens be subjected to stresses that modify its shape and change the GRIN distribution and consequently the refractive power of the eye. It is known that during this process the lens volume V remains constant [34,51,52]. This process can be modelled by manipulating the geometrical parameters of the PG refractive index while maintaining constancy of lens volume with accommodation (the volume of D defined by Eq. (9)).

Volume of the lens
The volume V of the 3D domain D is the volume of the solid of revolution obtained by rotating the function y(z; h) given in Eq. (10): where Eq. (12) defines the superior limit; setting z p + z a = d as the thickness of the lens determines the inferior limit. The sum of the terms within the parenthesis defines a length that will be called L 0 . Therefore the last equation can be written as V = a 2 r πL 0 , that is the volume of a cylinder of radius a r and length L 0 . The last equation can be rewritten as Since L 0 does not depend on a r , the latter modifies according to Eq. (15) to maintain constant volume during accommodation. This provides an analogy for accommodation with a cylinder of length L 0 and radius a r ; if the length L 0 of the cylinder increases the radius must decrease following Eq. (15) to keep the volume constant.

Setting the PG GRIN parameters from the physiology of the lens
To represent accommodation realistically with the PG refractive index (13) the radius of curvature of both the anterior and posterior surfaces of the lens must decrease while the lens thickness increases [9]. Hence, it is necessary to obtain the radii of curvature of both the anterior and posterior surfaces of the lens. The radius of curvature R of the iso-indicial contours along the z-axis can be obtained, as shown in Appendix B, by setting y = 0 in Eq. (6); that is Since z>0, R(z) is always finite. The radii of curvature of the anterior and posterior surfaces of the lens can be written as R a = R(z e − z a ) and R p = R(z e + z p ), respectively. Hence, the gradient of the curvature radius dR(z)/dz, of the anterior part of the lens is −(a r /a z ) 2 − ma 2 r /(2z 2 ), and for the posterior part is (a r /a z ) 2 + ma 2 r /(2z 2 ). Since this gradient differs from −1, the isoindicial contours are not concentric (see Fig. 3), as is required for a realistic lens model [53].
By increasing Poisson parameter m (see Fig. 1(a)) the accommodative process can be modelled since the PG GRIN function becomes wider and the corresponding curvatures change, Eq. (16). Below we show that this is indeed the case. Now, it is necessary to find the PG GRIN lens parameters from experimental data. Shao et al. report in vivo measurements of the radii of curvature of the anterior R a and posterior R p surfaces as well as of the central thickness d of the lens under relaxed (0 D) and accommodated (4.00 D) states [21], see Table 1. Through an optimization process, we obtain the parameters of the PG GRIN lens, m, b, a z , a r , to reproduce the relaxed state measurements of R a , R p and d, very close to the experimental error tolerance; R a and R p , deviate from the maximum of their corresponding experimental values by 2.5% and by 0.5%, respectively, while d deviate from the minimum experimental value by 0.83%, as presented in Table 1. Once the PG GRIN lens parameters are determined, z a , z p and R can be obtained.
The accommodated state will be simulated by increasing the Poisson parameter m to reproduce the experimental data within a very good approximation; the PG GRIN physiological parameters R a , R p and d, deviate from the maximum of their corresponding experimental values by 1.4%, 13.8% and 4.0%, respectively. The volume V of the lens with the parameters shown in Table 1 is V = 106.5 mm 3 . From Fig. 3 we can see that the PG function is always cut at height h, were R c is the axial radius of curvature and Q is the conic constant that describes the level of asphericity. As can be observed in Eq. (6), the logarithmic term does not allow the description of our lens model in terms of a conic section. However, we can expand the logarithmic term in Eq. (6) in the neighborhood of a given z 0 . After some algebra the following equation is obtained z 0 is the position where the curve under study crosses the z−axis, namely, y(z 0 ) = 0 (if z 0 = z e − z a (z 0 = z e + z p ), the curve under study corresponds to the external anterior (posterior) surface of the lens). Therefore, the expansion in Eq. (18) is written in the form of Eq. (17). Note that comparing Eqs. (17) and (18), the expression of the radius of curvature rigorously obtained in Eq. (16) is exactly reproduced in the expansion. Therefore, the conic constant Q l associated to the PG lens reads It is important to see that even though we have been able to relate the conic representation Eq. (17) with Eq. (6) using the expansion of the logarithm function up to the second order, we remark that the actual lens surfaces are not necessarily described by conic curves. The contribution of the logarithm function is relevant; when using its expansion it should consider a larger number of power terms since the series of the logarithm has a very slow convergence. In the next section we will show that having the whole expression allows for the building up of an efficient schematic eye.

Schematic eye and imaging properties
We constructed a schematic eye model to test the imaging capabilities of the proposed PG GRIN lens. We take the conic representation of the cornea described by Eq. (17). As in Ref. [46], the biometrical parameters of the cornea were obtained from Ref. [55], and those of the rest of the eye from Ref. [21] and are shown in Table 2.
We studied several situations, two of them are presented in Table 1. The first case corresponds to the relaxed state of the lens. We consider light coming from a point source at 6 metres from the external surface of the cornea and traversing the unaccommodated lens, Fig. 5(a). By a finite Table 2. Biometric data: the radius of curvature of the external R c a and internal R c p surfaces of the cornea and their respective conic constants Q a and Q p , the thickness of the cornea d c , the distance between cornea and lens d c−c , in mm. The refractive indices of the cornea n cor and of the aqueous humor n a . The central and peripheral refractive indices of the lens, n c and n s , resp., are taken from Ref. [ ray-tracing procedure [9] we found that the paraxial image position is at 23.44 mm (on the retina).
In the second case, the accommodative state of 4 D is simulated using the same position of the light source. To do this we have set m = 7.4 and observed that the paraxial image position is at 22.31 mm (prior to the position of the retina shown in Fig. 5(a)). This is equivalent to a myopic eye of 4.03 D. In order to focus the image on the retina, the light source needs to move closer to the eye. We placed the point source at 270 mm from the cornea and observed that the paraxial rays are focussed at 23.44 mm, that is, on the retina (see Fig. 5(a)). [It should be noted that given that the distance between the object and the eye is 270 mm and the axial length of the eye is under 25 mm, rays diverging from the object appear to be parallel at incidence on the cornea].
As shown at the end of Section 3, the accommodation process can be simulated by increasing the Poisson parameter m using fixed values of b, a z , n c and n s . Although the experimental data reported in Ref. [21] are given for two accommodative states, our model allows calculation of a continuous of accommodative states by simply increasing the parameter m, with the corresponding a r given by Eq. (15) to maintain constancy in the lens volume. The particular cases ranging from 0 to 6 diopters are shown in Table 3. Fig. 6 shows the axial radius of curvature of the anterior and posterior apices of the lens as a function of m, i.e., as a function of accommodation. As can be observed, the posterior radius of curvature changes more slowly than the anterior with accommodation, in agreement with the expected physical behavior of the biological lens; the anterior segment of the lens widens more than the posterior segment [9]. The PG GRIN model agrees with this behavior, as shown in Figure 4(b). For this reason we propose that the dynamics of the ciliary body can be represented by the changes of the parameter m, which parametrizes the function of the ciliary body.  Table 3). Table 3. PG GRIN parameters. The parameters are written in mm, with exception of m which is dimensionless and A that is the accommodation written in Diopters. As in Table 1, b = 0.67 and a z = 3.15. The accommodation process is shown in Fig. 4 Fig. 7. Analogous to an elastic medium, with accommodation the GRIN profile could be expected to flatten in the axial direction and steepen in the equatorial direction from the centre to the periphery of the lens. Fig. 7 shows that this is indeed the case for our model. This phenomenon was confirmed experimentally by using advanced 3T clinical magnetic resonance imaging techniques [47].
To test the imaging capabilities of our schematic eye, finite ray tracing was performed. Figure 8 shows, in blue, the resulting longitudinal spherical aberration (LSA) as a function of the height of the ray for the unaccommodated and the accommodative states (4.03 D) while in red is shown the polynomial fit of the LSA of order 8th, with only even order terms. In both cases, the LSA remains below 1 D for a 6 mm pupil diameter. This level of aberration lies within the experimental measurements, as reported previously [9], for the unaccommodated state which ranges from 0 to 1 D for ray height between 0 and 3 mm. Also, in the same ray height range, the LSA of the PG GRIN model is lower than predicted in the single function model of Popiolek-Masajada and Kasprzak [45], which ranges from −1 to 0.5 D, the model of Navarro et al. [41], which ranges approximately from −1.5 to 0 D and Gullstrand's model [13] from −4 to 0 D. In normal lighting conditions, the pupil diameter is around 4 mm diameter. In this case, for the unaccommodated state, the LSA lies between −0.1 and 0 D, and between −0.5 and 0 D for the accommodative state. The change in sign of the LSA as a function of ray height is also observed in the Popiolek-Masajada and Kasprzak model [45]. The fourth Zernike polynomial coefficient Z 0 4 can be calculated with the coefficient of second order of the polynomial fit of the LSA shown in Figure 8, since that coefficient times 1/4 gives the wavefront aberration coefficient W 4,0 [9]. In this case, Z 0 4 is related to W 4,0 through Z 0 4 = W 4,0 /(6 √ 5). Therefore, for a pupil diameter of 6 mm, we obtain Z 0 4 = −0.0014 ± 0.0008 µm, for the unaccommodated state, whilst Z 0 4 = −0.0055 ± 0.0007 µm, for the 4.03 D accommodated state. These values lie within the experimental data reported [56][57][58]. The fact that the coefficients are negative is likely to be because the model is based on a relatively young eye [57]. The fact that the increase of a single parameter permits a departure from the unaccommodated state to the 4 D accommodative state obtaining values close to the ones measured in Ref. [21], with values of spherical aberration within experimental ranges, makes it possible to predict further accommodative states, as shown in Table 3.

Discussion
The model described in this paper can simulate accommodation with a single parameter inducing concomitant changes in the the GRIN profile and surface curvatures that are biologically feasible. We showed in Section 3 that the values of R a , R p and d of our model lie close to the error tolerance of the measurements reported in Ref. [21]. These values are also close to the measurements of radius of curvature performed by Ortiz et al. [5]: for one subject they find R a = 12.48 ± 0.20 mm, R p = 7.25 ± 0.25 mm and d = 3.18 ± 0.02 mm, while for another subject R a = 11.43 ± 0.04 mm, R p = 6.12 ± 0.15 mm and d = 3.36 ± 0.02 mm; both were for the unaccommodated state, in which a conic external shape of the lens was assumed. Conversely, Khan et al. [47] report a decrease of equatorial diameter of −0.30 ± 0.23 mm and an increase of the axial thickness of 0.34 ± 0.16 mm with accommodation from 0 D to 5 D. The changes in the equatorial diameter and axial thickness of our model with the same level of accommodation are, as can be seen in Table 3, −0.4 mm and 0.42 mm, respectively. Hence, the prediction of our model lies within the experimental range of accommodative changes.
As regards the asphericity of our model, we performed an approximation of our Poisson-Gauss function by expanding the logarithmic term, with which we have obtained an analytical approximation to the conic constant for any iso-indicial surface of the lens shown in Eq. (18). An approximation of this type is not reported in the Kasprzak model [43] since this model does not possess intrinsically a conic description. The conic constant associated to our model is presented in Eq. (19). With this we have obtained positive values of the conic constant for our model. It is important to note that we have constructed our lens model using parameters with in vivo high experimental accuracy, such as axial lens thickness and axial radii of curvature [59]. It is also important to acknowledge existing relationships between the optical and biometric properties of the lens [60], particularly, between the radius of curvature and the asphericity of the lens [61]. It has been demonstrated that this a feature of surface measurements when fitting to conics [62]. Taking these correlations into account would lead to creation of lens models with greater precision. It is worth noting that whilst the equatorial region of the lens is not directly relevant to vision, recent literature describing imaging methods that can show the entire lens shape [63] offer the prospect of further fine refinement of the model for future applications for which the equatorial region is of interest.
The parameters z a , z p and R of our model were obtained as a result of reproducing the in vivo measurements of the radii of curvature and the thickness of a 26-year-old lens. Nonetheless, there are in vitro age-dependent measurements of z a , z p and R [20], that, by performing scaling of approximately 18% for each of our parameters, can be reproduced with excellent approximation.
This result notwithstanding, the ultimate value of a model is its relevance and adaptability to the individual eye, thereby allowing for more personalized predictions. The shape parameters of the Poisson-Gauss lens are designed to be fitted to experimental data from individual eyes and to identify the ideal GRIN variation for that particular eye, reproducing the varying aberrations reported in the eye [64] and in the lens [65].

Conclusions
In this work we have proposed and demonstrated that the human crystalline lens can be modelled with a single function that can, simultaneously, describe the behaviour of the GRIN distribution and shape of its surfaces, and with it we built an schematic eye. By varying one single parameter, that can be related to the ciliary body function in providing the force needed for acccommodation, our model allows the reproduction of in vivo biometrical experimental data of relaxed and accommodated states of a human lens and at the same time can find any intermediate state as well as extrapolated accommodative states capable of focusing on the retina. The accommodative imaging abilities of our model have been tested by studying ray propagation through a schematic eye constructed with a conic cornea and the lens model presented in this work, showing image formation with aberrations within experimental range.
The proposed schematic eye can be used in investigations that consider variations in refractive error and age, incorporating different ocular biometries and extending to pathological conditions such as keratoconus. This will enable significant progress to be made towards the understanding of whole eye imaging and ultimately in the design of biologically relevant intraocular lenses.
Γ denotes the Gamma function and 1 F 1 the confluent hypergeometric function [67]. In Fig. 9 it is shown the behavior of h in Eq. (5) as function of m through the expressions for σ z and µ z above. The values of h, corresponding to each accommodative state shown in Table 3 The equation of the external iso-indicial contour of the lens can be written as y(z) = a r ln (1/h) + m ln (bz) − bz − z 2 /a 2 z , where h is given by Eq. (5). Now consider the inverse function z(y). Given that we seek the axial radius of curvature, each derivative in Eq. (29) must be evaluated at y = 0, point in which z has a local maximum as function of y, that is, dz/dy| y=0 = 0. Therefore Since z(y) cannot be inverted analytically we employ the inverse function theorem to write