Journal of Optometry (2012) 5, 110-120

Journal

Optometry

www.journalofoptometry.org

ORIGINAL ARTICLE

Hyperelastic modelling of the crystalline lens: Accommodation and presbyopia

Elena Lancharesab*, Rafael Navarroc, Begona Calvo

a Aragón Institute of Engineering Research (I3A), University of Zaragoza, Zaragoza, Spain b Centro de Investigación Biomédica en Red en Bioingeniería, Biomateriales y Nanomedicina (CIBER-BBN), Spain c ICMA, Consejo Superior de Investigaciones Científicas (CSIC) and University of Zaragoza, Zaragoza, Spain

Received 20 February 2012; accepted 25 April 2012 Available online 12 June 2012

KEYWORDS

Finite element method;

Biological tissues;

Transversely isotropic

hyperelastic

behaviour;

Accommodation; Presbyopia

Abstract

Purpose: The modification of the mechanical properties of the human crystalline lens with age can be a major cause of presbyopia. Since these properties cannot be measured in vivo, numerical simulation can be used to estimate them. We propose an inverse method to determine age-dependent change in the material properties of the tissues composing the human crystalline lens.

Methods: A finite element model of a 30-year-old lens in the accommodated state was developed. The force necessary to achieve full accommodation in a 30-year-old lens of known external geometry was computed using this model. Two additional numerical models of the lens corresponding to the ages of 40 and 50 years were then built. Assuming that the accommodative force applied to the lens remains constant with age, the material properties of nucleus and cortex were estimated by inverse analysis.

Results: The zonularforce necessary to reshape the model of a 30-year-old lens from the accommodated to the unaccommodated geometry was 0.078newton (N). Both nucleus and cortex became stiffer with age. The stiffness of the nucleus increased with age at a higher rate than the cortex.

Conclusions: In agreement with the classical theory of Helmholtz, on which we based our model, our results indicate that a major cause of presbyopia is that both nucleus and cortex become stiffer with age; therefore, a constant value of the zonular forces with aging does not achieve full accommodation, that is, the accommodation capability decreases. © 2012 Spanish General Council of Optometry. Published by Elsevier España, S.L. All rights reserved.

* Corresponding author at: Mechanical Engineering Department, University of Zaragoza, María de Luna 3, E-50018 Zaragoza, Spain. E-mail address: elanchar@unizar.es (E. Lanchares).

1888-4296/$ - see front matter © 2012 Spanish General Council of Optometry. Published by Elsevier España, S.L. All rights reserved. http://dx.doi.org/10.1016/j.optom.2012.05.006

Modelado hiperelástico del cristalino: acomodación y presbicia Resumen

Objetivo: La modificación de las propiedades mecánicas del cristalino humano con la edad puede constituir una causa principal de la presbicia. Como dichas propiedades no pueden medirse ''in vivo'', se puede utilizar la simulación numérica para su cálculo. Proponemos un método inverso para la determinación del cambio con la edad de las propiedades materiales de los tejidos que componen el cristalino humano.

Métodos: Se desarrolló un modelo de elementos finitos de cristalino de 30 anos de edad en estado acomodado. Se calculó la fuerza necesaria para lograr la acomodación plena en un cristalino de 30 años de edad utilizando este modelo. A continuación se construyeron dos modelos numéricos adicionales de cristalino para edades de 40 y 50 anos. Suponiendo que la fuerza acomodativa del cristalino permanece constante con la edad, se calcularon las propiedades de material del núcleo y la corteza mediante un análisis inverso.

Resultados: La fuerza zonular necesaria para reconstruir el modelo de un cristalino de 30 años de edad, partiendo de la geometría acomodada hasta alcanzar la no acomodada, era de 0,078 Newton (N). Tanto el núcleo como la corteza adquirieron más rigidez con la edad. La rigidez del núcleo se incrementaba con la edad a un porcentaje superior a la de la corteza. Conclusiones: De acuerdo con la teoría clásica de Helmholtz, en la que nos basamos, nuestros resultados indican que una de las principales causas de la presbicia es que tanto el núcleo como la corteza adquieren más rigidez con la edad y que, por tanto, el valor constante de las fuerzas zonulares no logra una acomodación plena con el envejecimiento. En consecuencia, se produce una disminución de la capacidad de acomodación.

© 2012 Spanish General Council of Optometry. Publicado por Elsevier España, S.L. Todos los derechos reservados.

PALABRAS CLAVE

Método de elementos finitos;

Tejidos biológicos;

Comportamiento

hiperelástico

transversalmente

isótropo;

Cristalino;

Acomodación;

Presbicia

Introduction

The mechanism of accommodation and its gradual decline with age is a subject of growing interest due to the high prevalence of presbyopia. One of the most significant problems and current limitations in analyzing accommodation and presbyopia is the lack of a single universally accepted theory on the subject. Ex vivo1,2 and in vivo3-5 experiments seem to suggest that the Helmholtz theory is probably the most appropriate. Its main assumption is that full accommodation of the lens corresponds to its 'natural shape' that is maximum axial thickness and maximum surface curvatures, when no external forces are applied. The lens flattens under the relaxation of the ciliary muscle, then the axial thickness and both anterior and posterior curvatures decrease. As a result, the power of the lens gradually decreases until the unaccommodated state is reached.6,7

In order to develop a plausible theory or model of accommodation and its gradual decline with age we develop a finite element (FE) model of the human lens with a twofold purpose. The first goal is to estimate the zonular forces acting during the process of accommodation using a FE model of a 30-year-old lens. The second goal is to estimate the material properties at any age, knowing the geometry in both accommodated and unaccommodated states. To implement changes for the geometry of the lens with age and accommodation, we followed the empirical studies of Dubbelman and co-workers.5,8 The capsule shows a nonlinear behaviour over finite strains.9 Following this finding, the tissue was modelled using a non-linear constitutive model.

