## Abstract

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 different 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.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. 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–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–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–5,18–21]. Optical coherence tomography (OCT) is a leading modality in imaging ocular parameters [22–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–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–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–36,40–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–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.

## 2. Construction of the Poisson-Gauss function as a dynamic 3D model for the human crystalline lens

#### 2.1 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\in \mathbb {N}$, where $\mathbb {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=\sqrt {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 $\sigma$ 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 $\sigma _z$ from an average value located at $z = \mu _z$, $y = \mu _y = 0$ as shown in Appendix A. Both $\mu _z$ and $\sigma _z$ are functions of the parameters of the PG function, that is, $\mu _z(m, b, a_z)$ and $\sigma _z=\sigma _z(m, b, a_z)$. It should be noted that neither $\mu _z$ nor $\sigma _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}=\mu _y=0$, $z=\mu _z+\sigma _z$, namely The corresponding intersection results in a level curve given by $PG(y, z)=h$, that can be rewritten as $\frac {z^2}{a_z^2}+\frac {{y}^2}{a_r^2}+bz-m\ln (bz)=\ln {(1/h)}$. 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 $\mathcal {D}_0$ of the lens in the $z-y$ plane is limited by Eq. (6) and is defined byThe 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 $\mathcal {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.#### 2.2 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=\mu _z + \sigma _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 $\mu _z$ and $\sigma _z$ are calculated as described in Apendix A.Taking the above equations into account, the 3D GRIN profile of the crystalline lens can be described as:

for $z>0$, $r\geq 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 $\mathcal {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.## 3. 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 $\mathcal {D}$ defined by Eq. (9)).

#### 3.1 Volume of the lens

The volume $V$ of the 3D domain $\mathcal {D}$ is the volume of the solid of revolution obtained by rotating the function $y(z;h)$ given in Eq. (10):

Therefore the last equation can be written as $V=a_r^2\pi 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.#### 3.2 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 $\mathcal {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$, $\mathcal {R}(z)$ is always finite. The radii of curvature of the anterior and posterior surfaces of the lens can be written as $\mathcal {R}_a=\mathcal {R}(z_e-z_a)$ and $\mathcal {R}_p=\mathcal {R}(z_e+z_p)$, respectively. Hence, the gradient of the curvature radius $d\mathcal {R}(z)/dz$, of the anterior part of the lens is $-(a_r/a_z)^2-ma_r^2/(2z^2)$, and for the posterior part is $(a_r/a_z)^2+ma_r^2/(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 $\mathcal {R}_a$ and posterior $\mathcal {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 $\mathcal {R}_a$, $\mathcal {R}_p$ and $d$, very close to the experimental error tolerance; $\mathcal {R}_a$ and $\mathcal {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 $\mathcal {R}_a$, $\mathcal {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$, determined by the variance, Eq. (5). This, and the fact that the volume $\mathcal {D}$ of the lens remains constant by means of Eq. (15), ensures the optical integrity of the GRIN structure during the accommodative process.

On the other hand, a conic representation can be described by the expression

were $\mathcal {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## 4. 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 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)).

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.

Changes in the GRIN profile with accommodation are shown in 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_4^0$ 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_4^0$ is related to $W_{4,0}$ through $Z_4^0=W_{4,0}/(6\sqrt {5})$. Therefore, for a pupil diameter of 6 mm, we obtain $Z_4^0=-0.0014\pm 0.0008\,\mu$m, for the unaccommodated state, whilst $Z_4^0=-0.0055\pm 0.0007\,\mu$m, for the 4.03 D accommodated state. These values lie within the experimental data reported [56–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.

## 5. 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 $\mathcal {R}_a$, $\mathcal {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 $\mathcal {R}_a=12.48\pm 0.20$ mm, $\mathcal {R}_p=7.25\pm 0.25$ mm and $d=3.18\pm 0.02$ mm, while for another subject $\mathcal {R}_a=11.43\pm 0.04$ mm, $\mathcal {R}_p=6.12\pm 0.15$ mm and $d=3.36\pm 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\pm 0.23$ mm and an increase of the axial thickness of $0.34\pm 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].

## 6. 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.

## Appendix A: The variance of the Poisson-Gauss function

In this work we have used the concept of variance to measure the spatial extent of the Poisson-Gauss function. In statistics and probability theory the variance is the squared deviation of a probability distribution from its expected value [48]. In other words, the variance measures the spread of a probability distribution.

However, this statistical notion can be interpreted geometrically, since the spread of a function is a measure of how much a function is spread or extended in a particular direction. For instance, in optics the variance has been used to measure the spot size of laser optical beams [50]. This spot size is defined as the transverse illuminated area of a laser optical beam, that is, the transverse area in which most of its energy lies. The square root of twice the variance was used to measure the radius of the spot size of a cylindrical laser optical beam [66]. Consider the Gaussian $e^{-x^2/a^2}$. Its variance is $a^2/2$. The square root of twice the variance of the above Gaussian is therefore $a$, and at the same time $a$ is sometimes considered as the “width” of the Gaussian. Hence, we can say that the square root of twice the variance measures the “width” of a function, so we use the same notion to find the “width” of the Poisson-Gauss function in the radial and longitudinal directions.

Twice the variance of a function of two variables $f(y,z)$, whose domain is $z>0$ and $-\infty <y<\infty$, can be written for $z$ as

## Appendix B: Radius of curvature of the Poisson-Gauss lens

The radius of curvature $R$ of a function $y(z)$ is written as

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, $\left .dz/dy\right |_{y=0}=0$. Therefore

Since $z(y)$ cannot be inverted analytically we employ the inverse function theorem to write After some algebra it is obtained where## Funding

CHRISTUS-LATAM HUB CENTER OF EXCELLENCE AND INNOVATION, S.C..

## Acknowledgments

We thank CHRISTUS-LATAM HUB CENTER OF EXCELLENCE AND INNOVATION, S.C. for supporting the open access article publication charge. The first author acknowledges support from INAOE, Mexico.

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## References

**1. **S. Nakao, T. Ono, R. Nagata, and K. Iwata, “Model of refractive indices in the human crystalline lens,” Jpn. J. Clin. Ophthalmol. **23**, 903–906 (1969).

**2. **B. K. Pierscionek and D. Y. C. Chan, “Refractive index gradient of human lenses,” Optom. Vis. Sci. **66**(12), 822–829 (1989). [CrossRef]

**3. **C. E. Jones, D. A. Atchison, R. Meder, and J. M. Pope, “Refractive index distribution and optical properties of the isolated human lens measured using magnetic resonance imaging (MRI),” Vision Res. **45**(18), 2352–2366 (2005). [CrossRef]

**4. **S. R. Uhlhorn, D. Borja, F. Manns, and J. M. Parel, “Refractive index measurement of the isolated crystalline lens using optical coherence tomography,” Vision Res. **48**(27), 2732–2738 (2008). [CrossRef]

**5. **S. Ortiz, P. Pérez-Merino, E. Gambra, A. de Castro, and S. Marcos, “In vivo human crystalline lens topography,” Biomed. Opt. Express **3**(10), 2471–2488 (2012). [CrossRef]

**6. **J. J. Esteve-Taboada, R. Montés-Micó, and T. Ferrer-Blasco, “Schematic eye models to mimic the behavior of the accommodating human eye,” J. Cataract Refractive Surg. **44**(5), 627–641 (2018). [CrossRef]

**7. **B. K. Pierscionek and J. W. Regini, “The gradient index lens of the eye: an opto-biological synchrony,” Prog. Retinal Eye Res. **31**(4), 332–349 (2012). [CrossRef]

**8. **G. Smith and D. A. Atchison, * The Eye and Visual Optical Instruments* (Cambridge University Press, 1997).

**9. **D. A. Atchison and G. Smith, * Optics of the Human Eye* (Butterworth-Heinemann, 2000).

**10. **H. H. Emsley, * Visual Optics* (Butterworth, 1952).

**11. **Y. Le Grand and S. G. El Hage, * Physiological Optics* (Springer-Verlag, 1980).

**12. **L. Matthiesen, “Ueber Begriff und Answerthung des sogenannten Totalindex der Krystalllinse,” Pfluegers Arch. **36**(1), 72–100 (1885). [CrossRef]

**13. **H. Helmholtz, * Helmholtz’s Treatise on Physiological Optics* (Dover, 1962), Vol. 1. Appendix II.

**14. **D. A. Atchison and L. N. Thibos, “Optical models of the human eye,” Clin. Exp. Optom. **99**(2), 99–106 (2016). [CrossRef]

**15. **Y. Wu, A. Liu, H. Lv, X. Yi, Q. Li, X. Wang, Y. Ding, and J. Tong, “Finite Schematic Eye Model with Maxwell Fish-eye Spherical lens,” in * 2010 Symposium on Photonics and Optoelectronics* (IEEE, 2010), pp. 1–4.

**16. **R. G. Zainullin, A. B. Kravtsov, and E. P. Shaitor, “The crystalline lens as a Luneburg lens,” Biofizica **19**, 913–915 (1974).

**17. **T. Liu and L. N. Thibos, “Customized models of ocular aberrations across the visual field during accommodation.,” J. Vis. **19**(9), 13 (2019). [CrossRef]

**18. **B. A. Moffat, D. A. Atchison, and J. M. Pope, “Age-related changes in refractive index distribution and power of the human lens as measured by magnetic resonance micro-imaging in vitro,” Vision Res. **42**(13), 1683–1693 (2002). [CrossRef]

**19. **R. A. Schachar, “Growth patterns of fresh human crystalline lenses measured by in vitro photographic biometry,” J. Anat. **206**, 575–580 (2005). [CrossRef]

**20. **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]

**21. **Y. Shao, A. Tao, H. Jiang, M. Shen, J. Zhong, F. Lu, and J. Wang, “Simultaneous real-time imaging of the ocular anterior segment including the ciliary muscle during accommodation,” Biomed. Opt. Express **4**(3), 466–480 (2013). [CrossRef]

**22. **M. Shen, L. Cui, M. Li, D. Zhu, M. R. Wang, and J. Wang, “Extended scan depth optical coherence tomography for evaluating ocular surface shape,” J. Biomed. Opt. **16**(5), 056007 (2011). [CrossRef]

**23. **C. Du, D. Zhu, M. Shen, M. Li, M. R. Wang, and J. Wang, “Novel optical coherence tomography for imaging the entire anterior segment of the eye,” Invest. Ophthalmol. Visual Sci. **52**(2), 987 (2011). [CrossRef]

**24. **L. A. Lossing, L. T. Sinnott, C. Y. Kao, K. Richdale, and M. D. Bailey, “Measuring changes in ciliary muscle thickness with accommodation in young adults,” Optom. Vis. Sci. **89**(5), 719–726 (2012). [CrossRef]

**25. **S. Ortiz, D. Siedlecki, I. Grulkowski, L. Remon, D. Pascual, M. Wojtkowski, and S. Marcos, “Optical distortion correction in optical coherence tomography for quantitative ocular anterior segment by three-dimensional imaging,” Opt. Express **18**(3), 2782–2796 (2010). [CrossRef]

**26. **D. Siedlecki, A. de Castro, E. Gambra, S. Ortiz, D. Borja, S. Uhlhorn, F. Manns, S. Marcos, and J. M. Parel, “Distortion correction of OCT images of the crystalline lens: gradient index approach,” Optom. Vis. Sci. **89**(5), E709–E718 (2012). [CrossRef]

**27. **A. de Castro, J. Birkenfeld, B. M. Heilman, M. Ruggeri, E. Arrieta, J. M. Parel, F. Manns, and S. Marcos, “Off-axis optical coherence tomography imaging of the crystalline lens to reconstruct the gradient refractive index using optical methods,” Biomed. Opt. Express **10**(7), 3622–3634 (2019). [CrossRef]

**28. **J. Yao, J. Huang, P. Meemon, M. Ponting, and J. P. Rolland, “Simultaneous estimation of thickness and refractive index of layered gradient refractive index optics using a hybrid confocal-scan swept-source optical coherence tomography system,” Opt. Express **23**(23), 30149–30164 (2015). [CrossRef]

**29. **H. L. Liou and N. A. Brennan, “Anatomically accurate, finite model eye for optical modeling,” J. Opt. Soc. Am. A **14**(8), 1684–1695 (1997). [CrossRef]

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

**31. **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]

**32. **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]

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

**34. **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**(5), 1649–1663 (2014). [CrossRef]

**35. **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]

