## Abstract

We present in this work a numerical model for characterizing the scattering properties of the human lens. After analyzing the scattering properties of two main scattering particles actually described in the literature through FEM (finite element method) simulations, we have modified a Monte Carlo’s bulk scattering algorithm for computing ray scattering in non-sequential ray tracing. We have implemented this ray scattering algorithm in a layered model of the human lens in order to calculate the scattering properties of the whole lens. We have tested our algorithm by simulating the classic experiment carried out by Van der Berg *et al* for measuring “in vitro” the angular distribution of forward scattered light by the human lens. The results show the ability of our model to simulate accurately the scattering properties of the human lens.

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

## 1. Introduction

The human lens is a remarkable biological structure. Its mission is to focus the light refracted at the cornea into the retina. The human lens (sometimes we will simply refer to it as the *lens* for the sake of clarity) is a layered structure formed by cells which contain an intracellular medium composed by a solution of proteins (mainly of the crystallin protein family) in water [1]. The concentration of the proteins is higher in the lens core than in the cortex. This explains why the refractive index is greater at the core than at the periphery, as it depends on the concentration of crystallin proteins in water. Therefore, the human lens acts as a gradient index (GRIN) lens. This fact has been demonstrated both in animals [2] and humans [3]. The fiber cells provide the mechanical properties of the lens, which is rigid enough to maintain its biconvex lens-like shape, but also flexible to allow accommodation. On the other hand, the intracellular media provides the lens transparency and refractive index inhomogeneity of the human eye lens [1]. It is important to point out that, contrary to what happens with other tissues, the normal human lens does not present structures that might produce a high amount of scattered light, such as blood vessels, cell nuclei, organelles, or other potential scattering structures. Moreover, in normal conditions, there is an important index matching between the cell walls and the intracellular medium, particularly at the lens core, which reduces light scattering and/or diffraction at these boundaries [1].

Despite the mechanisms described before, some scattering is always present in the human lens. In normal eyes, scattering only manifests itself when the eye is illuminated with high amounts of energy. Even for normal eyes, ocular scattering gives rise to the so-called “straylight” defined by Van der Berg [4] which is related to the peripheral part of the ocular PSF and to the physiological phenomenon of glare [4]. The amount of light scattered by the eye lens is greater for aging individuals, particularly when cataracts are present. In this case, the individual will experience some of the different symptoms associated with this condition, such as loss of visual acuity and contrast, more frequent appearing of glare, etc. being all of them related to an increasing light scattering within the lens.

Light scattering at the human lens has been thoroughly studied both theoretical- and experimentally by many authors. Bettelheim *et al* [5] measured “in vitro” the angular distribution of scattered light by lens sections (of around 10 μm thick and perpendicular to the optical axis) of different donors. They fitted their results to a theoretical model based on random fluctuations of density and orientation of scattering particles, and they concluded that the lens scattering was due to a mixture of two scattering particles: spherical protein clusters with radius comprised between 130 and 450 nm and cylindrical cytoskeletal bodies with a length comprised between 600 and 2600 nm [5]. In any case, a strong dependence between the model parameters and the spatial location of the section within the human lens was observed. Subsequent works expanded this model by studying the dependence with age [6], and eyes presenting different kinds and degrees of cataracts [7]. Notice that all of these experiments have been performed “ex vivo” that is by extracting the lens after the decease of the donor, and we have to keep in mind that scattering properties of extracted lenses (“ex vivo”) may not be the same as that of lenses in a living eye (“in vivo”).

Van der Berg and Ijspeert measured “in vitro” the light scattering of whole lenses [8] from several donors of different ages, some of them clear and some of them presenting cataracts. In most cases, the fitting of the resulting scattered intensity against the polar angle θ shown an inverse squared dependence. In a later work [9], Van der Berg and Spekreijse were able to measure “in vitro” the light scattered by different parts of the human lens using slit illumination, avoiding, in this way, the need for slicing the lens. The result of this analysis was a phenomenological model which explains the lens scattering through the contribution of three sources. The first source is the Rayleigh scattering produced by small particles with a radius of around 10 nm, closely related to that of the α-crystallin protein aggregates [9]. The second source is the non-Rayleigh scattering produced by relatively big (700 nm radius) “effective” particles. Finally, the third source described is the surface scattering at the zones of discontinuity (“Wasserspalten”) corresponding to the different layers of the human lens [4,9]. It is noteworthy that the phenomenological model presented by Van der Berg and Spekreijse shows a spatial dependence due to the inhomogeneity of the lens.

The exact nature of the large scattering particles within the human lens as multilamellar bodies was reported by Gilliland *et al* [10] and thoroughly analyzed by Costello *et al* [11,12]. These multilamellar bodies, or vesicles, as we will call them in the following, are structures formed by a lipid shell that surrounds a body with the same composition of water and proteins than that of the lens intracellular medium. The mean size of these particles is around 2.13 microns in diameter. As it is shown by Costello *et al* [11] and Mendez-Aguilar *et al* [13], the size of these particles, together with the index jump between the vesicle shell, the inner body, and the external medium makes these vesicles highly efficient scatterers with a relatively high cross-section. In all of these works, the concentration of the vesicles within the human lens is deemed as constant for a given eye, but it can vary between eyes. In particular, cataractous eyes present a greater vesicle concentration than normal ones. On its part, the size of vesicles is distributed randomly along the lens following a Gaussian distribution with a mean diameter of 2.13 μm and a standard deviation of 0.64 μm [11].

