A validated finite element model to reproduce Helmholtz’s theory of accommodation: a powerful tool to investigate presbyopia

To reproduce human in vivo accommodation numerically. For that purpose, a finite element model specific for a 29‐year‐old subject was designed. Once the proposed numerical model was validated, the decrease in accommodative amplitude with age was simulated according to data available in the literature.

division is somewhat arbitrary. There is a gradual transition from the outermost longitudinal muscle fibres through the radial fibres to the innermost circular muscle fibres, with some intermingling of the different fibre types. [4][5][6] Contraction of the entire ciliary muscle as a whole pulls the anterior choroid forward, moving the apex of the ciliary processes towards the lens equator, and serves the primary function of releasing resting zonular tension at the lens equator to allow accommodation. 4,5 During contraction of the ciliary muscle, the circular portion of the ciliary muscle tends to increase in thickness, whilst the radial and longitudinal portions decrease in thickness.
The crystalline lens is the ocular element that alters the focal length of the eye. During accommodation, the human lens undergoes several changes: its diameter decreases, its thickness increases, the anterior and posterior surfaces of the lens move anteriorly and posteriorly, respectively and the curvatures of the anterior and posterior surfaces of the lens increase. Moreover, the thickness of the lens nucleus increases, but without a significant change in thickness of the cortex. 4,7 Several authors have designed finite element (FE) models to explain and understand accommodation. [8][9][10][11][12][13][14][15] However, most of these biomechanical models are based on accommodated geometry, and reproduced the accommodation process in reverse by applying the forces exerted by the relaxation of the ciliary muscle to reach the non-accommodated state. These FE models are a powerful tool to understand the opto-mechanical crystalline lens changes. However, in vivo accommodation cannot be reproduced faithfully through these tests, and thus some key aspects of accommodation which involves other tissues (ciliary muscle, sclera) cannot be analysed properly. In contrast to these previous studies, the reference configuration of our FE model is the non-accommodated state, which is the resting state of the physiological system. We propose a 3-dimensional FE model that incorporates contraction of the ciliary muscle to reproduce in vivo human accommodation. This FE model is standardised for a 29-year-old subject as there is sufficient comparative and physiological data in the literature to develop the model and evaluate it with experimental data. Ciliary muscle contraction is simulated with an electro-chemical-mechanical continuum model. 16,17 Moreover, to allow the lens to change its shape after release of zonular tension, the initial stresses of the lens capsule are included in the model. The proposed FE model was validated by comparison with the experimental optical and main biometric changes during human in vivo accommodation provided by Ramasubramanian and Glasser. 18,19 To achieve this aim, the ciliary muscle contraction was calibrated by comparing the morphological changes with the recent work of Wagner et al., 20 who evaluated ciliary muscle thickness (CMT) profiles derived from optical coherence tomography (OCT) images during in vivo human accommodation.
With the proposed FE model, the physiological decrease in accommodation amplitude with age was evaluated and compared with experimental data. [21][22][23] A sensitivity analysis of the mechanical properties of the crystalline lens as a function of age was performed to observe the most influential mechanical factors in the development of presbyopia. [24][25][26] Justification of the geometry and the 3-dimensional FE model are presented in the Method section, followed by a description of the constitutive laws of all materials involved, including the formulation of the active behaviour of the ciliary muscle and the mechanical properties used to analyse the decrease in accommodation amplitude with age. Methodology is provided to analyse the optical and biometric changes during accommodation. The results of these changes are presented together with some calibration parameters of the muscle activation. Once the ciliary muscle is engaged, it is coupled with the accommodation system to simulate the change in shape of the lens, and thus the refractive power change. An analysis of the decrease in accommodation amplitude with age is shown.

