Alteration in refractive index proﬁle during accommodation based on mechanical modelling

: The lens of the eye has a gradient refractive index (GRIN). Ocular accommodation, which alters the shape of the lens in response to visual demand, causes a redistribution of the internal structure of the lens leading to a change in the GRIN proﬁle. The nature of this redistribution and the consequence of change in the GRIN proﬁle are not understood. A modelling approach that considers how the GRIN proﬁle may change with accommodation needs to take into account optical and mechanical parameters and be cognisant of individual variability in the shape and size of lenses. This study models the normalised axial GRIN proﬁle during accommodation using reduced modelling and incorporating ﬁnite element analysis to connect inhomogenous mechanical characteristics of the lens to optical performance. The results show that simulated stretching changes the length of the plateau but does not alter the cortical gradient, which supports clinical ﬁndings. There is a very small change to the accommodated and non-accommodated proﬁles when normalised, yet this yields measurable changes in aberrations with around 11% and almost 13% difference in spherical aberration and astigmatism respectively. The results can be used in reconstruction of the refractive index and for investigating gradual changes with age.


Introduction
The lens of the eye is a gradient refractive index (GRIN) structure, a characteristic that attenuates the monochromatic and chromatic aberrations [1][2][3][4][5][6][7]. In younger lenses, the GRIN profile will alter as the lens changes its shape and external geometry, during the process of accommodation. The mechanical properties of the lens will have an effect on the shape alterations that occur in accommodation and will also influence the way in which the cells respond to forces of stretching and relaxation. The nature of this internal response, which may involve a redistribution of cytoplasmic proteins, the crystallins, remains an unknown.
Biomicroscopic studies give a useful perspective of geometrical changes and some indication of how structural ring-like features seen in the living lens as sharp demarcations and termed, the zones of discontinuity, alter with accommodation [8][9][10]. A study by Bahrami et al. [11] suggests a direct connection between the zones of discontinuity and the iso-indicial contour distribution across the lens. Whilst biomicroscopic and Scheimpflug techniques can indicate changes in gross features, the effect on the refractive index from a change in accommodative state cannot yet be directly obtained from a clinical observation. It is not known how the refractive index profiles alters with accommodation and indeed whether the redistribution of tissue with cell layers is the same in all individuals. Early clinical studies suggest that accommodation thickens the nuclear portion of the lens [12,13] which corresponds to the flatter region of the refractive index profile. In vitro studies of the lens have found a difference in the moduli of elasticity of the nucleus and cortex [14]. How these differences may manifest in changes to the isoindicial contours during the accommodative process can be gleened from finite element modelling.
Classical experiments have suggested that for younger eyes the nucleus of the lens has a higher modulus of elasticity than the cortex [14]. This would appear counterintuitive and in particular for older lenses where an increase in stiffness of the lens core has been suggested as the main cause of presbyopia [15][16][17][18][19].
The magnitude of refractive index increases with progression from the periphery and is relatively constant in the centre of the human lens [20][21][22]. The GRIN is created from varying distributions and concentrations of lens proteins [22,23]. Such variations in protein density can be expected to produce a non-constant distribution of material properties. The function of accommodation is fundamentally an alteration in the optics of the lens and is affected by its mechanical and material properties. The lens material is constituted from the crystallin proteins but how these affect the mechanics of the lens is not known.
The aspect of the GRIN structure that is most pertinent to vision is the axial GRIN profile which plays an important role in optical power of the eye and has an effect on the third-order aberrations [6]. The steepness of this profile has been suggested as a measure of age-related change [20,[24][25][26]. This study examines a spatially normalised representation of this axial profile to explore its change with accommodation. Such a profile overcomes the need to consider absolute thickness changes as the GRIN profile alters thereby providing potential for comparing lenses of various age and shape as they undergo accommodation. An axial reduced model of accommodation has been derived to examine GRIN profile adjustments. This incorporates analysis from a finite element model that uses a distribution of modulus of elasticity based on the optical structure introduced by Bahrami and Goncharov [6]. The normalised axial profile of the lens can be used in comparative analyses of mechanical properties of the lens and in advanced modelling of accommodative behaviour.

Reduced mechanical modelling
As the thickness of the lens alters with accommodation, there is a change in the axial GRIN profile. The spatial coordinates of this profile are divided along the equator into anterior and posterior sections thereby producing a normalised GRIN profile that can be approximated by the power law equation [20,24], where r is the normalised distance from the centre of the lens, n C is the maximum refractive index at the centre of the lens, n S is the minimum refractive index at the periphery of the lens and the exponent p is a constant that adjusts for the steepness of the profile. Figure 1 shows a series of profiles for different p values; as p increases the profile central flat portion increases in thickness. r values of −1 and +1 correspond to the anterior and the posterior sections of the lens, respectively. Although parabolic profiles (p = 2) are common in animal eyes and are found in the porcine lens which is physiologically close to the human lens [22,27], the exponent p is greater than 2 for the human lens [11,20,22] and it may alter with age [11,22].
Given that the contours are iso-indicial and that the protein densities that give rise to the optical characteristics also influence the mechanical properties, it is feasible to assume that the material characteristics of the medium along any iso-indicial contour are constant and are described as iso-mechanical contours.
In devising a reduced (axial) model of accommodation, using iso-mechanical contours allows the lens to be treated as a combination of springs in series along its optic axis. The geometrical parameters of these springs are identical except for their initial thickness δt, which can be assumed to be equal to the axial thickness of the corresponding iso-indicial layer. In accordance with Newton's 3rd law of motion, the forces applied to both sides of these springs in series are identical to the force applied at the poles of the lens when it changes thickness. This force F produces a change in the thickness of each spring δ x, which is defined by where δ K is the spring constant and subscript r indicates the dependency of the coefficients to the normalised distance. The summation of all small springs δ x r gives the total change in the axial thickness of the lens during accommodation. The measured mechanical characteristics of the lens of the eye at each point can be quantified as Young's modulus E, which is a measure of the stiffness of a material defined as where A 0 is the cross-sectional area through which the force is applied and δt r is the initial thickness of a small spring at the normalised distance r.
Assuming that the maximal and minimal refractive index values n C and n S are not altered during accommodation, then the normalised refractive index profile will be maintained (i.e. the exponent p will not change) if the proportional change in thickness of each iso-indicial, and hence iso-mechanical, contour (δt r /δ x r ) remains the same along the optic axis (for r ∈ {−1, +1}) under accommodation. For two specific springs located at the centre and surface of the lens this can be described as where the subscripts C and S stand for the central regions (nucleus) and the surface layers (cortex) of the lens, respectively. Given that the cross section for all small springs is paraxially the same and that magnitudes of all applied forces are also identical, substituting Eq. (3) in Eq. (4) gives Equation (5) does not concur with experimental measurements, which suggest differences in mechanical characteristics of the nucleus and cortex [14][15][16][17][18][19] and this challenges the assumption that normalised profile does not alter with accommodation. The majority of experimental data for younger eyes indicate E C < E S [14][15][16][17][18][19], which results . Since the nucleus of the human lens is axially thicker than the cortex [11,20,22] (δt C > δt S ), it can be concluded that as the axial thickness of the lens increases, the nucleus will bear a larger part of the deformation, (i.e. δ x C > δ x S ). This result predicts a larger p value (see Eq. (1)) for the profile at its accommodated compared to its unaccommodated state, yet, any discussion on the magnitude and importance of this difference requires a numerical modelling based on experimental data.

Finite element analysis
The refractive index profile can be translated to a geometrical model [6] as shown in Fig. 2 where R is the radius of curvature of the lens, K the conic constant, and T the lens thickness; subscripts a and p refer to the anterior and posterior sections of the lens, respectively. Taking an initial lens shape equivalent to 10 dioptres change in lens power [28] gives: R a = 6.893mm, R p = 5.450mm, K a = −8.870, K p = −3.810, T a = 2.118mm, and T p = 1.412mm. Taking p = 4.00, n S = 1.370, and n C = 1.420 [29]. The GRIN is distributed over 10 layers of the lens with an incremental difference in refractive index magnitude between each adjacent layer of 0.005 (Fig. 2).
The GRIN distribution presented as shells of incrementally increasing refractive index from periphery to lens centre is shown in Fig. 2. This model also incorporates the distribution of Young's modulus following similar increments and taking E C = 600 Pa and E S = 3000 Pa for the values of Young's modulus of the centre and the surface of the lens, respectively. The distribution of Young's moduli were based on the values from measurements for central and peripheral parts of a 27 year old lens undertaken by spinning lenses of varying ages [14] and have been applied to earlier FEA models that were constructed based on data from human lenses that had undergone stretching in a manner mimicking the action of the ciliary muscle. To The geometry of the lens model in its 10 dioptre accommodation using the GIGL model [6], the study by Sheil et al. [28] and mechanical measurements by Fisher [14].
extend this measurement to a lamellar structure, the incremental change in the modulus in the modelling procedure is considered to be 267 Pa per layer. These values are used in a 10 shell finite element model (Fig. 3). Figure 3(a) depicts the cross-section of the meshed geometry which is comprised of 45184 mapped face meshing elements in total. The model is rotationally symmetric around the optical axis z (Fig. 2). To increase the accuracy of analyses for the fine structure of the capsule, the size of elements used for this thin layer is reduced by 100x compared with those used for the lens (Fig. 3(b)). The thickness of the capsule varies from 0.01 mm at the poles to 0.02 mm at the equator [30]. The capsular modulus = 5 × 10 6 Pa [31] and the cross section of the zonule is represented by a continuous plane with a Young's modulus of 10 7 Pa [32]. Poisson's ratio of 0.47 was used for the lens capsule [33] and 0.49 for the lens material. The non-linear equations of the finite element analysis were solved using ANSYS ® (Academic Research, Release 15.0, ANSYS Inc.). The zonule is represented as a single entity covering a range of insertion points to the lens capsule and to the ciliary body. This is modelled as a curved structure because this accords better than a single point to the in vivo situation in which there is no single insertion point for the different parts of the zonule but a spread of insertions. As the origin of all zonular insertions is not definitively known and anatomical studies have shown new insertions recently [34], an exact replication of the in vivo situation is not yet possible.
Displacing the ciliary body by 0.1 mm in the radial direction, the change in the shape of the  lens and the relevant elastic strain applied to the lens structure are shown in Fig. 4. The 0.1 mm displacement decreases the axial thickness from 3.53 mm to 3.18 mm, the latter represents the thickness of the unaccommodated lens [28]. The change in contour is most marked in the centre with clear boundaries depicting the greater degrees of strain compared to the outer layers. For the innermost contour, the strain around the poles is higher and that in the equatorial zone is lower than for any other contour. With progression towards the periphery, the strain increases in the equatorial region and decreases at the poles. The distribution of high strain in the equatorial zone shows a bifurcation in the outer-third of the lens and this continues to the surface. The outer polar regions have very low strain levels.

Axial refractive index profile
The axial positions of different lens layers for the accommodated and unaccommodated Finite Element Model taken from Fig. 4 are shown in Fig. 5 as blue points. Mathematica (Academic Research, Release 9.0, Wolfram Research Inc.) was used to solve the least squares equations to fit power law curves to the points. The increase in the width of the axial profile is caused by the expansion of the flatter, central section. This concurs with the conclusion, discussed in Section 2, based on experimental findings of a lower value of Young's modulus in the core than in the cortex of the lens [14]. The normalised distributions of axial profiles are shown in Fig. 6. Whilst, the profiles appear to be superimposable, there is a slight difference ie the p exponent equals 4.0 and 3.85 for accommodated and unaccommodated profiles respectively. Hence, the profile becomes flatter with accommodation. The small difference between the exponent coefficients of the profiles indicates a potential for the normalised profile to be treated as an invariant characteristic of accommodating lens, yet its effect on the optics and imaging performance is the main criteria.   In optical modelling of the accommodative lens, ignoring the small geometrical alteration in the normalised profile (i.e. using a constant p value for all states of accommodation) may cause larger errors in aberration calculations. Table 1 presents the Seidel aberration coefficients of a lens with the external geometry as seen in Fig. 2 for the exponents p = 4.00 and p = 3.85. The marginal ray angle is taken as 0.036352 rad with a height of 1 mm, the chief ray angle is 0.015000 rad and the refractive index of the media surrounding the lens is equivalent to aqueous n = 1.336. The pupil is located at the anterior surface of the lens.
The small difference in the value of p between the accommodated and unaccommodated states results in greater differences in spherical aberration, coma and astigmatism ( Table 1). The change in the optical power is less than 1 percent and there is very little change in field curvature. For the lens in Table 1 the average refractive index calculated as in [7] is 1.4100 and 1.4097 for p = 4.00 and 3.85, respectively.

Discussion
The lamellar structure of cells filled with proteins of different proportions across the lens tissue [35], forms the refractive index gradient [21,22]. The non-homogeneity that is caused by variations in refractive index and density across the lens would be expected to create anisotropy in mechanical properties. The nature of this has not been measured. Whilst the GRIN nature of the lens has been investigated in many studies [22], there is comparatively less work on the mechanics of the lens and indeed on the opto-mechanical changes with accommodation. How the lens adjusts its cellular layers and the subsequent alterations to the GRIN are not known because it is difficult to measure this in vivo [10]. Koretz and colleagues [8,9] undertook seminal studies that showed how the zones of discontinuity: demarcations seen in the living lens that have been linked to the refractive index gradient [11] alter with age and accommodation. This offered a perspective on where movement of internal material occurs when the lens changes shape but could not indicate how the GRIN may alter. Using the optics of the GRIN to model a lens with layers that vary gradually in material properties has provided a means of estimating alterations to the GRIN and how these alterations could affect the aberrations of the lens.
Early experimental work on human lenses subjected to spinning has suggested that E C < E S [14]. Subsequent work has not been conclusive. Mathematical interpretations that have recalculated the early experimental work report considerable scatter particularly for the cortical modulus [16]. An improved replication of the spinning experiment has indicated the need for refinement of the original work and has concluded that nuclear stiffness is a less accurate value to measure than cortical stiffness [19]. Dynamic mechanical analysis of human lenses ranging in age from 14 to 78 years found that in younger lenses up to age 30, the nucleus was softer than the cortex but that over the age of 50 this reversed and the nucleus was more rigid than the cortex [15]. This was supported by findings of a greater stiffness in the cortex in lenses capable of accommodation and a greater stiffness in the nucleus in lenses above the age of 50 [18]. An acoustic radiation technique also found a greater nuclear than cortical rigidity but did not test younger lenses capable of high accommodative capacity [17]. The model in this paper considered 10D of accommodation in order to compare changes in optical properties across a wide accommodative difference. Since only young lenses aged within the first two decades are capable of such a level of accommodation, the model incorporated a lower value of Young's modulus in the nucleus than in the cortex and used the values originally found by Fisher 1971 [14]. The finding of a softer nucleus compared to the cortex has been confirmed in subsequent studies [15,18] that measured shear modulus across lenses of different ages. Gradients of moduli that varied across the nuclear region were found in older presbyopic lenses [15,18] but were not so definitive in younger lenses capable of accommodation. It should also be noted that these studies had been conducted on lenses that has been frozen and thawed [15,18]. How this may translate to the fresh accommodating lens is not known. The model presented in this paper serves to show how the nature of the GRIN may change with accommodation and the effect that this may have on image quality. It is not intended to provide a definitive mechanical model. The variations in modulus were modelled as shells or contours and in any given shell, Young's modulus is an average over that section of the gradient. The central shell in the model which approximates around half of the width of the lens and provides a constant central region is not inconsistent with the trend seen for a young lens [15].
The FEA model simulated 10D of accommodation with a displacement of 0.1mm, which is an 0.2mm change in equatorial diameter and which resulted in an 0.35mm change in lens thickness. The latter is within the range of values reported in the in vivo study of Dubbelman et al. [36] who found a mean change in lens thickness of 0.045 ± 0.012mm/D. This equates to a mean of 0.45mm and a range of 0.33 to 0.57mm for a 10D change in accommodation to the simulation modelled in this study. In vivo measurements by Strenk et al. 1999 [37] show a change in lens thickness of −0.15 to 0.63mm for lenses aged between 22 and 50 years for a change in accommodation of between 0.1 and 8.0D. This study also reported that the alteration in lens equatorial diameter for the same accommodative change ranged from 0.93mm to 0mm. The 0.2mm change in lens equatorial diameter modelled in this study is within the lower end of the range of values reported by Strenk et al. 1999 [37]. Figure 4 shows that with stretching in the radial direction, the strain which represents geometrical change or displacement is greater in the polar regions of the inner layers than in the same region of the outer layers. The greater deviation, in the central layers compared to those in the periphery, when comparing stretched and unstretched forms of the FEA model, translates to a widening of the central plateau part of the GRIN with simulated accommodation (Fig. 5).This is dependent on the forces applied and also on the material properties of the lens. It is notable that the central flat section of GRIN profile widens with accommodation, consistent with early clinical reports of an widening of the nuclear region when the lens accommodates [12,13] and with greater changes seen in internal zones of discontinuity as the lens accommodates [8].
The very high levels of strain in the outer layer equatorial regions reflect, to some extent, the proximity of these parts of the lens to the sources of force application. This is the region pertinent to vision and the changes in lens shape with accommodation need to be such as to control aberrations and maintain image quality for all focussing distances. The thickening of the nucleus, which corresponds to the plateau region of the refractive index ensures that aberrations are controlled. Further investigations on the opto-mechanical properties of the lens are needed to understand changes with age as well as with accommodation.
Modelling provides a means of depicting the effect of mechanical parameters on the optical performance of the lens by examining how external shape and internal structure may alter with accommodation. Advanced models have considered the iso-indicial contours and whether these adhere to the external shape of the lens applying a contour distribution that accords with the normalised axial refractive index profile.
This study, linking optical to mechanical characteristics, shows that even if the change in the normalised GRIN profile, with accommodation, is small, it has a measurable effect on the aberrations, and consequently on the emergent ray path. These results are very useful for examining reconstructions of refractive index based on in vivo measurements such as those applied using Optical Coherence Tomography [38]. The normalisation method presented in this study offers the potential of setting a standard for verification of accuracy in the reconstruction. If the reconstructions of a lens in different states of accommodation show more than a given difference in p values, the method and parameters used need to be re-examined. Optomechanical modelling can provide a further means of probing the relationship between optical and mechanical properties of the lens and how these affect lenticular function and its alteration with age. The ultimate aim of the modelling approach is to reconstruct a system that is as close as possible to the biological tissue. Future models of the eye lens need to incorporate material and optical parameters that mimic the constituents of the lens and strive for a continuous GRIN profile that is similar to that found in the human lens [22]. Such methods can be developed to study gradual changes in the lens that affect the process of accommodative loss with age and present a novel means of studying presbyopia.