The variation of the capsule thickness with the radial position10 and with age11,12 was also included in the model. We assume that the maximum zonular forces obtained in the first step do not change substantially with age13 since the ciliary muscle seems to remain functional throughout the lifespan.14,15 This method has been applied to obtain the changes of the material properties from 30- to 40- and 50-year-old (y.o.) lenses.

It is worth remarking that numerical models and computer simulation are often the only way to obtain quantitative information about forces or material properties that are difficult to measure experimentally. Mechanical measurements can be taken in vitro, but they are extremely difficult to obtain in vivo since even minimally invasive techniques can potentially modify the material or the equilibrium of forces. The non-physiological conditions of measurements in vitro make the extrapolation of data to the in vivo lens uncertain. Here again, models and simulations are necessary to evaluate the feasibility of these extrapolations. Given their crucial role, biomechanical models of the lens have been previously developed to quantify the stress during the process of accommodation and to examine contributions of individual constituents,16"22 most of them assuming a linear elastic behaviour of the tissues.

The optimization-based methodology presented in this work allows estimation of material constants of complex constitutive laws frequently used in biomechanics. This technique may be also used in patient-specific cases for pre-surgical planning as well, which would suppose a crucial advance in the customized modelling of biological tissues.

Methods

The method developed by Burd et al.17 was applied here, after several improvements and updates, to create a finite element model of a 30y.o. fully accommodated up to 7.5 D, human lens. Our main assumption is that the zonular force acting on the crystalline lens in this fully accommodated state (near vision) is negligible, according to Helmholtz's theory. Conversely, zonular forces are required to flatten the lens in such a way that it reaches the unaccommodated state (far vision). Numerical simulation allows estimation of these forces based on an optimization procedure.

A detailed description of the lens geometry used to develop the FE model is presented in Appendix. Fig. 1 depicts the geometrical parameters of the two-dimensional model.

Fig. 2 shows the process to obtain the FE model of the crystalline lens from the parametric geometric model described in Appendix. Both accommodated and unaccommodated geometries can be established for a given age (Fig. 2(a)). The 3D FE model depicted in Fig. 2(b) is composed of three parts: nucleus (red), cortex (cyan) and capsule (purple).

To compute the change of power (D) with force, or accommodation, here we assumed a simple optical model with a homogeneous equivalent refractive index. In principle, we could use an adaptive GRIN lens model,23 or estimate the relative contribution of the GRIN structure to accommodation (by extrapolating in vitro data obtained from primate lenses24). However, there are several reasons for choosing the homogeneous index model. One reason is that empirical GRIN data are still scarce in humans, but the main reason is that our geometrical model is based on the Dubbelman et al. empirical data,5,8 which were obtained also assuming a homogeneous equivalent index for the lens. Therefore that choice was made for the sake of simplicity and (what is more important) for consistency. The optical properties of the lens depend on both the inner distribution of refractive index and the curvatures of the external and internal surfaces. Thus the FE model must faithfully reproduce the in vivo lens geometry. As a result of the trade-off between minimum complexity and maximum rigor, different material behaviour for nucleus, cortex and capsule are considered.

To model the material behaviour of the three tissues, we used a quasi-incompressible anisotropic hyperelastic constitutive model.25 Quasi-incompressibility was assumed because of the high level of water of the lens26 so the lens volume does not vary during accommodation.27,28 A strain energy function V written in a decoupled volumetric-isochoric form (Vvol, Viso) is used29,30:

V = &voi + Vf«, = DLnJ2 + y (71 - 3)

C3 - 2

+ ^fexp[C4(/4 -1) ] -1},

where 1/D is a penalty coefficient to quasi-enforce null volumetric change, J = detF, being F the deformation gradient tensor, 7| is the first modified strain invariant of the symmetric modified right Cauchy-Green tensor C, l4 is the square of the stretch along the fibre preferential direction, C1 is the Neo-Hookean constant and the parameters C3 and C4

characterize the stiffness of the preferential direction of deformation.

The main difficulty was to estimate the values for the material properties. Due to the scarce data in the literature of properties of nucleus and cortex, the values reported by Fisher11 were used in this work. In contrast, measurements reported in literature20,21 seem sufficient to completely characterize the capsular properties. To obtain these parameters for the capsular tissue, we used the uniaxial tests performed by Krag and Andreassen19,20 as a basis, fitting the stress versus strain curves by means of least squares method. These values are compiled in Table 1. The increase in stiffness of the capsular tissue with age is negligible,2,19,21 so equal values were assigned to the material parameters of the capsule in the three numerical models (aged 30, 40 and 50). Since the behaviour of the capsular tissue is increasingly stiffer circumferentially towards the equator,21 a preferential circumferentia^direction of deformation was considered and introduced by /4 (Fig. 2(c)).

The fibres in the nucleus are not clearly arranged, therefore the nucleus is considered as an isotropic material and modelled by a Neo-Hookean model. In spite of the clear arrangement of the fibres in the cortex31 we modelled this tissue as an isotropic material due to the lack of data for the specific material parameters of these fibres.

To establish the boundary conditions, some assumptions were made. First of all, the Helmholtz theory was followed assuming that the force delivered by the ciliary muscle through the zonular fibres flattens the lens until the unaccommodated state configuration is reached. The rigid body motion of the model is not relevant for the purpose of our work, which is focused on the deformation of the crystalline lens. Thus, the movement of the lens along the axial direction was neglected.