Regarding the nature of the small size scatterers, we will not take into account in our model the soluble alpha-crystallin aggregates reported by Van den Berg and Spekreise [9]. Due to the small size of these aggregates and their close packing (particularly at the lens core), we will not consider them as independent scatterers, and we will suppose that their contribution is to raise the refractive index of the lens [14]. Instead, we will focus our attention on the insoluble high molecular weight (HMW) aggregates that are formed by gamma-crystallin denaturalization [15,16]. These aggregates were first reported in the 1970s and they have been thoroughly studied ever since [17,18]. In this context, it was Benedek [19] who first predicted theoretically that light scattering in cataractous eyes could be due to protein aggregates with a molecular weight of around $50\times {10}^{6}$ g/mol (with a corresponding estimated diameter of 50 nm). This prediction was confirmed experimentally by several authors [20,21] using different techniques. Nowadays there is a good knowledge of the biochemical processes that lead to the formation of HMW aggregates [22,23], basically through the inhibition of the chaperone activity of the alpha-crystallin and denaturalization of gamma-crystallin which becomes susceptible to aggregate forming structures such as amyloid fibers [24] or protein condensates [25]. These structures have been measured with different techniques such as electron microscopy or magnetic resonance [24,25]. From these measurements, diameters ranging from 20 nm to 100 nm [24,25] have been reported. An additional cause of light scattering in the eye are the zones of discontinuity but their effect is greater for the back-scattered light, so it will not be considered in our model.

In this work, we have developed a model for the numerical computation of the scattering of the eye lens. According to the literature, we have considered that light scattering in the lens is produced by a mixture of small (clusters) and large (vesicles) particles and we have performed a comprehensive study of the scattering properties of each kind of particle. We have also taken into account the inhomogeneous nature of the human lens which is manifested in the spatial variation of the protein concentration within the lens. As a consequence of this variation of the protein concentration, both the refractive index and some scattering parameters also vary spatially. In order to test our model, we have made a simulation of the human lens scattering experiment described in reference [8] and we have performed a non-linear fitting with the experimental results described in this paper in order to find the optimum values of the parameters of our model for both a normal and a cataractous eye.

The paper is organized as follows: in Section 2 we describe the basis of our model, including a thorough FEM electromagnetic study of the scattering properties of clusters and vesicles, and the implementation of a mixture model of those scattering particles (clusters and vesicles) for computing the bulk scattering of a ray using the Monte-Carlo technique. In Section 3, we give the results of the numerical simulation of the experiment reported in reference [8], in order to show the ability of our model to describe properly the light scattering by the human lens. Finally, conclusions are drawn to end the paper.

## 2. Inhomogeneous scattering model of the human eye lens

#### 2.1 Theoretical foundations

A healthy human lens acts like a high transmittance GRIN lens presenting a spatial variation of the refractive index. This spatial variation is due to the different concentration of protein within the human lens [1]. We can establish a relationship between this variable protein concentration and the refractive index of the eye, *n*, using a simple mixing model

Typically, a mix of heterogeneous particles is found across the intracellular medium. If the size of these particles remains small and their concentration low, the resulting scattering is negligible [26,27]. However, some factors, such as age, UV radiation, and local thermal variations, trigger the increase in particle size and number, generating a notable enhancement of scattering effects [5].

The scattering of a particle is strongly related to its geometry, the light wavelength, and the relationship between the refractive index of the particle and its surrounding medium. In this study, we neglect the absorption of the scatterers within the lens. Within this approach, two scattering magnitudes become relevant when analyzing the light-particle interaction. In our model, these magnitudes depend on the local concentration of protein within the lens. The first one is the scattering cross-section,${\sigma}_{s}\left(c\right)$, which describes the probability of incident electromagnetic radiation being scattered by a particle, calculated as [26,28–34]

**n**the normal vector, and

*I*

_{0}is the incoming irradiance. The integration domain corresponds to the area of an imaginary sphere which surrounds completely the scattering particle and the normal vector points outwards of this sphere. The second parameter of interest, the coefficient of asymmetry, describes the anisotropy of the scattered radiation, and it is computed as [32]

**n**the normal vector at the same point, and

*θ*the polar angle defined as the angle formed by the direction vector of the incoming radiation and the normal vector

**n**. In our study, we have used the far-field pattern for computing both the cross-section and the asymmetry coefficient according to [28,32,36].

We assume that the concentration of the scattering particles is low enough to disregard interference effects between scatterers. Moreover, as a first approximation, we also consider a medium made of identical particles (we will consider later the mixture of two scattering particles). Then, the scattering coefficient is given as