Finite element model
To obtain the most realistic representation possible, this model of the accommodation system comprises the most relevant components in the accommodation process: the shell of the eye composed of the cornea and sclera; iris; lens nucleus, cortex and capsule; zonules of Zinn divided into anterior and equatorial zonules and the ciliary body, which is composed of the ciliary muscle and processes. Figure 1 shows the dimensions of the different components of the FE model in the non-accommodated state standardised for a 29-year-old subject. Lens dimensions were obtained from the OCT study of Chang et al., 27 with a radius of curvature of the anterior (r a ) and posterior (r p ) lens surfaces of 11.00 and 6.10 mm, respectively, sagittal thickness of the lens (T L ) of 3.44 mm and an anterior chamber depth (ACD) of 3.03 mm. As the lens equatorial diameter cannot be obtained by OCT, it was extracted from the MRI work of Kasthurirangan et al. 29 with the value of 9.04 mm being obtained. The anterior and posterior thickness of the lens cortex were 0.80 mm and 0.50 mm, respectively, with

Key Points
• A novel finite element model that reproduces ocular accommodation is presented. The decrease in accommodative amplitude with age was reproduced using the proposed model. • The main cause of presbyopia is the stiffening of the lens nucleus with age. • The numerical model can help to understand human accommodation.
an anterior and posterior surface radius of curvature of 4.50 and 4.00 mm, respectively. 30,31 Zonules were anchored 0.50 mm in the anterior capsule based on Bernal et al., 32 and the ciliary body ring diameter was 11.90 mm. The averaged ciliary muscle geometry of emmetropic eyes reported by Wagner et al. 20 was used in our model. In their study, the geometry was provided by a CMT profile, which was obtained from OCT images following the methodology of Straßer et al. 33 Therefore, an in-house code in MATLAB 2020b (MathWorks, mathworks.com) with the inverse methodology 33 was created and implemented to obtain the geometry from the CMT profile.
The ciliary muscle was arranged into three families of muscle fibres, i.e., circular, radial and longitudinal, according to the distribution reported by Pardue and Sivak. 6 The volume percentage of the circular, radial and longitudinal muscle fibres is 12%, 33% and 55% respectively, and their distribution is shown in Figure 2. F I G U R E 1 Dimensions of the finite element (FE) model. It includes the lens, the zonules of Zinn (blue), the ciliary muscle (red), the iris and ciliary processes (light pink) as well as the sclera and cornea (white). A detailed view reflecting capsule thickness is shown on the left. Due to the lack of available data in the literature, the thickness distribution of the lens capsule used was taken from a 36-year-old subject 28 F I G U R E 2 3D finite element (FE) model of the accommodative system: lens, zonules, ciliary muscle, iris, sclera and cornea. The arrangement and orientation of the ciliary muscle fibres, longitudinal (yellow), radial (orange) and circumferential (red), is shown on the right. The FE model contained 402,373 elements Considering that the material model allows only one preferential direction, the orientation of the radial fibres was defined as a transition from the circular to the longitudinal fibres; see Figure 2.
Regarding the iris, a trabecular-iris angle of 40 • and a thickness at 750 µm from the scleral spur of 0.40 mm were considered for visual effects since its contraction was not considered. 34,35 Finally, the shell of the eye was designed according to the Le Grand 36 model. The thickness of the sclera and cornea were readjusted for a 29 year old subject 37 ; see Figure 1.
Considering all the above-mentioned factors regarding the geometry, a 3-dimensional FE mesh of the anterior half of the ocular globe was designed. Due to the symmetry, only a quarter of the geometry was considered, ensuring the corresponding symmetry boundary conditions. Axial displacements were restrained at the bottom scleral surface. Moreover, the exterior surfaces of the sclera were fixed by the extraocular muscles. The general-purpose finite element code Abaqus FEA version 14.1 (Simulia, 3ds.com) was used to build the FE model. A mesh sensitivity analysis was performed in order to establish the final mesh size; see Figure 2. The lens nucleus and cortex were considered solid bodies and were meshed with 96,640 hybrid linear hexahedral (C3D8H) elements; the capsule was meshed with 5290 4-node membrane (M4) elements; 546 zonules were modelled using Abaqus connector elements, whose nonlinear behaviour is explained below; the ciliary muscle was meshed with 94,809 linear hexahedral (C3D8) elements and the ciliary body was meshed with 201,110 C3D8 elements and 3978 linear wedge (C3D6) elements.
The fluids filling the eye, namely the aqueous and vitreous humours, were not included in the model as we considered that there is a homogeneous pressure inside the eye during accommodation. Due to the boundary conditions applied in the sclera, the eyeball structure is stable.