**36. **J. E. Gómez-Correa, S. E. Balderas-Mata, B. K. Pierscionek, and S. Chávez-Cerda, “Composite modified Luneburg model of human eye lens,” Opt. Lett. **40**(17), 3990–3993 (2015). [CrossRef]

**37. **J. E. Gómez-Correa, V. Coello, A. Garza-Rivera, N. P. Puente, and S. Chávez-Cerda, “Three-dimensional ray tracing in spherical and elliptical generalized Luneburg lenses for application in the human eye lens,” Appl. Opt. **55**(8), 2002–2010 (2016). [CrossRef]

**38. **L. Moser, “Über das Auge,” Dove’s Rep. Phys. **5**, 337–349 (1844).

**39. **H. Helmholtz, * Helmholtz’s Treatise on Physiological Optics* (Dover, 1962), Vol. 1.

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

**41. **R. Navarro, J. Santamaría, and J. Bescos, “Accommodation-dependent model of the human eye with aspherics,” J. Opt. Soc. Am. A **2**(8), 1273–1281 (1985). [CrossRef]

**42. **Y. Huang and D. T. Moore, “Human eye modeling using a single equation of gradient index crystalline lens for relaxed and accommodated states,” Proc. SPIE **6342**, 634201 (2006). [CrossRef]