being*N*the volumetric concentration (particles per unit of volume) of the scattering particles. We can define a scattering mean-free-path, ${l}_{s}\left(c\right)={\mu}_{s}{\left(c\right)}^{-1}$, as the distance traveled by light between two scattering events. In principle, the scattering phenomenon is described by the scattering coefficient and mean free path. However, a non-absorbing turbid media with particles presenting high asymmetry, such as the human lens, is better described [31] by defining the reduced scattering coefficient: $\mu {\text{'}}_{s}\left(c\right)={\mu}_{s}\left(c\right)\left[1-g\left(c\right)\right]$ and the effective mean-free-path as: $l{\text{'}}_{s}\left(c\right)=\mu {\text{'}}_{s}{\left(c\right)}^{-1}$.

Turbid media, particularly, biological tissues, show scattering probability distributions strongly related to *g*. The Henyey-Greenstein scattering phase function includes the anisotropy parameter to calculate the phase function,$p\left(\theta \right)$as [27,31]

*θ*is the polar angle. Notice that the phase function depends indirectly on the concentration through the asymmetry coefficient

*g*. Next, we will show the results obtained after computing the scattering parameters of different scatterers present in the human lens using finite element methods (FEM).

#### 2.2 Scatterers in the human lens

All FEM simulations considered in this paper have been evaluated using COMSOL^{TM} Multiphysics 5.2. In every simulation, we use electromagnetic radiation with a wavelength of 550 nm, according to the maximum sensitivity of the human eye. Also, we have taken in all cases the amplitude of the incident electric field as ${E}_{0}=1\text{\hspace{0.17em}}\text{V}/\text{m}$, equivalent to ${I}_{0}=1.33\times {10}^{-4}\text{\hspace{0.17em}}\text{W}/{\text{cm}}^{\text{2}}$.

### Protein clusters

We can consider the crystallin protein as a spherical particle with a radius of about 2 nm. When these proteins agglomerate, they form protein clusters that can be considered as a greater particle which radius up to 10 nm for normal eyes [24]. According to the optimal close sphere packing fraction whose value is around 75%, these new dispersers are not composed of pure proteins. Thus, the refractive index of these clusters decreases until a value of ${n}_{c}\approx 1.55$from the value of${n}_{p}=1.62$which presents the bulk protein. For normal eyes, the scattering effect of these clusters is small enough to not impair the transparency of the human lens [34]. However, in certain conditions (lens aging, temperature increment, etc.), the cluster radius increases up to 20 nm [24] or greater (around 50 nm) [25,35]. In this case, the radiation scattered by these clusters rises significantly according to Mie theory.

Figure 1 shows the near- and far-field patterns of the scattered radiation, the scattering cross-section, and the asymmetry coefficient, as a function of the radius of the cluster and the protein concentration of the surrounding medium. We have set the maximum value of the radius at 200 nm, in order to consider very large clusters, although, according to the sizes reported in references [24,25] the cluster radius varies between 10 and 50 nm. Protein cluster shows a typical dipolar near field pattern for a particle with 25 nm radius. However, we get a micro-jet response for a cluster radius of 150 nm. The scattering cross-section increases when the radius becomes larger as shown in Fig. 1(d), and it also depends on the protein concentration of the surrounding medium. When the protein concentration is low, the refractive index jump between the surrounding medium and the protein cluster increases, raising the value of the scattering cross-section. The anisotropy of protein clusters can be studied considering the far-field pattern of the scattered radiation. Figure 1(a), Fig. 1(b) and Fig. 1(c) show this far-field pattern (see insets). For low cluster radius, we see a dipolar behavior, being the anisotropy of cluster particles close to$g=0.$However, when the radius increases, the anisotropy grows up towards$g=\mathrm{0.8.}$On the other hand, anisotropy does not depend strongly with protein concentration, see Fig. 1(e). By interpolating the data depicted in Fig. 1 we can obtain the value of the scattering cross-section and asymmetry coefficient for any value of the protein concentration and cluster radius.

### Vesicles or multilamellar bodies

The fiber cells that structure the human lens show different shapes depending on when they were generated. Thus, the shape of a fiber cell grown up during the human fetal period shows a greater uniformity than a structure generated at a larger period of life [10]. Throughout the lens, breaks at the cells’ walls [11], generate a type of bulges known as multilamellar bodies or vesicles [11,13]. To model these vesicles, we consider a particle, with an average diameter of around 2.3 μm, having a lipid shell with a thickness of 50 nm that shows a refractive index around 1.55 [11,13]. Following Costello *et al* [11], we have considered that the refractive index of the inner medium of the vesicle is similar to the surrounding one but presenting a 6% increment in the refractive index. Figure 2 shows the modulus of the electric vector corresponding to the near- and far field patterns, the scattering cross-section, and the asymmetry coefficient for different protein concentration of the surrounding medium and particle size.

As shown in Fig. 2, vesicles display an electromagnetic near-field pattern that can be compared with that of a micro-lens. Due to their greater size, compared with clusters, the scattering cross-section of vesicles is larger than that of protein clusters, making them more efficient scatterers. Indeed, there is a strong dependence of the scattering parameters with the diameter as displayed in Fig. 2. Notice that for a vesicle with radius 1.4 μm, we get a scattering cross-section of around 1.75 × 10^{−11} m^{−2}. Dividing this cross-section by the area of the vesicle we get a scattering efficiency of 2.84, which is of the same magnitude order than that obtained by Costello et al [11]. However, as the ratio between the inner refractive index and that of the external medium is fixed, there is almost no dependence with the protein concentration as shown in Fig. 2(d).

The electromagnetic far-field pattern of the vesicles shows a strong forward scattering distribution -see red box insets of Figs. 2(a) and 2(b)-. This effect corresponds with an anisotropy value around$g=0.9$, a value that it is only obtained for large protein clusters (${r}_{c}=200$nm). The anisotropy grows when the protein concentration becomes larger.

As with the case of clusters, we can compute the values of the scattering cross-section and asymmetry coefficient by interpolating the data resulting from the FEM simulation.

#### 2.3 Monte Carlo’s bulk scattering algorithm for human lens particles

Once we have computed the properties of the different scattering particles within the human lens (clusters and vesicles) we can now proceed with the implementation of a computer bulk scattering algorithm. In this work, we have modified the Monte Carlo’s bulk scattering algorithm provided by Zemax OpticSudio (ZOS in the following), although our method can be easily extrapolated to other ray tracing programs or programming languages. More details on the ZOS bulk scattering algorithms can be found in reference [37], while its implementation in an ocular model is reported in reference [38]. Our departing point is the algorithm that works with the Henyey-Greenstein probability distribution (5), described in reference [39].

In order to modify this algorithm for the scattering particles of the human eye, we must take into account some considerations. First, we are dealing with an inhomogeneous medium whose refractive index depends on the local concentration of protein. Second, the scattering parameters of both clusters and vesicles depend on the external medium refractive index (through the protein concentration). Finally, the size (and volumetric concentration in the case of clusters) of both types of scattering particles varies also within the lens. In the case of clusters, we will model this variation in size and concentration after Siew *et al* [6], while for vesicles we will consider a random variation of size with a constant concentration following Costello *et al* [11].

As we will see with more detail in the next section, we have modeled the lens as a medium composed of several layers of different refractive index, but each of these layers is a homogeneous one. In this way, we may keep the same homogeneous Monte Carlo’s bulk scattering algorithm for all the layers of the lens.

The parameters of the bulk scattering algorithm are the cluster radius${r}_{c}$expressed in nm, the mean distance between clusters *d*, also in nanometers, and the vesicle concentration *N*_{v}, measured in vesicles per mm^{3}. As said before, each layer of the lens may have different values of these parameters in order to model the non-homogeneous scattering properties of the human lens. The bulk scattering algorithm can be described by the following steps:

- 1. For each ray passing through the bulk medium of the layer, the algorithm computes a random vesicle diameter using a Gaussian distribution with mean 2.13 μm and a standard deviation of 0.64 μm [11]. This is an approximation to the fact that the ray encounters vesicles with random sizes when it propagates through the lens. We also compute the volumetric cluster concentration as ${N}_{c}=\sqrt{2}\text{}\text{\hspace{0.05em}}{d}^{-3}$, being
*d*the mean distance between clusters. - 2. Next, the local protein concentration is determined from the value of the local refractive index and the simple mixture model
with the same parameters as discussed before, see Eq. (1). From the protein concentration and the size of the scattering particles, we compute the scattering coefficients${\sigma}_{s}$and

*g*for both clusters and vesicles using the results of the electromagnetic simulation depicted in Figs. 1 and 2. - 3. Afterward, we compute, for each scattering particle, the mean free path,$l{\text{'}}_{s}$, as the inverse of the reduced scattering coefficient$\mu {\text{'}}_{s}$ which, in turn, depends on the asymmetry coefficient,
*g*, the scattering cross-section, ${\sigma}_{s}$, and the volumetric concentration of the scatterers,*N,*see Eq. (4). - 4. Next, we have to manage the mixture of two different scattering particles in Monte Carlo’s algorithm. To do so, we have followed the procedure described in reference [39]. According to this work, when we have a mixture of two kinds of non-interacting scattering particles, clusters, and vesicles in our case, characterized by their free mean paths$l{\text{'}}_{sc}$ and $l{\text{'}}_{sv}$, respectively, the probability that the scattering particle is a cluster is ${p}_{c}=l{\text{'}}_{sc}/\left(l{\text{'}}_{sc}+l{\text{'}}_{sv}\right)$. Conversely, the probability that the scattering particle is a vesicle is ${p}_{v}\equiv 1-{p}_{c}=l{\text{'}}_{sv}/\left(l{\text{'}}_{sc}+l{\text{'}}_{sv}\right)$. In these conditions, for each ray that passes through the medium, we can determine whether a cluster or a vesicle might scatter the ray by extracting a random number $\xi \in \left[0,1\right]$. If $\xi \le {p}_{c}$then the ray might be scattered by a cluster, otherwise, it might be scattered by a vesicle.
- 5. Finally, the bulk scattering algorithm will determine whether or not the ray is actually scattered, and, if so, it computes the direction of the scattered ray using the Henyey-Greenstein probability distribution (5). To do this, the algorithm uses the values of the mean free path and the asymmetry coefficient of the scattering particle determined in the previous step. See references [27,37] for further details.
- 6. The procedure is repeated for each ray traced through the medium to get the final distribution of the scattered light. In our simulations, we have used a number of rays comprised between 10
^{7}and 10^{8}in order to obtain a significant number of scattered rays.

To implement this algorithm we have programmed a Dynamic-link Library (or DLL) in C language which can be used by the bulk scattering feature of ZOS.

#### 2.4 Modeling human lens as a layered medium

The final step of our model is the modeling of the eye lens as a layered GRIN medium. To do so we have used the capabilities of the non-sequential ray tracing mode of ZOS. In this mode it is possible to define nested volumes, so we can simulate a layered GRIN medium by defining a set of nested volumes each one with a different refractive index. The reader can find further information and examples in reference [40]. Obviously, in a GRIN layered medium defined as a set of nested volumes, the surfaces that delimit these volumes are iso-indical surfaces.

Therefore, we need a model for the refractive index distribution of the lens so we can find the iso-indical surfaces. In the literature, there are several refractive index models available and we have selected the model proposed by Navarro et al [41]. We have done so because this model takes into account the age as a parameter, which makes it very well suited for studying cataracts associated with age, and because of its great optical and morphological accuracy. In Navarro’s model, the refractive index of the human lens is given by the following function

*z*the axial coordinate, ${n}_{0}$ the maximum central index value and ${\delta}_{n}$, ${f}_{ant}$, ${\Delta}_{ant}$, ${a}_{ant}$, ${b}_{ant}$, ${\in}_{ant}$, ${t}_{ant}$, ${f}_{pos}$, ${\Delta}_{pos}$, ${a}_{pos}$, ${b}_{pos}$, ${\in}_{pos}$and $p$ are constants which depend on the age and accommodation. These constants can be computed from the data tabulated in Table 1 of reference [41]. In this work, we have used the constants corresponding to a young subject (24 year old) and an elderly one (80 year old) with no accommodation. In Eq. (7) $\left({z}_{ant},{\omega}_{ant}\right)$,$\left({z}_{i},{\omega}_{i}\right)$, and$\left({z}_{pos},{\omega}_{pos}\right)$are the coordinates that describes the anterior, intermediate and posterior surface of the lens (see details on reference [41]).

By imposing the condition ${n}_{ant}\left(z,\omega \right)-n=0$, where *n* is a given refractive index, we get an iso-indical surface located in the anterior portion of the lens. Similarly, the solution of the equation ${n}_{pos}\left(z,\omega \right)-n=0$ defines an iso-indical surface located at the rear of the lens. Solving both of these equations for different values of *n* results in a family of iso-indical surfaces$z=f\left(\omega ;n\right)$ that define the different layers of the human lens. According to Navarro’s model, these iso-indical surfaces are conicoids described by the general equation

*Q*the conic constant. All of these parameters depend on whether the surface belongs to the anterior or posterior zone of the lens and on the value of the refractive index [41].

In this work, we have modeled the human lens with 11 layers corresponding to a set of 11 equispaced values of the refractive index between the minimum and maximum refractive index of the lens, that are given by Eq. (7). Thus, each of these layers is defined as a lenticular volume formed by two conicoid surfaces described by Eq. (8), so our lens model consists of a set of nested aspheric lenses. In Fig. 3 we can see the layered model of the human lens for a 24 and an 80-year-old subject. Notice the significant differences in the lens morphology and refractive index distribution for these two ages.

One of the basic assumptions of our model is that each of the lenticular volumes that constitutes the lens, contains scatterers (clusters and vesicles) with different size and concentration. This is in agreement with the experiments of Siew *et al* [6] and Van der Berg and Speckrijse [9], which show that the scattering particles are not homogeneously distributed through the human lens so that their mean size and concentration vary along the axial direction *z*, particularly for clusters [6]. For vesicles, it is usually assumed a uniform concentration (which varies from eye to eye) along the eye lens [11] with a randomly distributed diameter. In our layered model of the human lens, we can change independently for each layer the parameters that define the size of clusters and the concentration of both clusters and vesicles (note that, as stated before, the vesicle size is selected randomly by the Monte Carlo algorithm), although we keep this later parameter equal for all the layers in accordance with [11].

In order to compute the cluster radius and distance for each layer, we have fitted the data provided by Siew *et al* [6], obtaining the following spatial dependence for the cluster radius

*z*is the axial distance measured from the apex of the human lens, and the remaining parameters are ${\beta}_{1}=-0.0059$,${\beta}_{2}=0.0865$and$w=2.74$ rad/mm. To implement this function in our model we select first a value for ${r}_{0}$ and then we compute the cluster radius for each layer by evaluating Eq. (9) at the locations given by the values of${z}_{0}$corresponding to the vertex of the iso-index surfaces that define the layered model of the lens.

We have proceeded in a similar way for the cluster distance *d*. Indeed we have selected this parameter for characterizing the cluster concentration because its spatial dependence is given by Siew *et al* [6]. After fitting the experimental data reported by these authors, we have found that the spatial dependence of the cluster distance is given by the following function

^{−1},${\alpha}_{2}=0.7$mm

^{−2}, and${\alpha}_{3}=-0.225$mm

^{−3}. Therefore by selecting a value for${d}_{0}$and substituting the values of${z}_{0}$, in Eq. (10) we get the cluster distance for each layer.

Using Eqs. (9) and (10) we have effectively reduced the parameters of the scattering particles in our model to three:${r}_{0}$, ${d}_{0}$and${N}_{v}$, being this latter parameter being the volumetric concentration of vesicles.

To summarize, we have modeled the eye lens as a set of nested volumes, each one of them an aspheric lens whose anterior and posterior surfaces are iso-indical surfaces of Navarro’s gradient index model. For each layer, the concentration and size of clusters vary according to Eqs. (9) and (10) but the vesicle concentration is constant throughout the lens. When a ray is traced within this structure, the modified Monte Carlo’s scattering algorithm determines, for each layer, first which particle (cluster or vesicle) might scatter the ray and then, whether the ray is scattered at all and, if so, the direction of the scattered ray. The ray is then propagated to the next layer where the process continues. Tracing a great number (10^{7}-10^{8}) rays through the system we may have a reliable depiction of the forward scattering of light by our lens model.

## 3. Numerical simulation

In order to test our model, we have performed a numerical simulation in ZOS based on the actual experimental set-up used by Van der Berg and Ijspeert [8] for measuring “in vitro” the angular distribution of the scattered light by a whole human lens. A drawing of the simulated set-up can be seen in Fig. 4(a). Following Van der Berg and Ijspeert, we have located the human lens within a cylindrical holder containing a solution with a refractive index close to that of the water. The holder is illuminated by a point source with a divergence of around${3}^{\circ}$radius. We have simulated the mobile detector used by Van der Berg *et al*, with a ZOS polar detector. In the non-sequential mode of ZOS, a detector is a surface, which may be divided in pixels, used to compute radiometric magnitudes such as flux, intensity, etc, by computing the energy of the rays that pass through this surface. A polar detector is then a semispherical surface which can compute the angular distribution of the light scattered by the lens, more specifically the radiant intensity. The curvature radius of the polar detector is 550 mm, the same as the distance between the lens holder and mobile photodetector used in Van der Berg and Ijspeert experiment [8]. A magnified image of the lens and the cylindrical vessel can be seen in Fig. 4(b). All the simulations have been computed by tracing 10^{8} rays in non-sequential mode using an AMD Ryzen Threadripper 3.4 GHz computer with 32 Gb of RAM memory and the version 16.5 of ZOS.

We have analyzed first the influence of the algorithm parameters (cluster radius, inter-cluster distance, and vesicle concentration) on the angular distribution of the scattered light given by the function${\mathrm{log}}_{10}I\left(\theta \right)$, where $I\left(\theta \right)$is the average radiant intensity measured by the polar detector. In order to compute the transmittance of the lens, we have repeated the simulation disabling the scattering calculation, so ZOS traces the rays without using our bulk scattering DLL. In these conditions, an estimation of the transmittance is the quotient [8]

In Fig. 5 we have plotted the angular distribution of the scattered light computed for three different concentrations of vesicles:${N}_{v}=500$mm^{−3},${N}_{v}=5000$mm^{−3} and${N}_{v}=50000$mm^{−3} (these concentrations should be read as, for example, 500 vesicles per mm^{3} but we have written 500 mm^{−3} for simplicity) for three different values of${r}_{0}$: 20, 60 and 100 nm, respectively. For all these plots, the cluster distance was ${d}_{0}=1200$ nm and we selected the 80 years old model. Vesicle concentration are similar to the ones given in Costello’s paper [11] for a normal (556 mm^{−3}) and cataractous (4071 mm^{−3}) lens and the third concentration (50000 mm^{−3}) would represent a highly cataractous lens.

In Fig. 5(a) we can see that, for a small cluster radius, all the scattering is due to the vesicles so the dispersion is greater for lower angles (which is characteristic of the scattering produced by particles whose asymmetry coefficient g is close to 1). When the cluster radius is 20 nm, the transmittance values obtained for the low and medium concentrations of vesicles are 98.6% and 85.6%, respectively. Costello *et al* [11] reported a relative scattering (fraction of the incident energy which is scattered by the lens) of 2.1% (vesicle concentration of 556 mm^{−3}) and a 14.6% (vesicle concentration 4071 mm^{−3}). Therefore the values of the transmittance we have obtained are compatible with the values of the relative scattering for similar vesicle concentrations. In Figs. 5(b) and 5(c) we show that, when the cluster radius becomes greater, light is scattered more evenly in all directions, particularly for small vesicle concentration. It is interesting the drop of transmittance that occurs for high vesicle concentration (50000 mm^{−3}) or for large cluster radius (100 nm), indicating that both mechanisms could explain the greater opacity of the cataractous lens. However, notice that, when we have, at the same time, high vesicle concentration and large cluster radius, the scattering is mainly due to vesicles as shown by the bell shape of the curve in Fig. 5(c). This is due because we have kept a low cluster concentration (${d}_{0}=1200$ nm) in order to accomplish our algorithm hypothesis and disregard interference effects between scatterers.

Finally, we have performed a non-linear least squares fitting of our model to the data obtained by Van der Berg and Ijspeert for a young (24-year-old) subject with clear lens and an elderly (80-year-old) subject with cataracts [8]. To fit these data we have used the layered models of these lenses depicted in Fig. 3 following the gradient index model of Navarro *et al* [41]. The fitting parameters are average cluster radius, ${r}_{0}$, inter cluster distance at $z=0$, ${d}_{0}$, vesicle concentration, ${N}_{v}$, and a normalization parameter, ${I}_{N}$, which has been introduced to take into account the fact that Van der Berg data were normalized [8]. We have employed a global optimization algorithm programed in Matlab who can control ZOS files through an API (application programming interface). Figure 6 shows the results obtained for both a normal and a cataractous eye.

The values of the fitting parameters are${r}_{0}=31.4$ nm,${d}_{0}=1560$ nm,${N}_{v}=300$mm^{−3}, and${I}_{N}=65$for the normal eye, and${r}_{0}=29$ nm,${d}_{0}=1040$ nm,${N}_{v}=1025$mm^{−3}, and${I}_{N}=31.6$for the cataractous one. The value of cluster diameter, for both normal and cataractous lens, is close to the 50 nm diameter estimated from the molecular weight ($50\times {10}^{6}\text{\hspace{0.17em}}\text{g}/\text{mol}$) of the aggregates given by Benedek [19] from a purely theoretical calculation, and it is compatible with the size of the amorphous and fibril aggregates shown in [25]. However, the cluster radius obtained is larger than the size of cluster reported by [24] from “in-vitro” UV induced protein aggregation. The difference between the clusters of normal and cataractous eyes is not their size (it is slightly larger for a normal eye), but their concentration, as clusters are more densely packed for the cataractous eye than for the normal one. Besides, the vesicle concentration for cataractous eye is three times that of the normal eye with values of the same order of magnitude of that of Costello *et al*. Thus, for a cataractous eye, we have found a concentration of 1025 mm^{−3} against 4071 mm^{−3} reported by Costello *et al* [11] and, for a normal eye, we get 300 mm^{−3} against 556 mm^{−3}. It is clear that in our fit, both clusters and vesicles contributed to the scattering and this may explain the reduction in vesicle concentration.

## 4. Conclusions

In this work, we present an inhomogeneous model for characterizing numerically the scattering properties of the human eye lens. Our model is based on three points: 1) a rigorous FEM simulation of the properties of particles (protein clusters and vesicles) with the actual dimensions and refractive index as reported by the literature, 2) a two-particle Monte Carlo bulk scattering ray tracing, with the proper modifications for taking into account the dependence of the scattering parameters with the local protein concentration, and 3) a layered model of the human lens which takes into account the inhomogeneous distribution of the particles within the human lens. Notice that our model considers two sources of inhomogeneity in the human lens: the different protein concentration (and, thus, refractive index) within the human lens, and the variation in size and number of particles along different zones of the lens.