Ciliary muscle
The process of smooth muscle contraction can be initiated by mechanical, electrical or chemical stimuli. The contractile mechanism involves several signal transduction pathways, all of which converge to increase the intracellular calcium concentration, resulting in phosphorylation of myosin. 38 Regarding the ciliary muscle, the action of the neurotransmitter acetylcholine on post-synaptic muscarinic receptors controls muscle activation. This system forms part of the parasympathetic branch of the autonomic nervous system. On the other hand, an inhibitory effect is driven by the sympathetic system, which is mediated principally by the action of the neurotransmitter noradrenaline on β 2 receptors. 39 In this work, the passive and active finite strain response of the ciliary muscle was simulated within the framework of continuum mechanics using a common methodology based on postulating the existence of a strain energy density function (SEF). 16,17 This function depends on a series of state variables related to the deformation of the active and passive elements, the contraction level and the muscle fibre arrangement: C = F T F is the right Cauchy-Green deformation tensor with F the deformation gradient. C e is an equivalent deformation tensor associated with the deformation of the elastic components in the muscle, and it arises from a two-step formulation of the contraction process. 16 From this formulation, a decomposition of the strain deformation gradient is proposed as F = F e F a . Therefore, the elastic deformation gradient can be defined as C e = F T e F e . a is the contraction or stretch of the active component of the muscle fibres that forms the active strain deformation gradient F a assuming N).
F I G U R E 3 (a) Activation function of the ciliary muscle and (b) mechanical properties of the crystalline lens used in this study. There is a circle in the lens nucleus and cortex plot indicating the values reported by Wilde et al. 26 The mechanical properties of the lens capsule were characterised by Krag and Andreassen 41,42 incompressibility. 16 The structural tensor N = n 0 ⊗ n 0 defines the anisotropy of the muscle due to the direction of muscular fibres n 0 .
The SEF is decoupled into volume-changing and volume-preserving parts in order to handle the quasiincompressibility constraint. Furthermore, the deviatoric part is divided into a passive contribution due to the collagen and elastin, Ψ p , and an active contribution associated with the muscular fibres, Ψ a . Thus, the total strain energy function Ψ can be expressed as follows: a have been applied with J = det(F) being the Jacobian of the transformation. Equation 2 can be particularised for the behaviour of the ciliary muscle and formulated in a more proper way for computational purposes as: where f is the force-stretch relationship that account for the actin and myosin overlap, f act is the activation function and Ψ a is the active SEF. The passive SEF, Ψ p , is defined as a function of the invariants of the isochoric right Cauchy-Green tensor: where I 1 and I 2 are the first and second modified strain invariants and I 4 is the invariant related to the anisotropy of the passive response. According to the SEF proposed by Calvo et al., 40 the passive response Ψ p can be written: The active SEF, Ψ a , associated with the active response, and consequently with the actin-myosin interaction, is expressed as the product of two functions that scale the energy stored due to contraction in the elastic element Ψ � a . 16 The influence of filament overlap on the active response of the muscle f λ is formulated in terms of the muscle fibre stretch a , the optimum stretch λ opt and the coefficient ξ: The activation function f act is assumed to represent the excitation input signal that triggers the contraction independently of its origin. In this work, the product of two squared hyperbolic tangents has been selected as: where s 1 and s 2 regulate the initial and final slope of the function, respectively, t is the time variable, and t i and t s define the start of the activation and the stimulus duration, respectively. Figure 3a shows the activation function used in the study. The energy stored in the cross-bridges is expressed in terms of the invariant associated to C e in the direction of the muscle fibres n 0 : where P 0 is a proportionality factor related to the maximum active stress due to muscle contraction. 16 Due to the non-recoverable deformation that occurs in the active part of the muscle fibres (its initial shape is only recovered with the help of the elastic components), the model incorporates a dissipation energy rate associated with active stresses and strains. 16 Using the second law of thermodynamics in the shape of the Clausius-Planck inequality and neglecting the thermal dissipation rate, some of the power produced internally is stored while another portion is dissipated. 16 From this expression, two constitutive relationships can be derived. The first one is used to obtain the classical relation between the derivative of the so-called strain energy density function with respect to the elastic strain tensor and the second the Piola-Kirchhoff stress tensor. The second relationship allows establishment of the evolution of the contraction velocity a as a function of the same strain energy and stresses 16 : where P a = −̃ P 0 f f act is the active stress, with ̃ being a friction parameter that considers the relative sliding speed between actin and myosin. The parameter C is defined as: where v 0 is associated with the initial contraction velocity. A set of preliminary analyses, explained in the Results section, was performed to obtain the parameters P 0 and v 0 of the active behaviour of the ciliary muscle; see Tables 1 and 2.