**43. **H. T. Kasprzak, “New approximation for the whole profile of the human crystalline lens,” Oph. Phys. Optics **20**(1), 31–43 (2000). [CrossRef]

**44. **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**(1), 250–261 (2008). [CrossRef]

**45. **A. Popiolek-Masajada and H. Kasprzak, “Model of the optical system of the human eye during accommodation,” Oph. Phys. Optics **22**(3), 201–208 (2002). [CrossRef]

**46. **A. Popiolek-Masajada and H. T. Kasprzak, “A new schematic eye model incorporating accommodation,” Optom. Vis. Sci. **76**(10), 720–727 (1999). [CrossRef]

**47. **A. Khan, J. M. Pope, P. K. Verkicharla, M. Suheimat, and D. A. Atchison, “Change in human lens dimensions, lens refractive index distribution and ciliary body ring diameter with accommodation,” Biomed. Opt. Express **9**(3), 1272–1282 (2018). [CrossRef]

**48. **T. Yamane, * Statistics An Introductory Analysis* (Harper & Row, 1967).

**49. **F. A. Haight, * Handbook of the Poisson Distribution* (John Wiley, 1967).

**50. **W. H. Carter, “Spot size and divergence for Hermite Gaussian beams of any order,” Appl. Opt. **19**(7), 1027–1029 (1980). [CrossRef]