We have tested our model by simulating numerically a classical experiment [8] in which the angular distribution of the radiant intensity scattered by the human lens was measured. By varying the model parameters: radius of protein cluster, distance between clusters, and vesicle concentration along the different layers of the human lens, we have been able to calculate the relative effect of these parameters in the amount of light scattered and the shape of the scattering distribution.

Finally, we have fitted our model to the experimental results reported in reference [8] showing a good agreement between them for both a normal and a cataractous eye. The values of the fitting parameters obtained are feasible and they show how the cataractous eye present a small increase in the cluster protein concentration while tripling the concentration of vesicles.

Given the characteristics of our model, it might be used as a workbench for analyzing the effect of a different configuration of the particles responsible for the ocular scattering in a realistic way. It also may be helpful in establishing new diagnostic techniques based on the scattering of light by the human lens particles.

## Funding

European Fund for Regional Development (EFRD-FEDER, EU), and the Spanish Government’s Agencia Estatal de Investigación (AEI) belonging to Ministerio de Economia y Competitividad (MINECO) (grant DPI2016-75272-R).

## Acknowledgements

We wish to thank Dr. Javier Alda for his useful comments and his help in the elaboration of this manuscript.

## Disclosures

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

## References

**1. **P. J. Donaldson, A. C. Grey, B. Maceo Heilman, J. C. Lim, and E. Vaghefi, “The physiological optics of the lens,” Prog. Retin. Eye Res. **56**, e1–e24 (2017). [CrossRef] [PubMed]