The nodes placed along the optical axis were only allowed to move along this direction. The nodes located in the equator are not allowed to move along the optical axis. To achieve an equilibrium solution, the sum of axial components of the zonular forces must be zero.

Although the zonular fibres themselves were not modelled, an assumption has to be made about the equatorial region where they are attached to the lens. In general, the insertion region of the zonular fibres is divided into three ring bands: anterior, central (equatorial) and posterior32,33 provides data concerning the location of the insertion regions: anterior 1.5 mm and posterior 1.25 mm from the equator. The width of the anterior and the posterior regions (0.4 mm and 0.5 mm, respectively) was obtained from Ludwig.34 Because the width of the central region is unknown, the same width as for the posterior insertion region of zonular fibres was used (0.5 mm). Fig. 3 shows the directions of application of the ciliary forces and the location of the bands where ciliary processes are inserted in the capsule. The values of the distances of these bands from the equator of the lens are shown in Table 2. The nodes of the FE mesh situated in these bands represent each fibre insertion hence the points of application of the zonular forces.

A set of zonular forces was applied on the FE model of the lens and then the unaccommodated configuration was computed. To determine the zonular forces acting on the human crystalline lens we considered four parameters: total thickness TT, anterior curvature Cant, anterior conic constant Qant

Figure 1 Geometrical parameters of the two-dimensional model.

•JSHEtfJ■■■

Figure 2 Process to obtain the FE model of the crystalline lens from the parametric geometric model described in Appendix.

Table 1 Values of the material parameters of nucleus, cortex and capsule at the age of 30.

Parameters C1 (MPa) D (MPa-1 ) C3 (MPa) C4

Nucleus11 9.3667 x 10-5 214.96 0.0 0.0

Cortex11 5.8295 x 10-4 34.54 0.0 0.0

Capsule19-20 0.2160 0.2835 0.0339 9.7406

Figure 3 Directions of application of the ciliary forces and the location of the bands where zonular fibres are inserted in the capsule.

and posterior curvature Cpos. The corresponding squared relative errors (ECant, EKant, ECpos and ETT) are summed to build the cost function W35:

W = Ecant + 10 EQant + Ecpos + Ejj . (2)

The cost function minimization gave the optimal set of forces.

Then, two models of the lens at the ages of 40 and 50 years were developed, following the same process as described previously. The zonular forces obtained previously were now an input to the process, and further simulations were carried out. Now, the parameters of the material properties varied and the cost function (2) corresponding to each model was minimized.

Results

Fig. 4 shows the unaccommodated and accommodated geometries of the 30-year FE model after application of zonular forces. When the cost function reaches a minimum (0.015 and 0.020 for the 40 and 50-year model, respectively), the unaccommodated geometries match the previously established anterior and posterior curvatures, corresponding to the fully unaccommodated state of the lens. Fig. 5(a) and (b) shows the conicoids (green curves)

Table 2 Distances of insertion bands from the equator (Xant, Yant, Xpos and Ypos in Fig. 3).32"34 The right column depicts the width of the band.

Band X (mm) Y (mm) Width (mm)

Anterior 1.50 1.04 0.4

Posterior 1.25 1.29 0.5

Equatorial 0.00 0.00 0.5

Figure 4 Unaccommodated and accommodated geometries of the 30-year FE model after application of zonular forces.

that fit the nodes of the anterior and posterior surface in the unaccommodated configuration superimposed on the reference curves (dashed curves) which are the previously described anterior and posterior surfaces for the 30-year model, in the fully unaccommodated state. A central area of radius 2.5 mm and 2 mm for the anterior and posterior parts, respectively, was used to compare the reference surface with that obtained by numerical simulation. The posterior curvature of the lens matches the reference curve, and the agreement of the anterior curvatures is reasonably good in a central area of radius 2 mm.

Fig. 6(a) shows the maximum principal strain distribution of the capsule in the unaccommodated state of the 30-year FE model. Since the higher logarithmic strain (LE) in the capsule is 0.0703, the strain is Xprindpal = eLE = e0 0703 = 1.0728 i.e. 7.28%. This value supports the initial assumption of using a hyperelastic constitutive model for the crystalline lens tissues, since they show a finite strain behaviour.21

Fig. 6(b) shows the maximum principal stress distribution in the capsule corresponding to the unaccommodated state of the 30-year numerical model. The maximum value corresponds to the posterior pole area, where the capsular tissue is thinner. The stress decreases towards the equator (160.3-82.84kilopascal (kPa)) corresponding to increases in thickness.

Fig. 6(c) shows the maximum principal stress distribution in cortex and nucleus corresponding to the unaccommodated geometry of the 30-year FE model. The maximum value (0.5526 kPa) corresponds to the red spot at the location of the zonular insertion band, where the forces act. The maximum principal stress in the optical area of the cortex is a tensile stress of about 0.2 kPa, meanwhile a compressive stress of 0.25 kPa turns up in the nuclear tissue. These values are 2-3 orders of magnitude lower than those obtained for the capsule.

The maximum principal stress distribution resulting from numerical simulation for the 40- and 50-year FE models was similar to that of the 30-year model. The maximum values were 96 and 61.97kPa, respectively, for the capsule, and

Figure 5 Conicoids (green curves) that fit the nodes of the anterior (a) and posterior (b) surface in the unaccommodated configuration superimposed on the reference curves (dashed curves).

1.364 and 2.050kPa for nucleus and cortex. The observed tendency of increasing stress in the interior of the lens and decreasing stress in the capsule with age during the process of accommodation is due to increases in capsular thickness with age.

