Simulation of Active Eye Motion Using Finite Element Modelling

Dynamic simulation of human eye movements can significantly improve our understanding about the visuomotor system and also may be used for clinical applications. In this study we have developed a 3D finite element model of the eye orbit. The major novel issues which are considered in this work are using dynamic finite element method for properly mass distribution modelling during simulation and active motion simulation by employing muscle activation parameter. In this modelling, extraocular muscles are modelled using a continuum constitutive hyper‐elastic energy function. Simulation results of force duction and active motion are presented and compared with the experimental results to verify the accuracy and reliability of the proposed modelling. Then, orbital fat deformation during eye motion is studied through various simulations and compared with experimental results. Finally, the proposed eye model and available experimental data are exploited to investigate more details about eye suspension deformation and extraocular muscles characteristics.


Introduction
Many parameters such as muscle neural excitation, muscle dimension and muscle insertion position to eye ball greatly influence eye motion. Hence, curing visuomotor disorders such as strabismus by surgery is complicated. An accurate and efficient eye modeling may facilitate curing eye disorders. Surgeons can use the model to simulate how to cure these disorders before surgery and reduce surgical inaccuracies and risks. Furthermore, it can provide us a better understanding of the eye motion. Meanwhile, the model should include the main eye tissues affecting eye movements in order to avoid high computational cost.
Eye modelling is an interesting problem for researchers. Descartes pioneering work for eye modelling in 1630 see Descartes 1972 has been continued by Listing, Ruete, Donder, Helmholtz and many others until now details can be found in works by Simonsz & Spekreijse 1990, Donders 1864and Fetter et al. 1997. Simonsz & Spekreijse 1996, Clement 1982and Miller & Robinson 1984 suggested some computerized model for studying eye behaviour. Recently Ghosh et al. 2009 employed quaternions andBayro-Corrochano 2003 used geometric algebra for eye modelling with lumped parameters. By using lumped parameters, mass distribution and some other details in eye modelling are neglected. Hence, it is hard or even impossible to use these models for accurate investigation of mechanical characteristics affecting eye motion.
"Strands" are used by Wei et al. 2010 as a new modelling approach for eye. By using this method, more emphasis is devoted on muscles model by more accurate mass distribution modelling. Therefore, it is possible to obtain more realistic behaviour of both the extraocular muscles EOMs and eye motion. Pai 2010 showed that using the "strands" and "continuum mechanics" are two reliable solutions for properly representing parts' inertia and also can account for variations of inertia distribution during the motion. Continuum-mechanics method used for modelling eye orbit by Schutte et al. 2006 . They used finite element FE model to have a more accurate modelling of eye. However, the model presented by Schutte et al. 2006 has some simplification which affects the model performance. On previous FE model muscles were assumed to be orthotropic homogenous elastic tissues and other eye parts were homogenous isotropic elastic tissues. Muscle activation was not considered in muscle modeling and the simulation was static without encountering parts' mass . Eye motion simulation was done by considering thermo elastic behavior for the muscle and varying the muscles' temperature see more in Schutte et al. 2006 .
Abbas Karami a * Mohammad Eghtesad a ** In the present work, eye motion is simulated by using dynamic FE approach and a hyper-elastic model is used for muscle.
Muscle is a fundamental part affecting the eye motion. Eye displacement magnitude, velocity and direction are functions of muscles activation variation. To the best of our knowledge, this is the first work that employs an activepassive hyper-elastic model for EOM. Many different hyper-elastic constitutive equations have been proposed and used for modelling muscle behaviour such as Paetsch et al. 2012, Blemker et al. 2005, Grasa et al. 2011, Tang et al. 2009and Ehret et al. 2011 . Most of these models are proper for static or quasi-static simulations. In this work, we have exploited a hyper-elastic constitutive model which is proper for dynamic simulation and it is compatible with EOMs behaviour proposed by Ehret et al. 2011 .
The purpose of this study is simulating eye active motion with realistic properties and investigate various details in eye orbit suspension and muscle. To this end, we take into account as many parts of eye orbit as required in our FE model. Orbit geometry, eye-ball position in orbital fat and also mechanical characteristics of the eye orbital tissues may have different effects on the resulted eye motion. Many of these features are not studied till now. Meanwhile reliable available experimental data about the eye parts mechanical characteristics are limited details can be found in works by Schoemaker et al. 2006, Robinson et al. 1969and Yoo et al. 2011 . We are intended to obtain some information which are not convenient through the conventional experiments through our model. Previously Sigal & Ethier 2009, Blemker et al. 2005, Stavness et al. 2012and Stitzel et al. 2002 used FE models as a complementary to experimental data in biomechanical problems.
This article is organized as follows. In sections 2 and 3, we will take a look on the key elements which are needed for eye modelling; eye orbit geometry, material properties and mechanical characteristics. The FE modelling simulations will be reported in section 4. In this section two common tests for eye are simulated and their results are compared and confirmed with the available experimental data. Then, orbital fat deformation in our model is compared with experimental data. Finally, orbital fat elastic modulus and muscle pre-tension effect on eye motion are investigated through various simulations. This article concludes in section 5.