**2. **B. Philipson, “Distribution of protein within the normal rat lens,” Invest. Ophthalmol. **8**(3), 258–270 (1969). [PubMed]

**3. **P. P. Fagerholm, B. T. Philipson, and B. Lindström, “Normal human lens - the distribution of protein,” Exp. Eye Res. **33**(6), 615–620 (1981). [CrossRef] [PubMed]

**4. **T. J. T. P. van den Berg, “Intraocular light scatter, reflections, fluorescence and absorption: what we see in the slit lamp,” Ophthalmic Physiol. Opt. **38**(1), 6–25 (2018). [CrossRef] [PubMed]

**5. **F. A. Bettelheim and M. Paunovic, “Light scattering of normal human lens I: Application of random fluctuation and orientation theory,” Biophys. J. **26**, 85–100 (1979). [CrossRef] [PubMed]

**6. **E. L. Siew, D. Opalecky, and F. A. Bettelheim, “Light scattering of normal human lens. II. Age dependence of the light scattering parameters,” Exp. Eye Res. **33**(6), 603–614 (1981). [CrossRef] [PubMed]

**7. **E. L. Siew, F. A. Bettelheim, L. T. Chylack Jr., and W. H. Tung, “Studies on human cataracts. II. Correlation between the clinical description and the light-scattering parameters of human cataracts,” Invest. Ophthalmol. Vis. Sci. **20**(3), 334–347 (1981). [PubMed]