**51. **G. Smith, P. Bedggood, R. Ashman, M. Daaboul, and A. Metha, “Exploring ocular aberrations with a schematic human eye model,” Optom. Vis. Sci. **85**(5), 330–340 (2008). [CrossRef]

**52. **E. A. Hermans, P. J. Pouwels, M. Dubbelman, J. P. Kuijer, R. G. van der Heijde, and R. M. Heethaar, “Constant volume of the human lens and decrease in surface area of the capsular bag during accommodation: an MRI and Scheimpflug study,” Invest. Ophthalmol. Visual Sci. **50**(1), 281–289 (2009). [CrossRef]

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

**54. **D. Siedlecki, H. Kasprzak, and B. K. Pierscionek, “Schematic eye with a gradient-index lens and aspheric surfaces,” Opt. Lett. **29**(11), 1197–1199 (2004). [CrossRef]

**55. **M. Guillon, D. P. Lydon, and C. Wilson, “Corneal topography: a clinical model,” Oph. Phys. Optics **6**(1), 47–56 (1986). [CrossRef]

**56. **M. T. Sheehan, A. V. Goncharov, V. M. O’Dwyer, V. Toal, and C. Dainty, “Population study of the variation in monochromatic aberrations of the normal human eye over the central visual field,” Opt. Express **15**(12), 7367–7380 (2007). [CrossRef]