The change in thickness of anterior and posterior cortex (ACT, PCT) and nucleus (NT) caused by the process of accommodation measured by Dubbelman et al.8 is compared to the values obtained by numerical simulation (Table 3). The FE simulation reproduces the empirical observation that the variation of the total thickness of the lens with the accommodation is mainly due to the change in thickness of the nucleus: 0.3 mm changes NT, TTchanges0.33mm, according to Dubbelman, and 0.234 mm versus 0.387 mm, according to the numerical simulation. In both cases, Dubbelman and FE simulation, the values of the variation of the nucleus and cortex thicknesses were of the same order of magnitude. There is a small discrepancy between the experimental results by Dubbelmann, where the thickness

of the cortex was found to be nearly constant, and the FE simulation that predicts a non-negligible change with accommodation.

The lens radius Rlens increases from 7.5 diopters (D) to 0 D. For the 30-year FE model, the simulation provided an increment of 0.32 mm (from 4.56 mm to 4.88 mm), which is somewhat higher than the 0.28 mm measured by Strenk et al.36 The change in Rlens obtained for 40 and 50 years old were 0.24 mm and 0.17 mm, respectively.

Fig. 7(b) shows a test of consistency between the empirical change of lens radii with accommodation (Dubbelman et al.8) and our finite element model predictions for the age of 30. Dubbelman and co-authors provided the expressions of the lens radii R as a function (f) of accommodation A(R = f(A)). On the other hand our model predicts the change of these radii as a function of applied force F(R = g(F)) where g is a linear function (see Fig. 7(a)). One simple and direct way to cross check these two models is to combine both functions (models) f and g. If both models are

LE. Max. In-Plane Principal (Ave. Crit: 75%) +7.030e-02 +6.794e-02 +6.557e-02 +6.321e-02 +6.085e-02 +5.848e-02 +5.612e-02 +5.376e-02 +5.139e-02 +4.903e-02 +4.667e-02 +4.431e-02 +4.194e-02

S. Max. Principal (Ave. Crit.: 75%) +1.603e-01 +1.538e-01 +1.474e-01 +1.409e-01 +1.345e-01 +1.280e-01 +1.216e-01 +1.151e-01 +1.087e-01

g+1.022e-01 +9.575e-02 +8.930e-02 +8.284e-02

S. Max. Principal

(Ave. Crit.: 75%)

r +5.526e-04

- +4.860e-04

- +4.193e-04

- +3.527e-04

- +2.861e-04

- +2.194e-04

- +1.528e-04

- +8.612e-05

- +1.948e-05

- -4.717e-05

1 --1.1386-04 - -1.805e-04 L-2.471e-04

Figure 6 Maximum principal strain distribution of the capsule with the FE model (a), numerical model (b) and cortex and nucleus (c) in the unaccommodated state of the 30-year FE model.

Table 3 Values at the age of 30 of anterior cortex (ACT), nucleus (NT), posterior cortex (PCT) and cortex (CT), total thickness

of the lens (TT) and lens radius (Riens ;) for the 7.5 D and 0D states, and increment (A) caused by accommodation comparing

values obtained by FE simulation to those by Dubbelman et al. 46

7.5 D 0 D (FEM) 0 D (Dubbelman) 4 (FEM) 4 (Dubbeman)

ACT (mm) 0.89 0.804 0.87 0.086 0.02

NT (mm) 2.49 2.256 2.19 0.234 0.30

PCT (mm) 0.58 0.513 0.57 0.067 0.01

CT (mm) 1.47 1.317 1.44 0.153 0.03

TT (mm) 3.96 3.573 3.63 0.387 0.33

Riens (mm) 4.56 4.88 4.84 0.32 0.28

consistent they must correspond to the same value of the radius so that f = R = g, or simply f=g. This provides a new expression with only two variables, accommodation (A) and force (F). Then we solve the equation for A (twice, for

both anterior and posterior radii) to obtain the curves of Fig. 7(b). In case when both models were totally consistent (compatible), we should expect both curves to be identical (red line in Fig. 7(b)). The difference observed between the

Anterior

Posterior

-1-1-1-1 1

0,04 0,06

Zonular force (N)

-Anterior - Posterior -Mean

0,02 0,04 0,06

Zonular force (N)

a4 DC 4

Anterior Posterior ____ ■ h-■

1-----

0,04 0,06 0,08 Zonular force (N)

0,04 0,06 0,08 Zonular force (N)

Anterior Posterior m---—■ --«

-------

0,04 0,06 0,08 Zonular force (N)

— Anterior

-m- Posterior

_ Mean

0,04 0,06 0,08 Zonular force (N)

Figure 7 Consistency between the empirical change of lens radii with accommodation.

two curves is a measure of the discrepancy between empirical data and model predictions. The physical meaning of Fig. 7 is that compared to our FE simulation, in the Dubbel-man's expressions the posterior radius changes too much (and the anterior too little) with ciliary force. The discrepancy is of the order of 1 D in a range of 8D (12.5%). Conversely, our model predicts a smaller change of the posterior radius (and greater for the anterior radius) with accommodation than that observed empirically.

The empirical expressions of both radii Rant, Rpos as a function of the applied force F for the 30y.o. model of the lens are Rant(F) =44.72F +6.71 (R2 =0.9983) and Rpos(F) = 11.21F+4.71 (R2 =0.9903). The same analysis was performed for the ages of 40 and 50 years old (see Fig. 7(c)-(f)) and the corresponding linear regressions were obtained: for the 40-year FE model, Rant (F) = 31.813F +7.3206 (R2 =0.9925) and Rpos(F) = 11.254F +4.8426 (R2 =0.9936); for the 50-year FE model, Rant (F) = 20.459F +8.1851 (R2 =0.9964) and Rpos(F) = 7.3749F +5.0518 (R2 =0.999).