**8. **T. J. Van den Berg and J. K. Ijspeert, “Light scattering in donor lenses,” Vision Res. **35**(1), 169–177 (1995). [CrossRef] [PubMed]

**9. **T. J. van den Berg and H. Spekreijse, “Light scattering model for donor lenses as a function of depth,” Vision Res. **39**(8), 1437–1445 (1999). [CrossRef] [PubMed]

**10. **K. O. Gilliland, C. D. Freel, C. W. Lane, W. C. Fowler, and M. J. Costello, “Multilamellar bodies as potential scattering particles in human age-related nuclear cataracts,” Mol. Vis. **7**, 120–130 (2001). [PubMed]

**11. **M. J. Costello, S. Johnsen, K. O. Gilliland, C. D. Freel, and W. C. Fowler, “Predicted light scattering from particles observed in human age-related nuclear cataracts using Mie scattering theory,” Invest. Ophthalmol. Vis. Sci. **48**(1), 303–312 (2007). [CrossRef] [PubMed]

**12. **K. O. Gilliland, S. Johnsen, S. Metlapally, M. J. Costello, B. Ramamurthy, P. V. Krishna, and D. Balasubramanian, “Mie light scattering calculations for an Indian age-related nuclear cataract with a high density of multilamellar bodies,” Mol. Vis. **14**, 572–582 (2008). [PubMed]