Orbit geometry
Eyeball and other parts of human eye orbit have a complex arrangement very close to each other. This intricate situation has compelled many researchers to ignore some parts that seem to be less effective on eye motion. Hence, models may obtained which are efficient and convenient for analysis; also accurate and useful in results see Figure  1 .
For modelling the eye geometry, we will use published data in the recent works such as ones by Schutte et al. 2006and Raab et al. 2010 In our model, the eye ball, its suspension and four rectus muscles are considered see Figure 2 . Optic nerve will be ignored since in comparison to orbital fat it seems to have little effect on eye's motion more on optic nerve mechanical characteristics is reported by Sigal & Ethier 2009 . In order to accurately analyse the eye motion, we must take into account the necessary details about muscle direction. EOMs tangency to the eye ball after the insertion position is considered based on the arc contacts that are reported in major anatomy references see the book by Chalam et al., 2010Chalam et al., -2011 . Various details are extracted from MRI reports and employed for preparing eye geometry. The resulted geometry with tetrahedral elements can be seen in Figure 2. Contacting surfaces in the model are simulated by using thin low elasticity layers between the contacting parts. Maxilla which is shown in Figure 1 is the fixed part in our FE model and the outer surface of the yellow elements on part "e" and "f" of Figure 2 are in contact with that.