C1 is a material parameter that indicates the stiffness, and the estimated values of this parameter for the ages of 40 and 50 years are shown in Table 4. Both cortex and nucleus become stiffer with age. The cortex at the age of 40 is 1.12fold stiffer than at the age of 30, and at 50y.o. becomes 1.03-fold stiffer than at 40. The nucleus at 40y.o. is 1.40fold stiffer than at 30y.o., and at 50y.o. is 1.64-fold stiffer than at 40.

Discussion

In this paper a parametric 3D FE model of the human crystalline lens was developed to estimate both the zonular forces acting during the process of disaccommodation and the change of the material properties of the lens with aging. A non-linear behaviour was considered for the three tissues and the capsular tissue was modelled as a non-isotropic material.

Both the geometry and the method are based on previous studies.13,17,37 As proposed by some authors17,23 the equatorial surface of the lens is non-planar, and so it was modelled. To model the outermost geometry of the lens, a more compact expression could have been used to describe the lens geometry, for instance that proposed by Kasprzak,28 which consists of a single continuous function capable of reproducing the whole lens surface. We developed a parametric geometry, compiling the work of several authors, to make it dependent on age and state of accommodation.

Most finite element models in the literature assumed axial symmetry9,16,17 and the tissues were modelled as linear elastic materials.13,16,17 We created a three-dimensional model which allows us to introduce the anisotropy of the capsular tissue and a hyperelastic quasiincompressible model was used16,28,38 which is justified by the empirical observation of a nonlinear pseudoelastic behaviour over finite strains.18,20,21,39

An analysis of the position of the attachment of the zonular fibres to the lens was also performed. Bands of insertion wider than those used in the presented work were considered. Each one joined the one placed next, forming a single continuous band at the equator. The outcomes of numerical

simulation showed a slight variation of the total error with respect to the presented work.

The forces required to achieve the change of shape of the model from the accommodated to the unaccommodated state is 0.078 N, about five times greater than the 0.015 N estimated by Fisher.40 However, our result agrees with previous findings by finite element simulation of the accommodative process. Burd et al.17 obtained values of zonular forces ranging from 0.08 N to 0.1 N; Hermans et al.35 obtained a value of 0.081 N. These values are slightly greater than ours, but these authors did not include the variation of the capsular thickness with the radial position, considering this value as a constant for anterior and posterior surface of the lens and therefore being less stiff.

The estimated value of the zonular force is consistent with previous empirical in vitro measurements for the whole lens41 or for the lens capsule2,39 who obtained values of the same order of magnitude.

We assumed that the zonular forces remain constant with age. Some authors support this assumption13 since the ciliary muscle seems to remain functional with age. Presbyopia can be attributed, besides increase in lens stiffness, to continually decreasing zonular tension secondary to life-long increases in lens thickness, making accommodative ciliary muscle movement irrelevant (modified geometric theory, Strenk et al.42). Nevertheless, the age-related changes in zonular forces are not yet known and the assumption of their preservation is not widely accepted. In order to evaluate the effect of varying the zonular forces, 1.5-fold the estimated value was applied in the model. A linear increase of the anterior and posterior radii was observed (Fig. 7(a)), leading to negative accommodation at values of zonular force higher than 0.078N. Nevertheless, this effect slows down at greater values of the force, which is also observed at the ages of 40 and 50 years old. Perhaps increasing zonular forces, causing negative accommodation, could help to compensate the ''lens paradox''3 (decrease of dioptric power of the eye despite the increase of lens surface curvatures with age).

The discrepancy found in the predicted disaccommodation caused by increasing force between the anterior and posterior surfaces (Fig. 7(b)) is in part a consequence of the outcome of the cost function. It also can be due to the mismatch between the rates of change of anterior and posterior curvatures. Therefore, to obtain a finer tuning this residual inconsistency must be minimized to achieve the coincidence of both curves. If we trust our FE simulation, then the parameters of the experimental data should be modified: increasing the rate of change of the anterior curvature and conversely decreasing the rate of change of the posterior curvature with accommodation. The opposite strategy would be to modify other parameters of the model. Since the resultant force acts radially, any modification of the FE model would require modification of the geometry (no better agreement can be obtained by modification of the material properties) but in all cases some modification of the geometric parameters should be needed to obtain a perfect match between model and data.

According to the images obtained by a Scheimpflug camera of the inner structure of the lens8 the thickness of the cortex appears to remain nearly constant during the process of accommodation, the variation of the total thickness

Table 4 Material properties obtained for nucleus and cortex corresponding to 40 and 50-year FE models, compared to the values of the 30-year FE model. Ci in MPa and D in MPa-1.

Parameters Cortex C1 Cortex D Nucleus C1 Nucleus D

30y.o. 0.58295 x10-3 40y.o. 1.3423 x 10-3 50y.o. 2.5168 x10-3 34.54 15.00 8.00 0.93667x 10-4 2.8523 x 10-4 6.3255 x 10-4 214.96 70.59 31.83

being mainly due to the change in the nucleus thickness. This effect is also obtained by the FE simulation. Nevertheless, Dubbelman et al.8 proposed that ACT and PCT remain constant during accommodation, which was not observed in our simulation. The values of the variation of the nucleus and cortex thicknesses given by numerical simulation and measured by Dubbelman et al.8 were of the same order of magnitude.