**13. **E. M. Méndez-Aguilar, I. Kelly-Pérez, L. R. Berriel-Valdos, and J. A. Delgado-Atencio, “Simulation and analysis of light scattering by multilamellar bodies present in the human eye,” Biomed. Opt. Express **8**(6), 3029–3044 (2017). [CrossRef] [PubMed]

**14. **S. Bassnett, Y. Shi, and G. F. J. M. Vrensen, “Biological glass: structural determinants of eye lens transparency,” Philos. Trans. R. Soc. Lond. B Biol. Sci. **366**(1568), 1250–1264 (2011). [CrossRef] [PubMed]

**15. **R. J. W. Truscott, “Age-related nuclear cataract-oxidation is the key,” Exp. Eye Res. **80**(5), 709–725 (2005). [CrossRef] [PubMed]

**16. **L. Acosta-Sampson and J. King, “Partially folded aggregation intermediates of human γD-, γC-, and γS-Crystallin are recognized and bound by human αB-Crystallin chaperone,” J. Mol. Biol. **401**(1), 134–152 (2010). [CrossRef] [PubMed]

**17. **A. Spector and D. Roy, “Disulfide-linked high molecular weight protein associated with human cataract,” Proc. Natl. Acad. Sci. U.S.A. **75**(7), 3244–3248 (1978). [CrossRef] [PubMed]