Material properties and mechanical characteristics
Precise material description is required for accurate FE model. There are a few works on mechanical properties of different eye parts experimentally or numerically . Most of the available experiments are limited to the cases of the dead animal eyes see for example the work by Schoemaker et al. 2006 . Only a few of the experiments are done on live human eye such as Robinson et al. 1969 . The dead organ characteristics vary based on how they are preserved in laboratory till the experiment is done see Schoemaker et al. 2006 . On the other hand, for parts such as orbital fat, because of its position and also its delicacy, it is too hard to carry out the experiments and finding mechanical characteristics. These issues make the eye modelling not only demanding but also challenging.
For the eye modelling, one of the major issues is the muscle model which must appropriately simulate both active and passive behaviours of the EOMs. Meanwhile, other parts such as eye ball and eye orbit suspension should be considered accurately according to reliable experiments and observations. The major part of the eye orbit is filled with orbital fat and most of the eye ball surface is in contact with it. In order to have a more realistic model, two parts in eye suspension with different elastic modulus, but with the same density of 1000 kg/m 3 , is considered. This assumption is done since orbit surgeons reported different fat cell sizes for the intraconal retrobulbar fat the fat inside the muscle cone and the extraconal retrobulbar fat the fat outside the muscle cone Schutte et al. 2006 . Viscosity of the intraconal and extraconal retrobulbar fats is ignored as a result of low viscous shear modulus which is reported through the experiments by Schoemaker et al. 2006 . The mechanical properties assumption for these two parts are reported in Table 1. These assumptions are employed from Schutte et al. 2006 and seems to be proper according to experimental reports by Yoo et al. 2011 . Moreover, eye ball Young modulus and Poisson's ratio are also assumed based on the data reported by Schutte et al. 2006 and its density is calculated to be 1270 kg/m 3 according to the data stated by Robinson 1981 . As stated before, Schutte et al. 2006 considered EOMs as a simple orthogonal elastic material which is always passive. EOMs are modelled in the present work by using the active-passive hyper-elastic model. Among the models proposed so far as a constitutive model that simulates both active and passive behaviours of the muscles, only a few of them can model dynamic characteristics of the muscles Tang et al. 2009, Ehret et al. 2011, Lu et al. 2010 . The constitutive model presented by Ehret et al. 2011 considers the muscle as a continuous biological material whose properties depends mainly on the activation parameter. This model is employed in this work. Strain energy function is proposed as where I  and K  denote generalized invariants and 0, 0, 0       are material constant parameters used in defining these invariants. Additionally, C is Cauchy-Green tensor and    M m m is structural tensor where m denote the preferred orientation unit vector and  is the dyadic product symbol. In these equations, a w is used to define muscle activation and 0 w is used as a weight factor to denote the isotropic tissue constituent. More details can be seen in Ehert &Itskov 2007 andEhret et al. 2011 .
in equation 1 should be computed according to experimental data for simulating EOMs. According to the test procedure described by Robinson et al. 1969 , the reported muscle force almost in fibre direction. One can compute the Piola-Kirchhoff stress tensor in the fibre direction as can be found in the following paragraphs.   Figure 3: Comparing exploited muscle model passive response and experimental data reported in Robinson et al. 1969 .
Finally, nonlinear least square method exploited to fit equation 2  . Figure 3 displays the comparison between muscle model response in fibre direction and experimental data. The muscle density is assumed to be 1870 kg/m 3 as stated by Lockwood-Cooke et al. 1999 . Additionally, since available data for EOM pretension is limited, we assume pretension in rectus muscles to be 80mN. This assumption for medial and lateral rectus muscles is adopted from the work by Simonsz & Spekreijse 1996 . However, for superior and inferior rectus muscles, due to lack of data, we will assume similar properties for the lateral and medial rectus muscles. The result of this assumption on eye motion will be investigated later on.

Finite element simulations
In modelling of biomechanical systems, validating the proposed model by the available experimental data is an important step to make sure that the model is properly and accurately constructed. In order to do this, we will simulate force duction test and 14 degree saccade motion. Force duction test simulation results can be used to check the geometry and material assumptions considered for the modelling of eye orbit section 4.2 . Furthermore, in section 4.3, 14 degree saccade simulation is performed as a forward dynamic simulations in order to check active eye motion properties. 14 degree saccade motion is chosen to simulate because experimental data of the EOM force during this saccade motion is available.
After ensuring that the eye model force and motion are in accordance with experimental data, it is possible to use the modelling for various analysis. Orbital fat deformation is checked with experimental data in three different cases in subsection 4.4. In subsections 4.5 and 4.6 different simulations are reported to enhance our understanding about visuomotor dependencies on eye's part mechanical characteristics.

Sketch of finite element modelling
According to the complexity of the system, dynamic FE simulation of the system is too time consuming. However, this issue is inevitable to obtain a proper representation of mass and inertia besides accurate results. In this work, the eye orbit geometry is discretized to 31158 tetrahedral elements. Assumptions mentioned in sections 2 and 3 are considered to construct following matrices: where n M denotes the mass matrix and n i and n p are the internal force and the external force vectors for the model. In equations 4 -6 , "n" is the element number n 1, …, 31158 , N is the element shape function and B denotes the element strain tensor. In 5 , σ is the Cauchy stress tensor and in 6 , g represents the gravitational vector and t is the surface traction vector. Element matrices will be assembled to construct; M , i and p in 7 .
More details about finite element procedure can be found in various nonlinear finite element book such as the one by Borst et al. 2012 .