The increment of Rlens obtained by the simulation, for the 30-year FE model, was of 0.32 mm, versus the 0.28 mm estimated by Strenk et al.36 The lens equatorial radius was not included as a restriction in the cost function, thus this could have caused the mismatch between the two values. This minor discrepancy could be related with the small volume change allowed in the FE simulation, but again the uncertainty in the experimental estimation could also explain it. The changes in Rlens for 40 and 50 years old were 0.24mm and 0.17 mm, respectively.

As was expected11,37 our results showed an increase in stiffness with age for both nucleus and cortex. The stiffness of the nucleus increased with age at a higher rate than the cortex, which is in agreement with other works.43,44 In spite of the assumptions made (zonular forces constant with age, no motion along the optical axis during the accommodative process, etc.), our results support classic theories that the increasing stiffness of the tissues of the lens can be a major cause of presbyopia.

The methodology developed so far is intended to be a useful tool to estimate parameters that are non-measurable in vivo, such as the zonular forces and the material properties of the tissues at a given age. Moreover, in this way we hope it could be applied not only to study presbyopia but also to find ways to restore accommodation. In the last few years, with the development of cataract surgery techniques, intraocular lenses have been continuously improved. Numerical simulation of lens and accommodation can help in design of intraocular lenses. The optics of the intraocular lens could be determined for a given case if a patient specific numerical model is developed. Then, the estimation of parameters such as the zonular forces or the antero-posterior displacement of the lens will help in determining the best optic design of an accommodative lens for that patient.

Financial disclosure

The authors have no financial interest in this work. Acknowledgements

The authors gratefully acknowledge the Instituto de Salud Carlos III (ISCIII) ahd the CIBER-BBN (Centro de

Investigación Biomédica En Red en Bioingeniería, Bioma-teriales y Nanomedicina) initiative, and also the research support of the Spanish Ministry of Education and Science through the research projects DPI2008-02335 and FIS2008-00697. We also wish to acknowledge the collaboration of the student M. Blasco.

Appendix.

(i) External Geometry

Fig. 1 shows a schematic diagram of the basic lens geometry assumed here. The way to build the lens geometry is based on the method proposed by.13 The main difference is that here the central surface dividing the anterior and posterior part of the lens is a conicoid surface instead of the equatorial plane. As a result, the equators of nucleus and cortex lie on different planes in agreement with recent experimental findings.45

Another important difference is that the external lens surfaces are divided into three parts in order to join the two curves, corresponding to the central46 and equatorial37 edges, based on experimental data. The intermediate part is computed to guarantee continuity in the surface and in its first derivative.

The central part of the anterior and posterior surfaces of the lens (red in Fig. 1) are conicoid surfaces with revolution symmetry (see Eq. (A1)) where C is the curvature and Q is the conic constant.

z1 = z0 +--= (A1)

1 + V1 -(Q + 1)C2x2

According to,37 the equatorial edge z3 (green curve in Fig. 1) is assumed to be an arc of circumference with radius re (Eq. (A2)). It connects the anterior and posterior surfaces and is positioned at the equator at the previously established lens radius Rlens.

z3 = \Jr¡ - (x - Rlens + re)2 (A2)

Following,37 the ratio between re and Rlens is constant (re/Rlens = 0.1208), as well as the values of the anterior and posterior angles which determine the length of the equatorial arc (9ant = 63°, 6pos = 37°).

The central part of the lens z1 (red curves) and the equatorial edge z3 (green curves) must be connected through an intermediate surface (blue in Fig. 1) to complete the outer geometry of the lens. The blue curves, z2, are conic surfaces, designed to guarantee continuity with the red and green curves. The connecting points are named Pu between the red and blue curves, and P2, between the blue

Table 5 Geometrical parameters of the lens as a function of age A (y.o.) and vergence of the stimulus S in diopters

(D) 3,5,8,42,47

Regression function

Cant (mm 1 ) (12.7-0.058A)-1 +0.0077S

Qant -4-0.5S

Cpos (mm-1 ) (5.90 - 0.013A)-1 +0.0043S

Qpos -3

ACT (mm) 0.51(±0.04) + 0.0116(±0.007)A + 0.0040S

NT (mm) 2.11(±0.04) + 0.0030(±0.001)A + 0.0400S

PCT (mm) 0.33(±0.04) + 0.0082(±0.007)A + 0.0006S

TT (mm) 2.95 + 0.0228A + 0.045S

Rnucleus (mm) 3.05+ (0.5542 - 0.0091A)(1 -S/(15 -0.25 S))

Rlens (mm) 0.0069(±0.001)A + 4.35(±0.07)

ARlens (mm) 0.5542 -0.0091A

For the amplitude of accommodation (AA) we used an empirical law48 which gives the maximum accommodation response in dioptres as a function of age (A).

AA = 15 -0.25A (A3)

This expression states that the amplitude of accommodation declines by 0.25 D per year, predicting an amplitude of 0D for 60y.o. lenses. This expression has been further verified by other authors.1

According to the Helmholtz theory, the in vitro forcefree lens geometry corresponds to the fully accommodated state. We used an empirical linear equation obtained by Rosen et al.47 (see Table 5) to compute Rlens at each age of interest. To estimate the equatorial diameter of the lens in the unaccommodated state, we used an expression for the ARlens36 (see Table 5).

and green curves. Following,13 P1 is set at x = 2.5mm for the anterior part and x = 2 mm for the posterior surface. As previously explained, P2 is fixed by 9ant and 0pos, using the cited values.37 The computed values for C2,ant and Q2,ant were 0.1021 mm-1 and 1.1302, respectively. The expressions for z2ant, Q2,ant and C2,ant can be obtained similarly to for z1ant, Q1,ant and C1,ant. The computed values for C2,pos and Q2,pos were 0.0036 mm-1 and 53.5284, respectively.