Constitutive laws of the remaining tissues
The whole lens was modelled with an elastic behaviour. The lens capsule was modelled via membrane elements, assuming pure traction behaviour. In turn, the lens capsule was divided into anterior and posterior regions, characterising their mechanical properties differently based on the work of Krag and Andreassen, 41,42 see Table 3.
In contrast with those of the lens capsule, the mechanical properties of the lens nucleus and cortex are unclear and somewhat controversial. Burd et al. 43 reported that Fisher's 24 spinning lens measurements might not be (C e , a , N), (7) f act = (tanh 2 (s 1 (t − t i ))tanh 2 (s 2 (t − (t i + t s ))), reliable, and relatively few subsequent studies regarding the mechanical properties of the internal lens have been performed, 25,26,44 with slight differences being observed between them. We determined the values reported by Wilde et al., 26 taking into account the stress-strain relationship E = 3G, with G being the shear modulus, can be considered for isotropic elastic linear and incompressible materials. Due to the lack of data, we extrapolated values from a 33-years-old lens to a 29-years-old one; see Figure 3b. The cornea, sclera and ciliary processes were modelled through an isotropic, quasi-incompressible Neo-Hookean hyperelastic model. The material properties for the 29-year-old subject were chosen from the data available in the literature, as displayed in Table 3.
The zonules were modelled with linear connector elements as Helmholtz's accommodation theory postulates, 1,2 i.e., they are stretched in the non-accommodated state, which is the reference configuration. This was considered as an initial force value F 0 in the force-displacement equation that defines their behaviour: with F 0 = 50 mN, k = 100 mN mm −1 and u being the displacement in the direction of the zonule (element). 45 No compression was considered.

Analysis procedure
The most common theory of accommodation, proposed by Helmholtz, 1 was considered, and the contraction and relaxation of the ciliary muscle were reproduced in order to simulate the distant (non-accommodated) and near vision (accommodated) states.
The reference configuration of the FE model corresponds to the non-accommodated geometry, where the lens capsule and the Zinn zonules are not stress free. The decapsulated lens substance assumes a maximally unaccommodated form after removing the capsule, and thus the decapsulated lens was initially stress free. 49 The initial stress distribution ( 0 ij ) was obtained with an extra analysis, stretching the lens from the accommodated to the non-accommodated state; see Figure 4a. The stress results at the integration points were mapped to the new lens capsule geometry.
All dynamic optical and biometric lens measurements of the accommodation mechanism were evaluated for each increment of the evolution process modelled in the simulation. For that purpose, displacement of the lens nodes and nucleus contour were registered by an URDFIL Abaqus subroutine and post-processed with MATLAB R2020b.
To measure the actual accommodative change (difference between the non-accommodated and corresponding accommodated states), the central optical power of the eye (COP Eye ) was evaluated throughout the simulation using Gullstrand's eye formula: with the central optical powers of the lens (COP L ) and the cornea (COP C ) 8 : where n a = 1.336 is the refractive index of the aqueous humour, n L = 1.42 is the estimated overall refractive index of the lens, n c = 1.376 is the refractive index of the cornea and n air = 1.00 is the refractive index of air. The remaining biometric terms, r a , r p , T L , ACD, rc a , rc p and T c , defined in Figure 1, vary throughout the simulation. The power of the cornea was COP L = n L − n a r a + n L − n a r p − T L (n L − n a ) 2 r p r a n L ,  kept constant throughout the simulation. The radii of curvature were calculated throughout the stretching process through an equation for a conic section having an apex at the origin and tangent to the y axis.
To obtain the corresponding radius of curvature (r) and the conic constant (K), a non-linear system of equations formed by the coordinates of the nodes (x, y) was solved. The goodness of the fit was R 2 > 99%. The standardised design of the eye had 58 dioptres (D) in the farvision state.
Lastly, the deformed geometry of the ciliary muscle was also measured during accommodation for each increment of the evolution process modelled in the simulation.
Straßer's 33 methodology was used to measure the CMT profile along the overall length of the muscle, which consists of measuring the distance between the sclera and the outer boundary of the ciliary muscle from perpendicular rays originating in the exterior surface of the sclera; see Figure 4b.