Modelling of force duction test
At the first step, in order to validate our model, a prevalent clinical approach to observe the behaviour of eye orbit is required. The force duction test is done in order to see whether patients' eye ball can move in eye orbit regardless of the muscle excitation. Simulation is done by applying 200mN force as shown in Figure 4. The experimental test done by Collins et al. 1981 have shown the rotational stiffness of 0.09 to 0.19 mNm/deg when exerting the force nasally and 0.08 to 0.13 when exerting the force temporally. We have simulated force duction test both temporally and nasally. The stiffness is 0.12 mN/deg nasally and 0.16 mN/deg temporally. These results are acceptable according to the experimental results by Collins et al. 1981 . Consequently, the model geometry seems to be consistent with available real human eye data.

Modelling of eye active motion
For the first time, we use a constitutive hyper-elastic model whose activation varies with time In this subsection, we will see if the characteristics of our active motion simulation are similar to the available experimental results.

Fourteen degree eye saccade motion
In order to validate the EOM force in eye motions, we will use the available saccade data reported by Robinson 1964 andRobinson et al. 1969 . Neurologists often propose open loop control for saccade, which will be considered in the following simulation one can see for example the book by Chalam et al. 2010Chalam et al. -2011 . Herein, the medial rectus activation will be increased from 0.1129 to 0.5 in 20 milliseconds. Activation time period is assumed to be 20 milliseconds according to the force profile reported in Robinson et al. 1969 . Initial and final value for muscle activation are chosen by try and error in order to acquire the desired motion.
The forward dynamics simulation resulted displacement is shown in Figure 5. Eye has rotated 14 degree horizontally in approximately 0.025 sec. The eye trajectory is in consistent with the experimental results reported in Robinson 1964 .
The corresponding medial rectus force for this motion is shown in Figure 6. Robinson et al. 1969 have shown experimentally that in 15 degree saccade, muscle tension is about 450mN. This shows that the value obtained in our simulation is sufficiently comparable with the available experimental reports. The chattering in Figure 6 is due Latin American Journal of Solids and Structures, 2018, 15 3 , e24 7/12 to the nonlinearities and implicit integration. Thus, we can realize that the modelling gives reliable muscle force and eye motion.