**57. **H. Cheng, J. K. Barnett, A. S. Vilupuru, J. D. Marsack, S. Kasthurirangan, R. A. Applegate, and A. Roorda, “A population study on changes in wave aberrations with accomodation,” J. Vis. **4**(8), 272–280 (2004). [CrossRef]

**58. **E. Tepichín-Rodríguez, A. S. Cruz Felix, E. López-Olazagasti, and S. Balderas-Mata, “Emmetropic eyes: objective performance and clinical reference,” Proc. SPIE **8785**, 878551 (2013). [CrossRef]

**59. **P. Artal, ed., * Handbook of Visual Optics: Fundamentals and Eye Optics* (CRC Press, 2017).

**60. **A. S. Vilupuru and A. Glasser, “Optical and biometric relationships of the isolated pig crystalline lens,” Oph. Phys. Optics **21**(4), 296–311 (2001). [CrossRef]

**61. **F. Manns, V. Fernandez, S. Zipper, S. Sandadi, M. Hamaoui, A. Ho, and J. M. Parel, “Radius of curvature and asphericity of the anterior and posterior surface of human cadaver crystalline lenses,” Exp. Eye Res. **78**(1), 39–51 (2004). [CrossRef]

**62. **A. Perez-Escudero, C. Dorronsoro, and S. Marcos, “Correlation between radius and asphericity in surfaces fitted by conics,” J. Opt. Soc. Am. A **27**(7), 1541–1548 (2010). [CrossRef]

**63. **E. Martinez-Enriquez, M. Sun, M. Velasco-Ocana, J. Birkenfeld, P. Pérez-Merino, and S. Marcos, “Optical coherence tomography based estimates of crystalline lens volume, equatorial diameter, and plane position,” Invest. Ophthalmol. Visual Sci. **57**(9), OCT600 (2016). [CrossRef]

**64. **R. Navarro, L. González, and J. L. Hernández-Matamoros, “On the prediction of optical aberrations by personalized eye models,” Optom. Vis. Sci. **83**(6), 371–381 (2006). [CrossRef]

**65. **A. de Castro, J. Birkenfeld, B. Maceo, F. Manns, E. Arrieta, J. M. Parel, and S. Marcos, “Influence of shape and gradient refractive index in the accommodative changes of spherical aberration in nonhuman primate crystalline lenses,” Invest. Ophthalmol. Visual Sci. **54**(9), 6197–6207 (2013). [CrossRef]

**66. **R. L. Phillips and L. C. Andrews, “Spot size and divergence for Laguerre Gaussian beams of any order,” Appl. Opt. **22**(5), 643–644 (1983). [CrossRef]

**67. **M. Abramowitz and I. A. Stegun, * Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables* (National Bureau of Standards, 1972).