R ESULTS
During ciliary body contraction, the zonules release tension to allow the lens to round up to its accommodated state, see Figure 5. The FE solution, corresponding to an actual accommodation amplitude of 5.82 D, is shown on the left whilst the reference configuration, the nonaccommodated state is depicted on the right to show (15) y 2 − 2rx + (K + 1)x 2 = 0. To evaluate the accuracy of the proposed FE model of accommodation, we compared the dynamic changes produced in the ciliary muscle and lens during accommodation with experimental values collected from the literature. Figure 6 shows the morphological comparison between the numerical and experimental ciliary muscle contraction. 20 The experimental ciliary muscle contraction corresponds to an accommodation stimulus of 4 D. 20 The actual accommodative response is lower than the accommodation stimulus. 50,51 Figure 6a presents the contour of both the non-deformed and deformed geometries, while Figure 6b depicts comparison between the numerical and experimental CMT profiles. 20 Once the numerical and experimental contractions of the ciliary muscle were compared, we calibrated the main parameters of active behaviour of the ciliary muscle, maximum active stress (P 0 ) and the initial contraction velocity (v 0 ) in order to reproduce the contraction and relaxation of the ciliary muscle. The maximum active stress is related to the maximum change in the ciliary muscle ring diameter (ϕ CM ), whereas the initial contraction velocity (v 0 ) is related to the temporal performance of the muscle contraction. Due to the lack of data in the literature, we calibrated the muscle properties based on a virtual case. 52 Glasser et al. 52 provided electrical innervation to quantify the changes in lens diameter in rhesus monkeys. The virtual case attempted to replicate this situation but by quantifying the change in the ciliary body ring diameter. An accommodation square wave was implemented for 4 s by simulating electrical innervation. For the maximum electrical innervation, a maximum change in Δϕ CM of 0.70 mm was presumably obtained. The methodology proposed would be the same for an actual case. These properties were calibrated by means of a sensitivity analysis observing Δϕ CM . First, the maximum active stress (P 0 ) was obtained at the maximum response amplitude, 0.70 mm. 52 To do so, P 0 ranged from 0.50 to 2.00 MPa, with v 0 = 3.00 based on a previous study. 17 Figure 7a shows that a Δϕ CM of 0.70 mm is obtained for P 0 = 1.50 MPa. Once the maximum active stress was achieved, sensitivity analysis of the initial contraction velocity, v 0 , was performed to calibrate the temporal performance (activation and deactivation) of the ciliary muscle; see Figure 7b.

Coupling of the ciliary muscle and lens
Before analysing the dynamic biometric and optical measurements of the lens, we will describe the joint performance of the ciliary muscle and lens. Figure 8 presents three conditions of the ciliary muscle and lens based upon different accommodated states, namely, at 0.00, 0.75 and 1.50 s, which correspond to different ciliary body ring diameters (Δϕ CM ). The numerical FE model presented a maximum accommodation of 5.82 D.
The accommodative response followed a linear relationship with the Δϕ CM , with an increase of 5.82 D and an increase in lens thickness of 0.44 mm for a reduction of 0.66 mm in Δϕ CM or 67.54 µm in ∆CMT max . Due to the way in which the CMT profiles are obtained with the perpendicular rays from the sclera, the changes observed in the CMT profiles are of lower magnitude.