Our model introduces an improvement with respect to the geometry proposed by Hermans. The interface between the anterior and posterior halves of the lens, instead being an equatorial plane, is modelled as a conicoid surface of revolution.45 We used the conicoid proposed by Navarro et al.23 From this curve a quadratic surface of revolution is obtained.

The interfaces between nucleus and cortex, both anterior and posterior, are assumed to be concentric with the outer surfaces of the cortex. No equatorial circumference is used for the equator of the nucleus. The total thickness (TT) of the lens is described by the anterior half lens thickness (tant) and the posterior half lens thickness (tpos). Fig. 1 shows these parameters and also the anterior cortex thickness (ACT), posterior cortex thickness (PCT) and nucleus thickness (NT) as the sum of anterior nucleus thickness (ANT) and posterior nucleus thickness (PNT).

Rosen et al.47 found a nearly constant ratio of 0.7 between anterior and posterior thicknesses, measured as the distances from the anterior and posterior poles to the equator. This factor, along with the expression for the lens thicknesses TT in Table 5 was used here to compute the thickness of nucleus and cortex for different ages and accommodations.

(ii) Age and accommodation

(iii) Geometry of the nucleus

To design the geometry of the nucleus we used in vivo experimental data obtained by Dubbelman et al.8 They observed that the anterior and posterior thicknesses of the cortex do not change appreciably with accommodation. Moreover, the growth rate of the cortex with aging is greater than that of the nucleus, which is consistent with the fact that new fibres are continuously formed in the outermost layer of the cortex. The expressions for ACT, NT, PCT and TT at a given age and state of accommodation are shown in Table 5.

The equatorial radius of the nucleus was estimated by the expression proposed by Brown3 for Rnucleus at one age and state of accommodation (see Table 5). Since the anterior and posterior cortex -nucleus interface and the outer surface of the lens are coaxial in the optical zone, the geometry of the nucleus is completed with two anterior and posterior curves which extend the central parts to the point of the medium surface situated at Rnucleus from the optical axis (see Fig. 1).

(iv) Thickness of the capsule

The capsular thickness varies with the radial position and grows with age.49 We followed the equations proposed by Burd et al.17 based on Fisher's data49 to estimate the capsular thickness as a function of the radial position, for different ages. The anterior capsular thickness is greater than the posterior. For the 30-year numerical model, the capsular thickness was assigned a value of 10.87 ^m at the anterior pole, increasing towards the equator (20 ^m), and then decreasing towards the posterior pole, up to a value of 3.159 ^m. For the 40-year numerical model, the values for the anterior and posterior pole were 13.37^m and 3.454^m, respectively, and for the 50-year numerical model, 15.88 ^m and 3.75 ^m.

The model described above is a parametric general model of the lens (nucleus and cortex) geometry. To obtain the shape of the lens at a given age and at an accommodative state, we need to particularize the values of all parameters involved. For the external geometry, we used the empirical expressions obtained by Dubbelman et al.5 (Table 5), where A is the age in years and S is the vergence of the stimulus.

References

1. Glasser A, Campbell M. Presbyopia and the optical changes in the human crystalline lens with age. Vis Res. 1998;38:209-229.

2. Ziebarth NM, Borja D, Arrieta E, et al. Role of the lens capsule on the mechanical accommodative response in a lens stretcher. Invest Ophthalmol Vis Sci. 2008;49:4490-4496.

3. Brown N. Change in shape and internal form of the lens of the eye on accommodation. Exp Eye Res. 1973;15:441-459.

4. Koretz J, Cook C, Kaufman P. Aging of the human lens: changes in lens shape upon accommodation and with accommodative loss. J Opt SocAm A-Opt Image Sci. 2002;19:144-151.

5. Dubbelman M, der Heijde GV, Weeber H. Change in shape of the aging human crystalline lens with accommodation. Vis Res. 2005;45:117-132.

6. Belaidi A, Pierscionek B. Modeling internal stress distributions in the human lens: can opponent theories coexist? J Vis. 2007;7:1-12.

7. Von Helmholtz H. Helmholtz's treatise on physiological optics. In: James PC, Southall, editors. Mechanism of accommodation, Vols 1, 2, Translation of 3rd German Voss: Hamburg; 1909.

8. Dubbelman M, der Heijde GV, Weeber H, Vrensen G. Changes in the internal structure of the human crystalline lens with age and accommodation. Vis Res. 2003;43:2363-2375.

9. David G, Pedrigi R, Heistand M, Humphrey J. Regional multiaxial mechanical properties of the porcine anterior lens capsule. J Biomech Eng. 2007;129:97-104.

10. Chien C, Huang T, Schachar R. Analysis of human crystalline lens accommodation. J Biomech. 2006;39:672-680.

11. Fisher R. The elastic constants of the human lens. J Physiol (Lond). 1971;212:147-180.

12. Fincham E. The mechanism of accommodation. Br J Ophthalmol. 1937;8(Suppl.):5-80.

13. Hermans E, Dubbelman M, der Heijde GV, Heethaar R. Change in the accommodative force on the lens of the human eye with age. Vision Res. 2008;48:119-126.

14. Pardue MT, Sivak J. Age-related changes in human ciliary muscle. Opt Vis Sci. 2000;77:202-210.

15. Stachs O, Martin H, Kirchhoff A, Stave J, Terwee T, Guthoff R. Monitoring accommodative ciliary muscle function using three-dimensional ultrasound. Graefe's Arch Clin Exp Ophthalmol. 2002;240:906-912.

16. Ljubimova D, Eriksson A, Bauer S. Aspects of eye accommodation evaluated by finite elements. Biomech Model Mechanobiol. 2007;7:139-150.