**18. **R. J. W. Truscott and R. C. Augusteyn, “Changes in human lens proteins during nuclear cataract formation,” Exp. Eye Res. **24**(2), 159–170 (1977). [CrossRef] [PubMed]

**19. **G. B. Benedek, “Theory of transparency of the eye,” Appl. Opt. **10**(3), 459–473 (1971). [CrossRef] [PubMed]

**20. **G. B. Benedek, “Cataract as a protein condensation disease: the Proctor Lecture,” Invest. Ophthalmol. Vis. Sci. **38**(10), 1911–1921 (1997). [PubMed]

**21. **A. P. Bruckner, “Picosecond light scattering measurements of cataract microstructure,” Appl. Opt. **17**(19), 3177–3183 (1978). [CrossRef] [PubMed]

**22. **H. Bloemendal, W. de Jong, R. Jaenicke, N. H. Lubsen, C. Slingsby, and A. Tardieu, “Ageing and vision: structure, stability and function of lens crystallins,” Prog. Biophys. Mol. Biol. **86**(3), 407–485 (2004). [CrossRef] [PubMed]

**23. **K. K. Sharma and P. Santhoshkumar, “Lens aging: Effects of Crystallins,” Biochim. Biophys. Acta **1790**(10), 1095–1108 (2009). [CrossRef] [PubMed]

**24. **S. D. Moran, T. O. Zhang, S. M. Decatur, and M. T. Zanni, “Amyloid Fiber Formation in Human γD-Crystallin Induced by UV-B Photodamage,” Biochemistry **52**(36), 6169–6181 (2013). [CrossRef] [PubMed]

**25. **J. C. Boatz, M. J. Whitley, M. Li, A. M. Gronenborn, and P. C. A. van der Wel, “Cataract-associated P23T γD-crystallin retains a native-like fold in amorphous-looking aggregates formed at physiological pH,” Nat. Commun. **8**, 15137 (2017). [CrossRef] [PubMed]

**26. **H. C. van de Hulst, *Light Scattering by Small Particles* (Dover, 1957).

**27. **L. V. Wang and H. Wu, *Biomedical Optics: Principles and Imaging* (Wiley, 2007).

**28. **C. F. Bohren and D. R. Huffman, *Absorption and Scattering of Light by Small Particles* (Wiley, 1983).

**29. **A. Cuadrado, J. Toudert, and R. Serna, “Polaritonic-to-Plasmonic Transition in Optically resonant Bismuth nanospheres for High-Contrast Switchable Ultraviolet Meta-Filters,” IEEE Photonics J. **8**(3), 4801811 (2016). [CrossRef]

**30. **A. Cuadrado, J. Toudert, B. García-Cámara, R. Vergaz, F. J. González, J. Alda, and R. Serna, “Optical tuning of nanospheres through phase transition: An optical nanocircuit analysis,” IEEE Photonics Technol. Lett. **28**(24), 2878–2881 (2016). [CrossRef]

**31. **F. Martelli, S. Del Bianco, A. Ismaelli, and G. Zaccanti, *Light Propagation through Biological Tissue and other Diffusive Media* (SPIE, 2010).

**32. **A. Doicu, Th. Wriedt, and Y. A. Eremein, *Light Scattering by Systems of Particles, Null-Field Method with Discrete Sources: Theory and Programs* (Springer, 2006).

**33. **V. Valery Tuchin, *Tissue Optics, Light Scattering Methods and Instruments for Medical Diagnostics* (SPIE, 2015)

**34. **C. John, Stover *Optical Scattering: Measurement and Analysis* (SPIE, 2012)

**35. **F. A. Bettelheim and S. Ali, “Light scattering of normal human lens. III. Relationship between forward and back scatter of whole excised lenses,” Exp. Eye Res. **41**(1), 1–9 (1985). [CrossRef] [PubMed]

**36. **B. García-Cámara, F. Algorri, A. Cuadrado, V. Urruchi, J. M. Sánchez-Pena, R. Serna, and R. Vergaz, “All-Optical nanometric switch based on the directional scattering of semiconductor nanoparticle,” J. Phys. Chem. C **199**(33), 19558–19564 (2015). [CrossRef]

**37. **S. Gandadhara, “How to Create a User-Defined Scattering Function,” (Zemax Knowledge base, 2008) http://customers.zemax.com/os/resources/learn/knowledgebase/how-to-create-a-user-defined-scattering-function

**38. **J. M. P. Coelho, J. Freitas, and C. A. Williamson, “Optical eye simulator for laser dazzle events,” Appl. Opt. **55**(9), 2240–2251 (2016). [CrossRef] [PubMed]

**39. **N. Pfeiffer and G. H. Chapman, “Successive order, multiple scattering of two-term Henyey-Greenstein phase functions,” Opt. Express **16**(18), 13637–13642 (2008). [CrossRef] [PubMed]

**40. **K. Nam-Hyong, “How to create complex non-sequential objects,” (Zemax Knowledge base, 2005) http://customers.zemax.com/os/resources/learn/knowledgebase/how-to-create-complex-non-sequential-objects?lang=en-US

**41. **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] [PubMed]