Dynamic optical and biometric lens measurements
The experimental in vivo results in humans obtained by Ramasubramanian and Glasser 18,19 were used for comparison with the numerical dynamic optical and biometric lens measurements. Figure 9a shows the changes in the anterior and posterior lens radius of curvature against the actual accommodative change of the eye. This graph shows where the optical change is happening. Both surface curvatures fit perfectly according to the experimental data. However, the anterior radius of curvature does not fit as well at the end of accommodation. This result may be due to a slight inconsistency in the definition of the initial stresses of the lens capsule compared with the physiological values. Figure 9b shows the changes in lens thickness (∆T L ) and ACD as a function of the accommodative change. There is a ∆T L of 0.44 mm and ∆ACD of 0.31 mm for an accommodative change of 5.82 D. The orange line shows the change in the numerical nucleus thickness, which influences approximately 80% of the lens thickness, as previously reported. 7,31 As a consequence of the increase in lens thickness, the lens moves anteriorly and posteriorly; see Figure 9c. The anterior lens movement is notably higher, approaching 0.31 mm anteriorly and 0.13 mm posteriorly for an accommodative change of 5.82 D.
Lens diameter (see Figure 9d) was compared with ex vivo experimental tests due to the difficulty of measuring the lens diameter in vivo. Data for a 29-year-old individual were obtained by interpolation from human groups aged 14 to 39.5 years as reported by Manns et al. 53 The numerical data fit perfectly with the interpolated data for a 29-year-old. All optical and biometric eye measurements were within the deviation reported by Ramasubramanian and Glasser. 18,19 Analysing the influence of lens mechanical properties in presbyopia The mechanical properties of the lens are fundamentally related to how it changes its shape, 8,10,26 above other factors such as the lens geometry. 15 Therefore, we analysed the decrease in accommodation amplitude by varying the mechanical properties, i.e., the lens capsule, nucleus and cortex, while we kept the remaining tissues of the FE model proposed for the 29-year-old subject constant. Thus, in this hypothetical case, we kept the same lens geometry and the initial stresses of the lens capsule calculated for the validated FE model, as reflected in Figure 4a. We analysed three different cases: Case A, corresponding to all mechanical properties being fixed for a 29-year-old subject, except for the mechanical properties of the lens nucleus and cortex; Case B is the same as A, but with the change in the anterior and posterior lens capsule and Case C is a combination of Cases A and B. Both the mechanical properties of the lens capsule, nucleus and cortex were modified according to aging experimental data, 25,26 see Figure 3b. Figure 10 shows the accommodation amplitude data as a function of age for three different numerical cases against the experimental data. [21][22][23] Case B, changing the mechanical properties of the lens capsule, presents a small decrease in accommodation amplitude, from 5.82 D for a 29-year-old subject to 4.24 D for a 51-year-old. By contrast, Cases A and C, where the mechanical properties of the lens nucleus and cortex change with age, present decreasing amplitudes similar to the experimental data. [21][22][23]