17. Burd H, Judge S, Cross J. Numerical modelling of the accommodating lens. Vis Res. 2002;42:2235-2251.

18. Fisher R. Elastic constants of the human lens capsule. J Physiol (Lond). 1969;201:1-19.

19. Krag S, Andreassen T. Biomechanical characteristics of the human anterior lens capsule in relation to age. Invest Ophthalmol Vis Sci. 1997;38:357-363.

20. Krag S, Andreassen T. Mechanical properties of the human lens capsule. Prog Retin Eye Res. 2003;22:749-767.

21. Pedrigi R, David G, Dziezye J, Humphrey J. Regional mechanical properties and stress analysis of the human anterior lens capsule. Vis Res. 2007;47:1781-1789.

22. Burd HJ. A structural constitutive model for the human lens capsule. Biomech Model Mechanobiol. 2009;8:217-231.

23. Navarro R, Palos F, Gonzalez LM. Adaptive model of the gradient index of the human lens. II. Optics of the accommodating aging lens. J Opt Soc Am A Opt Image Sci Vis. 2007;24: 2911-2920.

24. Maceo BM, Manns F, Borja D, et al. Contribution of the crystalline lens gradient refractive index to the accommodation amplitude in non-human primates: in vitro studies. J Vis. 2011;11:23, 1-13.

25. Lanchares E, Calvo B, Cristóbal JA, Doblaré M. Finite element simulation of arcuates for astigmatism correction. J Biomech. 2008;41:797-805.

26. Burd HJ, Wilde GS, Judge SJ. An improved spinning lens test to determine the stiffness of the human lens. Exp Eye Res. 2011;92:28-39.

27. Hermans EA, Pouwels PJW, Dubbelman M, Kuijer JPA, van der Heijde RGL, Heethar RM. Constant volume of the human lens and decrease in surface area of the capsular bag during accommodation: an MRI and Scheimpflug study. Invest Ophthalmol Vis Sci. 2009;50:281-289.

28. Kasprzak H. New approximation for the whole profile of the human crystalline lens. Ophthalmic Physiol Opt. 2000;20:31-43.

29. Flory PJ. Thermodynamic relations for high elastic materials. Trans Faraday Soc. 1961;57:829-838.

30. Holzapfel GA. Nonlinear solid mechanics. Chichester: (John Wiley & Sons Ltd.), Wiley; 2000.

31. Kuszak J, Zoltoski R, Tiedemann C. Development of lens sutures. IntJDev Biol. 2004;48:889-902.

32. Glasser A, Kaufman P. The mechanism of accommodation in primates. Ophthalmology. 1999;106:863-872.

33. Streeten B. Anatomy of the zonular apparatus. In: Tasman W, Jaeger W, eds. Duane's clinical ophthalmology. Hangerstown: Lippincott and Wilkins; 2003.

34. Ludwig K. Zonular apparatus, anatomy, biomechanichs and coupling to the lens. In: Guthoff R, Ludwig K, eds. Current aspects of human accommodation. Kaden: Heidelberg; 2001:71-92.

35. Hermans E, Dubbelman M, van der Heijde GV, Heethaar R. Estimating the external force acting on the human eye lens during accommodation by finite element modelling. Vis Res. 2006;46:3642-3650.

36. Strenk S, Semmlow J, Strenk L, et al. Age-related changes in human cliliary muscle and lens: a magnetic resonance imaging study. Invest Ophthalmol Vis Sci. 1999;40:1162-1169.

37. Schachar R, Huang T, Huang X. Mathematic proof of Schachar's hypothesis of accommodation. Ann Ophthalmol Clin. 1993;33:103-112.

38. Burd H, Wilde G, Judge S. Can reliable values of young's modulus be deduced from Fisher's spinning lens measurements? Vis Res. 2006;46:2235-2251.

39. Krag S, Andreassen T. Biomechanical measurements of the porcine lens capsule. Exp Eye Res. 1996;62:253-260.

40. Fisher R. The force of contraction of human ciliary muscle during accommodation. J Physiol (Lond). 1977;270:51-74.

41. Manns F, Parel JM, Denham D, et al. Optomechanical response of human and monkey lenses in a lens stretcher. Invest Ophthalmol Vis Sci. 2007;48:3260-3268.

42. Strenk SA, Strenk LM, Guo S. Magnetic resonance imaging of the anteroposterior position and thickness of the aging, accommodating, phakic and pseudophakic ciliary muscle. J Cataract Refract Surg. 2010;36:235-241.

43. Heys KR, Leigh-Cram S, Willis-Truscott RJ. Massive increase in the stiffness of the human lens nucleus with age: the basis for presbyopia? Mol Vis. 2004;10:956-963.

44. Weeber HA, Eckert G, Pechhold W, der Heijde RGLV. Stiffness gradient in the crystalline lens. Graefe's Arch Clin Exp Ophthalmol. 2007;245:1357-1366.

45. Jones C, Atchison D, Meder R, Pope J. Refractive index distribution and optical properties of the isolated human lens measured using magnetic resonance imaging (MRI). Vis Res. 2005;45:2352-2366.

46. Dubbelman M, der Heijde GV. The shape of the aging human lens: curvature, equivalent refractive index and the lens paradox. Vis Res. 2001;41:1867-1877.

47. Rosen A, Denham DB, Fernandez V, et al. In vitro dimensions and curvatures of human lenses. Vis Res. 2006;46:1002-1009.

48. Weale R. The senescence of human vision. Oxford: Oxford University Press; 1992.

49. Fisher R, Pettet B. The postnatal growth of the capsule of the human crystalline lens. J Anat. 1972;112:207-214.