Eye suspension deformation
In this subsection, we will analyze the eye suspension deformation to evaluate the validity of the assumption for the elastic modulus in the modelling of eye's motion. To this end, the findings reported based on MRI scans by Schoemaker, et al. 2006 will be employed. They have reported that in the distance range of 1mm to 5mm from sclera, the studied points in the orbital fat have rotated about 53% to 37% of the total eye rotation and the studied points posterior in orbit which are 15mm to 22mm from sclera have rotated about 18% to 7% of the total eye rotation. As shown in Figure 7, the simulated deformation in the orbital fat with a 300Pa Young modulus assumption is in consistent with the aforementioned experimental data. Hence, we can conclude that the assumed mechanical properties are suitable for modelling of eye motion.
Since the data reported by Schoemaker, et al. 2006 is based on 2D MRI images and their calculated rotations ratio are based on a planar motion assumption, to have a deeper look at the problem, we will investigate both 2D calculated based on a planar motion and 3D based on a spatial motion rotational motions. For this purpose, we will compare both 2D and 3D rotational motion of some of the points studied in Figure 7. The acquired data are reported in Table 2. Since the differences between 2D and 3D angles in Table 2   What we know about the orbital fat mechanical properties is limited. In this subsection, we are intended to see the effect of one of these properties, such as its bulk elastic modulus, on the eye's motion. As we have shown in the previous subsection 4.4, the intraconal retrobulbar fat deformation with the assumed elastic modulus is acceptable according to the experimental results by Schoemaker, et al. 2006 . We simulated three cases with 200Pa, 300Pa and 400Pa to see the effect of the bulk elastic modulus on eye's motion and orbital fat's deformation. In these simulations, other assumptions are similar to those of the previous subsections for forward dynamics simulations.
As shown in Figure 8, the muscle force during the simulations for different cases are in acceptable range according to the experimental results reported in Schoemaker, et al. 2006 andCollins et al. 1981 . Figure 9 shows the resulted motion, both horizontally and vertically, in these three situations. Motion paths are completely similar and show that the orbital fat Young modulus has minor effect on eye motion. Figure 9: Eye motion in "x" and "y" directions for different orbital fat elastic modulus vs. time.
As shown in Table 3, variations of the elastic modulus have approximately 0.0001 rad effect on cornea rotation, 0.007 rad effect on deformation of points in 15 to 22mm range behind sclera and maximum 0.002 rad effect on rotation of points in 1 to 5mm range behind sclera. We should note that the major disorders corresponding to the orbital fat such as "exophthalmos" are mostly related to the orbital fat geometry and the muscle thickness. This issue confirms the above simulation result which states that the orbital fat mechanical properties does not have significant effect on eye ball motion. In this subsection, EOM pretension is studied through the proposed model. To the best of our knowledge, muscle pretension magnitude is not obtained experimentally so far. As mentioned earlier, it is common to assume medial and lateral rectus muscles pretension equivalent to 80Nm Spekreijse 1996 andSchutte et al. 2006 . In this subsection, we are trying to study the effects of this assumption. While the medial and lateral rectus muscles pretensions are the same as previous simulations 80mN , superior and inferior rectus muscles pretension are increased to 100mN to see their influence on the motion. Figure 10 and Figure 11 show that by increasing pretension in the superior and inferior rectus muscles, the magnitude of the vertical displacement is increased while they do not have any effect on the eye horizontal motion. Consequently, the model can be employed along with accurate experimental data of eye trajectory to estimate muscle pretension. One may record various eye motions of a human and then uses an identification method to compute the corresponding EOMs pretension according to the experimental data. Figure 10: Influence of the superior and inferior rectus muscles pretension on eye vertical displacement. Figure 11: Influence of the superior and inferior rectus muscles pretension on eye horizontal displacement.

Summary of FE simulations
In this section, we reviewed the computation of the FE model matrices subsection 4.1 . Then the eye orbit geometry, EOMs force and eye motion validated by simulating the experiments which were previously done in real human eye. In subsection 4.2, force duction test is reported to check the model stiffness and geometry and in subsection 4.3, saccadic motion simulation results are discussed which include muscle force and eye motion. Then, orbital fat deformation in multiple active motion simulation was compared with experimental results subsection 4.4 . In all three subsections, the results are completely consistent with the real human eye motion behaviour. Finally, orbital fat and EOMs effect on eye motion was studied. These two subsections show the model applications to obtain a better comprehension of different eye parts effect on the visuomotor function.

Conclusions
Although utilizing continuum-mechanics methods such as FE for human eye tissues with nonlinear properties such as hyper-elastic active-passive material assumption is complicated and time consuming, it helps us to obtain more accurate results and better comprehension. In this work, multiple dynamic simulations of eye FE model were reported. Validation of our model was done by simulating force duction test and saccadic motion. Our assumptions seem to be appropriate for eye modelling since the corresponding stiffness, resulted by force duction test simulation, and EOMs force obtained by saccadic motion simulation was completely valid according to the experimental data reported in previous works. Moreover, the model employed for studying the eye parts such as orbital fat in eye orbit. The resulted orbital fat deformation in various points were compared with those of the MRI scan from a real case and the computed deformations were in an acceptable range. Orbital fat effect on eye motion was also investigated through multiple simulations. We observed that the eye motion, with different intraconal retrobulbar fat elasticity, were almost the same. It is noteworthy that the orbital fat is not convenient for experimental studies. Finally, muscle pretension effect on eye motion was studied and a new method was proposed in order to find out the muscle pretension through FE modelling.
Our simulation clearly showed the FE modelling usefulness for computing different mechanical characteristics of the eye orbit according to the convenient experimental data by using the model in inverse problem.