DISCUSSION
The main goal of this study was to develop an FE model that reproduced the in vivo accommodation process in humans. The proposed model reproduced Helmholtz's theory of accommodation and compared successfully with quality experimental data. [18][19][20] We accomplished this goal by quantifying contraction of the ciliary muscle and the optical and biometric measurements of accommodation, see Figures 6 and 9.
The resulting movement of the ciliary muscle comes from simultaneous contraction of the circular, radial and longitudinal muscle fibres. We faithfully reproduced the morphological changes in the ciliary muscle contraction of an emmetropic eye 20 by means of a computational eye model, which includes continuum and structural mechanics. 16,17 The morphological changes in the numerical and experimental ciliary muscle were similar, contracting towards the lens equator to release the resting tension on the zonular fibres, see Figure 6. This action was produced mainly by the circular muscle fibres, which governed ciliary muscle contraction. The numerical ciliary muscle became F I G U R E 9 (a) Change in the anterior (r a ) and posterior (r p ) lens radius of curvature, (b) lens thickness (T L ) and anterior chamber distance (ACD), (c) anterior and posterior lens movement and (d) lens diameter, as functions of the change in accommodation slightly thinner below the apex of the ciliary muscle, in contrast with the experimental values. 20 However, detecting small changes is highly difficult experimentally because the OCT measurement technique depends upon the patient's sclera radius of curvature and on the operator. 33 The longitudinal muscle fibres appear to have little influence on ciliary muscle contraction. It is true that there is a slight decrease in the length of the ciliary muscle as a result of pulling the posterior part of the ciliary muscle forward to the scleral spur. However, the circular fibres contract slightly upwards intrinsically due to the geometry of the muscle itself, and not due to the contraction of the longitudinal fibres. This might be due to the boundary condition applied in the sclera, which made the interaction between the ciliary muscle and the sclera too stiff. On the other hand, Glasser reported 4 that the longitudinal fibres seem to brace the system rapidly so that the contraction of the inner portion is most effective; this is consistent with the numerical results obtained in the present study.
One weakness of the FE model is that we modelled the radial muscle fibres as a mixture of circular and longitudinal muscle fibres due to limitations in numerically designing the radial fibres. The radial muscle fibres are supposed to perform a combined effect of the circular and longitudinal fibres. 4 Another weakness is that the iris constriction, which increases the depth of focus of the eye, was not simulated because it is not relevant from a mechanical point of view.
Two of the limitations of the FE model arose from a lack of data in the literature. First, we considered that all muscular fibres were activated equally, although the force exerted, which depends on the force-stretch relationship (f λ ), was higher for the circular fibres. Second, the muscle properties were calibrated according to a virtual case based on the electrical innervation provided to rhesus monkeys by Glasser et al. 52 Nevertheless, the ciliary muscle could be updated according to new investigations with the proposed methodology.
The joint performance of the ciliary muscle and lens with respect to accommodation was analysed numerically. We obtained a linear response with an increase in accommodation of 5.82 D and an increase in T L of 0.44 mm for a reduction of 0.66 mm in the ciliary muscle ring diameter (Δϕ CM ), or 67.54 µm in the maximum ciliary muscle thickness (∆CMT max ). These values were within the range reported in the literature. Strenk et al. 54 reported Δϕ CM values ranging from 0.50 to 1.00 mm for an 8 D accommodative stimulus. Ruggeri et al. 55  Human accommodation is similar to that of other hominids. However, depending on the species and age, accommodation occurs in slightly different ways; for example, lens thickness and diameter change in humans is reduced compared with rhesus monkeys for the same accommodation amplitude. 3,19,53,56 For this reason, validating the main optical and biometric lens measurements was complicated and meaningful.
In this study, we present a numerical approach standardised for a 29-year-old individual, accurately reproducing human accommodation. The optical and biometric measurements of the lens were compared with in vivo accommodation dynamic measurements by Ramasubramanian and Glasser,18,19 and all the numerical data presented were within the range of the experimental measurements. To achieve this aim, we previously reported that the mechanical properties of the crystalline lens are key to how the lens shape changes. 15 The presented model indicates that several factors such as lens mechanical properties, the zonular arrangement and capsule thickness distribution must be included to simulate the accommodation of humans.
We hypothesise that the stresses in the non-accommodated lens capsule should vary little over time for two reasons. Were they to increase, then accommodation would become higher over time while if they were to decrease, then the stiffness added with age would cause an earlier appearance of presbyopia. With this hypothesis, we determined that stiffening of the lens nucleus was mainly responsible for the decrease in accommodation amplitude (see Figure 10), rather than the increased stiffness in the lens capsule with increasing age. This conclusion is supported by our previous work, 15 where we observed that the stiffness ratio between the lens nucleus and cortex is key to how the lens changes its shape. In contrast to Wilkes et al., 12 our FE model can explain some causes of presbyopia over a wider age range. Additional studies considering more factors involved in accommodation as well as mechanical data regarding the internal properties of the lens are needed to understand presbyopia further. In contrast to other FE models that stretched the lens, simulating ex vivo conditions, [8][9][10][11][12][13][14] we designed our reference configuration using non-accommodated geometry. Designing the 3D FE model in this way allowed us to take lens movement into consideration, which is essential because changes in lens curvature do not reflect all significant accommodative changes. 57 In the numerical validation, 0.97 D (16%) of the 5.82 D of the accommodation amplitude was produced by lens movement. Furthermore, this model could be compared with novel research from the literature [18][19][20] and provides the necessary background to analyse and evaluate fundamental aspects of accommodation, as well as the detailed geometry of a 29-year-old individual for other research studies.